Broad-band study of gamma-ray blazars at redshifts z = 2.0–2.5 Open Access

ABSTRACT

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 γ-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 γ-ray band, the spectral analysis revealed a wide range of flux and photon indices, from 5.32 × 10⁻¹⁰ to 3.40 × 10⁻⁷ photon cm⁻² s⁻¹ 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 γ-ray luminosity and flux variations, reaching peak luminosity exceeding 10⁵⁰ erg s⁻¹. The temporal analysis extended to X-ray and optical/ultraviolet (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 γ-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 disc luminosity derived from fitting the excess in the UV band.

1 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-time-scale 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 non-thermal 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 very high energy (VHE; >100 MeV) γ-ray bands. Their broad-band spectral energy distribution (SED) 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¹⁴ Hz, 10¹⁴ Hz < ν_p < 10¹⁵ Hz, and ν_p > 10¹⁵ 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, Maraschi & Treves 1985; Maraschi, Ghisellini & Celotti 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, Begelman & Rees 1994) suggest that the photons originate outside the jet, coming directly from the accretion disc (Dermer, Schlickeiser & Mastichiadis 1992; Dermer & Schlickeiser 1994), or they may be reprocessed by the broad-line region (BLR; Sikora, Begelman & Rees 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 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 & Biermann 1989; Mannheim 1993; Mücke & Protheroe 2001; Mücke et al. 2003; Böttcher et al. 2013; Petropoulou & Mastichiadis 2015; Gasparyan, Bégué & Sahakyan 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 (Padovani et al. 2018; IceCube Collaboration 2018a, b), and the detection of multiple neutrino events concurrent with the active phase of PKS 0735 + 178 in optical/ultraviolet (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, Oikonomou & Petropoulou 2018; Padovani et al. 2018; Sahakyan 2018; Cerruti et al. 2019; Gao et al. 2019; Righi, Tavecchio & Pacciani 2019; Sahakyan 2019; Gasparyan, Bégué & Sahakyan 2022).

The emission from blazars is highly beamed, and their bolometric luminosity can exceed 10⁴⁸ erg s⁻¹, 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 discs 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 Pacciani et al. 2012; Paliya 2015; Paliya et al. 2015, 2016, 2017, 2019; D’Ammando & Orienti 2016; Orienti et al. 2016; Marcotulli et al. 2017; Liao et al. 2019; Li et al. 2020; Sahakyan et al. 2020; Sahakyan, Harutyunyan & Israyelyan 2023b). In Sahakyan et al. (2020), the origin of emission from the most distant blazars detected in the HE γ-ray band (z > 2.5) was investigated analysing 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 broad-band 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 z = 2.0–2.5. A prolonged observation period of ∼14 yr 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 analysed data broad-band 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 broad-band SEDs; and Section 6 outlines the modeling results. Finally, the conclusions are summarized in Section 7.

2 SOURCE SAMPLE

The fourth catalogue of AGNs detected by the Fermi Large Area Telescope (Fermi-LAT; Ajello et al. 2022) contains 3814 blazars, among which 792 are FSRQs, 1458 are BL Lacs, and 1493 are blazar candidates of uncertain type (BCUs). The most distant blazar, GB 1508 + 5714, is at z = 4.31. A small fraction of the blazars, 79 (2.07  per cent), have a redshift between 2.0 < z < 2.5, and these are the ones selected for the current study (the source sample). This group includes 64 FSRQs, 9 BL Lacs, and 6 BCUs. The BL Lacs generally have redshifts lower than z = 2.1, with the exception of SDSS J145059.99 + 520111.7, which is at z = 2.47. Among the BL Lacs, there are three LBLs, four IBLs, and only two HBLs. The FSRQs have a more homogeneous redshift distribution, being observed at almost all redshifts; the most distant FSRQ in the sample is S5 1053 + 70 at z = 2.49. BCUs, which exhibit characteristics similar to blazars but lack reliable optical associations, have been observed across a range of redshifts. For instance, 4FGL J1139.0+4033 (CRATES J113903 + 403303) has been observed at z = 2.36, while 4FGL J1003.4+0205 (SDSS J100326.63 + 020455.6) is located at z = 2.08. In Sahakyan, Vardanyan & Khachatryan (2022a) BCUs were classified by training machine-learning algorithms on the γ-ray properties of FSRQs and BL Lacs. According to those criteria, four BCUs from our sample show probability similar to FSRQs (CRATES J113903+403303, MG4 J162750+4802, TXS 2315+189, and SDSS J120542.82+332146.9), one to BL Lacs (SDSS J100326.63+020455.6) and one (SDSS J120542.82 + 332146.9) can not be classified.

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 z > 2.5, as reported by Sahakyan et al. (2020).

3 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 analysed.

3.1 Fermi-LAT data

The 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 2008 August 4 to 2022 December 4 (∼14.5 yr), were considered to study the properties of all 79 blazars under consideration. The standard data-reduction procedure was performed following the recommendations from the Fermi-LAT science team.¹ 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-centred on the γ-ray position of the sources, were downloaded and analysed. 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 |${\rm (DATA\_QUAL \gt 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 TS = 2 × (log L − log L₀), where L is the likelihood with the source at the position of interest included, and L₀ 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-yr period was divided into equal intervals (e.g. 5, 7, 10 d, depending on the source’s overall detection significance) 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-time-scale flux variations in blazar emissions (e.g. see Rani et al. 2013; Britto et al. 2016; Baghmanyan, Gasparyan & Sahakyan 2017; Sahakyan & Gasparyan 2017; Zargaryan et al. 2017; Gasparyan et al. 2018; Sahakyan, Baghmanyan & Zargaryan 2018; Sahakyan 2021; Sahakyan & Giommi 2022; Sahakyan et al. 2022b).

3.2 NuSTAR

NuSTAR (Harrison et al. 2013) is a hard X-ray telescope operating in the 3–79 keV range, equipped with two focal plane modules: FPMA and FPMB. Among the sources considered, PKS 0528+134, S3 0458–02, PKS 0446+11, PKS 1329–049, 87GB 080551.6+535010, and TXS 0322+222 were observed by NuSTAR once; PKS 2149–306, S5 0212+73, and PKS 0227–369 were observed twice; and 4C + 71.07 was observed three times. In total, these amount to 15 observations that provide critical information on the hard X-ray band emission of the sources.

The analysis of all NuSTAR data was performed using the NuSTAR_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).

3.3 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 analysed 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 centred 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 analysed using the tool provided by the UK Swift Science Data Centre (Evans et al. 2009).

3.4 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 analysed.

The data were analysed 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 arcsec centred on the source, while background counts were obtained from a larger circular region with a radius of 20 arcsec, 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 E(B − V) obtained from the Infrared Science Archive.² The corrected 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 Ly α absorption systems. These effects were corrected for in the SEDs during theoretical modeling, following the procedures described in Ghisellini et al. (2011, 2010).

3.5 Archival data

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

4 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 low-parabola model in 4FGL an additional analysis is performed using power-law model.

The analysis results are presented in Fig. 2(a), 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⁻¹⁰ to (3.40 ± 0.02) × 10⁻⁷ photon cm⁻² s⁻¹, 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⁻⁸ photon cm⁻² s⁻¹. 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⁻⁸ photon cm⁻² s⁻¹ for flux and 1.91–2.62 for photon index, while BCUs between (0.05–1.60) × 10⁻⁸ photon cm⁻² s⁻¹ 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 grey in Fig. 2(a). 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 F_γ plane.

Fig. 2(b) presents the γ-ray luminosity L_γ versus the photon index. This contrasts with the flux, luminosity, on the other hand, accounts for the total energy emitted by the source per unit time, so showing the intrinsic power of the sources. The luminosity of the sources under consideration spans from (3.67 ± 1.37) × 10⁴⁶ erg s⁻¹ to (6.62 ± 0.05) × 10⁴⁸ erg s⁻¹, with the lowest estimated for SDSS J105707.47+551032.2 and the highest for 4C +01.02. Notably, the luminosity of 4C +01.02 also exceeds that of B3 1343 + 451, which is the most luminous in the sample of sources with a redshift beyond 2.5. The range of luminosities for the new sample is slightly shifted towards a higher luminosity range. For instance, the luminosities of PKS 0226–559 (2.21 × 10⁴⁸ erg s⁻¹), PKS 0601–70 (1.12 × 10⁴⁸ erg s⁻¹), PKS 2149–306 (1.09 × 10⁴⁸ erg s⁻¹), S3 0458–02 (2.76 × 10⁴⁸ erg s⁻¹), 4C + 71.07 (1.48 × 10⁴⁸ erg s⁻¹), S4 0917 + 44 (1.80 × 10⁴⁸ erg s⁻¹), PKS 0446 + 11 (1.21 × 10⁴⁸ erg s⁻¹), PKS 1329–049 (1.16 × 10⁴⁸ erg s⁻¹), and 4C + 01.02 (6.62 × 10⁴⁸ erg s⁻¹) are exceeding luminosity of 10⁴⁸ erg s⁻¹, making them among the luminous blazars detected in the γ-ray band.

Fig. 2(c) shows the relationship between the X-ray photon index and the flux of the considered sources. The X-ray flux for the sources studied range from (1.06 ± 0.32) × 10⁻¹³ erg cm⁻² s⁻¹ for PMN J1344-1723 to (2.96 ± 0.02) × 10⁻¹¹ erg cm⁻² s⁻¹ for 4C +71.07. The X-ray photon index is predominantly soft (less than 2.0) for the majority of the sources, suggesting that the X-ray emissions are likely dominated by the rising part of inverse Compton component. Notably, the brightest sources in the sample, such as 4C + 71.07 and PKS 2149–306, have fluxes of (2.96 ± 0.02) × 10⁻¹¹ erg cm⁻² s⁻¹ and (2.01 ± 0.04) × 10⁻¹¹ erg cm⁻² s⁻¹, with photon indices of 1.28 ± 0.01 and 1.19 ± 0.03, respectively. These indices indicate their X-ray spectra are particularly hard compared to the rest of the sample. Correspondingly, these two sources also show the highest X-ray luminosities, being (1.13 ± 0.11) × 10⁴⁸ erg s⁻¹ for 4C + 71.07, and (8.79 ± 1.83) × 10⁴⁷ erg s⁻¹ for PKS 2149–306.

In Fig. 2(d), 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⁻⁷ photon cm⁻² s⁻¹ and (7.30 ± 0.20) × 10⁻⁸ photon cm⁻² s⁻¹, respectively. Conversely, the γ-ray bright source 4C + 01.02, has only a moderate X-ray flux of (2.18 ± 0.11) × 10⁻¹² erg cm⁻² s⁻¹, 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, 60 002 045 002 for 4C + 71.07, and 60 367 002 002 for PKS 0227–369) were relatively short (on the order of a few hundred seconds), and hence no spectral analysis was conducted. In the 3–10 keV range, the highest flux of (2.02 ± 0.01) × 10⁻¹¹ erg cm⁻² s⁻¹ was observed for 4C + 71.07 on MJD 56675.22, while the lowest flux of (1.95 ± 0.27) × 10⁻¹³ erg cm⁻² s⁻¹ 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.

4.1 γ-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 z = 2.464 showing variability, computed for energies above 239.66 MeV, is displayed in Fig. 3(a). From the start of Fermi-LAT observations until MJD 55 504 (2010 November 4), the γ-ray emission of this source was in a low state being (4.89 ± 0.71) × 10⁻⁹ photon cm⁻² s⁻¹. Subsequent periods of flaring activity occurred between MJD 55504.30–55584.50, 56153.18–57472.15, and MJD 58084.93–58664.52. During its most active state, between MJD 58084.93 and 58664.52, the maximum flux was (3.99 ± 0.67) × 10⁻⁷ photon cm⁻² s⁻¹ 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⁴⁸ erg s⁻¹, the luminosity during flares exceeded 10⁴⁹ erg s⁻¹. For example, the source luminosity was >10⁴⁹ erg s⁻¹ 42 times, with the peak luminosity of (5.33 ± 0.92) × 10⁴⁹ erg s⁻¹ observed on MJD 58167.46.

Among the flaring sources, significant increases in γ-ray emission have been observed for PKS 2149–306, 4C+71.07, PKS 1329–049, and 4C + 01.02. For example, the time-averaged γ-ray flux of PKS 2149–306 (z = 2.345) is (2.01 ± 0.47) × 10⁻⁷ photon cm⁻² s⁻¹, but there was a notable period of elevated γ-ray emission between MJD 56288.92 and 56444.75 (Fig. 3b). During this interval, the flux above 144.85 MeV exceeded 10⁻⁶ photon cm⁻² s⁻¹ on eight time intervals. The peak flux during this period was (1.68 ± 0.33) × 10⁻⁶ photon cm⁻² s⁻¹ observed on MJD 56302.61. Additionally, the source exhibited another active γ-ray emission state between MJD 55554.39 and 55650.67. In this period, the γ-ray flux consistently exceeded 10⁻⁷ photon cm⁻² s⁻¹, with a maximum of (7.08 ± 1.22) × 10⁻⁷ photon cm⁻² s⁻¹ observed on MJD 55613.98. During multiple flaring periods of 4C + 71.07 (Fig. 3d), which is at a redshift of z = 2.218, the γ-ray emission notably exceeded its time-averaged flux of (1.06 ± 0.01) × 10⁻⁷ photon cm⁻² s⁻¹. Specifically, in three distinct flaring periods – MJD 55862.01–55932.97, MJD 57227–57245.79, and MJD 57296.92–57338.47 – the flux above 138.16 MeV exceeded 10⁻⁶ photon cm⁻² s⁻¹. The peak flux during these flaring events was at (8.02 ± 1.40) × 10⁻⁶ photon cm⁻² s⁻¹, and was observed on MJD 57335.21. The γ-ray emission of PKS 1329–049, with a redshift of z = 2.15, is predominantly in its average emission state, as depicted in Fig. 3(f). However, it exhibited an elevated emission state during the period from MJD 55442.09 to 55507.87 when the peak flux reached (2.71 ± 0.52) × 10⁻⁶ photon cm⁻² s⁻¹ observed on MJD 55445.32. In contrast, 4C + 01.02, with a redshift of z = 2.099, experiences alternating periods of flaring activity, with the most intense γ-ray flares observed after MJD 57000. The source entered an active γ-ray emission state starting from MJD 59615.91, when the peak γ-ray flux, measured above 171.79 MeV, reached (2.92 ± 0.50) × 10⁻⁶ photon cm⁻² s⁻¹ on MJD 59663.19. The other two sources, depicted in Fig. 3(c) and (e), also exhibit multiple flaring periods but with only modest increases in their γ-ray flux. Specifically, the peak γ-ray flux of S3 0458–02, measured above 202.55 MeV, was (8.32 ± 1.41) × 10⁻⁷ photon cm⁻² s⁻¹ observed on MJD 57258.92. Similarly, for S4 0917 + 44, it was (4.27 ± 0.75) × 10⁻⁷ photon cm⁻² s⁻¹ observed on MJD 57946.70.

Especially profound is the luminosity increase in 4C+71.07, PKS 1329–049, and 4C + 01.02 (see Fig. 3). For 4C + 71.07, during the flaring periods in MJD 57227–57245.79 and MJD 57296.92–57338.47, the source luminosity exceeded 10⁵⁰ erg s⁻¹ 11 times. The highest luminosity of (2.03 ± 0.36) × 10⁵⁰ erg s⁻¹ was observed on MJD 57335.21. PKS 1329–049 was in an extreme bright state on MJD 55445.23 and MJD 55468.18, with luminosities of (1.07 ± 0.21) × 10⁵⁰ erg s⁻¹ and (1.39 ± 0.26) × 10⁵⁰ erg s⁻¹, respectively. Similarly, 4C + 01.02 was in an extreme bright state between MJD 59662.03–59663.55, during which in 6 consecutive bins the flux exceeded 10⁵⁰ erg s⁻¹, with the highest flux of (2.08 ± 0.48) × 10⁵⁰ erg s⁻¹ observed on MJD 59662.62. Because of such elevated luminosity, 4C+71.07, PKS 1329–049, and 4C + 01.02 rank among the sources with the highest luminosity in the γ-ray band.

In Fig. 4, the γ-ray light curves for 4C+13.14, 87GB 080551.6+535010, B2 1436+37B, PKS 0227–369, PKS 0528+134, PKS 0446+11, PKS 0601–70, S5 1053+70, and TXS 0322 + 222 are presented. For the majority of the time, the emission from these objects remains in a low to average state. However, during certain flaring periods, their γ-ray emission shows modest increases. These sources are typically weak, with γ-ray fluxes generally on the order of 10⁻⁸ photon cm⁻² s⁻¹, but there are instances when the flux exceeds 10⁻⁷ photon cm⁻² s⁻¹. For example, during an elevated γ-ray emission state, the peak flux of S5 1053 + 70 above 179.20 MeV was (1.42 ± 0.24) × 10⁻⁷ photon cm⁻² s⁻¹ observed on MJD 57684.01. Similarly, for PKS 0601–70, the highest γ-ray flux reached (1.64 ± 0.30) × 10⁻⁷ photon cm⁻² s⁻¹ on MJD 55136.37. For PKS 0446 + 11, the peak γ-ray flux was (1.77 ± 0.34) × 10⁻⁷ photon cm⁻² s⁻¹ observed on MJD 55304.46, and so on.

4.2 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.

The X-ray flux variations for 4C + 71.07 (blue) and PKS 2149–306 (green) are shown in the upper panel of Fig. 5. For 4C + 71.07, the lowest observed X-ray flux was (1.00 ± 0.12) × 10⁻¹¹ erg cm⁻² s⁻¹ on MJD 55247.68, while the highest was (5.63 ± 1.02) × 10⁻¹¹ erg cm⁻² s⁻¹ on MJD 56832.97. Elevated X-ray emission states for this source were observed around MJD 56 000 and MJD 57000, during which most of the X-ray observations were performed. Conversely, PKS 2149–306 exhibited lower amplitude flux changes, with its highest X-ray flux being (2.92 ± 0.58) × 10⁻¹¹ erg cm⁻² s⁻¹, observed on MJD 56643.02. The lower panel of Fig. 5 shows the X-ray flux variations for PKS 1329–048 (red) and S3 0458–02 (blue). Initially, the X-ray emission for S3 0458–02 was on the order of ∼10⁻¹² erg cm⁻² s⁻¹, but it increased to approximately ∼5 × 10⁻¹² erg cm⁻² s⁻¹ around MJD 56200. For PKS 1329–048, the observed highest X-ray flux was (4.28 ± 0.68) × 10⁻¹² erg cm⁻² s⁻¹ on MJD 55390.19, whereas it decreased to (8.99 ± 4.03) × 10⁻¹³ erg cm⁻² s⁻¹ on MJD 57902.45.

4.3 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 V, B, and U filters are approximately |$10^{-12} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$|⁠, and approximately |$2 \times 10^{-13} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$| in the W1, W2, and M2 filters. After MJD 56500, the source exhibits increased flux across all filters during several observations. The highest observed flux in the B filter was |$(3.91 \pm 0.33) \times 10^{-12} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$| on MJD 57363.45. For PKS 0226–559, the mean flux is |$\sim 5 \times 10^{-13} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$|⁠, but the flux in V filter increased up to |$(3.1-0.46) \times 10^{-12} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$| during the flare on MJD 58161.21. The UV emission from PKS 2149–306, in the W1, W2, and M2 filters, remains relatively stable across various observations, whereas the optical emission (V, B, and U) shows variability. Specifically, on MJD 53717.92, MJD 54980.86, and MJD 58586.29, the flux increased to approximately |$3 \times 10^{-12} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$|⁠, with the lowest observed flux being around |$10^{-13} \, \text{erg} \, \text{cm}^{-2} \, \text{s}^{-1}$|⁠.

5 MODELING OF MULTIWAVELENGTH SEDS

The data analysed 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 analysed in this work are highlighted in different colors (see the legend), while archival data retrieved using the VOU-Blazar tool are depicted in grey. All the 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 R_dis from the central black hole and is assumed to have a spherical geometry with a radius R. 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 non-thermal particles (electrons and positrons), whose energy distribution follows a power law with an exponential cutoff described as

where p 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, N₀ denotes the normalization constant of the electron distribution, which defines the total energy content of the electrons as given by U_e = mc²∫γN(γ)dγ.

Within the emitting region (blob), electrons lose energy through synchrotron emission under the magnetic field B, 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, Maraschi & Treves 1985; Maraschi, Ghisellini & Celotti 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 disc (Dermer, Schlickeiser & Mastichiadis 1992; Dermer & Schlickeiser 1994), photons reprocessed by the BLR (Sikora, Begelman & Rees 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 (R_dis < R), and that accretion disc 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 disc luminosity (L_disc) 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 disc’s emission does not outshine the observed non-thermal emission from the jet. With L_disc, the radius of the BLR is calculated using the relation R_BLR = 10¹⁷L_{45, disc} (Ghisellini & Tavecchio 2015). Subsequently, the BLR is represented as a spherical shell, with an inner boundary of R_{in, BLR} = 0.9 × R_BLR and an outer boundary of R_{out, BLR} = 1.2 × R_BLR.

The SEDs were fitted using the publicly available jetset code (Massaro et al. 2006; Tramacere et al. 2009; Tramacere, Massaro & Taylor 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, disc, SSC, and EIC-BLR components, which are optimized during the fitting process. The free parameters in the model are the spectral index (p), the cutoff Lorentz factor (γ_cut), the minimum Lorentz factor (γ_min), the magnetic field (B), the Doppler factor (δ), the size of the emitting region (R), and the energy density of the electrons (U_e). Optimization is performed in two stages: initial fitting is done with the Minuit optimizer, followed by refinement using the Markov Chain Monte Carlo method. We note that R 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, R was treated as a parameter dependent on the variability time-scale (t_var) and δ, via the relation R = δct_var/(1 + z), with t_var assumed to be within a 0.1–7-d range. However, we refrain from commenting on t_var (and consequently R) 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, Rodighiero & Vaccari (2008).

6 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, disc, 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 disc components, while at HEs, the SSC and EIC-BLR components are dominant. The modeling facilitates the estimation of parameters describing the jet (B, δ) as well as those characterizing the emitting electrons (p, γ_min, and γ_cut). From the modeling of the SEDs of the FSRQs, the power-law index of the electrons varies within a range from p = 1.56 to p = 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, p = 1.56 ± 0.064 for PKS 0226–559, p = 1.60 ± 0.12 for OX 131, and p = 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.01 and 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³, determined by the optical/UV and/or X-ray data. A high value of (3.55 ± 0.49) × 10⁴ 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³, was estimated for PMN J1344-1723. In the modeling of BL Lacs, the cutoff Lorentz factor (γ_cut) is (1.53–25.41) × 10⁴ 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⁵, 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.71 for 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–0.04 G. For the FSRQs, except for the bordering cases of PKS B1412-096, PMN J2135-5006, and PMN J1344-1723, where δ of 5.24, 5.45, and 47.18 were estimated respectively, δ varies within a narrow range of 12.52–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–49.84 were estimated. This is because the synchrotron component of these objects peaks around 10¹³–10¹⁴ 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 |$L_{e}=\pi c R_b^2 \Gamma ^2 U_{e}$| and |$L_{B}=\pi c R_b^2 \Gamma ^2 U_{B}$|⁠, respectively. Assuming one proton per relativistic emitting electron and that protons are ‘cold’ in the comoving frame, the total jet luminosity is calculated as L_jet = L_e + L_B + L_p. The corresponding L_e and L_B luminosities are presented in Table 3. For the FSRQs, L_e varies within the range from 7.24 × 10⁴³ erg s⁻¹ to 2.94 × 10⁴⁵ erg s⁻¹, while L_B ranges from 2.72 × 10⁴² erg s⁻¹ to 4.06 × 10⁴⁵ erg s⁻¹. Among the considered sources, 33 exhibit an L_e/L_B ratio in the range of 0.1 to 10, suggesting that the emitting region is close to equipartition. For PKS B1412-096 the ratio L_e/L_B = 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 L_e/L_B > 100 is estimated, due to the dominance of the SSC component over the synchrotron component. In the modeling of BL Lacs, L_e is in the range of (0.04–3.85) × 10⁴⁷ erg s⁻¹, which is higher than the L_B, ranging from (0.10–3.92) × 10⁴³ erg s⁻¹. The relationship between the total jet luminosity (L_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; grey), 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 z > 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 z > 2.5. For the FSRQs the highest luminosity of 6.51 × 10⁴⁵ erg s⁻¹ was observed for 4C + 71.07 at a redshift of z = 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 z = 3.72 and z = 4.31, respectively.

The fit presented in Fig. 6 also enables the estimation of disc luminosity for cases where a thermal blue-bump component is observable. This excess is evident in a total of 37 sources. The disc luminosity, L_D, varies from 5.60 × 10⁴⁴ erg s⁻¹ to 3.97 × 10⁴⁷ erg s⁻¹. The relationship between disc luminosity and redshift for γ-ray blazars with z > 2.0 is depicted in the middle panel of Fig. 7. High disc luminosities of 1.25 × 10⁴⁷ erg s⁻¹, 2.32 × 10⁴⁷ erg s⁻¹, 3.35 × 10⁴⁷ erg s⁻¹, and 3.97 × 10⁴⁷ erg s⁻¹ were estimated for PKS 2149–306, 4C +71.07, 2MASS J16561677-3302127, and B2 0552 + 39A, respectively. Furthermore, the relationship between disc luminosity and jet luminosity is examined in the lower panel of Fig. 7. This relationship can be described by the equation log₁₀(L_jet) = 0.22 × log₁₀(L_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 disc luminosities exceeding 10⁴⁷ erg s⁻¹ and jet luminosities exceeding 10⁴⁵ erg s⁻¹, highlighting their significant energetic output.

7 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 z = 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⁻¹⁰ to 3.40 × 10⁻⁷ photon cm⁻² s⁻¹ 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⁴⁶ to (6.62 ± 0.05) × 10⁴⁸ erg s⁻¹, 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⁻¹³ erg cm⁻²s⁻¹, while the brighter ones reach (2.96 ± 0.02) × 10⁻¹¹ erg cm⁻²s⁻¹, 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.

Using the adaptive binning methods for light curve, we found flux variations in 31 sources, with the most pronounced variability in 4C+01.02, 4C+71.07, PKS 0226–559, PKS 1329–049, PKS 2149–306, S3 0458–02, and S4 0917+44. The sources like PKS 2149–306, 4C+71.07, PKS 1329–049 and 4C + 01.02 showed significant increases in γ-ray emission, with peak fluxes occasionally exceeding 10⁻⁶ photon cm⁻² s⁻¹ while for the other sources show more modest increases in γ-ray flux: the γ-ray fluxes exceed 10⁻⁷ photon cm⁻² s⁻¹ during flaring periods. Noticeably, the γ-ray luminosity of 4C+71.07, PKS 1329–049, and 4C + 01.02 occasionally was above 10⁵⁰ erg s⁻¹, putting them among the most luminous sources in the γ-ray band.

In addition to γ-ray variability, our analysis extended to X-ray and optical/UV bands. Despite the limited data, X-ray flux variability was evident for 4C +71.07, PKS 0226–559, PKS 1329–049, PKS 2149–306, S3 0458–02 and S4 0917+44, while in optical/UV band clear flux variation is found for 4C + 01.02, PKS 0226–559, and PKS 2149–306.

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³, reflecting the diverse characteristics of the sources studied. The total jet luminosity varied within a relatively narrow range from 3.20 × 10⁴⁴ to 6.51 × 10⁴⁵ erg s⁻¹. The disc luminosity, estimated for sources with discernible thermal components in their SEDs, provided insights into the accretion processes. These disc luminosities, ranging dramatically from 4.15 × 10⁴⁴ to 3.97 × 10⁴⁷ erg s⁻¹, 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.

ACKNOWLEDGEMENTS

We acknowledge the use of analysis tools and services from the Markarian Multiwavelength Data Center (https://mmdc.am), the Astrophysics Science Archive Research Center (HEASARC) platforms, as well as data from the Fermi-LAT, Swift, and NuSTAR telescopes. This work was support by the Higher Education and Science Committee of the Republic of Armenia, in the frames of the research project No. 21T-1C260.

DATA AVAILABILITY

The data underlying this article will be shared on reasonable request to the corresponding author.

Footnotes

References