Broadband Study of Gamma-Ray Blazars at Redshifts $z=2.0-2.5$

High redshift blazars are among the most powerful non-explosive sources in the Universe and play a crucial role in understanding the evolution of relativistic jets. To understand these bright objects, we performed a detailed investigation of the multiwavelength properties of 79 $\gamma$-ray blazars with redshifts ranging from z = 2.0 to 2.5, using data from Fermi LAT, Swift XRT/UVOT, and NuSTAR observations. In the $\gamma$-ray band, the spectral analysis revealed a wide range of flux and photon indices, from $5.32 \times 10^{-10}$ to $3.40 \times 10^{-7}$ photons cm$^{-2}$ s$^{-1}$ and from 1.66 to 3.15, respectively, highlighting the diverse nature of these sources. The detailed temporal analysis showed that flaring activities were observed in 31 sources. Sources such as 4C+71.07, PKS 1329-049, and 4C+01.02, demonstrated significant increase in the $\gamma$-ray luminosity and flux variations, reaching peak luminosity exceeding $10^{50}$ erg s$^{-1}$. The temporal analysis extended to X-ray and optical/UV bands, showed clear flux changes in some sources in different observations. The time-averaged properties of high redshift blazars were derived through modeling the spectral energy distributions with a one-zone leptonic scenario, assuming the emission region is within the broad-line region (BLR) and the X-ray and $\gamma$-ray emissions are due to inverse Compton scattering of synchrotron and BLR-reflected photons. This modeling allowed us to constrain the emitting particle distribution, estimate the magnetic field inside the jet, and evaluate the jet luminosity, which is discussed in comparison with the disk luminosity derived from fitting the excess in the UV band.


INTRODUCTION
Blazars are jetted active galactic nuclei (AGN) and are among the most powerful persistent sources of electromagnetic radiation in the universe.Current unification theories assumes that blazars are a subtype of AGNs with jets oriented at a small angle relative to the observer's line of sight (Urry & Padovani 1995).Their emission is believed to be primarily powered by the accretion of matter into supermassive black holes (Blandford & Znajek 1977) and exhibit extreme characteristics, such as high-amplitude and short-timescale variability, core dominance, superluminal motion, and significant optical polarization.Blazars are typically grouped into two main groups based on the presence of emission lines in their spectra (Urry & Padovani 1995): BL Lacertae objects (BL Lacs), which have weak or absent emission lines, and flat spectrum radio quasars (FSRQs), which display strong emission lines.
The nonthermal emission from the jets of blazars is observable across almost all accessible bands of the electromagnetic spectrum (Padovani et al. 2017) up to high energy (HE; > 100 MeV) and (VHE; > 100 MeV) -ray bands.Their broadband spectral energy distribution typically shows two broad humps.The lower-energy component, from the radio to the optical/X-ray band, is generally attributed to the synchrotron emission of electrons.The frequency of the synchrotron peak ( p ) serves as a criterion to further classify blazars: as low synchrotron peaked sources (LSPs or LBLs), intermediate synchrotron peaked sources (ISPs or IBLs), or high synchrotron peaked sources (HSPs or HBLs) when  p < 10 14 Hz, 10 14 Hz <  p < 10 15 Hz, and ★ E-mail: narek@icra.it p > 10 15 Hz, respectively (Padovani & Giommi 1995;Abdo et al. 2010).The origin of the second component, however, remains a subject of debate.In leptonic scenarios, this second peak is interpreted as the result of inverse Compton scattering of low-energy photons, which may be of internal or external origin.In the synchrotron self-Compton (SSC) models (Ghisellini et al. 1985;Maraschi et al. 1992;Bloom & Marscher 1996), it is the synchrotron photons that are upscattered through inverse Compton processes.On the other hand, external inverse Compton (EIC) scenarios (e.g., Sikora et al. 1994) suggest that the photons originate outside the jet, coming directly from the accretion disk (Dermer et al. 1992;Dermer & Schlickeiser 1994), or they may be reprocessed by the broad-line region (BLR) (Sikora et al. 1994), or emitted from the torus (Błażejowski et al. 2000).In Bégué et al. (2023), a novel approach for fitting the blazar Spectral Energy Distribution (SED) utilizing convolutional neural networks is introduced which allows self-consistent modeling, enabling a more detailed interpretation of the observed results.
In alternative models such as the hadronic or lepto-hadronic scenarios, the second spectral bump is assumed to be from the synchrotron emission of ultra-HE protons or from the decay of secondary particles produced during hadronic interactions (Mannheim 1993;Mannheim & Biermann 1989;Mücke & Protheroe 2001;Mücke et al. 2003;Böttcher et al. 2013;Petropoulou & Mastichiadis 2015;Gasparyan et al. 2022).Interest in these models has increased, particularly after establishing a potential link between blazars and VHE neutrinos, following the association of TXS 0506+056 with the IceCube-170922A neutrino event (IceCube Collaboration et al. 2018a,b;Padovani et al. 2018), and the detection of multiple neutrino events concurrent with the active phase of PKS 0735+178 in optical/UV, X-ray, and -ray bands (Sahakyan et al. 2023a).These multimessenger observations have started extensive discussions, with various models being applied to explain the multimessenger observations of blazars (Ansoldi et al. 2018;Keivani et al. 2018;Murase et al. 2018;Padovani et al. 2018;Sahakyan 2018;Righi et al. 2019;Cerruti et al. 2019;Sahakyan 2019;Gao et al. 2019;Gasparyan et al. 2022).
The emission from blazars is highly beamed, and their bolometric luminosity can exceed 10 48 erg s −1 , allowing them to be observed even at very high redshifts (e.g., see Ackermann et al. 2017;Sahakyan et al. 2020).These distant blazars are particularly interesting, as their study offers insights into the formation and evolution of supermassive black holes, relativistic jets, and the connections between accretion disks and jets.Moreover, their -ray emission is important for probing the early universe; -ray emission from distant blazars undergoes attenuation via  absorption when interacting with extragalactic background light (EBL) photons, thereby enabling observations that can constrain the EBL's density.The multiwavelength properties of distant blazars have been extensively studied in a number of publications (e.g., see Sahakyan et al. 2023bSahakyan et al. , 2020;;Li et al. 2020;Liao et al. 2019;Marcotulli et al. 2017;Orienti et al. 2016;D'Ammando & Orienti 2016;Paliya et al. 2015;Pacciani et al. 2012;Paliya 2015;Paliya et al. 2016Paliya et al. , 2017Paliya et al. , 2019)).In Sahakyan et al. (2020), the origin of emission from the most distant blazars detected in the HE -ray band ( > 2.5) was investigated analyzing data in the optical/UV, X-ray, and -ray bands.From the temporal evolution of emission in these bands, flaring periods were identified when the luminosity substantially increased.Their broadband emission was modeled using one-zone SSC and EIC models, assuming that the external photons are infrared (IR) photons from the dusty torus (see Arsioli & Chang (2018) for a discussion on the contribution of different external fields in shaping the HE emission of LSPs).As a result, the parameters that characterize the particle emission and jet in these distant blazars were estimated.
In this work, we expand the study by Sahakyan et al. (2020) to -ray emitting blazars which have estimated redshifts that are in the range of  = 2.0 − 2.5.A prolonged observation period of ∼ 14 years enables a comprehensive investigation of emissions from these sources across optical/UV, X-ray, and -ray bands as well as to perform detailed spectral and temporal analysis.Using the analyzed data broadband SEDs for a substantial number of sources were constrained and modeled which allows for a systematic comparison of blazar emission parameters at varying distances, which, in turn, could enhance our understanding of these brightest objects.
The structure of the paper is as follows: Section 2 introduces the sample of sources under consideration; Section 3 details the data analysis methodology; Section 4 presents the results of the data analysis; Section 5 discusses the modeling of broadband SEDs; and Section 6 outlines the modeling results.Finally, the conclusions are summarized in Section 7.
In summary, the source sample considered in this study comprises 79 objects.In Table 1, the three leftmost columns list these sources along with their 4FGL names and classes.The spatial distribution of these blazars in Galactic coordinates is depicted in the Hammer-Aitoff projection shown in Fig. 1.For completeness, this figure also includes blazars with redshifts greater than  > 2.5, as reported by Sahakyan et al. (2020).

MULTIWAVELENGTH OBSERVATIONS OF CONSIDERED BLAZARS
In order to investigate the multiwavelength properties of the selected sources, data collected by the Fermi-LAT, Nuclear Spectroscopic Telescope Array (NuSTAR), Neil Gehrels Swift Observatory (hereafter Swift) X-Ray Telescope (XRT), and Ultraviolet/Optical Telescope (Swift UVOT) were downloaded and analyzed.

Fermi-LAT data
The Large Area Telescope (LAT), onboard the Fermi Gamma-ray Space Telescope, is a HE instrument that uses the pair-production technique to detect -rays in the energy range between 20 MeV and > 300 GeV.By default, it operates in scanning mode, continuously monitoring the -ray sky since its launch in 2008 (Atwood et al. 2009).
The Fermi-LAT PASS8 data collected from August 4, 2008, to December 4, 2022 (∼ 14.5 years), were considered to study the properties of all 79 blazars under consideration.The standard datareduction procedure was performed following the recommendations from the Fermi-LAT science team 1 .For each blazar, events in the energy range between 100 MeV and 500 GeV from a region of interest (ROI) of 12 • -reduced to 10 • for several sources to better represent the ROI-centered on the -ray position of the sources, were downloaded and analyzed.The Fermi ScienceTools version 2.0.8 and the P8R3_SOURCE_V3 instrument response function were used.To reduce contamination from the Earth's limb, a zenith angle cut of 90 • was applied.Events with a higher probability of being photons were selected using the filter evclass = 128 and evtype = 3, and the good time intervals were chosen with the expression (DATA_QUAL > 0)&&(LAT_CONFIG == 1).The analysis model file was created based on the Fermi Fourth Source Catalog (4FGL-DR3 Abdollahi et al. 2022) and includes all sources within an ROI radius plus an additional 5 • .During the likelihood fitting, the spectral parameters of all sources within the ROI were allowed to vary, while those of sources outside the ROI were fixed to their 4FGL values.The model file also includes the Galactic diffuse emission model gll_iem_v07 and the isotropic component iso_P8R3_SOURCE_V3_v1.The spectral parameters of the sources, along with the normalization of both background models, were optimized by applying a binned likelihood analysis using the fermiPy tool (Wood et al. 2017).The detection significance of the sources is quantified by the likelihood test statistic (TS), defined as   = 2 × (log  − log  0 ), where  is the likelihood with the source at the position of interest included, and  0 is the likelihood without the source.
The -ray variability of the considered sources was investigated by generating light curves in two distinct manners.Initially, for all sources, the 14-year period was divided into equal intervals (e.g., 5, 7, 10 days, depending on the source's overall detection significance) 1 https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/ to ensure that the light curves did not contain a significant number of upper limits.Within these intervals, the flux and photon index were estimated by applying the unbinned likelihood analysis method.This approach provides a general overview of the flux changes over time, but since the fluxes are averaged over several days, any potential short-scale flux variations are likely to be smoothed out.For a more detailed examination of flux evolution over time, light curves were also generated using the adaptive binning method (Lott et al. 2012), which is applied when photon statistics are sufficient.This technique allows for flexible time bin widths that are determined by assuming a constant uncertainty in the flux estimation.Consequently, brighter source states yield shorter bins, whereas longer bins are used during lower and/or average source states.Light curves produced by this method have been extensively utilized to study short-timescale flux variations in blazar emissions (e.g., see Rani et al. 2013;Britto et al. 2016;Sahakyan & Gasparyan 2017;Zargaryan et al. 2017;Baghmanyan et al. 2017;Gasparyan et al. 2018;Sahakyan et al. 2018;Sahakyan 2021;Sahakyan et al. 2022b;Sahakyan & Giommi 2022).
The analysis of all NuSTAR data was performed using the NuS-TAR_Spectra pipeline, a shell script built upon the NuSTAR Data Analysis Software (Middei et al. 2022).This script streamlines the process by autonomously retrieving calibrated and filtered event files, employing Ximage for accurate source positioning, and utilizing the nuproducts command to retrive science ready products.Source counts are extracted from a predefined circular region, while background counts are from an annular region, with inner and outer radii sizes dynamically determined by the source's count rate.Post data selection, the spectra are binned to ensure a minimum of one count per bin, and spectral fitting is executed within the XSPEC framework (Arnaud 1996), applying Cash statistics (Cash 1979) for the energy range from 3 keV to the upper energy limit of detectable signal, which varies between 20 and 79 keV.The fitting is performed using both power-law and log-parabola models, the observed flux in the 3-10 and 10-30 keV bands is estimated and the corresponding SED is computed using the best-fit parameters.Details describing NuSTAR_Spectra pipeline can be found in Middei et al. (2022).

Swift XRT
The high-redshift blazars selected for this study were also frequently monitored by the Swift XRT in the 0.3-10 keV energy range.Of the 79 considered sources, 60 were observed at least once by the Swift telescope.The source PKS 0528+134 was the most frequently observed-141 times-while 4C +71.07,PKS 2149-306, PKS 0226-559, S3 0458-02, 4C +01.02, and PMN J1344-1723 each were observed more than 20 times.
All the Swift XRT data accumulated from the observations of selected sources was retrieved and analyzed using swift_xrtproc script.This automated tool accesses both PC and WT mode observations from Swift XRT and executes the standard analysis procedure which includes the creation of exposure maps, calibration of observational data.The source spectral files are obtained by estimating the source counts within a 20-pixel radius circular region and the background from a surrounding annular region centered on the source.Additionally, it performs corrections for pile-up effects and conducts spectral fitting employing both power-law and log-parabola models within the XSPEC framework.Subsequently, it computes the SED spectral points using the optimal spectral model, computes the flux across specified energy intervals, and estimates the photon index for the selected energy range.For more details on the swift_xrtproc tool see Giommi et al. (2021).
The majority of these sources considered here exhibit no significant variability in the X-ray band, as discussed in the next section, so to enhance the photon statistics and refine the estimation of the X-ray flux, if for a source multiple observations are available, they were combined and analyzed using the tool provided by the UK Swift Science Data Centre (Evans et al. 2009).

Swift UVOT
Together with the XRT, the sources were also observed by Swift UVOT, which provided data in three optical filters (V, B, and U) and three UV filters (W1, M2, and W2).All available Swift UVOT data from the observations of the considered sources were downloaded and analyzed.
The data were analyzed using the uvotsource task included in the HEAsoft package, version 6.29.Source counts were extracted from a circular region with a radius of 5" centered on the source, while background counts were obtained from a larger circular region with a radius of 20", located in a nearby source-free area.The observations of all sources were individually checked to ensure the accuracy of source and background region selection.For all sources, the magnitudes were derived using the uvotsource tool and corrected for reddening and Galactic extinction using the reddening coefficient  (−) obtained from the Infrared Science Archive 2 .The corrected 2 http://irsa.ipac.caltech.edu/applications/DUST/fluxes measured for each filter were used to construct the light curves and SEDs.As the considered sources are at high redshift, their optical/UV fluxes could be affected by absorption from neutral hydrogen in intervening  −  absorption systems.These effects were corrected for in the SEDs during theoretical modeling, following the procedures described in Ghisellini et al. (2011) and Ghisellini et al. (2010).

Archival data
To construct the most comprehensive multiwavelength SEDs possible for the considered sources, archival data were also extracted and analyzed in addition to the data discussed herein.This was accomplished using the VOU-Blazar tool (Chang et al. 2020) through Markarian Multiwavelength data center3 , which retrieves multiwavelength data from 71 catalogs and spectral databases through various online services.

RESULTS OF DATA ANALYSES
In this section, a comprehensive spectral and temporal analysis of the sources included in our sample are performed.The outcomes of the -ray data analysis are provided in Table 1.For each source, the power-law index ( p ), or () together with the curvature parameter () when the data are best described using a low-parabola model, as well as the flux and luminosity (computed from the power-law fit), both with their respective uncertainties and the redshift of each source is reported.When the source spectrum is described by lowparabola model in 4FGL an additional analysis is performed using power-law model.
The analysis results are presented in panel a of Fig. 2, which depicts the flux and photon index derived from a power-law fit.The FSRQs are shown by blue markers, BL Lacs by orange, and BCUs by green.The -ray flux of sources within our sample ranges from (5.32 ± 2.25) × 10 −10 to (3.40 ± 0.02) × 10 −7 photon cm −2 s −1 , with the lowest value corresponding to SDSS J100326.63+020455.6 and the highest to 4C +01.02.The mean flux is at 2.48×10 −8 photon cm −2 s −1 .The photon index span from 1.66 ± 0.12 to 3.15 ± 0.12, with the lowest and highest indices observed for SDSS J100326.63+020455.6 and PKS 1915-458, respectively.The FSRQs, which are the most numerous in our sample, largely define these observed ranges.On the other hand, BCUs and BL Lacs exhibit narrower distributions in both -ray flux and photon index.Specifically, BL Lacs range between (0.18−1.86)×10 −8 photon cm −2 s −1 for flux and 1.91−2.62for photon index, while BCUs between (0.05 − 1.60) × 10 −8 photon cm −2 s −1 for flux and 1.66 − 2.73 for photon index.For comparison, blazars with a redshift exceeding 2.5 from Sahakyan et al. (2020) are depicted in light gray in panel a of Fig. 2.This visualization shows that blazars with a redshift beyond 2.5, as well as those included in the current sample, exhibit similar features, occupying similar regions on the photon index versus   plane.
In panel d of Fig. 2, a comparison of the -ray and X-ray fluxes for the selected sources is shown.The wide spread observed in the data suggests that there is no direct or obvious correlation between the ray and X-ray fluxes when considering time-averaged measurements.It should be noted, however, that these are average values and that during shorter time-scale events, such as flares, a correlation may appear.Interestingly, the two bright in the X-ray band sources, 4C +71.07 and PKS 2149-306, also have notably high -ray fluxes of (1.11 ± 0.02) × 10 −7 photon cm −2 s −1 and (7.30 ± 0.20) × 10 −8 photon cm −2 s −1 , respectively.Conversely, the -ray bright source 4C +01.02, has only a moderate X-ray flux of (2.18 ± 0.11) × 10 −12 erg cm −2 s −1 , indicating that a high -ray flux does not necessarily imply a correspondingly high X-ray flux.This discrepancy shows the complexity of the emission mechanisms and the potential influence of other factors such as beaming, the environment of the source, or the presence of different emission processes at different wavelengths.
The NuSTAR analysis results are presented in Table 2, where for each source the observation sequence, observation time, and flux in the 3 − 10 keV and 10 − 30 keV ranges, along with the photon index are provided.Three NuSTAR observations (60160099002 for S5 0212+73, 60002045002 for 4C +71.07, and 60367002002 for PKS 0227-369) were relatively short (on the order of a few hundred sec-onds), and hence no spectral analysis was conducted.In the 3 − 10 keV range, the highest flux of (2.02 ± 0.01) × 10 −11 erg cm −2 s −1 was observed for 4C +71.07 on MJD 56675.22,while the lowest flux of (1.95 ± 0.27) × 10 −13 erg cm −2 s −1 was observed for 87GB 080551.6+535010.The photon index for all considered sources is hard, ranging from 1.09 to 1.67, suggesting that the hard X-ray component corresponds to the rising part of the inverse Compton component.The variability of the 3 − 10 keV and 10 − 30 keV fluxes could only be investigated for 4C +71.07 and PKS 2149-306, as multiple observations in different periods are available; however, the flux remained relatively stable.

𝛾-ray variability
The use of adaptive binning methods for computing light curves allowed detailed investigation of the -ray flux variation of considered sources.Flux variations, characterized by several times increases from average levels, have been observed in 31 sources.The most variable sources in the -ray band, where multiple flaring activities have been observed, are PKS 0226-559, PKS 2149-306, S3 0458-02, 4C+71.07,S4 0917+44, PKS 1329-049 and 4C+01.02.For example, the adaptively binned -ray light curve of PKS 0226-559, the most distant source in our sample at a redshift  = 2.464 showing variability, computed for energies above 239.66MeV, is displayed in panel a) of Fig. 3. From the start of Fermi-LAT observations until MJD 55504 (November 04, 2010), the -ray emission of this source was in a low state being (4.89 ± 0.71) × 10 −9 photon cm −2 s −1 .Subsequent periods of flaring activity occurred between MJD 55504. .During its most active state, between MJD 58084.93 and 58664.52, the maximum flux was (3.99 ± 0.67) × 10 −7 photon cm −2 s −1 observed on MJD 58507.69.Also, the luminosity of the source significantly increased during these flaring periods.While the long-term averaged luminosity was (1.76 ± 0.04) × 10 48 erg s −1 , the luminosity during flares exceeded 10 49 erg s −1 .For example, the source luminosity was > 10 49 erg s −1 42 times, with the peak luminosity of (5.33 ± 0.92) × 10 49 erg s −1 observed on MJD 58167.46.

X-ray variability
In the considered source sample, some of the sources were observed multiple times by the Swift telescope, permitting investigation of the temporal variability of X-ray flux across different years.Among the sources considered, variability in X-ray flux was observed in 4C +71.07,PKS 0226-559, PKS 1329-049, PKS 2149-306, S3 0458-02 and S4 0917+44.For PKS 0226-559 and S4 0917+44, a limited number of observations were available but the ratio between the maximum and minimum fluxes is approximately 3 − 4, indicating temporal flux variability.

Variability in optical/UV bands
The observations from the Swift UVOT of selected sources allows to study the flux variability within the optical/UV bands.Investigating variability is challenging when the number of observations is limited, as changes in flux up to a factor of 2 can be observed across different filters; however, this does not provide a comprehensive understanding of the flux changes in time.Notably, clear flux variability is evident in the emissions from 4C+01.02,PKS 0226-559, and PKS 2149-306.For 4C+01.02, the initial flux measurements in the , , and  filters are approximately 10 −12 erg cm −2 s −1 , and approximately 2 × 10 −13 erg cm −2 s −1 in the 1, 2, and 2 fil-ters.After MJD 56500, the source exhibits increased flux across all filters during several observations.The highest observed flux in the  filter was (3.91 ± 0.33) × 10 −12 erg cm −2 s −1 on MJD 57363.45.For PKS 0226-559, the mean flux is ∼ 5 × 10 −13 erg cm −2 s −1 , but the flux in  filter increased up to (3.1 − 0.46) × 10 −12 erg cm −2 s −1 during the flare on MJD 58161.21.The UV emission from PKS 2149-306, in the 1, 2, and 2 filters, remains relatively stable across various observations, whereas the optical emission (, , and ) shows variability.Specifically,on MJD 53717.92,MJD 54980.86,and MJD 58586.29, the flux increased to approximately 3 × 10 −12 erg cm −2 s −1 , with the lowest observed flux being around 10 −13 erg cm −2 s −1 .

MODELING OF MULTIWAVELENGTH SEDS
The data analyzed in this study, combined with those collected from the archives, enabled to construct the SEDs for the selected sources from the radio to the HE -ray bands.The multiwavelength SEDs are presented in Fig. 6 where the data analyzed in this work are highlighted in different colors (see the legend), while archival data retrieved using the VOU-Blazar tool are depicted in gray.All the 10 10 10 12 10 14 10 16 10 18 10 20 10 22 10 24 [Hz] 10 14    SEDs exhibit the traditional double-humped structure with a clear Compton dominance (i.e., the Compton peak luminosity is significantly higher than the synchrotron peak luminosity).
The selected SEDs were modeled using a leptonic one-zone synchrotron and inverse Compton model.In this scenario, the emission region is assumed to be located at a distance  dis from the central black hole and is assumed to have a spherical geometry with a radius .This region moves along the jet with a bulk Lorentz factor Γ and is observed at a small viewing angle, resulting in the amplification of radiation by a factor of  ≈ Γ.The emitting region is populated with nonthermal particles (electrons and positrons), whose energy distribution follows a power law with an exponential cutoff described as: where  is the spectral index of the power-law distribution for the emitting electrons,  cut is the cutoff Lorentz factor, and  min represents the minimum energy.Here,  0 denotes the normalization constant of the electron distribution, which defines the total energy content of the electrons as given by  e =  2 ∫  ().Within the emitting region (blob), electrons lose energy through synchrotron emission under the magnetic field , resulting in the first bump in the multiwavelength SED.The second bump is attributed to inverse Compton scattering of low-energy photons, which may originate either internally or externally to the jet.Synchrotron photons can be inverse Compton scattered to higher energies via SSC radiation (Ghisellini et al. 1985;Maraschi et al. 1992;Bloom & Marscher 1996), and the inverse Compton scattering of external photons may also contribute to the formation of the second component.Depending on the emitting region's proximity to the central black hole, different low-energy photon fields can be inverse Compton up-scattered: photons directly from the accretion disk (Dermer et al. 1992;Dermer & Schlickeiser 1994), photons reprocessed by the BLR (Sikora et al. 1994), or photons from the dusty torus (Błażejowski et al. 2000).
In the current study, we assume that the emitting region is within the BLR ( dis < ), and that accretion disk photons reprocessed by the BLR (EIC-BLR) alongside synchrotron photons are the primary targets for inverse Compton scattering.However, we note that alternative scenarios, wherein the emitting region lies significantly closer to the central black hole or outside the BLR, cannot be excluded.Constraining the location of the emitting region requires high-quality and high-resolution data at low energy bands, which is unavailable for distant blazars and falls outside the scope of this study.
During the modeling, the accretion disk luminosity ( disk ) and the temperature (energy) of the photons were estimated by fitting the UV band's excess emission with a blackbody component.If the thermal component could not be distinguished, an upper limit was established by ensuring that the disk's emission does not outshine the observed nonthermal emission from the jet.With  disk , the radius of the BLR is calculated using the relation  BLR = 10 17  45,disk (Ghisellini & Tavecchio 2015).Subsequently, the BLR is represented as a spherical shell, with an inner boundary of  in,BLR = 0.9 ×  BLR and an outer boundary of  out,BLR = 1.2 ×  BLR .
The SEDs were fitted using the publicly available JetSet code (Massaro et al. 2006;Tramacere et al. 2009Tramacere et al. , 2011;;Tramacere 2020), which facilitates the comparison of numerical models with observed data.For the FSRQs, the synchrotron/EIC-BLR model was used, while for the BL Lacs, the synchrotron/SSC model.Only one BCU has sufficient data for the modeling, SDSS J100326.63+020455 which was also modeling using synchrotron/SSC-BLR modeling considering the shape of the SED.The model consists of a combined contribution from synchrotron, disk, SSC, and EIC-BLR components, which are optimized during the fitting process.The free parameters in the model are the spectral index (), the cutoff Lorentz factor ( cut ), the minimum Lorentz factor ( min ), the magnetic field (), the Doppler factor (), the size of the emitting region (), and the energy density of the electrons ( e ).Optimization is performed in two stages: initial fitting is done with the Minuit optimizer, fol-lowed by refinement using the Markov Chain Monte Carlo method.We note that  can also be constrained by variability considerations.However, due to insufficient data for all sources, variability analysis could not be conducted.Therefore, during the fitting process,  was treated as a parameter dependent on the variability timescale ( var ) and , via the relation  =  var /(1 + ), with  var assumed to be within a 0.1 to 7-day range.However, we refrain from commenting on  var (and consequently ) as it can only be accurately determined with high-quality data.The pair-production absorption effect due to interaction with EBL photons was incorporated using the model of Franceschini et al. (2008).

RESULTS
The modeling results are shown in Fig. 6 and the results are given in Table 3.For this modeling, we selected only those sources with sufficient data, particularly those having the X-ray data, as the HE component cannot be constrained without them.There is significant amplitude variability in some sources across optical, IR, and X-ray bands (e.g., PKS 1329-049, PKS 0446+11, 4C +01.02, etc.), which our current modeling does not adequately interpret.To accurately model these flaring periods, a careful selection of (quasi) contemporaneous data is required.However, this is beyond the scope of the current paper, whose primary aim is to estimate the time-averaged properties of the selected sources.
The EIC-BLR model was used to model all the SED of all FSRQs (46 sources).In Fig. 6, the sum of all components -the combined contribution of the synchrotron, disk, SSC, and EIC-BLR components -is depicted in blue.The emission at lower energy bands (up to X-rays) is dominated by the synchrotron and disk components, while at HEs, the SSC and EIC-BLR components are dominant.The modeling facilitates the estimation of parameters describing the jet (, ) as well as those characterizing the emitting electrons (,  min , and  cut ).From the modeling of the SEDs of the FSRQs, the powerlaw index of the electrons varies within a range from  = 1.56 to  = 2.85, predominantly constrained by the SSC fit of the X-ray data.The hardest spectral indices are observed in sources lacking high-quality X-ray data.For instance,  = 1.56 ± 0.064 for PKS 0226-559,  = 1.60 ± 0.12 for OX 131, and  = 1.70 ± 0.21 for PKS 1348+007, etc.The power-law index of the electrons, inferred from the SSC fitting of BL Lacs (7 sources), is harder than that for FSRQs, ranging between 1.71 and 2.25.Relatively hard indices of 1.71±0.01and 1.79 ± 0.01 were estimated for SDSS J145059.99+520111.7 and NVSS J090226+205045, respectively.
For the FSRQs, the cutoff Lorentz factor ( cut ) for the considered sources varies within a range from 237.6 to 3.53×10 3 , determined by the optical/UV and/or X-ray data.A high value of (3.55±0.49)×10 4 was estimated for PMN J2135-5006, which has sparse X-ray data and is therefore not considered in the discussions.The lowest value of 237.6 ± 21.8 was estimated for TXS 0322+222, as the shape of the low-energy data (below optical) and the optical/UV data impose a strong constraint on this parameter.A comparably high cutoff Lorentz factor,  cut = (3.53± 0.31) × 10 3 , was estimated for PMN J1344-1723.In the modeling of BL Lacs, the cutoff Lorentz factor ( cut ) is (1.53 − 25.41) × 10 4 which is larger than those typically estimated for EIC-BLR model, as the average energy of the synchrotron photons is lower than that of BLR photons.The highest cutoff Lorentz factor,  cut = (2.54 ± 0.05) × 10 5 , was estimated for SDSS J105707.47+551032.2, which is due to the fact that the synchrotron emission of this source extends up to the X-ray band.The minimum Lorentz factor ranges from 1.08 to 207.0 and is given in Table 3, constrained by the requirement that the emission at radio bands does not overshoot the observed data.
For the FSRQs the estimated magnetic field ranges between 1.30 and 15.17 G.The highest value, 15.17 ± 0.90 G, was estimated for PMN J0134-3843 which is attributed to the fact that the synchrotron component of this source exhibits a comparable flux compared to its inverse Compton component.Consequently, when the external photons are boosted into the jet frame (with a  = 24.71for this object), their density increases.This requires the presence of lower density electrons, thereby requiring a higher magnetic field to account for the observed synchrotron emission.In contrast to the modeling of FSRQs, when considering BL Lacs, the magnetic field is found to be in the range of 0.001 to 0.04 G.For the FSRQs, except for the bordering cases of PKS B1412-096, PMN J2135-5006, and PMN J1344-1723 were estimated respectively,  varies within a narrow range of 12.52 to 33.82.This range is typical for blazars.For PMN J1344-1723, a  of 47.18 was estimated due to the pronounced Compton dominance observed in this object.In contrast, a lower  of 5.24 for PKS B1412-096 was estimated because its Compton dominance is nearly 1.Meanwhile,  = 5.45 for PMN J2135-5006 was estimated due to its relatively high cutoff Lorentz factor ( cut ), which is unconstrained by the data.For the BL Lacs, larger  in the range of 28.94 to 49.84 were estimated.This is because the synchrotron component of these objects peaks around 10 13 − 10 14 Hz, so a higher  is required to explain the second peak, which is in the MeV band.
Another key parameter to assess following the modeling is the jet luminosity.The jet power, in the form of electrons and magnetic field, is computed as   =  2  Γ 2   and   =  2  Γ 2   , respectively.Assuming one proton per relativistic emitting electron and that protons are 'cold' in the comoving frame, the total jet luminosity is calculated as  jet =   +   +   .The corresponding   and   luminosities are presented in Table 3.For the FSRQs,   varies within the range from 7.24×10 43 erg s −1 to 2.94×10 45 erg s −1 , while   ranges from 2.72 × 10 42 erg s −1 to 4.06 × 10 45 erg s −1 .Among the considered sources, 33 exhibit an   /  ratio in the range of 0.1 to 10, suggesting that the emitting region is close to equipartition.For PKS B1412-096 the ratio   /  = 0.05 is estimated which is because of the synchrotron component is at the same level as the inverse Compton component.Conversely, for OX 131, 4C +71.07, and PKS 2149-306, a significantly higher ratio with   /  > 100 is estimated, due to the dominance of the SSC component over the synchrotron component.In the modeling of BL Lacs,   is in the range of (0.04 − 3.85) × 10 47 erg s −1 , which is higher than the   , ranging from (0.10 − 3.92) × 10 43 erg s −1 .The relationship between the total jet luminosity ( jet ) and redshift is illustrated in Fig. 7 upper panel, where our current source sample (FSRQs in blue and BL Lacs in orange) is compared with sources from Sahakyan et al. (2020) (gray), all of which are FSRQs.Clearly, BL Lacs exhibit different properties and tend to have higher luminosities.However, discussions about the distribution of BL Lacs remain inconclusive.This is because BL Lacs are rarely identified at distances  > 2.0 but a significant fraction of BL Lacs lack robust redshift determination, suggesting that BL Lacs might indeed be present at larger redshifts.In contrast, FSRQs show a similar luminosity range to those at redshifts  > 2.5.For the FSRQs the highest luminosity of 6.51 × 10 45 erg s −1 was observed for 4C +71.07 at a redshift of  = 2.218.Upon comparing both samples, it is evident that the jet luminosity tends to increase with source distance, with two notable exceptions: PKS 1351-018 and GB 1508+5714, which have the highest redshifts of  = 3.72 and  = 4.31, respectively.
The fit presented in Fig. 6 also enables the estimation of disk lumi- nosity for cases where a thermal blue-bump component is observable.This excess is evident in a total of 37 sources.The disk luminosity,   , varies from 5.60 × 10 44 erg s −1 to 3.97 × 10 47 erg s −1 .The relationship between disk luminosity and redshift for -ray blazars with  > 2.0 is depicted in the middle panel of Fig. 7. High disk luminosities of 1.25 × 10 47 erg s −1 , 2.32 × 10 47 erg s −1 , 3.35 × 10 47 erg s −1 , and 3.97×10 47 erg s −1 were estimated for PKS 2149-306, 4C +71.07, 2MASS J16561677-3302127, and B2 0552+39A, respectively.Furthermore, the relationship between disk luminosity and jet luminosity is examined in the lower panel of Fig. 7.This relationship can be described by the equation log 10 ( jet ) = 0.22 × log 10 ( D ) + 34.60, represented by the green line in Fig. 7, with the 2 uncertainty region shown in light green.Notably, sources such as PKS 2149-306 and 4C +71.07 exhibit disk luminosities exceeding 10 47 erg s −1 and jet luminosities exceeding 10 45 erg s −1 , highlighting their significant energetic output.

CONCLUSION
High redshift blazars have consistently been the subject of intense investigation due to their crucial role in understanding the connection between jet formation, accretion processes, and black hole development in the early stages of black hole formation.In this paper, we investigate the multiwavelength properties of 79 blazars with redshifts ranging from  = 2.0 to 2.5 using data accumulated by Fermi, Swift XRT/UVOT, and NuSTAR.
In the -ray band, the flux and photon index of the selected sources span from 5.32×10 −10 to 3.40×10 −7 photon cm −2 s −1 and between 1.66 and 3.15, respectively, illustrating the diverse characteristics of the sources under consideration.Meanwhile, the luminosity ranges from (3.67±1.37)×10 46 to (6.62±0.05)×10 48 erg s −1 , positioning some of these sources among the brightest blazars detected in the -ray band.Similarly, in the 0.3-10 keV band, the weakest sources exhibit fluxes on the order of (1.06 ± 0.32) × 10 −13 erg cm −2 s −1 , while the brighter ones reach (2.96±0.02)×10−11 erg cm −2 s −1 , with the majority of the sources exhibiting a soft photon index (less than 2.0).This soft photon index extends into the hard X-ray band (3 − 30 keV), where the index ranges from 1.09 to 1.67.The comparison of -ray and X-ray fluxes reveals no clear correlation in time-averaged measurements.This is exemplified by sources such as 4C +71.07 and PKS 2149-306, which are bright in both -ray and X-ray bands, and 4C +01.02, which, despite its high -ray flux, exhibits only a moderate X-ray flux.
In this study, we employed a leptonic one-zone synchrotron and inverse Compton model, considering both internal and external photons, to interpret the multiwavelength SEDs of selected FSRQs and BL Lacs.The range of parameters estimated from the modeling provides a general view of the sources' emissions in an average state.The power-law index of the electrons, indicative of particle acceleration mechanisms, varied from 1.56 to 2.85.For FSRQs the highest energy emissions, as characterized by the cutoff Lorentz factor, varied broadly from 237.6 to 3.53 × 10 3 , reflecting the diverse characteristics of the sources studied.The total jet luminosity varied within a relatively narrow range from 3.20 × 10 44 to 6.51 × 10 45 erg s −1 .The disk luminosity, estimated for sources with discernible thermal components in their SEDs, provided insights into the accretion processes.These disk luminosities, ranging dramatically from 4.15 × 10 44 to 3.97 × 10 47 erg s −1 , highlighted the varying scales of accretion efficiency and energy conversion in these sources.
In summary, this comprehensive study of 79 high-redshift blazars illuminates the intricate dynamics and energetic processes at play in these distant and powerful objects.Our findings, encompassing a broad spectrum of multiwavelength observations, highlight the diverse characteristics and behaviors of blazars, underscoring their value in probing the physics of jet formation and accretion processes, and black hole development in the early Universe.The observed variations in flux and luminosity across different bands provide key insights into the particle acceleration and emission mechanisms.Future investigations, when a larger number of distant blazars are known, will allow a more systematic and statistical comparison of their emission properties, thereby unveiling the complexities of these energetic sources and providing a window into the HE universe and its evolution over cosmic time.

Figure 1 .
Figure 1.Hammer-Aitoff projection in Galactic coordinates showing the distribution of -ray blazars with redshifts above  > 2.0.BL Lacs within the redshift range 2.0 ≤  ≤ 2.5 are represented by orange symbols, FSRQs are depicted in blue, BCUs in green, and blazars with redshifts  > 2.5 are shown as grey diamonds.
Photon index

Figure 2 .
Figure 2. Panels (a) and (b): The -ray flux (> 100 MeV) and luminosity versus the photon index.Panel (c) The X-ray flux versus the photon index.Panel (d) The X-ray flux versus the -ray flux.

Figure 3 .Figure 4 .
Figure3.The -ray light curve of the sources considered in this study, which show a high amplitude -ray flux increase during the flares.The left axis shows the variation of flux over time, while the right axis displays the luminosity.MNRAS 000,1-20 (2022)

Figure 6 .
Figure 6.The broadband SED modeling results for the sources under consideration.The data analyzed from Swift UVOT, XRT, NuSTAR, and Fermi-LAT observation are depicted in cyan, blue, orange, and red, respectively, while archival data are shown in gray.The blue line represents the combined contribution of the synchrotron, disk, SSC, and EIC-BLR components.The blue shaded area indicates the region of uncertainty associated with the model.MNRAS 000, 1-20 (2022) Figure 6.(continued)

Figure 7 .
Figure 7. Upper panel: The dependence of the total jet luminosity ( jet ) on redshift.Middle panel: The dependence of the disc luminosity ( D ) on redshift.Lower panel: The correlation between jet luminosity and disc luminosity.The sources from Sahakyan et al. (2020) are shown in gray.

Table 1 .
The table shows the source sample, detailing the outcomes of the -ray analysis.For each source, the name, associated 4FGL name, and class are provided. p represents the photon index when the source spectrum is best modeled with a power-law model, whereas  and  denote the slope and curvature, respectively, when the spectrum is modeled with a log-parabola.The flux is reported in units of 10 −8 photons cm −2 s −1 , and the luminosity is expressed in 10 47 ergs s −1 .The redshift  is given in last column.The result for this object is listed from the catalogue, as its ROI contains 10 extended sources, which complicates the analysis.

Table 2 .
NuSTAR analysis results.The source name, NuSTAR sequence ID, observation time in MJD, the logarithm of the 3-10 keV (log  3−10 ) and 10-30 keV (log  10−30 ) bands flux and the photon index for each observation are given.

Table 3 .
Broadband SED Modeling of the Considered Sources.[1] Source Name.[2] Doppler Factor.[3] Power-Law index of emitting electrons.[4] Minimum Lorentz factor of the electrons and [5] cutoff Lorentz factor in units of 100.[6] Magnetic field, expressed in units of .[7] Accretion Disk Luminosity, in units of 10 46 erg s −1 .[8] and [9] Power of the jet in the form of relativistic electrons (  ) and Magnetic field (  ), both in units of 10 44 erg s −1 .[10] Luminosity of the jet in the form of protons assuming one proton per relativistic emitting electron in units of 10 45 erg s −1 .