###### Abstract

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 _Fermi_ satellite. We
also study the properties of their emission regions and the central engines
and discuss cosmological and astrophysical implications.

###### keywords:

active galactic nuclei; supermassive black holes; non-thermal emission; high
redshift Universe

1 1 0 11 December 2023 1 March 2024 7 March 2024 https://doi.org/ Central
Engine and Spectral Energy Distribution Properties of High Redshift Gamma Ray
Blazars Central Engine and Spectral Energy Distribution Properties of High
Redshift Gamma Ray Blazars Anilkumar Tolamatti 1,2,*, Krishna Kumar Singh 1,2,
and Kuldeep Kumar Yadav 1,2 Anilkumar Tolamatti, Krishna Kumar Singh and
Kuldeep Kumar Yadav Tolamatti, A.; Singh, K.K.; Yadav, K.K.
<EMAIL_ADDRESS>(A.T.)

##  1 Introduction

Blazars are bright beacons in the Universe. They represent a subclass of the
radio-loud active galactic nuclei (AGNs) which are characterized as elliptical
galaxies hosting supermassive black holes
(SMBHs;$M_{BH}:10^{6}-10^{10}M_{\odot}$) at their centers and oppositely
oriented powerful jets of plasma ejected from the poles of SMBH Urry &
Padovani (1995); Padovani et al. (2017). 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 Blandford & Znajek (1977). 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
Blandford & Payne (1982). 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 Maraschi & Tavecchio (2003); Ghisellini et al. (2011); Chen
et al. (2015); Zhang et al. (2022). The hybrid model, a combination of the
above two mechanisms, speculates about the observed differences in different
types of AGNs with jets Meier (1999). 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) $\gamma$-rays. They are also thought to be the apparent
manifestations of AGN activity during their early evolutionary stages
Komissarov (2021). 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 Urry & Padovani (1995); Padovani et
al. (2017). 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) $\gamma$-ray emissions, and apparent superluminal motion Wagner (1995);
Giommi et al. (2021); Lister et al. (2019); Weaver (2020); Shukla (2020);
Singh (2020a, b); Tolamatti et al. (2022). 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 Singh et al. (2022).

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 $\gamma$-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 Marscher (1985)) 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 Sikora (1994); Ghisellini (2009). 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 Aharonian (2000);
Bottcher (2013); Petropoulou (2015).

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 Stickel (1991). 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 Abdo (2010);
Padovani (2012); Giommi (2012). For LSP blazars, the synchrotron peak
frequency of $\leq 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
Ghisellini et al. (2011); Padovani et al. (2017). 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
$\gamma$-ray luminosity and the synchrotron peak frequency Fossati (1998);
Prandini (2022). 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
$\gamma$-ray luminosity Ghisellini (2016, 2017); Singh (2019). 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 Ghisellini (1998); Fan (2017).

In recent years, a class of high-$z$ blazars has acquired particular
importance in the field of blazar research Paliya (2020); Sahakyan (2020). 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 $\gamma-\gamma$ pair production Stecker (2006);
Singh (2020). 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 $\gamma$-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.

##  2 High-$z$ Blazars

We use a purely arbitrary definition of $z>0.4$, in the present work, for
characterizing the high redshift $\gamma$-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
$\gamma$-ray spectra measured by the Large Area Telescope on the Fermi
spacecraft (_Fermi_ -LAT) in the energy range above few tens of GeV Tolamatti
(2022). At these redshifts, the effects of cosmological evolution on structure
formation and HE $\gamma$-ray propagation also becomes important. Thus,
high-$z$ $\gamma$-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 Ghisellini (2013); Ackermann et al.
(2015); Paliya (2020); Sahakyan (2020); Inayoshi (2020); Belladitta (2022);
Ajello (2014). The host galaxies of these sources are generally observed to
host extremely massive black holes ($M_{BH}>10^{9}M_{\odot}$) 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 Johnson (2016). Thus, on account of being a
bright HE $\gamma$-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
($\Gamma~{}\approx$10–20) of the plasma within the jet, there must exist
$2\Gamma^{2}$ sources sharing the same properties but appearing different due
to their jet pointing elsewhere Ghisellini (2013). 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~{}\approx$ 6.1 when the age of the Universe was less than $\approx$ 1
billion years Belladitta (2020). Its central engine is found to be a black
hole of mass 1.45${}_{-0.89}^{+1.89}\times 10^{9}$ M⊙ with a bolometric
luminosity of $\approx~{}8\times 10^{46}$ erg s-1. This corresponds to a seed
black hole mass of 10${}^{6}M_{\odot}$ at $z~{}=~{}30$ (when the first stars
and galaxies are thought to have formed $\approx$ 15 million years after the
Big Bang) Belladitta (2022). 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 Ajello (2009); Ghisellini (2010); Volonteri
(2010, 2011); Ghisellini (2013); Sbarrato (2015).

##  3 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:

  * •

4LAC-DR3 Blazars Catalog: The fourth catalog of AGNs detected by the _Fermi_
-LAT (4LAC), derived from the third data release (DR3) of the 4FGL catalog
based on 12 years of $\gamma$-ray data in the energy range above 50 MeV,
provides measurements of spectral parameters, SEDs, yearly light curves, and
associations for all sources Ajello (2022); Abdollahi (2022). We consider both
low and high Galactic latitude sources of blazar type (i.e., BL Lacs, FSRQs,
and Blazar Candidates of Uncertain types (BCUs)). Out of 3743 blazar-type
sources, only 1811 have redshift measurements. Among 1811 blazars with known
$z$, there are 875 BL Lacs, 792 FSRQs and 144 BCUs.

  * •

Blazar Central Engines Catalog : The central engine properties like black hole
mass (MBH) and accretion disk luminosity (Ldisk) have been reported in
literature Paliya (2021) for 674 emission line blazars and 346 absorption line
blazars observed with the _Fermi_ -LAT. For these 1030 blazars, the values of
MBH and Ldisk are derived using optical spectroscopic line features reported
either in the 16th data release of the Sloan Digital Sky Survey (SDSS-DR16)
Ahumada (2020)) or from the NASA/IPAC Extragalactic Database (NED)
(https://ned.ipac.caltech.edu/) and SIMBAD Astronomical Databases
(http://simbad.cfa.harvard.edu/simbad/). It is important to note that since
MBH values are derived from the scaling relations with line features, high
correlations are expected between MBH and Ldisk. In our present study, we have
ignored 4 blazars whose Ldisk values are unusually low ($<10^{10}$ erg s-1).

  * •

Blazar Emission Regions Catalog : Recently, Fan et al. (2023) Fan et al.
(2023) have studied the properties of emission zones using the SED fitting
results of more than 2700 blazars reported in Yang et al. (2023). 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 (Ue) and magnetic field energy
density (UB) from the SED fitting. For such blazars, we have collected the
luminosities of $\gamma$-ray at 1 GeV (log Lγ), X-ray at 1 keV (log LX),
optical at 2.43 $\times 10^{14}$ Hz (log LO) and radio at 1.4 GHz (log LR)
from Yang et al. (2022).

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) $\subseteq$ blazar emission
regions catalog (1791 sources) $\subseteq$ 4LAC-DR3 blazars catalog (1811
sources).

##  4 Results and Discussion

###  4.1 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\approx 0.31$ and is relatively
narrow with a standard deviation of $\sigma\approx 0.39$. Whereas, the
distribution of FSRQs peaks around a median of $z\approx 1.11$ and is broader
with $\sigma=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 Tolamatti et al.
(2023). 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 ($\sim 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 Abdo et al. (2009); Giommi (2012). 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 Bottcher and Dermer (2002); Ghisellini et al. (2011);
Ajello (2014). 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~{}\approx$ 3.59 and z$\approx$ 4.31,
respectively Paliya (2020b); Ackermann et al. (2017). They correspond to a
population of blazars existing just a few billion years after the Big Bang.

Figure 1: (Left) Histograms of the redshifts for FSRQs, BL Lacs and BCUs taken
from 4LAC-DR3 catalog Ajello (2022); Abdollahi (2022). (Right) The number
distributions of FSRQs, BL Lacs and BCUs in the high redshift ($z>0.4$) and
low redshift ($z\leq 0.4$) regimes.

###  4.2 $\gamma$-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$\sigma$ 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 $\gamma$-ray flux points (i.e., without taking EBL attenuation into
account) have been fitted with both Power-Law (PL)

$\frac{dN}{dE}=K{\Bigl{(}\frac{E}{Eo}\Bigr{)}}^{-\Gamma}$ (1)

and Log Parabola (LP)

$\frac{dN}{dE}=K{\Bigl{(}\frac{E}{Eo}\Bigr{)}}^{-\alpha-\beta\text{~{}ln}(E/Eo)}$
(2)

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,
$\Gamma$ is PL spectral index, $\alpha$ is LP index and $\beta$ is the LP
curvature parameter. For a fair comparison, we have calculated the modified LP
index ($\alpha^{\prime}$) 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, $\alpha^{\prime}$ 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 $\Gamma$, $\alpha^{\prime}$ and
$\beta$ for high and low-z blazars. The two sample Kolmogorov–Smirnov (K-S)
test results in _p_ -values of $\approx$0 for all three parameters. This
suggests that distributions of $\Gamma$, $\alpha^{\prime}$, and $\beta$
parameters are intrinsically different for high and low-z blazars. A Gaussian
fit to these distributions gives mean and standard deviation, $\mu\pm\sigma$,
values of 2.33 $\pm$ 0.25 (2.02 $\pm$ 0.23), 2.10 $\pm$ 0.39 (1.84 $\pm$
0.38), and 0.11 $\pm$ 0.07 (0.07 $\pm$ 0.7) for $\Gamma$, $\alpha^{\prime}$
and $\beta$, 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 $\Gamma$, $\alpha^{\prime}$ and
$\beta$ parameters are positively correlated with (log (1 + z)) with the
correlation coefficient ($\rho$) values of 0.6, 0.33, and 0.34, respectively,
and _p_ -values of $\approx$0\. Therefore, blazars at higher redshifts will
have softer $\gamma$-ray spectra (higher values of $\Gamma$ and
$\alpha^{\prime}$) and with more curvature (higher values of $\beta$) than
low-z blazars. Similar correlations have also been obtained between BL Lacs
and FSRQs in the literature Iyida et al. (2022); Ajello (2022). 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 $\gamma$-ray
photons by EBL via $\gamma$-$\gamma$ pair production during their propagation
from source to the observer Gould and Schréder (1966, 1967). The absorption
probability of HE photons traveling over cosmological distances is
characterized by optical depth ($\tau$) or opacity, which strongly depends on
$z$, $\gamma$-ray photon energy (E) and density of EBL photons Stecker (2006);
Singh et al. (2014). The EBL absorption effects start dominating at energies
above 30 GeV for redshifts beyond $z~{}>~{}0.4$ Tolamatti (2022). Therefore,
the intrinsic (emitted) $\gamma$-ray spectrum of the source gets modified and
the observed spectrum is given as

${\Bigl{(}\frac{dN}{dE}\Bigr{)}}_{obs}={\Bigl{(}\frac{dN}{dE}\Bigr{)}}_{int}e^{-\tau(E,z)}.$
(3)

For high-z blazars, $\tau$ will be large due to strong EBL absorption.
Therefore, the observed photon spectrum will be soft with more curvature Finke
et al. (2010); Dominguez et al. (2011). 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 Stecker (2006); Mankuzhiyil2023 .

Figure 2: (Left) Histograms of the Power-Law (PL) spectral index ($\Gamma$),
Log Parabola (LP) index ($\alpha^{\prime}$) and curvature parameter ($\beta$)
of blazars measured by the _Fermi_ -LAT. (Right) Variations of $\Gamma$,
$\alpha^{\prime}$ and $\beta$ parameters as a function of redshift (log (1 +
z)).

###  4.3 Luminosity Features

For luminosity studies, we have selected sources from the blazar emission
regions catalog. Out of 1791 sources in the catalog, 1789 have $\gamma$-ray
luminosities (Lγ) at 1 GeV, 1030 have X-ray luminosities (LX) at 1 keV, 1059
have optical luminosities (LO) at 2.43 $\times$ 1014 Hz and 1581 have radio
luminosities (LR) at 1.4 GHz. The left panels in Figure 3 show distributions
of log Lγ, log LX, log LO and log LR for high and low-z blazars Yang et al.
(2022). The right panels represent variations of these parameters as a
function of luminosity distance (DL) measured in Mega-parsec (Mpc). Note that,
DL is a function of z. A Gaussian fit to the distributions of log Lγ, log LX,
log LO, and log LR gives $\mu\pm\sigma$ values of 45.80 $\pm$ 0.89 (43.78
$\pm$ 0.86), 45.20 $\pm$ 0.72 (44.01 $\pm$ 0.88), 45.65 $\pm$ 0.67 (44.41
$\pm$ 0.52), and 42.91$\pm$ 0.89 (41.0 $\pm$ 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 LX, log LO, and log LR show strong positive
correlations with log DL having corresponding $\rho$ 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 DL. The slopes of linear
regression indicate that observed Lγ, LX, LO and LR vary as $\sim$
D${}_{L}^{2.27}$, D${}_{L}^{1.47}$, D${}_{L}^{1.51}$ and D${}_{L}^{2.14}$,
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.

###  4.4 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 ($\nu_{syn}$), inverse Compton (IC) peak
frequency ($\nu_{IC}$), energy fluxes at synchrotron ($\nu F_{\nu,syn}$) and
IC ($\nu F_{\nu,IC}$) peak frequencies and Compton-Dominance parameter
($CD=\nu F_{\nu,IC}/\nu F_{\nu,syn}$). For fair comparisons, we estimated
$\nu_{syn}$ and $\nu_{IC}$ in the source frame by multiplying them with a
factor (1 + z). The left panels in Figure 4 show the histogram distributions
of $\nu_{syn}$, $\nu_{IC}$, $\nu F_{\nu,syn}$, $\nu F_{\nu,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 $\approx$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 $\nu_{syn}$, $\nu_{IC}$, $\nu F_{\nu,syn}$, $\nu
F_{\nu,IC}$ and $CD$ gives $\mu\pm\sigma$ values of 13.42 $\pm$ 1.04 (15.32
$\pm$ 1.46), 22.00 $\pm$ 1.07 (23.53 $\pm$ 1.68), $-$11.90 $\pm$ 0.45
($-$11.65 $\pm$ 0.46), $-$11.42 $\pm$ 0.55 ($-$12.05 $\pm$ 0.60) and 0.48
$\pm$ 0.50 ($-$0.4 $\pm$ 0.40) for high-z (low-z) blazars, respectively. The
Pearson correlation analysis indicates that $\nu_{syn}$ and $\nu_{IC}$ are
anti-correlated with log (1 + z). Whereas, $\nu F_{\nu,IC}$ is positively
correlated with log (1 + z) with $\rho$=0.46. However, a weak anti-correlation
is observed between $\nu F_{\nu,syn}$ and log (1 + z) with $\rho$ = $-$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_ Ghisellini et
al. (2010); Ghisellini (2010). 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 Shakura and Sunyaev (1973). It is interesting to notice the
variation of $\nu_{syn}$ with respect to log (1 + z) in the top right panel of
Figure 4. The low-z blazars show large variation in $\nu_{syn}$ (from 1012 Hz
to 1019 Hz) making a vertical branch, whereas high-z blazars show little
variation in $\nu_{syn}$ (from 1012 to 1014 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 Prandini
(2022) 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 $\nu_{syn}$ and two branch feature in the
variation of $\nu_{syn}$ with log (1 + z) in Figure 4.

Figure 3: (Left) Histograms of the luminosities of $\gamma$-ray at 1 GeV (log
Lγ), X-ray at 1 keV (log LX), optical at 2.43 $\times 10^{14}$ Hz (log Lo) and
radio at 1.4 GHz (log LR) for sample high and low-z blazars. (Right)
Variations of log Lγ, log LX, log Lo and log LR as a function of luminosity
distance (log $D_{L}$).

Figure 4: (Left) Histograms of synchrotron peak frequency ($\nu_{syn}$), IC
peak frequency ($\nu_{IC}$), energy fluxes at synchrotron ($\nu F_{\nu,syn}$)
and IC ($\nu F_{\nu,IC}$) peak frequencies, and Compton dominance parameter
($CD$). (Right) Variations of log $\nu_{syn}$, log $\nu_{IC}$, log $\nu
F_{\nu,syn}$, log $\nu F_{\nu,IC}$ and log $CD$ as a function of log (1 + z).

###  4.5 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 (MBH), accretion
efficiency ($\eta$), 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 MBH values for high and low-z blazars.
In both the cases, MBH values range from $\sim$2.5 $\times$ 106 M⊙ to
$\sim$1.75 $\times$ 1010 M⊙. The average MBH values for high and low-z blazars
are similar with $\langle$log M${}_{BH,\text{Low-z}}\rangle\approx\langle$log
M${}_{BH,\text{High-z}}\rangle\approx 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 MBH 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.

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

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 MBH and log (1 + z)
as indicated by $\rho=$ 0.16 with _p_ -value of 2.88 $\times 10^{-7}$.
However, by considering only high-z blazars, the correlation becomes
relatively strong with $\rho=$ 0.41 and _p_ -value of 3.89$\times$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 MBH and log (1 + z) with $\rho$ = $-$0.23 and _p_
-value of 1.75 $\times$ 10 -5. Although the K-S test suggests that there is a
high probability that MBH 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, MBH 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 LDisk for blazars included in
this study. _p_ -value of 2.08 $\times 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
LDisk ($\langle$ log L${}_{Disk}\rangle$) 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
($\langle$L${}_{Disk,\text{High-z}}\rangle\approx 5.6\times
10^{45}\text{erg}s^{-1}$) have one order of magnitude higher accretion disk
luminosity than low-z blazars ($\langle$L${}_{Disk,\text{Low z}}\rangle\approx
3.0\times 10^{44}\text{erg}s^{-1}$). The upper right panel in Figure 6 shows
the variation of log Ldisk as a function of log (1 + z). The Pearson
correlation analysis suggests a strong positive correlation between the two
parameters with $\rho=$ 0.67 and _p_ -value of 1.28 $\times 10^{-82}$. Also,
the linear regression analysis of the scatter plot gives the slope of line $m$
= 3.3 $\pm$ 0.15. Thus, the observed LDisk broadly varies as $\sim$ (1 +
z)3.3.

Figure 6: (Upper Left) Histograms of LDisk for high and low redshift blazars.
(Upper Right) Scatter plot of log LDisk vs log (1 + z). (Lower Left)
Histograms of the Eddington ratios ($\lambda_{Edd}$) for high and low redshift
blazars. (Lower Right) Scatter plot of log $\lambda_{Edd}$ as function of log
(1 + z).

We also calculate the Eddington ratio ($\lambda_{Edd}=L_{Disk}/L_{Edd}$) using
the definition of Eddington luminosity as L${}_{Edd}=1.26\times
10^{38}(M_{BH}/M_{\odot})\text{ erg }s^{-1}$. The corresponding distributions
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
$\lambda_{Edd}$ distributions for high and low-z blazars are intrinsically
different as indicated by the K-S test _p_ -value of 8.74 $\times 10^{-13}$.
The Gaussian function fitting to the distributions with mean and standard
deviation ($\mu,~{}\sigma$) 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 $\langle$
$\lambda_{Edd,\text{High-z}}\rangle=0.12$ compared to the low-z blazars with
$\langle$ $\lambda_{Edd,\text{Low-z}}\rangle=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
$\lambda_{Edd}$ and log (1 + z) with $\rho$ = 0.37 and _p_ -value of
$2.46\times 10^{-4}$. The linear regression line in the scatter plot with
m$=0.37\pm 0.10$ suggests that observed $\lambda_{Edd}$ varies as $\sim$(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 Bottcher and Dermer (2002); Ghisellini et al. (2011);
Ajello (2014). 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 $\lambda_{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 Narayan et al. (1997).
Alternatively, the observed features can also be attributed to the selection
bias effects in observations as it is difficult to detect blazars having low
$\lambda_{Edd}$ and low LDisk at high redshifts.

###  4.6 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 Volonteri (2010); Volonteri et al.
(2015). Following the detailed calculations given in Shapiro (2005);
Belladitta (2022), we estimate the mass of black hole seed (MBH,seed) required
to produce currently observed SMBHs at a given $z$, assuming that it evolves
with constant $\lambda_{Edd}$ and constant accretion efficiency defined as
$\eta=L_{Disk}/\dot{M}_{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 MBH growth is given by

$\frac{dM_{BH}}{dt}=(1-\eta)\dot{M}_{BH}=\frac{(1-\eta)\lambda_{Edd}L_{Edd}}{\eta
c^{2}}=\frac{(1-\eta)\lambda_{Edd}}{\eta}\frac{M_{BH}}{t_{s}}$ (4)

where ts is the characteristic accretion timescale (also known as Salpeter
time) independent of MBH. It is given as

$t_{s}=\frac{\sigma_{T}c}{4\pi Gm_{p}}\approx 0.45\text{ Gyr}$ (5)

Therefore, required black hole seed mass (MBH,seed) for the formation of SMBH
of mass MBH is given by

$M_{BH,seed}=M_{BH}\times\exp{\Bigl{(}\frac{-t}{\tau}\Bigr{)}}$ (6)

where $t$ is the total time available during which the seed black hole
accretes and grows, $\tau$ is the characteristic $e-folding$ timescale defined
as

$\tau=\frac{0.45\eta}{\lambda_{Edd}(1-\eta)}\text{ Gyr}$ (7)

The accretion efficiency parameter $\eta$, 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 Thorne (1974); Shapiro (2005). 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 MBH,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) Bromm and Larson (2004); Bromm and Yoshida (2011). In our
calculations, we use currently observed $\lambda_{Edd}$ throughout their
evolution to see if such a scenario explains the observed SMBHs. The
calculations are performed for two constant accretion efficiencies of $\eta$ =
0.1 and $\eta$ = 0.3 Shapiro (2005). For calculating t in Equation (6), we use
the standard flat $\Lambda$CDM model of the Universe with parameters H0 = 70
km s-1 Mpc-1, $\Omega_{M}$ = 0.30 and $\Omega_{\Lambda}$ = 0.70. The final
results depend on the assumptions made regarding various parameters like
redshift of the seed formation, $\eta$ and $\lambda_{Edd}$.

Figure 7 shows the variations of the estimated values of MBH,seed as a
function of log (1 + z) for $\eta=0.1$ (left panel) and $\eta=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_{\odot}$ (cyan color)
Abel et al. (2000); Bond et al. (1984), (ii) seeds from stellar dynamical
processes: $10^{3}-10^{4}M_{\odot}$ (gray color) Devecchi and Volonteri
(2009); Lupi et al. (2014); Alexander and Natarajan (2014); Volonteri et al.
(2016); Lupi et al. (2016), and (iii) direct collapse black holes:
$10^{4}$–$10^{6}M_{\odot}$ (green color) Haehnelt and Rees (1993); Begelman et
al. (2006). 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.

Figure 7:  Variations of the estimated seed MBH,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 $\lambda_{Edd,obs}$ and with
$\eta=$ 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}\text{
M}_{\odot}$; cyan color), seed black holes from stellar dynamical processes
($\geq 10^{3}\text{ M}_{\odot}$; gray color) and black holes from direct
collapse ($10^{4}$–$10^{6}\text{ M}_{\odot}$; green color).

The left panel in Figure 7 indicates that, for $\eta=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
$\sim$69% when we consider $\eta=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
Shapiro (2005). Alternatively, these blazars must accrete the matter with a
higher Eddington ratio than the currently observed $\lambda_{Edd}$ for a
significant time during their evolution to explain the growth of observed
SMBHs. Thus, the higher the $\eta$, the higher the $\lambda_{Edd}$ we need to
explain the growth of the SMBHs. A major fraction of the SMBHs in BL Lacs
($\sim$75% and $\sim$92% of the total BL Lacs corresponding to $\eta=$ 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 ($\lambda_{Edd}$). Therefore, these
objects must have accreted the matter with higher $\lambda_{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
$\lambda_{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
$\lambda_{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)
Onoue et al. (2019) and PSO J006+39 (z = 6.621) Tang et al. (2019). Although
there are possibilities for super Eddington accretion episodes in gas-rich or
dusty strong wind environments near the central engines Kim et al. (2015);
Kubota and Done (2019), 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 Blandford & Znajek (1977). Therefore, for a given
accretion efficiency, there will be fewer luminous disks and hence lower
($\lambda_{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 Jolley and Kuncic (2008); Jolley
et al. (2009); Ghisellini et al. (2010). 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 Alexander and Hickox
(2012); Feng et al. (2014). Similar calculations and conclusions have also
been reported in a recent study Belladitta (2022). 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
$\Lambda_{Edd}$, LDisk, and MBH 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 Shen et al. (2011). Therefore, the estimation of
MBH from different approaches like reverberation mapping, etc., are crucial
Peterson (1993); Singh et al. (2022).

###  4.7 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 $\gamma_{p}$. We use log B and log
$\gamma_{p}$ parameters reported in the blazar emission regions catalog for
this study. For some sources, ranges of log B and log $\gamma_{p}$ are
reported in the catalog, and we take the midpoint of the range as the value
for log B and log $\gamma_{p}$ for them. The left panels in Figure 8 show the
histograms of log B and log $\gamma_{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 $\gamma_{p}$ parameters are close to 0. This suggests that the
distributions of B and $\gamma_{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 $\mu\pm\sigma$ of 0.65 $\pm$
1.68 (B = 4.46 G) and 0.04 $\pm$ 1.53 (B = 1.09 G) for high and low-z blazars,
respectively. Also, there exists a mild positive correlation between B values
and log (1 + z) with $\rho$ = 0.26 and the P value of 1.24$\times 10^{-28}$.
Therefore, high-z blazars require emission zones with higher magnetic fields
compared to low-z blazars. The histogram distribution of log $\gamma_{p}$ for
high-z blazars shows a clear double hump structure with peaks around log
$\gamma_{p}$ values of 2 and 4. Similarly, the scatter plot of log
$\gamma_{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 $\gamma_{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 $\gamma_{p}$,
we calculated the arithmetic mean ($\mu^{\prime}$) and standard deviation
($\sigma^{\prime}$) for the distributions. The $\mu^{\prime}$ values for high
and low-z blazars are found to be 2.90 and 3.79 corresponding to peak Lorentz
factors of 7.94 $\times 10^{2}$ and 6.17 $\times 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 $\gamma_{p}$ and
log (1 + z) with $\rho$ = $-$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 $\gamma_{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 $\gamma_{p}$ in the absence of target photon fields.
Therefore, their SED peak lies at higher frequencies with very low Compton
dominance Bottcher and Dermer (2002); Ghisellini et al. (2011); Paliya (2021);
Fan et al. (2023). The variation of log $\gamma_{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 $\gamma_{p}$. Therefore, this can
also contribute to the observed soft photon spectral indices and higher
curvatures for high-z blazars.

Figure 8: (Left) Histogramsof log B and log $\gamma_{p}$ for high and low
redshift blazars. (Right) Variation of log B and log $\gamma_{p}$ values with
log (1 + z).

##  5 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.

  * •

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 ($\Gamma$), LP slope at 1 GeV/(1 + z)
($\alpha^{\prime}$) and LP curvature ($\beta$) roughly follow Gaussian
distributions with $\mu\pm\sigma$ values of $\Gamma=2.33\pm 0.25(2.02\pm
0.23)$, $\alpha^{\prime}=2.10\pm 0.39(1.84\pm 0.38)$ and $\beta=0.11\pm
0.07(0.07\pm 0.07)$ for high (low) redshift blazars. Also, $\Gamma$,
$\alpha^{\prime}$ and $\beta$ 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 $\gamma$ photons, leading to softer $\gamma$-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 $\gamma$-ray spectral features.

  * •

Observationally, high-z blazars are the most luminous objects in the
$\gamma$-ray Universe. The mean $\gamma-$ray, X-ray, optical and radio
luminosities (in $erg~{}s^{-1}$) of high (low) redshift blazars are estimated
as $\langle L_{\gamma}\rangle=6.31\times 10^{45}(6.03\times 10^{43})$,
$\langle L_{X}\rangle=1.58\times 10^{45}(1.02\times 10^{44})$, $\langle
L_{O}\rangle=4.46\times 10^{45}(2.57\times 10^{44})$ and $\langle
L_{R}\rangle=8.13\times 10^{42}(1.0\times 10^{41})$, respectively. These
luminosities have a strong positive correlation with the luminosity distance
(DL) with functional dependencies of $L_{\gamma}\propto D_{L}^{2.27}$,
$L_{X}\propto D_{L}^{1.47}$, $L_{O}\propto D_{L}^{1.51}$ and $L_{R}\propto
D_{L}^{2.14}$.

  * •

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 ($\nu_{syn}$), rest
frame IC peak frequency ($\nu_{IC}$), and Compton dominance (CD) roughly
follow a Gaussian profile with $\mu\pm\sigma$ values of log
$\nu_{syn}=13.42\pm 1.04(15.32\pm 1.46)$, log $\nu_{IC}=22.00\pm 1.07(23.53\pm
1.68)$ and log $CD=0.48\pm 0.50\pm(-0.4\pm 0.40)$ for high (low) redshift
blazars. Also, the Pearson correlation study indicates that log $\nu_{syn}$
and log $\nu_{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
$\langle\text{log }M_{BH}\rangle=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 $\langle
L_{Disk}\rangle=5.62\times 10^{45}\text{erg s}^{-1}(2.4\times 10^{44}\text{erg
s}^{-1})$, $\langle\lambda_{Edd}\rangle=0.12(0.03)$ with both the parameters
showing a positive correlation with log (1 + z).

  * •

We also estimated the masses of black hole seed (MBH,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 $\eta=0.1(0.3)$ and currently observed constant
$\lambda_{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 $\lambda_{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 $\eta$ = 0.1 and 0.3, respectively) are
required to have higher $\lambda_{Edd}$ than currently observed
$\lambda_{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 $\langle\text{B
}\rangle=4.46G(1.09)$ G and electron energy distributions with $\gamma_{p}$
values $\langle\gamma_{p}\rangle=7.94\times 10^{2}$ (6.17 $\times 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
$\gamma_{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 $\gamma_{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 $\gamma$-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 Merloni (2012)), the Large High Altitude Air Shower Observatory
(LHAASO Bai (2019)), the Southern Wide-field Gamma-ray Observatory (SWGO
Albert (2019)), and the Cherenkov Telescope Array (CTA CTA (2021)) are
expected to detect large number of high redshift blazars.

A.T.—Conceptualization,data curation, software and original draft preparation;
K.K.S.—Conceptualization, review and editing; K.K.Y.—review and editing. All
authors have read and agreed to the published version of the manuscript. This
research received no external funding. The data set used in this work is
derived from the literature and can be shared on reasonable request.

###### Acknowledgements.

Authors thank all three anonymous reviewers for their important comments and
suggestions. We also thank Hubing Xiao A. for sharing the blazar emission
region data used in this study. The authors declare no conflicts of interest.
References

## References

  * Urry & Padovani (1995) Urry, C.M.; Padovani, P. Unified Schemes for Radio-Loud Active Galactic Nuclei. _Publ. Astron. Soc. Pac._ 1995, _107_ , 803.
  * Padovani et al. (2017) Padovani, P.; Alexander, D.M.; Assef, R.J.; De Marco, B.; et al. Active galactic nuclei: What’s in a name? _Annu. Rev. Astron. Astrophys._ 2017, _25_ , 2.
  * Blandford & Znajek (1977) Blandford, R.D.; Znajek, R.L. Electromagnetic extraction of energy from Kerr black holes. _Mon. Not. R. Astron. Soc._ 1977, _179_ , 433–456.
  * Blandford & Payne (1982) Blandford, R.D.; Payne, D.G. Hydromagnetic flows from accretion disks and the production of radio jets. _Mon. Not. R. Astron. Soc._ 1982, _199_ , 883–903.
  * Maraschi & Tavecchio (2003) Maraschi, L.; Tavecchio, F. The Jet-Disk Connection and Blazar Unification. _Astrophys. J._ 2003, _594_ , 667–675.
  * Ghisellini et al. (2011) Ghisellini, G.; Tavecchio, F.; Foschini, L.; Ghirlanda, G. The transition between BL Lac objects and flat spectrum radio quasars. _Mon. Not. R. Astron. Soc._ 2011, _414_ , 2674–2689.
  * Chen et al. (2015) Chen, Y.Y.; Zhang, X.; Zhang, H.J.; Yu, X.L. Core-dominance parameter, black hole mass and jet-disc connection for Fermi blazars. _Mon. Not. R. Astron. Soc._ 2015, _451_ , 4193–4206.
  * Zhang et al. (2022) Zhang, L.; Liu, Y.; Fan, J. Classification and Jet Power of Fermi Blazars. _Astrophys. J._ 2003, _594_ , 667–675.
  * Meier (1999) Meier, D.L. A Magnetically Switched, Rotating Black Hole Model for the Production of Extragalactic Radio Jets and the Fanaroff and Riley Class Division. _Astrophys. J._ 1999, _522_ , 753–766.
  * Komissarov (2021) Komissarov, S.; Porth, O. Numerical simulations of jets. _New Astron. Rev._ 2021, _92_ , 101610.
  * Wagner (1995) Wagner, S.J.; Witzel, A. Intraday Variability In Quasars and BL Lac Objects. _Annu. Rev. Astron. Astrophys._ 1995, _33_ , 163–198.
  * Giommi et al. (2021) Giommi, P.; Perri, M.; Capalbi, M.; D’Elia, V.; Barres de Almeida, U.; et al. X-ray spectra, light curves and SEDs of blazars frequently observed by Swift. _Mon. Not. R. Astron. Soc._ 2021, _507_ , 5690–5702.
  * Lister et al. (2019) Lister, M.L.; Homan, D.C.; Hovatta, T.; Kellermann, K.I.; Kiehlmann, S.; et al. MOJAVE. XVII. Jet Kinematics and Parent Population Properties of Relativistically Beamed Radio-loud Blazars. _Astrophys. J._ 2019, _874_ , 43.
  * Weaver (2020) Weaver, Z.R.; Williamson, K.E.; Jorstad, S.G.; Marscher, A.P.; et al. Multiwavelength Variability of BL Lacertae Measured with High Time Resolution. _Astrophys. J._ 2020, _900_ , 137.
  * Shukla (2020) Shukla, A.; Mannheim, K. Gamma-ray flares from relativistic magnetic reconnection in the jet of the quasar 3C 279. _Nat. Commun._ 2020, _11_ , 4176.
  * Singh (2020a) Singh, K.K.; Meintjes, P.J. Characterization of variability in blazar light curves. _Astron. Not._ 2020, _341_ , 713–725.
  * Singh (2020b) Singh, K.K.; Meintjes, P.J.; Bisschoff, B.; Ramamonjisoa, F.A.; van Soelen, B. Gamma-ray and optical properties of the flat spectrum radio quasar 3C 279 flare in June 2015. _J. High Energy Astrophys._ 2020, _26_ , 65–76.
  * Tolamatti et al. (2022) Tolamatti, A.; Ghosal, B.; Singh, K.K.; Bhattacharyya, S.; Bhatt, N.; et al. Long-term multi-wavelength study of temporal and spectral properties of 3C 279. _Astropart. Phys._ 2022, _139_ , 102687.
  * Singh et al. (2022) Singh, K.K.; Tolamatti, A.; Godiyal, S.; Pathania, A.; Yadav, K.K. A Machine Learning Approach for Predicting Black Hole Mass in Blazars Using Broadband Emission Model Parameters. _Universe_ 2022, _8_ , 539.
  * Marscher (1985) Marscher, A.P.; Gear, W.K. Models for high-frequency radio outbursts in extragalactic sources, with application to the early 1983 millimeter-to-infrared flare of 3C 273. _Astrophys. J._ 1985, _298_ , 114–127.
  * Sikora (1994) Sikora, M.; Begelman, M.C.; Rees, M.J. Comptonization of Diffuse Ambient Radiation by a Relativistic Jet: The Source of Gamma Rays from Blazars? _Astrophys. J._ 1994, _421_ , 153.
  * Ghisellini (2009) Ghisellini, G.; Tavecchio, F. Canonical high-power blazars. _Mon. Not. R. Astron. Soc._ 2009, _397_ , 985–1002.
  * Aharonian (2000) Aharonian, F.A. TeV gamma rays from BL Lac objects due to synchrotron radiation of extremely high energy protons. _New Astron._ 2000, _5_ , 377–395.
  * Bottcher (2013) Bottcher, M.; Reimer, A.; Sweeney, K.; Prakash, A. Leptonic and Hadronic Modeling of Fermi-detected Blazars. _Astrophys. J._ 2013, _768_ , 54.
  * Petropoulou (2015) Petropoulou, M.; Dimitrakoudis, S. Constraints of flat spectrum radio quasars in the hadronic model: The case of 3C 273. _Mon. Not. R. Astron. Soc._ 2015, _452_ , 1303–1315.
  * Stickel (1991) Stickel, M.; Padovani, P.; Urry, C.M.; Fried, J.W.; Kuehr, H. The Complete Sample of 1 Jansky BL Lacertae Objects. I. Summary Properties. _Astrophys. J._ 1991, _374_ , 431.
  * Abdo (2010) Abdo, A.A.; Ackermann, M.; Agudo, I.; et al. The Spectral Energy Distribution of Fermi Bright Blazars. _Astrophys. J._ 2010,_716_ , 30–70.
  * Padovani (2012) Padovani, P.; Giommi, P.; Rau, A. The discovery of high-power high synchrotron peak blazars. _Mon. Not. R. Astron. Soc._ 2012, _422_ , L48–L52.
  * Giommi (2012) Giommi, P.; Padovani, P.; Polenta, G.; Turriziani, S.; et al. A simplified view of blazars: Clearing the fog around long-standing selection effects._Mon. Not. R. Astron. Soc._ 2012, _420_ , 2899–2911.
  * Fossati (1998) Fossati, G.; Maraschi, L.; Celotti, A.; Comastri, A.; Ghisellini, G. A unifying view of the spectral energy distributions of blazars. _Mon. Not. R. Astron. Soc._ 1998, _299_ , 433–448.
  * Prandini (2022) Prandini, E.; Ghisellini, G. The Blazar Sequence and Its Physical Understanding. _Galaxies_ 2022, _10_ , 35.
  * Ghisellini (2016) Ghisellini, G. The Blazar Sequence 2.0. _Galaxies_ 2016, _4_ , 36.
  * Ghisellini (2017) Ghisellini, G.; Righi, C.; Costamante, L.; Tavecchio, F. The Fermi blazar sequence. _Mon. Not. R. Astron. Soc._ 2017, _469_ , 255–266.
  * Singh (2019) Singh, K.K.; Meintjes, P.J.; Ramamonjisoa, F.A.; Tolamatti, A. Extremely High energy peaked BL Lac nature of the TeV blazar Mrk 501. _New Astron._ 2019, _73_ , 101278.
  * Ghisellini (1998) Ghisellini, G.; Celotti, A.; Fossati, G.; Maraschi, L.; Comastri, A. A theoretical unifying scheme for gamma-ray bright blazars. _Mon. Not. R. Astron. Soc._ 1998, _301_ , 451–468.
  * Fan (2017) Fan, J.H.; Yang, J.H.; Xiao, H.B.; et al. Intrinsic Correlations for Flaring Blazars Detected by Fermi. _Astrophys. J. Lett._ 2017, _835_ , L38.
  * Paliya (2020) Paliya, V.S.; Ajello, M.; Cao, H.M.; et al. Blazars at the Cosmic Dawn. _Astrophys. J._ 2020, _897_ , 177.
  * Sahakyan (2020) Sahakyan, N.; Israyelyan, D.; Harutyunyan, G.; Khachatryan, M.; Gasparyan, S. Multiwavelength study of high-redshift blazars. _Mon. Not. R. Astron. Soc._ 2020, _498_ , 2594–2613.
  * Stecker (2006) Stecker, F.W.; Malkan, M.A.; Scully, S.T. Intergalactic Photon Spectra from the Far-IR to the UV Lyman Limit for $0<z<6$ and the Optical Depth of the Universe to High-Energy Gamma Rays. _Astrophys. J._ 2006, _648_ , 774–783.
  * Singh (2020) Singh, K.K.; Meintjes, P.J. Extragalactic background light models and GeV-TeV observation of blazars. _NRIAG J. Astron. Geophys._ 2020, _9_ , 309–320.
  * Tolamatti (2022) Tolamatti, A.; Singh, K.K.; Yadav, K.K. Feasibility study of observing $\gamma$-ray emission from high redshift blazars using the MACE telescope. _J. Astrophys. Astr._ 2022, _43_ , 49.
  * Ajello (2014) Ajello, M.; Romani, R.W.; Gasparrini, D.; et al. The Cosmic Evolution of Fermi BL Lacertae Objects. _Astrophys. J._ 2014, _780_ , 73.
  * Ghisellini (2013) Ghisellini, G., High redshift blazars. _Mem. Soc. Astron. Ital._ 2013, _84_ , 719.
  * Ackermann et al. (2015) Ackermann, M.; Ajello, M.; Atwood, W. B.; et al. The Third Catalog of Active Galactic Nuclei Detected by the Fermi Large Area Telescope. _Astrophys. J._ 2015, _810_ , 14.
  * Inayoshi (2020) Inayoshi, K.; Visbal, E.; Haiman, Z. The Assembly of the First Massive Black Holes. _Annu. Rev. Astron. Astrophys._ 2020, _58_ , 27–97.
  * Belladitta (2022) Belladitta, S.; Caccianiga, A.; Diana, A.; et al. Central engine of the highest redshift blazar. _Astron. Astrophys._ 2022, _660_ , A74.
  * Johnson (2016) Johnson, J.L.; Haardt, F. The Early Growth of the First Black Holes. _Publ. Astron. Soc. Astra._ 2016, _33_ , e007.
  * Belladitta (2020) Belladitta, S.; Moretti, A.; Caccianiga, A.; et al. The first blazar observed at $z>6$. _Astron. Astrophys._ 2020, _635_ , L7.
  * Ajello (2009) Ajello, M.; Costamante, L.; Sambruna, R. M.; et al. The Evolution of Swift/BAT Blazars and the Origin of the MeV Background. _Astrophys. J._ 2009, _699_ , 603–625.
  * Ghisellini (2010) Ghisellini, G.; Della Ceca, R.; Volonteri, M.; et al. Chasing the heaviest black holes of jetted active galactic nuclei. _Mon. Not. R. Astron. Soc._ 2010, _405_ , 387–400.
  * Volonteri (2010) Volonteri, M., Formation of supermassive black holes. _Annu. Rev. Astron. Astrophys._ 2010, _18_ , 279–315.
  * Volonteri (2011) Volonteri, M.; Haardt, F.; Ghisellini, G.; Della Ceca, R. Blazars in the early Universe. _Mon. Not. R. Astron. Soc._ 2011, _416_ , 216–224.
  * Sbarrato (2015) Sbarrato, T.; Ghisellini, G.; Tagliaferri, G.; et al. Blazar candidates beyond redshift 4 observed by Swift. _Mon. Not. R. Astron. Soc._ 2015, _446_ , 2483–2489.
  * Ajello (2022) Ajello, M.; Baldini, L.; Ballet, J.; et al. The Fourth Catalog of Active Galactic Nuclei Detected by the Fermi Large Area Telescope: Data Release 3. _Astrophys. J. Suppl._ 2022, _263_ , 24.
  * Abdollahi (2022) Abdollahi, S.; Acero, F.; Baldini, L.; et al. Incremental Fermi Large Area Telescope Fourth Source Catalog. _Astrophys. J. Suppl._ 2022, _260_ , 53.
  * Paliya (2021) Paliya, V.S.; Domínguez, A.; Ajello, M.; Olmo-García, A.; Hartmann, D. The Central Engines of Fermi Blazars. _Astrophys. J. Suppl._ 2021, _253_ , 46.
  * Ahumada (2020) Ahumada, R.; Allende, P.C.; Almeida, A.; et al. The 16th Data Release of the Sloan Digital Sky Surveys: First Release from the APOGEE-2 Southern Survey and Full Release of eBOSS Spectra. _Astrophys. J. Suppl._ 2020, _249_ , 3.
  * Fan et al. (2023) Fan, J.; Xiao, H.; Yang, W.; Zhang, L.; et al. Characterizing the Emission Region Properties of Blazars, _Astrophys. J. Suppl._ 2023, _268_ , 23.
  * Yang et al. (2023) Yang, J.; Fan, J.; Liu, Y.; et al. Estimation of inverse Compton peak frequency for 4FGL Blazars. _Sci. China Phys. Mech. Astron._ 2023, _66_ , 249511.
  * Yang et al. (2022) Yang, J.H.; Fan, J.H.; Liu, Y.et al. The Spectral Energy Distributions for 4FGL Blazars. _Astrophys. J. Suppl._ 2022, _262_ , 18.
  * Tolamatti et al. (2023) Tolamatti, A.; Singh, K.K.; Yadav, K.K.; Classification of the Fermi-LAT blazar candidates of uncertain type using extreme gradient boosting. _Mon. Not. R. Astron. Soc._ 2023, _523_ , 5341–5352.
  * Abdo et al. (2009) Abdo, A.A.; Ackermann, M.; Ajello, M.; et al. Bright Active Galactic Nuclei Source List from the First Three Months of the Fermi Large Area Telescope All-Sky Survey. _Astrophys. J._ 2009, _700_ , 597–622.
  * Bottcher and Dermer (2002) Bottcher, M.; Dermer, C.D.; An Evolutionary Scenario for Blazar Unification. _Astrophys. J._ 2002, _564_ , 86–91.
  * Paliya (2020b) Paliya, V.S.; Domínguez, A.; Cabello, C.; et al. The First Gamma-Ray Emitting BL Lacertae Object at the Cosmic Dawn. _Astrophys. J. Lett._ 2020, _903_ , L8.
  * Ackermann et al. (2017) Ackermann, M.; Ajello, M.; Baldini, L.et al. Gamma-Ray Blazars within the First 2 Billion Years. _Astrophys. J. Lett._ 2017, _837_ , L5.
  * Iyida et al. (2022) Iyida, E.U.; Odo, F.C.; Chukwude, A.E.; Ubachukwu, A.A. Multi-frequency study of spectral indices of BL Lacertae objects and flat-spectrum radio quasars._New Astron._ 2022, _90_ , 101666.
  * Gould and Schréder (1966) Gould, R.J.; & Schréder, G. Opacity of the Universe to High-Energy Photons. _Phys. Rev. Lett._ 1966, _16_ , 252–254.
  * Gould and Schréder (1967) Gould, R.J.; Schréder, G.P. Opacity of the Universe to High-Energy Photons. _Phys. Rev._ 1967, _155_ , 1408–1411.
  * Singh et al. (2014) Singh, K.K.; Sahayanathan, S.; Tickoo, A. K.; Bhatt, N. Intrinsic VHE gamma-ray spectra of blazars as a probe for extragalactic background light. _New Astron._ 2014, _27_ , 34–40.
  * Finke et al. (2010) Finke, J.D.; Razzaque, S.; Dermer, C.D. Modeling the Extragalactic Background Light from Stars and Dust. _Astrophys. J._ 2010, _712_ , 238–249.
  * Dominguez et al. (2011) Dominguez, A.; Primack, J.R.; Rosario, D.J.; et al. Extragalactic background light inferred from AEGIS galaxy-SED-type fractions. _Mon. Not. R. Astron. Soc._ 2011, _410_ , 2556–2578.
  * (72) Mankuzhiyil, N; Persic, M; Franceschini, A. Prospects of additional contribution at Optical-NIR band of EBL in the light of VHE spectra, _Mon. Not. R. Astron. Soc._ 2023, _524_ , 133–142
  * Ghisellini et al. (2010) Ghisellini, G.; Tavecchio, F.; Foschini, L.; et al. General physical properties of bright Fermi blazars. _Mon. Not. R. Astron. Soc._ 2010, _402_ , 497–518.
  * Shakura and Sunyaev (1973) Shakura, N.I.; Sunyaev, R.A . Black holes in binary systems. Observational appearance. _Astron. Astrophys_ , 1973, _24_ , 337–355.
  * Narayan et al. (1997) Narayan, R.; Garcia, M.R.; McClintock, J.E. Advection-dominated Accretion and Black Hole Event Horizons. _Astrophys. J. Lett._ 1997, _478_ , L79–L82.
  * Volonteri et al. (2015) Volonteri, M.; Silk, J.; Dubus, G. The Case for Supercritical Accretion onto Massive Black Holes at High Redshift. _Astrophys. J._ 2015, _804_ , 148.
  * Shapiro (2005) Shapiro, S.L. Spin, Accretion, and the Cosmological Growth of Supermassive Black Holes. _Astrophys. J._ 2005, _620_ , 59–68.
  * Thorne (1974) Thorne, K.S.; Disk-Accretion onto a Black Hole. II. Evolution of the Hole. _Astrophys. J._ 1974, _191_ , 507–520.
  * Bromm and Larson (2004) Bromm, V.; Larson, R.B. The First Stars. _Annu. Rev. Astron. Astrophys._ 2004, _42_ , 79–118.
  * Bromm and Yoshida (2011) Bromm, V.; Yoshida N. The First Galaxies. _Annu. Rev. Astron. Astrophys._ 2011, _49_ , 373–407.
  * Abel et al. (2000) Abel, T.; Bryan, G.L.; Norman, M.L. The Formation and Fragmentation of Primordial Molecular Clouds. _Astrophys. J._ 2000, _540_ , 39–44.
  * Bond et al. (1984) Bond, J.R.; Arnett, W.D.; Carr, B.J. The evolution and fate of Very Massive Objects. _Astrophys. J._ 1984, _280_ , 825–847.
  * Devecchi and Volonteri (2009) Devecchi, B.;Volonteri, M. Formation of the First Nuclear Clusters and Massive Black Holes at High Redshift. _Astrophys. J._ 2009, _694_ , 302–313.
  * Lupi et al. (2014) Lupi, A.; Colpi, M.; Devecchi, B.; Galanti, G.; Volonteri, M. Constraining the high-redshift formation of black hole seeds in nuclear star clusters with gas inflows. _Mon. Not. R. Astron. Soc._ 2014, _442_ , 3616–3626.
  * Alexander and Natarajan (2014) Alexander, T.; Natarajan, P. Rapid growth of seed black holes in the early universe by supra-exponential accretion. _Science_ 2014, _345_ , 1330–1333.
  * Volonteri et al. (2016) Volonteri, M.; Dubois, Y.; Pichon, C.; Devriendt, J. The cosmic evolution of massive black holes in the Horizon-AGN simulation. _Mon. Not. R. Astron. Soc._ 2016, _460_ , 2979–2996.
  * Lupi et al. (2016) Lupi, A.; Haardt, F.; Dotti, M.; et al. Growing massive black holes through supercritical accretion of stellar-mass seeds. _Mon. Not. R. Astron. Soc._ 2016, _456_ , 2993–3003.
  * Haehnelt and Rees (1993) Haehnelt, M.G.; Rees, M.J. The formation of nuclei in newly formed galaxies and the evolution of the quasar population. _Mon. Not. R. Astron. Soc._ 1993, _263_ , 168–178.
  * Begelman et al. (2006) Begelman, M.C.; Volonteri, M.; Rees, M.J. Formation of supermassive black holes by direct collapse in pre-galactic haloes. _Mon. Not. R. Astron. Soc._ 2006, 370, 289–298.
  * Onoue et al. (2019) Onoue, M.; Kashikawa, N; Matsuoka, Y.; et al. Subaru High-z Exploration of Low-luminosity Quasars (SHELLQs). VI. Black Hole Mass Measurements of Six Quasars at 6.1 $\leq$ z $\leq$ 6.7. _Astrophys. J._ 2019, _880_ , 77.
  * Tang et al. (2019) Tang, J.-J.; Goto, T; Ohyama, Y.; et al. Rapid black hole growth at the dawn of the Universe: A super-Eddington quasar at z = 6.6. _Mon. Not. R. Astron. Soc._ 2019, _484_ , 2575–2586.
  * Kim et al. (2015) Kim, Y.; Im, M; Jeon, Y.; et al. Discovery of a Faint Quasar at z $\sim$ 6 and Implications for Cosmic Reionization. _Astrophys. J. Lett._ , 2015, _813_ , L35.
  * Kubota and Done (2019) Kubota, A.; Done, C. Modelling the spectral energy distribution of super-Eddington quasars. _Mon. Not. R. Astron. Soc._ 2019, _489_ , 524–533.
  * Jolley and Kuncic (2008) Jolley, E.J.D.; Kuncic, Z. Jet-enhanced accretion growth of supermassive black holes. _Mon. Not. R. Astron. Soc._ 2008, _386_ , 989–994.
  * Jolley et al. (2009) Jolley, E.J.D.; Kuncic, Z.; Bicknell, G.V.; Wagner, S. Accretion discs in blazars. _Mon. Not. R. Astron. Soc._ 2009, _400_ , 1521–1526.
  * Alexander and Hickox (2012) Alexander, D.M.; Hickox, R.C. What drives the growth of black holes? _New Astron. Rev_ 2012, _56_ , 93–121.
  * Feng et al. (2014) Feng, Y.; Di Matteo, T.; Croft, R.; Khandai, N. High-redshift supermassive black holes: Accretion through cold flows. _Mon. Not. R. Astron. Soc._ 2014, _440_ , 1865–1879.
  * Shen et al. (2011) Shen Y.; Richards G.T.; Strauss M.A.; Hall, P.B.; Schneider, D.P.; Snedden, S.; Bizyaev, D.; et al. A Catalog of Quasar Properties from Sloan Digital Sky Survey Data Release 7. _Astrophys. J. Suppl. Ser._ 2011, _194_ , 45. https://doi.org/10.1088/0067-0049/194/2/45
  * Peterson (1993) Peterson, B.M. Reverberation Mapping of Active Galactic Nuclei. _Publ. Astron. Soc. Pac._ 1993, _105_ , 247.
  * Merloni (2012) Merloni, A.; Predehl, P.; Becker, W.; Böhringer, H.; Boller, T.; Brunner, H.; Brusa, M.; Dennerl, K.; Freyberg, M.; Friedrich, P.; et al. eROSITA Science Book: Mapping the Structure of the Energetic Universe. _arXiv_ 2012, arXiv:1209.3114.
  * Bai (2019) Bai, X.; Bi, B.Y.; Bi, X.J.; Cao, Z.; Chen, S.Z.; Chen, Y.; Chiavassa, A.; Cui, X.H.; Dai, Z.G.; della Volpe, D.; et al. The Large High Altitude Air Shower Observatory (LHAASO) Science White Paper—2021 Edition. _arXiv_ 2019, arXiv:1905.02773.
  * Albert (2019) Albert, A.; Alfaro, R.; Ashkar, H.; Alvarez, C.; Álvarez, J.; Arteaga-Velázquez, J.C.; Solares, H.A.; Arceo, R.; Bellido, J.A.; BenZvi, S.; et al. Science Case for a Wide Field-of-View Very-High-Energy Gamma-Ray Observatory in the Southern Hemisphere. _arXiv_ 2019, arXiv:1902.08429.
  * CTA (2021) Cherenkov Telescope Array Consortium. Science with the Cherenkov Telescope Array. _World Scientific_ 2019. Available online: https://www.worldscientific.com/doi/pdf/10.1142/10986 (accessed on 5 December 2021).