Central Engine and Spectral Energy Distribution Properties of High Redshift Gamma Ray Blazars

We report on the properties of central engines in the $\gamma$-ray blazars located at high redshifts beyond z~>~0.4, where the extra-galactic background light (EBL) starts affecting their $\gamma$-ray spectra. The physical engine that provides power to the blazars of very high bolometric luminosity is assumed to be a highly collimated jet of matter moving relativistically away from the supermassive black hole (SMBH), located in the central region of the host galaxy, in a direction aligned toward the Earth. Due to their peculiar geometry and special physical conditions, blazars at redshifts beyond z~>~0.4 are bright enough to be detected in the $\gamma$-ray energy band. In this work, we investigate the physical properties of high-$z$ $\gamma$-ray blazars detected by the Large Area Telescope (LAT) on board the \emph{Fermi} satellite. We also study the properties of their emission regions and the central engines and discuss cosmological and astrophysical implications.


Introduction
Blazars are bright beacons in the Universe.They represent a subclass of the radioloud active galactic nuclei (AGNs) which are characterized as elliptical galaxies hosting supermassive black holes (SMBHs;M BH : 10 6 − 10 10 M ⊙ ) at their centers and oppositely oriented powerful jets of plasma ejected from the poles of SMBH [1,2].The accretion of gaseous matter onto the spinning SMBH is considered to be the dominant physical process for powering these sources and launching the relativistic plasma jets.However, the exact mechanism of jet formation is not completely understood and three alternative theories have been proposed.In the so-called Blandford-Znajek (BZ) mechanism, the jet extracts energy from the rotation of SMBH, and the power of the jet depends on the black hole spin and square of magnetic flux threading its event horizon [3].In the second model, also known as the Blandford-Payne (BP) mechanism, the jet extracts energy only from the rotation of the accretion disk, and the role of SMBH is not important [4].In the case of these two theories, the formation of jets is sustained by the accretion of matter onto the SMBHs, and a relation between jet and accretion luminosity is expected and has been confirmed in different studies [5][6][7][8].The hybrid model, a combination of the above two mechanisms, speculates about the observed differences in different types of AGNs with jets [9].The jets in AGNs are observed to be ultra-luminous sources of non-thermal radiation spanning over the entire electromagnetic spectrum ranging from radio waves to very high energy (VHE; E > 20 GeV) γ-rays.They are also thought to be the apparent manifestations of AGN activity during their early evolutionary stages [10].An extreme class of such AGNs, in which one of the jets points close to the line of sight of the observer, is referred to as blazars.Due to the peculiar orientation of the blazar jets, their emission undergoes relativistic Doppler boosting and beaming [1,2].As a consequence, the intrinsic emission from blazars is strongly enhanced and the characteristic timescales are shortened for an observer on Earth.As a result, blazars are observed to exhibit distinct observational features such as rapid and large amplitude flux variability on most timescales across the electromagnetic spectrum, spectral changes (harder-when-brighter), strong and variable optical polarization, variable high energy (HE; E > 30 MeV) γ-ray emissions, and apparent superluminal motion [11][12][13][14][15][16][17][18].The intensity amplification effects make even faint blazars detectable in the multi-wavelength observations.Thus, blazars are the ideal laboratories to probe the physics properties of AGNs and study the growth and evolution of SMBHs using the measurements of their broadband emission [19].
The broadband spectral energy distribution (SED) of blazars typically exhibits two characteristic broad humps.The first hump, peaking at low energies in radio to UV/soft X-ray, is produced by the synchrotron emission from the population of relativistic electrons in the presence of jet magnetic field.The origin of the second hump, which peaks in the γ-ray band is not well understood.It is often attributed to inverse Compton (IC) up-scattering of low energy seed photons by the relativistic jet electrons in the emission region.The seed photons for the IC process may be the synchrotron photons produced within the jet (synchrotron self Compton or SSC [20]) or thermal photons (UV/IR) originating from various components such as accretion disk, broad line region (BLR) and dusty torus (DT) external to the jet.These processes are referred to as external Compton (EC) for HE and VHE emissions in blazars [21,22].Alternatively, hadronic models involving the ultra-relativistic protons and higher jet magnetic field are invoked to explain the observed VHE emission through proton synchrotron, photo-hadronic interactions, and subsequent meson decay [23][24][25].
Blazars are historically classified into two categories: Flat Spectrum Radio Quasars (FSRQs) and BL Lacertae Objects (BL Lacs), on the basis of strength, visibility, and rest-frame equivalent width (EW) of the emission/absorption lines in their optical continuum [26].FSRQs are characterized by strong and broad emission lines with EW > 5 Å while BL Lacs have typically weak or no emission lines.Another classification of blazars is also proposed on the basis of the known features of double-hump structures in their broadband SEDs.The position of the low energy hump or synchrotron peak distinguishes the blazars as Low Synchrotron Peak (LSP), Intermediate Synchrotron Peak (ISP), and High Synchrotron Peak (HSP) objects [27][28][29].For LSP blazars, the synchrotron peak frequency of ≤ 10 14 Hz, corresponds to the radio to infrared regime.ISP blazars have synchrotron peak frequency in the range 10 14 -10 15 Hz, which lies in the optical/UV regime.Finally, HSP blazars exhibit the synchrotron peak frequency above 10 15 Hz, corresponding to the X-ray regime.Further, the high-energy hump or Compton peak of BL Lacs is equally luminous as the low-energy synchrotron peak and is attributed to the SSC process.Whereas in the case of FSRQs, the Compton peak is more luminous than the synchrotron peak and is generally ascribed by the EC process.Also, the accretion disc in FSRQs is more radiatively efficient than that in BL Lacs [2,6].This provides a more physical distinction between the blazar subclasses.The population study of a large sample of blazar SEDs indicates a so-called blazar sequence, which characterizes an empirical inverse relationship between the bolometric γ-ray luminosity and the synchrotron peak frequency [30,31].This implies that the synchrotron peak frequency is anti-correlated with the ratio of Compton to synchrotron peak (Compton dominance), i.e., as the bolometric luminosity increases synchrotron peak frequency decreases.According to this, FSRQs exhibit the lowest synchrotron peak frequency and highest bolometric luminosity whereas HSP blazars have the lowest bolometric luminosity.It means FSRQs and BL Lacs can be quantified as high and low-power blazars, respectively.Also, this empirical observation feature suggests that accretion and blazar jet powers are related, and BL Lacs with heavier black holes and FSRQs with lighter black holes have the same γ-ray luminosity [32][33][34].The origin of the blazar sequence has not been fully understood yet and is generally attributed to the differences in the electron cooling efficiency or Doppler boosting [35,36].
In recent years, a class of high-z blazars has acquired particular importance in the field of blazar research [37,38].The GeV-TeV photons emitted from these sources suffer from partial attenuation during their propagation due to the interaction with low energy extragalactic background light (EBL) via γ − γ pair production [39,40].In this study, we mainly focus on the investigation of the physical properties of high-z blazars using publicly available data from broadband observations and SED modeling in the literature.The paper is organized as follows: A brief overview of the high redshift γ-ray blazars is given in Section 2. In Section 3, we describe the data sample used in this work and discuss the results in Section 4. A summary of this study and future prospects of high-z blazars are presented in Section 5.

High-z Blazars
We use a purely arbitrary definition of z > 0.4, in the present work, for characterizing the high redshift γ-ray blazars in the sky.However, this definition can be supported by the phenomenological evidence that beyond z > 0.4 the imprints of EBL absorption start becoming dominant in the HE γ-ray spectra measured by the Large Area Telescope on the Fermi spacecraft (Fermi-LAT) in the energy range above few tens of GeV [41].At these redshifts, the effects of cosmological evolution on structure formation and HE γ-ray propagation also becomes important.Thus, high-z γ-ray blazars detected by the Fermi-LAT are ideal tools to probe the early epochs in the evolution of the Universe as their bolometric luminosity is dominated by the HE emission [37,38,[42][43][44][45][46].The host galaxies of these sources are generally observed to host extremely massive black holes (M BH > 10 9 M ⊙ ) at their center.Therefore, high-z blazars can constrain the mass function at the high end of SMBHs and provide a theoretical understanding of the growth and evolution of black holes over cosmic time [47].Thus, on account of being a bright HE γ-ray source, blazars at high redshift are ideal objects for exploring the energetics and emission mechanisms under extreme conditions, setting strong constraints on the EBL intensity, and shedding light on the formation and growth of SMBHs in the early Universe and also their connection with the powerful relativistic jets.The detection and understanding of the physical properties of a few high-z blazars in a given redshift bin enable us to infer the total density and behavior of radio-loud AGNs in that redshift bin because for each blazar with the Bulk Lorentz factor (Γ ≈10-20) of the plasma within the jet, there must exist 2Γ 2 sources sharing the same properties but appearing different due to their jet pointing elsewhere [43].This complements the independent estimation based on AGNs with misaligned jets.Therefore, it is important to identify and study the high redshift blazars to understand the cosmic evolution of SMBHs.Recent radio observations have reported the discovery of the most distant blazar PSO J0309 + 27 located at a redshift of z ≈ 6.1 when the age of the Universe was less than ≈ 1 billion years [48].Its central engine is found to be a black hole of mass 1.45 +1.89  −0.89 × 10 9 M ⊙ with a bolometric luminosity of ≈ 8 × 10 46 erg s −1 .This corresponds to a seed black hole mass of 10 6 M ⊙ at z = 30 (when the first stars and galaxies are thought to have formed ≈ 15 million years after the Big Bang) [46].Such a rapid growth of SMBHs from the early stars has serious implications on various evolution models for the SMBHs in the jetted AGNs across cosmic time [43,[49][50][51][52][53].

Data Set
In this study, we have used publicly available data in the literature.A brief description of the data set, obtained mainly from 3 sources, is given below:  [58] have studied the properties of emission zones using the SED fitting results of more than 2700 blazars reported in [59].This large sample contains 1791 sources with known z (750 FSRQs, 843 BL Lacs and 198 BCUs).The catalog reports derived values of the magnetic field (B) in the emission region, energy density of electrons (U e ) and magnetic field energy density (U B ) from the SED fitting.For such blazars, we have collected the luminosities of γ-ray at 1 GeV (log L γ ), X-ray at 1 keV (log L X ), optical at 2.43 ×10 14 Hz (log L O ) and radio at 1.4 GHz (log L R ) from [60].
The blazar central engines catalog is a subset of the blazar emission regions catalog.And, both catalogs are subsets of the 4LAC-DR3 blazars catalog as they only include 4LAC-DR3 catalog sources with known redshift.Therefore, blazar central engines catalog (1026 sources) ⊆ blazar emission regions catalog (1791 sources) ⊆ 4LAC-DR3 blazars catalog (1811 sources).

Redshift Distribution
The measured z-distribution of blazars listed in the 4LAC-DR3 blazars catalog is shown in the left panel of Figure 1.The two blazar-subclasses, BL Lacs and FSRQs, show different z-distributions.The distribution of BL Lacs typically peaks around a median redshift value of z ≈ 0.31 and is relatively narrow with a standard deviation of σ ≈ 0.39.Whereas, the distribution of FSRQs peaks around a median of z ≈ 1.11 and is broader with σ = 0.66.The BCUs with relatively flatter distribution occupy the space between BL Lacs and FSRQs.The z-values play an important role in the classification of BCUs in the two subclasses of blazars [61].The right panel of Figure 1 shows the populations of BL Lacs, FSRQs, and BCUs in different redshift regions.It is observed that in the low-z region, the blazar population is mostly dominated by BL Lacs (∼ 80%) and few FSRQs (13%) and BCUs (7%).In the high-z regime, FSRQs dominate the blazar population with 62% contribution whereas BL Lacs (30%) and BCUs (8%) correspond to the remaining 38% population.Thus, BL Lacs are typically located in the nearby universe, within a physical distance of about 2-3 billion light-years, whereas FSRQs are found at longer distances than BL Lacs, and may be associated with more distant and younger galaxies.The exact cause of these redshift distributions is not well understood, yet many explanations have been proposed in the literature.The first explanation is the selection bias in observations, i.e., BL Lacs are low luminous sources and it is very difficult to detect them at high redshifts.Even if they are detected, measurement of their redshift is challenging due to their line-less optical spectra [29,62].Thus, the population of BL Lacs is less numerous at higher redshifts.As only strongly luminous sources can be detected in the high-z regime, FSRQs being intrinsically more luminous than BL Lacs dominate the high redshift Universe.The second explanation is that FSRQs evolve into BL Lacs over cosmological timescales.Therefore, high-z sources (early and young sources) are FSRQs and they evolve into BL Lacs, which are found at lower redshifts [6,42,63].The highest redshift BL Lacs and FSRQ, detected by the Fermi-LAT, are 4FGL J1219.0+3653 and 4FGL J1510.1+5702(with radio association of NVSS J151002+570243) located at z ≈ 3.59 and z≈ 4.31, respectively [64,65].They correspond to a population of blazars existing just a few billion years after the Big Bang.

γ-ray Spectral Features
For the spectral study, we consider sources reported in the 4LAC-DR3 blazars catalog with a detection significance ('Signif_Avg' in the Fermi catalog) of more than 10σ in the energy range 100 MeV-1 TeV.There are 1111 such sources comprising 535 BL Lacs, 521 FSRQs, and 55 BCUs.In this catalog, the observed γ-ray flux points (i.e., without taking EBL attenuation into account) have been fitted with both Power-Law (PL) and Log Parabola (LP) spectral models.Here, K is the flux density (erg cm −2 s −1 ), Eo is the pivot energy (in MeV) at which the error in differential flux is minimal, Γ is PL spectral index, α is LP index and β is the LP curvature parameter.For a fair comparison, we have calculated the modified LP index (α ′ ) by computing the slope of the LP spectrum at energy 1 GeV/(1 + z).This energy corresponds to the photon energy of 1 GeV in the source rest frame.Therefore, α ′ represents the slope of the observed spectrum at the energy corresponding to 1 GeV in the source frame.Left panels in Figure 2 show distributions of Γ, α ′ and β for high and low-z blazars.The two sample Kolmogorov-Smirnov (K-S) test results in p-values of ≈0 for all three parameters.This suggests that distributions of Γ, α ′ , and β parameters are intrinsically different for high and low-z blazars.A Gaussian fit to these distributions gives mean and standard deviation, µ ± σ, values of 2.33 ± 0.25 (2.02 ± 0.23), 2.10 ± 0.39 (1.84 ± 0.38), and 0.11 ± 0.07 (0.07 ± 0.7) for Γ, α ′ and β, respectively, in case of high-z (low-z) blazars.The right panels in Figure 2 depict variations of these parameters as a function of (log (1 + z)).Pearson correlation analysis shows that Γ, α ′ and β parameters are positively correlated with (log (1 + z)) with the correlation coefficient (ρ) values of 0.6, 0.33, and 0.34, respectively, and p-values of ≈0.Therefore, blazars at higher redshifts will have softer γ-ray spectra (higher values of Γ and α ′ ) and with more curvature (higher values of β) than low-z blazars.Similar correlations have also been obtained between BL Lacs and FSRQs in the literature [54,66].This may be the manifestation of the underlying particle spectrum involved in the HE emission (detail discussion in Section 4.7) and absorption of HE γ-ray photons by EBL via γ-γ pair production during their propagation from source to the observer [67,68].The absorption probability of HE photons traveling over cosmological distances is characterized by optical depth (τ) or opacity, which strongly depends on z, γ-ray photon energy (E) and density of EBL photons [39,69].The EBL absorption effects start dominating at energies above 30 GeV for redshifts beyond z > 0.4 [41].Therefore, the intrinsic (emitted) γ-ray spectrum of the source gets modified and the observed spectrum is given as For high-z blazars, τ will be large due to strong EBL absorption.Therefore, the observed photon spectrum will be soft with more curvature [70,71].Hence, the detection and accurate measurements of HE and VHE spectra of high-z blazars are important (as they will have an imprint of EBL absorption effects) in constraining the EBL energy density in extragalactic space [39,72].

Luminosity Features
For luminosity studies, we have selected sources from the blazar emission regions catalog.Out of 1791 sources in the catalog, 1789 have γ-ray luminosities (L γ ) at 1 GeV, 1030 have X-ray luminosities (L X ) at 1 keV, 1059 have optical luminosities (L O ) at 2.43 × 10 14 Hz and 1581 have radio luminosities (L R ) at 1.4 GHz.The left panels in Figure 3 show distributions of log L γ , log L X , log L O and log L R for high and low-z blazars [60].The right panels represent variations of these parameters as a function of luminosity distance (D L ) measured in Mega-parsec (Mpc).Note that, D L is a function of z.A Gaussian fit to the distributions of log L γ , log L X , log L O , and log L R gives µ ± σ values of 45.80 ± 0.89 (43.78 ± 0.86), 45.20 ± 0.72 (44.01 ± 0.88), 45.65 ± 0.67 (44.41 ± 0.52), and 42.91± 0.89 (41.0 ± 0.92) for high-z (low-z) blazars, respectively.We observe that, in all the wavebands considered above, high-z blazars are at least an order of magnitude more luminous than low-z blazars.All the luminosities log L γ , log L X , log L O , and log L R show strong positive correlations with log D L having corresponding ρ values of 0.90, 0.77, 0.85, and 0.85, respectively.We also perform a linear regression analysis of the scatter plots of these luminosities vs log D L .The slopes of linear regression indicate that observed L γ , L X , L O and L R vary as ∼ D 2.27  L , D 1.47  L , D 1.51   L and D 2.14 L , respectively.As mentioned in Section 4.1, these observed features are consistent with the fact that high-z sources with strong luminosity can only be detected due to selection bias.The high-z regime is dominated by FSRQs, which are intrinsically more luminous than BL Lacs, which are mostly observed at lower redshifts.Thus, high-z blazars (mostly FSRQs) are the most luminous sources observed in the sky across all wavebands.

Broadband SED Properties
We investigate the broadband SED properties using the sources in blazars central engines catalog.This catalog contains typical SED characteristics like synchrotron peak frequency (ν syn ), inverse Compton (IC) peak frequency (ν IC ), energy fluxes at synchrotron (νF ν,syn ) and IC (νF ν,IC ) peak frequencies and Compton-Dominance parameter (CD = νF ν,IC /νF ν,syn ).For fair comparisons, we estimated ν syn and ν IC in the source frame by multiplying them with a factor (1 + z).The left panels in Figure 4 show the histogram distributions of ν syn , ν IC , νF ν,syn , νF ν,IC and CD for high and low-z blazars.The right panels in Figure 4 depict variations of these parameters with redshift (log (1 + z)).The two sample K-S test between high and low-z blazars results in p-values of ≈0 for all the parameters.This suggests that the histogram distributions of these parameters are intrinsically different for high and low-z blazars.A Gaussian function fit to the distributions of ν syn , ν IC , νF ν,syn , νF ν,IC and CD gives µ ± σ values of 13.42 ± 1.04 (15.32 ± 1.46), 22.00 ± 1.07 (23.53 ± 1.68), −11.90 ± 0.45 (−11.65 ± 0.46), −11.42 ± 0.55 (−12.05 ± 0.60) and 0.48 ± 0.50 (−0.4 ± 0.40) for high-z (low-z) blazars, respectively.The Pearson correlation analysis indicates that ν syn and ν IC are anti-correlated with log (1 + z).Whereas, νF ν,IC is positively correlated with log (1 + z) with ρ=0.46.However, a weak anti-correlation is observed between νF ν,syn and log (1 + z) with ρ = −0.34.CD has a strong positive correlation with log (1 + z).This is broadly consistent with the blazar sequence since at high-z only the higher luminosity objects can be detected, and their peak frequencies are shifted to lower values with increasing CD.Therefore, high-z blazars are mostly LSP blazars, and their HE hump typically peaks at MeV energies just below the energy range covered by the Fermi-LAT.Their X-ray spectrum is harder and it usually corresponds to the rising part of the HE hump.Being strong and bright X-ray emitters, they are sometimes called MeV blazars [50,73].The observations in the hard-X-ray band are ideal to detect the high-z blazars.Also, due to the low synchrotron peak frequency (typically in radio/IR waveband), the accretion disk emission is visible in the optical-UV region for most of the high-z blazars.Therefore, disk luminosity and black hole mass can be estimated by modeling the accretion disk emission [74].It is interesting to notice the variation of ν syn with respect to log (1 + z) in the top right panel of Figure 4.The low-z blazars show large variation in ν syn (from 10 12 Hz to 10 19 Hz) making a vertical branch, whereas high-z blazars show little variation in ν syn (from 10 12 to 10 14 Hz) making a horizontal branch in the scatter plot.Consequently, we have a bimodality in the top left panel of Figure 4.This is in agreement with a recent blazar sequence study [31] in which it has been observed that the synchrotron peak variation is mostly observed in BL Lacs and is insignificant in FSRQs.Since we mostly have BL Lacs at low-z and FSRQs at high-z in the data sample, a bimodality is obtained in the distributions of ν syn and two branch feature in the variation of ν syn with log (1 + z) in Figure 4.

Central Engine Properties
The blazar central engine consists of an SMBH at the center and an accretion disk surrounding it.It is characterized by the mass of SMBH (M BH ), accretion efficiency (η), and accretion disk luminosity (L Disk ).The sources in the blazar central engines catalog are considered for this study.The left panel in Figure 5 shows distribution of M BH values for high and low-z blazars.In both the cases, M BH values range from ∼2.5 × 10 6 M ⊙ to ∼1.75 × 10 10 M ⊙ .The average M BH values for high and low-z blazars are similar with ⟨log M BH,Low-z ⟩ ≈ ⟨log M BH,High-z ⟩ ≈ 8.6 and standard deviations are 0.54 and 0.68, respectively, when fitted with a Gaussian function.A two sample Kolmogorov-Smirnov (K-S) test is performed to see whether the two populations of M BH values are derived from the same parent population.The obtained p-value of 0.11 suggests that there is a high probability both samples belong to the same parent source population.However, their correlations with z are different as discussed below.The right panel in Figure 5 shows the scatter plot of log M BH and log (1 + z) for high and low-z blazar populations.Overall, for the total blazar population, there is a very weak correlation between log M BH and log (1 + z) as indicated by ρ = 0.16 with p-value of 2.88 ×10 −7 .However, by considering only high-z blazars, the correlation becomes relatively strong with ρ = 0.41 and p-value of 3.89×10 −28 .It suggests that more massive central black holes are associated with high-z blazars.This is likely due to selection bias in observations, as at high redshift only the most luminous (blazars with more massive black holes) sources can be detected.However, for low-z blazars, there is very weak anti-correlation between log M BH and log (1 + z) with ρ = −0.23 and p-value of 1.75 × 10 −5 .Although the K-S test suggests that there is a high probability that M BH values of low and high-z blazars belong to the same parent population, their correlations with z are different in low and high-z regimes.Therefore, M BH values of low and high-z blazars may not belong to the same source population.
To study the accretion disk properties, we consider the sources in the blazar central engines catalog.We have excluded 342 sources from this catalog as only upper limit values for L Disk are reported.The upper left panel in Figure 6 shows the histogram distributions of L Disk for blazars included in this study.p-value of 2.08 ×10 −26 obtained from the two sample K-S test indicates that the distributions for low and high-z blazars are significantly different and belong to different parent distributions.A Gaussian function fitting to the distributions gives the mean values of log L Disk (⟨ log L Disk ⟩) as 45.75 and 44.48 with standard deviations of 0.67 and 0.86 for high and low-z blazars, respectively.Therefore, in general, high-z blazars (⟨L Disk,High-z ⟩ ≈ 5.6 × 10 45 ergs −1 ) have one order of magnitude higher accretion disk luminosity than lowz blazars (⟨L Disk,Low z ⟩ ≈ 3.0 × 10 44 ergs −1 ).The upper right panel in Figure 6 shows the variation of log L disk as a function of log (1 + z).The Pearson correlation analysis suggests a strong positive correlation between the two parameters with ρ = 0.67 and p-value of 1.28 We also calculate the Eddington ratio (λ Edd = L Disk /L Edd ) using the definition of Eddington luminosity as L Edd = 1.26 × 10 38 (M BH /M ⊙ ) erg s −1 .The corresponding distri- butions for high and low-z blazars and their variations as a function of redshift are shown in the lower left and lower right panels of Figure 6, respectively.The λ Edd distributions for high and low-z blazars are intrinsically different as indicated by the K-S test p-value of 8.74 ×10 −13 .The Gaussian function fitting to the distributions with mean and standard deviation (µ, σ) gives (−0.92,0.46)and (−1.53,0.66)for high and low-z blazars, respectively.This implies that high-z blazars accrete the surrounding matter with a high Eddington ratio of ⟨ λ Edd,High-z ⟩ = 0.12 compared to the low-z blazars with ⟨ λ Edd,Low-z ⟩ = 0.03.Therefore, high-z blazars have radiatively more efficient accretion disks than low-z blazars.This is also confirmed by the mild positive correlation obtained between λ Edd and log (1 + z) with ρ = 0.37 and p-value of 2.46 × 10 −4 .The linear regression line in the scatter plot with m= 0.37 ± 0.10 suggests that observed λ Edd varies as ∼(1 + z) ≈0. 37.
Observation of blazars with strong luminosity, radiatively efficient accretion disk, and located at high redshift is likely connected with the cosmic evolution of blazars [6,42,63].As discussed in Section 4.1, BL Lacs dominate the blazar population at low redshifts and FSRQs at high redshifts.According to the hypothesized evolutionary scenario, the FSRQs (likely high-z blazars) evolve into BL Lacs (likely low-z blazar) over cosmological time scales.At high-z, there is likely a gaseous matter-rich environment surrounding the central black hole for accretion, hence the higher values of λ Edd .As they evolve into low-z blazars, the environment gradually gets depleted by accretion onto a black hole, the accretion disk becomes radiantly inefficient and accretion proceeds via advection-dominated flow [75].Alternatively, the observed features can also be attributed to the selection bias effects in observations as it is difficult to detect blazars having low λ Edd and low L Disk at high redshifts.

Dynamics of Growth of SMBHs
High-z blazars provide a unique opportunity to explore the evolution of SMBHs and the role of relativistic jets in them [51,76].Following the detailed calculations given in [46,77], we estimate the mass of black hole seed (M BH,seed ) required to produce currently observed SMBHs at a given z, assuming that it evolves with constant λ Edd and constant accretion efficiency defined as η = L Disk / ṀBH c 2 during the entire evolutionary process.Accounting for the loss of accretion mass energy in the form of outgoing disk luminosity, the rate of M BH growth is given by where t s is the characteristic accretion timescale (also known as Salpeter time) independent of M BH .It is given as Therefore, required black hole seed mass (M BH,seed ) for the formation of SMBH of mass M BH is given by where t is the total time available during which the seed black hole accretes and grows, τ is the characteristic e − f olding timescale defined as The accretion efficiency parameter η, believed to depend on the black hole spin, can vary between 0.1-0.3 or 0.4 for slowly to rapidly spinning black holes [77,78].We consider the sources in blazar central engines catalog for whom the classes (i.e., FSRQs, BL Lacs, or BCUs) are defined in the 4LAC-DR3 blazar catalog.For such blazars (558 FSRQs,53 BL Lacs, and 50 BCUs), we have calculated M BH,seed using Equation ( 6) by evolving them back to z = 30 (corresponding to the Universe with age less than the first 100 million years when the first stars and galaxies are assumed to have formed) [79,80].In our calculations, we use currently observed λ Edd throughout their evolution to see if such a scenario explains the observed SMBHs.The calculations are performed for two constant accretion efficiencies of η = 0.1 and η = 0.3 [77].For calculating t in Equation ( 6), we use the standard flat ΛCDM model of the Universe with parameters H 0 = 70 km s −1 Mpc −1 , Ω M = 0.30 and Ω Λ = 0.70.The final results depend on the assumptions made regarding various parameters like redshift of the seed formation, η and λ Edd .
Figure 7 shows the variations of the estimated values of M BH,seed as a function of log (1 + z) for η = 0.1 (left panel) and η = 0.3 (right panel).In these plots, we also show the expected mass ranges of different astrophysical objects present at z = 30, which could have been the seeds for the formation of observed SMBHs in the centers of galaxies.They are (i) population III star black hole remnants: 10-10 3 M ⊙ (cyan color) [81,82], (ii) seeds from stellar dynamical processes: 10 3 − 10 4 M ⊙ (gray color) [83][84][85][86][87], and (iii) direct collapse black holes: 10 4 -10 6 M ⊙ (green color) [88,89].The presence of SMBHs in blazars corresponding to M BH,seed values lying below the horizontal line logM BH,seed = 6 in Figure 7 can be explained by the evolution of one of the 3 astrophysical seed objects mentioned above.However, a significant number of blazars lie above the line logM BH,seed = 6, and plausible explanations for this are given below.The left panel in Figure 7 indicates that, for η = 0.1, about 26% of the blazar population requires more massive progenitors than what direct collapse models or stellar dynamical processes predict.This fraction goes up to ∼69% when we consider η = 0.3 (right panel of Figure 7).Therefore, high values of accretion efficiencies, i.e., high conversion rate of rest mass energy to luminous energy, are not realistic to explain the observed SMBHs.This is consistent with the similar conclusions derived in the literature [77].Alternatively, these blazars must accrete the matter with a higher Eddington ratio than the currently observed λ Edd for a significant time during their evolution to explain the growth of observed SMBHs.Thus, the higher the η, the higher the λ Edd we need to explain the growth of the SMBHs.A major fraction of the SMBHs in BL Lacs (∼75% and ∼92% of the total BL Lacs corresponding to η = 0.1 and 0.3, respectively) require seed mass more than the expected seed mass from direct collapse models.This can be due to their observed low accretion disk luminosity and hence low Eddington ratio (λ Edd ).Therefore, these objects must have accreted the matter with higher λ Edd during their evolution to explain the observed SMBHs.This implies that there was a conducive gas-rich environment around the central engines in their earlier growth time (high-z).These objects would have radiatively efficient accretion disks in their adulthood, meaning they were in the FSRQ phase (high λ Edd phase) and as their gas-rich environment became depleted, their accretion disk became radiatively inefficient and they evolved into BL Lacs at relatively low-z.It is evident from the figure that a significant number of high-z blazars (or radio-loud AGN) require higher mass progenitors than the currently hypothesized seed populations.Therefore, we need higher λ Edd or sometimes the occurrence of super-Eddington episodes to explain the growth of SMBHs.Such scenarios are invoked to explain the growth of SMBHs in the recently discovered high redshift AGNs J1205-0000 (z = 6.699) [90] and PSO J006+39 (z = 6.621) [91].Although there are possibilities for super Eddington accretion episodes in gas-rich or dusty strong wind environments near the central engines [92,93], whether such episodes are sustainable or not for a significant time is an open question in the study of the growth of SMBHs at high redshifts.Another possible scenario has been proposed in the literature, where in the presence of a jet, not all the gravitation potential energy of in-falling matter is converted into luminous energy but some part of it can go to increase the magnetic field energy inside the jet, which is necessary for jet launching processes [3].Therefore, for a given accretion efficiency, there will be fewer luminous disks and hence lower (λ Edd ) compared to the normal accretion scenario.This implies that, in the presence of a jet, black holes grow faster with a large accretion rate without increasing the disk luminosity [73,94,95].Cosmological simulations were also been explored to investigate the growth of seed black holes through processes other than accretion such as merger and AGN feedback [96,97].Similar calculations and conclusions have also been reported in a recent study [46].Hence, detection and detailed study of high-z blazars are important to understand the growth of SMBHs and their connection with relativistic jets.It is important to note that the conclusions drawn here rely on the reported values of Λ Edd , L Disk , and M BH in the blazar central engines catalog.These values were derived from scaling relations using single-epoch optical spectra.They are indirect methods and have inherent assumptions and limitations which may lead to significant bias [98].Therefore, the estimation of M BH from different approaches like reverberation mapping, etc., are crucial [19,99].

Properties of Emission Region
In the one-zone leptonic emission model, the emission region is mainly characterized by its magnetic field B and relativistic electron energy distribution having peak Lorentz factor γ p .We use log B and log γ p parameters reported in the blazar emission regions catalog for this study.For some sources, ranges of log B and log γ p are reported in the catalog, and we take the midpoint of the range as the value for log B and log γ p for them.The left panels in Figure 8 show the histograms of log B and log γ p for high and low-z blazars.The right panels show the variations of these parameters as a function of log (1 + z).p-values of two sample K-S tests for high and low-z distributions of log B and log γ p parameters are close to 0. This suggests that the distributions of B and γ p for high and low-z blazars are intrinsically different and belong to the same parent distribution.A Gaussian function fitting to the histograms of log B gives µ ± σ of 0.65 ± 1.68 (B = 4.46 G) and 0.04 ± 1.53 (B = 1.09G) for high and low-z blazars, respectively.Also, there exists a mild positive correlation between B values and log (1 + z) with ρ = 0.26 and the P value of 1.24×10 −28 .Therefore, high-z blazars require emission zones with higher magnetic fields compared to low-z blazars.The histogram distribution of log γ p for high-z blazars shows a clear double hump structure with peaks around log γ p values of 2 and 4. Similarly, the scatter plot of log γ p and log (1 + z) shows two horizontal branches around these peak values.A careful analysis suggests that these peak values actually correspond to the peak values of distributions of log γ p for all FSRQs and BL Lacs.Overall, two peaks in the histogram distributions and two branches in the scatter plots actually correspond to the FSRQs and BL Lacs, respectively.Instead of fitting the Gaussian profile to the histograms of log γ p , we calculated the arithmetic mean (µ ′ ) and standard deviation (σ ′ ) for the distributions.The µ ′ values for high and low-z blazars are found to be 2.90 and 3.79 corresponding to peak Lorentz factors of 7.94 ×10 2 and 6.17 ×10 3 , respectively.The scatter plot shows two distinct branches corresponding to the FSRQs and BL Lacs, and there is also a mild anti-correlation between log γ p and log (1 + z) with ρ = −0.47.These results are consistent with the blazar emission models.At high redshifts, most of the sources are FSRQs and BL Lacs dominate the population at low redshifts.FSRQs contain photon-rich environments like BLR and DT owing to radiatively efficient accretion disk, and relativistic electrons quickly lose their energy by interacting with ambient photons via the IC process.Therefore, these blazars have high Compton dominance, and the electron population has low γ p .Consequently, we observe blazars with low synchrotron frequency peaked SEDs.On the other hand, BL Lacs have a photon-starved environment surrounding them.Relativistic electrons can attain high γ p in the absence of target photon fields.Therefore, their SED peak lies at higher frequencies with very low Compton dominance [6,56,58,63].The variation of log γ p as a function of log (1 + z) also contributes to the observed behavior of spectral parameters measured by the Fermi-LAT.The higher the z, the higher the efficiency of the disk, and the higher the density of external photon fields, resulting in the electron energy distributions with lower γ p .Therefore, this can also contribute to the observed soft photon spectral indices and higher curvatures for high-z blazars.

Summary
In this work, we studied the properties of high-z (z > 0.4) blazars using the data available in the literature.We compared the observed properties of SEDs and derived properties of the central engines and emission zones of high and low redshift blazars.A summary of the main results obtained in this study is given below:

•
Redshift distribution study of 4 LAC-DR3 blazars catalog suggests that the FSRQs dominate the high-z blazar population with 62% composition followed by BL Lacs and BCUs with 30% and 8% compositions, respectively.However, in the low-z regime, BL Lacs dominate the total blazar population with 80% composition and FSRQs (13%) along with BCUs (7%) being the minority.The reasons for these anomalies are likely due to observational biases and/or due to hypothesized cosmic blazar evolution theories.We also estimated the masses of black hole seed (M BH,seed ) required for the formation of currently observed SMBHs in blazars using simple exponential growth theory for the formation of SMBHs from the black hole seed via the accretion process.With η = 0.1(0.3)and currently observed constant λ Edd , it is found that 26% (69%) of the total blazar population requires more massive seed progenitors than what the direct collapse models or stellar dynamical processes predict.This implies that low accretion efficiencies and high λ Edd (sometimes super Eddington) episodes are required to explain the growth of observed SMBHs.It is also noticed that most of the BL Lacs (75% and 92% for η = 0.1 and 0.3, respectively) are required to have higher λ Edd than currently observed λ Edd to explain the growth of their SMBHs.Therefore, BL Lacs should have had radiatively efficient accretion disk (FSRQ phase) during their lifetime before evolving into BL Lacs having radiatively inefficient accretion disk.

•
According to the single zone leptonic emission model, emission zones of high (low) redshift blazars are required to have magnetic fields ⟨B ⟩ = 4.46G(1.09)G and electron energy distributions with γ p values ⟨γ p ⟩ = 7.94 × 10 2 (6.17 ×10 2 ) to explain their broadband SEDs.They are broadly in agreement with the widely used blazar emission models.High-z blazars (mostly FSRQs) contain low-energy photon-rich environments and relativistic electrons quickly lose their energies through IC scattering.Therefore, their distribution has low γ p .On the other hand, low-z blazars (mostly BL Lacs) have a low-energy photon-starved environment surrounding them and relativistic electrons can attain higher γ p .This also contributes to the observed distribution of synchrotron and IC peak frequencies in blazar SEDs.
Therefore, high-z blazars are a unique class of objects in the distant Universe.Their further detection and detailed broadband emission study will shed more light on radiative and particle acceleration mechanisms under extreme conditions, cosmic evolution of SMBHs and relativistic jets, and intensity of EBL.Deep X-ray and γ-ray surveys are required to have an unbiased large population of blazars for further study of high-z blazars.In the coming decade, the extended ROentgen Survey with an Imaging Telescope Array (eROSITA [100]), the Large High Altitude Air Shower Observatory (LHAASO [101]), the Southern Wide-field Gamma-ray Observatory (SWGO [102]), and the Cherenkov Telescope Array (CTA [103]) are expected to detect large number of high redshift blazars.

Figure 3 .
Figure 3. (Left) Histograms of the luminosities of γ-ray at 1 GeV (log L γ ), X-ray at 1 keV (log L X ), optical at 2.43 ×10 14 Hz (log L o ) and radio at 1.4 GHz (log L R ) for sample high and low-z blazars.(Right) Variations of log L γ , log L X , log L o and log L R as a function of luminosity distance (log D L ).

Figure 5 .
Figure 5. (Left) Histograms of log M BH values for high and low-z blazars.(Right) Scatterplot of log M BH and log (1 + z).

Figure 6 .
Figure 6.(Upper Left) Histograms of L Disk for high and low redshift blazars.(Upper Right) Scatter plot of log L Disk vs log (1 + z).(Lower Left) Histograms of the Eddington ratios (λ Edd ) for high and low redshift blazars.(Lower Right) Scatter plot of log λ Edd as function of log (1 + z).

Figure 7 .
Figure 7. Variations of the estimated seed M BH,seed values at z = 30, required for the formation of observed SMBHs in blazars, as a function of log (1 + z), assuming the BH accretes with observed λ Edd,obs and with η = 0.1 (left panel) and 0.3 (right panel).The shaded regions correspond to expected mass ranges of Population III star BH remnants (10-10 3 M ⊙ ; cyan color), seed black holes from stellar dynamical processes (≥ 10 3 M ⊙ ; gray color) and black holes from direct collapse (10 4 -10 6 M ⊙ ; green color).

Figure 8 .
Figure 8. (Left) Histogramsof log B and log γ p for high and low redshift blazars.(Right) Variation of log B and log γ p values with log (1 + z).
• Observed Fermi-LAT photon spectra for the high redshift blazars are found to be softer with higher curvature compared to low redshift blazars.The histograms of PL index (Γ), LP slope at 1 GeV/(1 + z) (α ′ ) and LP curvature (β) roughly follow Gaussian distributions with µ ± σ values of Γ = 2.33 ± 0.25(2.02± 0.23), α ′ = 2.10 ± 0.39(1.84± 0.38) and β = 0.11 ± 0.07(0.07± 0.07) for high (low) redshift blazars.Also, Γ, α ′ and β exhibit a positive correlation with log (1 + z).This is likely due to the fact that high-z blazars suffer more EBL attenuation of HE and VHE γ photons, leading to softer γ-ray spectra than low-z blazars.Also, different spectral shapes of the underlying energy distribution of emitting particles in low and high-z blazars contribute to the observed Fermi-LAT γ-ray spectral features.= 6.31 × 10 45 (6.03 × 10 43 ), ⟨L X ⟩ = 1.58 × 10 45 (1.02 × 10 44 ), ⟨L O ⟩ = 4.46 × 10 45 (2.57× 10 44 ) and ⟨L R ⟩ = 8.13 × 10 42 (1.0 × 10 41 ), respectively.These luminosities have a strong positive correlation with the luminosity distance (D L ) with functional dependencies of L γ ∝ D 2.27 L , L X ∝ D 1.47 L , L O ∝ D 1.51 L and L R ∝ D 2.14 L .• The investigation of broadband SED properties reveals that a major fraction of high-z blazars are LSP objects.Their HE hump peaks at MeV energies and they have higher Compton dominance than low-z blazars.The logarithmic histogram distributions of rest frame synchrotron peak frequency (ν syn ), rest frame IC peak frequency (ν IC ), and Compton dominance (CD) roughly follow a Gaussian profile with µ ± σ values of log ν syn = 13.42 ± 1.04(15.32± 1.46), log ν IC = 22.00 ± 1.07(23.53± 1.68) and log CD = 0.48 ± 0.50 ± (−0.4 ± 0.40) for high (low) redshift blazars.Also, the Pearson correlation study indicates that log ν syn and log ν IC are in anti-correlation with with log (1 + z).log CD shows a strong positive correlation with log (1 + z).This is in accordance with the blazar sequence.High-z blazars, mostly being FSRQs, are intrinsically more luminous and therefore have high Compton dominance.• The SMBH mass distribution study of the blazar central black holes suggests that the average mass of SMBHs for high and low-z blazars is the same with ⟨log M BH ⟩ = 8.6.Low-z blazars show a weak anti-correlation with log (1 + z) whereas high-z blazars have a mild positive correlation with log (1 + z).Also, the average accretion disk luminosity and Eddington ratio of high (low) redshift blazars are ⟨L Disk ⟩ = 5.62 × 10 45 erg s −1 (2.4 × 10 44 erg s −1 ), ⟨λ Edd ⟩ = 0.12(0.03)with both the parameters showing a positive correlation with log (1 + z).• • Observationally, high-z blazars are the most luminous objects in the γ-ray Universe.The mean γ−ray, X-ray, optical and radio luminosities (in erg s −1 ) of high (low) redshift blazars are estimated as ⟨L γ ⟩