Gemini Near-infrared Spectroscopy of High-Redshift Fermi Blazars: Jetted Black Holes in the Early Universe Were Overly Massive

Jetted active galactic nuclei (AGNs) are the principal extragalactic $\gamma$-ray sources. Fermi-detected high-redshift ($z>3$) blazars are jetted AGNs thought to be powered by massive, rapidly spinning supermassive black holes (SMBHs) in the early universe ($<2$ Gyr). They provide a laboratory to study early black hole (BH) growth and super-Eddington accretion -- possibly responsible for the more rapid formation of jetted BHs. However, previous virial BH masses of $z>3$ blazars were based on C IV in the observed optical, but C IV is known to be biased by strong outflows. We present new Gemini/GNIRS near-IR spectroscopy for a sample of nine $z>3$ Fermi $\gamma$-ray blazars with available multi-wavelength observations that maximally sample the spectral energy distributions (SEDs). We estimate virial BH masses based on the better calibrated broad H$\beta$ and/or Mg II . We compare the new virial BH masses against independent mass estimates from SED modeling. Our work represents the first step in campaigning for more robust virial BH masses and Eddington ratios for high-redshift Fermi blazars. Our new results confirm that high-redshift Fermi blazars indeed host overly massive SMBHs as suggested by previous work, which may pose a theoretical challenge for models of the rapid early growth of jetted SMBHs.


INTRODUCTION
Active supermassive black holes (SMBHs) with masses ≳ 10 9 M⊙ powering the most luminous quasars have formed when the universe was less than a Gyr after the Big Bang (e.g., Wang et al. 2021).Understanding how SMBHs formed so quickly is a major outstanding problem in modern astrophysics (Natarajan 2014;Reines & Comastri 2016;Inayoshi et al. 2020), because their presence in large number densities may pose a challenge to modeling their rapid formation and subsequent growth (Johnson & Haardt 2016).High-redshift (z > 3) blazars are ideal laboratories to study early SMBH growth within the first 2 Gyr of the universe.These sources are active galactic nuclei (AGNs) powered by billion solar mass black holes with our line of sight lying within an angle ∼ 1/Γ of the jet axis, where Γ is the jet bulk Lorentz factor (Γ ≈ 10-15; Ghisellini et al. 2014).
The Fermi Large Area Telescope (Fermi-LAT; Atwood et al. 2009) has detected high-redshift (z > 3) blazars, which are dominated by flat-spectrum radio quasars (Ackermann et al. 2017) with measured BH masses ≳ 10 9 M⊙ (Marcotulli et al. 2020).The observation of a single blazar implies the presence of ∼ 2Γ 2 misaligned, jetted systems with the same BH mass pointing at other directions, increasing the space density of SMBHs of jetted AGN by 10-20 percent depending on the redshift bin (Sbarrato et al. 2015;Ackermann et al. 2017).The space density of high-redshift Fermi blazars implies that radio-loud systems with SMBHs ≳ 10 9 M⊙ are as common as ⋆ E-mail: xinliuxl@illinois.eduradio-quiet systems, and they may be even more common at higher redshifts (Sbarrato et al. 2015).Therefore, the jetted phase is likely a key ingredient for understanding rapid black hole growth in the early universe.
While high-redshift Fermi blazars are not as distant as the highestredshift quasars known, they still pose a challenge for early SMBH growth models.Jetted BHs are thought to be less efficient accretors due to higher radiative loss (e.g., from having a larger spin, if jet activity is correlated with the BH spin).For example, a spinning BH accreting at Eddington rate would need ∼3.1 Gyr to grow from a seed of ∼ 100 M⊙ to ∼ 10 9 M⊙ (Ghisellini et al. 2013).This implies that such massive, jetted BHs would not have had formed at z > 2, whereas their observed space density peaks around z ∼ 4 (Sbarrato et al. 2015).Simulations suggest that BHs accreting at super-Eddington rates are characterized by strongly collimated outflows or jets (McKinney et al. 2014;Sa ¸dowski et al. 2014).Highredshift blazars provide a laboratory to study super-Eddington accretion which may have led to more rapid growth of the jetted SMBHs than non-jetted systems (Volonteri et al. 2015).
Accurate BH mass measurements are crucial in order to quantify blazar demographics and growth.However, the existing BH mass estimates of high-redshift blazars are highly uncertain.As listed in Table 1 and 2, the existing BH masses of these high-redshift blazars were estimated using two methods: SED modeling of the accretion disk, and virial mass using C IV λ1549 (Paliya et al. 2020).Both estimates may have significant problems.The SED masses are limited by: (1) the data quality to cleanly (with jet contamination) and fully (limited by the Lyα forest absorption at the higher frequency range) sample the accretion disk emission and (2) the assumption of the standard disk model (Shakura & Sunyaev 1973), which may break down for highly spinning BHs.
We present new Gemini/GNIRS NIR spectroscopy to obtain robust virial BH masses for nine Fermi blazars known at z > 3 with existing multi-wavelength observations that maximally sample the SEDs (Paliya et al. 2020).For the virial BH mass, C IV λ1549 is the only broad emission line covered by the existing optical spectra for z > 3 blazars (Ackermann et al. 2017).However, it is well known that C IV λ1549 is a poor virial mass estimator (Shen & Liu 2012;Jun et al. 2015).While a general agreement is found between the SED-based and virial BH masses in 116 Fermi blazars at z < 3.2 (Ghisellini & Tavecchio 2015), the virial masses of > 95% of the sample are based on Hβ and/or Mg II λ2800, which do not suffer from the biases of C IV λ1549.Indeed, the few objects with C IV λ1549based masses in low-redshift blazars exhibit the largest discrepancy with the SED-based masses.C IV λ1549 often shows significant blueshifts due to accretion disk winds and strong outflows, which may be particularly relevant for high-redshift blazars.Significant mass uncertainties induce errors in the mass function and Eddington ratio distribution, hampering a robust test of theories of early SMBH formation and growth.Therefore, we will use Gemini/GNIRS NIR spectroscopy to measure the broad Hβ and Mg II λ2800 to obtain robust virial BH masses.The sample selection, observations, and data analysis are presented in §2, our resulting virial BH masses, bolometric luminosities, and Eddington ratios, along with comparison with the SDSS spectra and SED fitting are presented in §3, and we conclude in §4.

Sample Selection
The targets are the farthest known Fermi-detected γ-ray blazars.A γ-ray detection is a definitive signature for the presence of a closely aligned relativistic jet.The targets are the nine γ-ray detected blazars presented in (Paliya et al. 2020).They were identified in a systematic search of γ-ray counterparts detected by the Large Area Telescope on board Fermi (Atwood et al. 2009) for all the z > 3 radio-loud (R > 10, where R is the ratio of the rest-frame 5 GHz to optical B-band flux density) quasars in the Million Quasar Catalog (Flesch 2015).Compared with the Fermi blazars at lower redshifts, the z > 3 Fermi blazars occupy the region of high γ-ray luminosities and soft photon indices, typical of powerful blazars; their Compton dominance (i.e., ratio of the inverse Compton to synchrotron peak luminosities) is large (> 20), placing them among the most extreme blazar population (Ackermann et al. 2017).
The targets are all flat spectrum radio quasars with prominent emission lines (Urry & Padovani 1995) (rest-frame EW> 5Å) suitable for virial BH mass measurements.They also have existing X-ray observations from the Chandra, XMM-Newton, and/or Swift-XRT data archives for maximally sampling the multi-wavelength SEDs.The targets include nine of the 10 γ-ray detected z > 3 blazars in the fourth catalog of the Fermi-LAT-detected AGNs (Ajello et al. 2020).The remaining object is not included because it lacks existing X-ray data.While there are X-ray selected blazars at even higher redshifts, they are less suitable for this program because their SEDs are not as fully sampled.

SDSS Spectroscopy
We searched the SDSS data release 18 database for optical spectra.Four of our blazars, J1354−0206, J1429+5406, J1510+5702, and J1635+3629 have SDSS spectra, from which we will obtain C IV λ1549 BH mass measurements.

Gemini/GNIRS Observations
We conducted Gemini-N/GNIRS NIR spectroscopy for the nine z > 3 Fermi blazars with the goal to measure Hβ (Mg II λ2800) in the K (J) band for the eight targets at redshift 3.0 < z < 3.7 and Mg II λ2800 in the H band for the z = 4.31 target in regions of high atmospheric transparency.We used the 32 line mm −1 grating with a cross-dispersion and a 0. ′′ 675 slit.This gives a spectral resolution of R ≈ 800 and covers the full J, H, and K bands in the observed 0.8-2.5 µm using nodded exposures along the slit for 120 s each.The observing conditions and resulting data quality varied.We were unable to detect broad emission lines in the NIR for three sources: J0337-1204, J1354-0206, and J1510+5702.

Data Reduction and Analysis
We used the pypeit version 1.10.0package (Prochaska et al. 2020a,b;Prochaska et al. 2020c) to reduce the GNIRS spectra.The reduction pipeline steps include sky subtraction using standard A−B mode subtraction and a B-spline fitting, wavelength calibration using nightsky OH lines, cosmic ray rejection with l.a.cosmic (van Dokkum et al. 2012), flux calibration using spectroscopic standard stars, and telluric correction derived from fitting telluric models from the Line-By-Line Radiative Transfer Model (lblrtm; Clough et al. 2005).The 1D spectra are extracted following the method of Horne (1986).
We fit the continuum and emission lines in each 1D spectrum using a modified version of the publicly-available PyQSOFit code (Guo et al. 2018;Shen et al. 2019).In this code, the continuum is modeled as a blue power-law plus a 3rd-order polynomial for reddening.No Fe II continuum emission was detected, and including it does not significantly improve the fits, so Fe II emission templates (Vestergaard & Wilkes 2001) were not used.The total model is a linear combination of the continuum and single or multiple Gaussians for the emission lines.Since uncertainties in the continuum model may induce subtle effects on measurements for weak emission lines, we first perform a global fit to the emission-line free region to better quantify the continuum.We then fit multiple Gaussian models to the continuum-subtracted spectra around the Hβ and Mg II λ2800 emission line complex regions locally.We define narrow Gaussians as having FWHM < 1800 km s −1 .The narrow line widths are tied together.We use 50 Monte Carlo simulations to estimate the uncertainty in the line measurements.Figure 1 shows an example GNIRS spectra of J0833−0454, the brightest source in our sample.
Given the significant uncertainty in our NIR spectra, we must determine whether the broad line detections are significant or not.We consider the broad emission line detection and associated flux measurement to be reliable using a signal-to-noise ratio criteria.We define the broad line signal-to-noise ratio as: where MAD is the median absolute deviation of the data minus our best-fit model residual (resid) in the window containing the line complex, and A br is the peak or amplitude of the broad emission line.
Exact definitions of the signal-to-noise ratio differ and depends on the details of the spectral modeling.We consider our broad line flux measurements to be reliable if S/N br > 2, which we justify using simulations.We simulate a single ∼ 4000 km s −1 broad Gaussian emission line (typical of our blazars) on top of a noisy power-law continuum, and repeat this procedure with varying emission line amplitudes.The recovered flux f br is then compared to the true input flux fi by measuring the S/N br from each simulated broad line.The S/N br > 2 threshold is determined as the approximate point where the error on the recovered flux is within ∼ 20, which is comparable to the systematic uncertainties from data reduction, flux calibration, and spectral modeling choices, as shown in Figure 2.

Virial black hole mass estimates
Following Shen et al. (2011), we estimate the BH masses using the single-epoch virial method.This method assumes that the broad-line region (BLR) is virialized and uses the continuum luminosity and broad-line FWHM as a proxy for the BLR radius and virial velocity respectively.Under these assumptions, the BH mass can be estimated by: where L br and FWHM br are the broad-line luminosity and fullwidth-at-half-maximum (FWHM) with an intrinsic scatter of ∼ 0.4 dex in BH mass.The coefficients a and b are empirically calibrated against local AGNs with BH masses measured from reverberation mapping.We adopt the calibrations (Vestergaard & Peterson 2006)  (5) We caution that the continuum luminosities we have measured may be over-estimated if they are strongly relativistically beamed by the jet (e.g., Wu et al. 2004;Shaw et al. 2012), resulting in an over-estimation of the BH mass.Previous work has shown systematically-larger BH masses of 0.14 dex on average when broad line luminosities are used instead of continuum luminosities in the BH mass estimation prescriptions (Shaw et al. 2012).To estimate BH masses using the broad line luminosities, we make use of the fact that the line luminosity correlates with the continuum luminosity for broad-line quasars.We substitute the broad line luminosity into Equation 2 and adopt the calibrations derived by Shen et al. (2011)  are ∼ 0.2 dex larger than the masses derived from the broad line luminosities, roughly consistent with Shaw et al. (2012).The BH mass comparisons are shown in Figure 3.

Eddington ratio estimates
We estimate the Eddington ratios by diving the AGN bolometric luminosity by the Eddington luminosity (λ Edd = L bol /L Edd ).Following e.g., Belladitta et al. (2022), we estimate the AGN bolometric luminosity using two methods.First, following Shen et al. (2011), we use the 5100 Å (or 3000 Å) continuum luminosities measured from our GNIRS spectra assuming bolometric corrections (BC) derived from a mean quasar SED (Richards et al. 2006).Those bolometric corrections are BC5100 = 9.26, BC3000 = 5.15.We ignore the small inclination factor of i = cos 0 • cos 30 • = 1.15 in the bolometric correction expected to arise due to differences the viewing angle for Type 1 radio quiet AGNs and blazars.Secondly, we use the disk bolometric luminosities given by Paliya et al. (2020) inferred from SED fitting.The total AGN bolometric luminosity is computed by assuming L bol = 2 L disk (e.g., Calderone et al. 2013).To compute the Eddington luminosities, we usplote the SED-based BH masses against (denoted λ Edd,SED ) our virial BH masses using the broad-line luminosites (denoted λ Edd,spec ).The results are shown in Figure 4.
In Figure 4, we show that the spectral continuum-based bolometric luminosities are, on average, larger than the SED fitting-based luminosities by ∼ 0.5 dex (or about a factor of 3 in luminosity) for our sample.One possibility is that the systematic uncertainty in the bolometric correction (∼ 0.2 dex; Runnoe et al. 2012) is larger for blazars, or that the bolometric corrections for blazars are different.2020) (x-axis) using the continuum (left) or line luminosity (right) and following the recipes in the text.The BH masses estimated from the continuum luminosity are slightly larger than those estimated from the broad line luminosity by ∼ 0.2 dex.The gray y = x line demonstrates there is no strong discrepancy between virial and SED -based BH mass estimates for our sample using either method.We assume typical uncertainties of 0.4 dex for the virial BH masses (Shen 2013) and 0.3 dex for the SED-based BH masses (Paliya et al. 2020).Comparison between bolometric luminosities measured from the continuum of our NIR spectra using the 5100 Å luminosity (or 3000 Å luminosity for J1510+5702) (y-axis) and from the SED fitting approach (x-axis).The gray y = x line demonstrates that the spectral continuum-based bolometric luminosities are, on average, larger than the SED fitting -based luminosities by ∼ 0.2 dex.Right: Comparison between Eddington luminosity ratios derived from the bolometric luminosity measured from the continuum and broad-line BH masses from the M BH, vir, line method from our NIR spectra (y-axis) and from the SED fitting approach (x-axis).We assume typical uncertainties of 0.2 dex for spectral continuum -based bolometric luminosities (L bol, spec = 9.26 L 5100 or 5.15 L 3000 ; Richards et al. 2006;Runnoe et al. 2012) and 0.4 dex for the SED fitting -based bolometric luminosities (L bol, SED = 2 L disk ; Paliya et al. 2020).The gray y = x line and the large uncertainties in the Eddington ratio estimates (summed-in-quadrature from the BH mass and bolometric luminosity uncertainties) demonstrates that they are broadly consistent with each other and that the Eddington ratios are λ Edd ∼ 0.1-1.
This could be explained by contamination of the quasar continuum emission by a relativistic jet.Taken at face value, the discrepancy of about a factor of ∼ 3 in luminosity is small given the Doppler factors of δ = 10 − 20 of our parent sample of high-redshift FSRQs (Paliya et al. 2020).It is likely that only a portion of the UV/optical continuum emission is contaminated by beamed emission toward the observer, in contrast to the jetted radio emission, which is dominated by beamed emission (where Doppler factors are typically measured).This result demonstrates the importance of considering the caveats with virial BH masses and Eddington ratio calculations for blazars, which may have implications for studies of how the Eddington luminosity ratios or accretion disk properties might relate to jet power (e.g., Celotti & Ghisellini 2008).

SUMMARY AND CONCLUSION
High-redshift Fermi blazars represent the most distant sources in the γ-ray sky, providing a laboratory to study early SMBH growth and super-Eddington accretion.They may form from the more rapid formation of jetted BHs, possibly by super-Eddington accretion, but their existing mass estimates are highly uncertain.These sources pose a challenge for models of early SMBH growth and formation.With Gemini/GNIRS spectroscopy of nine z > 3 Fermi blazars, we have: (i) Measured robust virial BH masses using broad Hβ and/or Mg II λ2800 for six of our nine target blazars with sufficient line signal-to-noise ratios.The Hβ and Mg II λ2800 virial masses have proven to be more robust than the existing C IV λ1549-based masses.C IV λ1549 is known to be biased and prone to systematic uncertainties due to strong outflows (Shen & Liu 2012;Jun et al. 2015).
(ii) Compared the improved virial BH masses (with a 0.4 dex systematic scatter) against SED masses independently derived from modeling the accretion disk (∼0.3 dex systematic uncertainty; Paliya et al. 2020).A significant discrepancy would imply deviation from the standard disk (e.g., due to high BH spins) in high-redshift blazars.We do not find a statistically significant discrepancy given the sample size.
(iii) Measured robust Eddington ratios and bolometric luminosities from the Hβ and/or Mg II λ2800 broad lines and continuum luminosities for high-redshift blazars.We do find any compelling evidence to suggest they are super-Eddington accretors (e.g., Begelman & Volonteri 2017;Yang et al. 2020).
Future work with a larger sample size and could help confirm or falsify our main conclusions.This work highlights the challenge of obtaining deep enough NIR spectra to detect and model the broad emission lines sufficiently to estimate BH masses using ground-based observatories.Future observations with JWST would be a compelling test for blazars at even higher redshifts.
Astropy, a community-developed core Python package for Astronomy (Astropy Collaboration et al. 2018, 2013).We thank Gabriele Ghisellini, Tullia Sbarrato, and Jianfeng Wu for useful discussion.We are grateful to Lea Marcotulli for comments which improved our manuscript.
Based in part on observations obtained at the international Gemini Observatory (Program ID GN-2022A-Q-138; PI: X. Liu), a program of NSF's NOIRLab, which is managed by the Association of Universities for Research in Astronomy (AURA) under a cooperative agreement with the National Science Foundation. on behalf of the Gemini Observatory partnership: the National Science Foundation (United States), National Research Council (Canada), Agencia Nacional de Investigación y Desarrollo (Chile), Ministerio de Ciencia, Tecnología e Innovación (Argentina), Ministério da Ciência, Tecnologia, Inovações e Comunicações (Brazil), and Korea Astronomy and Space Science Institute (Republic of Korea).This work was enabled by observations made from the Gemini North telescope, located within the Maunakea Science Reserve and adjacent to the summit of Maunakea.We are grateful for the privilege of observing the Universe from a place that is unique in both its astronomical quality and its cultural significance.

Figure 1 .Figure 2 .
Figure 1.Example GNIRS spectrum of the source J0833−0454, the brightest in our sample.A power-law plus 3rd-order polynomial and Gaussians are used to fit the continuum (yellow) and emission lines, respectively.The Fe II emission templates (teal) are set to zero, as including them does not significantly improve the fits.The data is shown in black and the best-fit model is overplotted in blue.The individual narrow line components are plotted in green, and the broad line components are plotted in red.The gray bands on the top are line-free windows selected to determine the continuum emission.The light gray shaded bands are masked regions effected by telluric absorption at 1.35−1.45µm and 1.8−1.92µm.The lower panels show the zoomed-in emission line regions of C III] λ1909, Mg II λ2800, and Hβ.The full figure set showing the fitting results for each source in our sample is shown in Appendix A.

Figure 3 .
Figure 3.Comparison between virial BH masses (y-axis) for the broad emission lines with S/N > 2 in our GNIRS spectra and SED masses from Paliya et al. (2020) (x-axis) using the continuum (left) or line luminosity (right) and following the recipes in the text.The BH masses estimated from the continuum luminosity are slightly larger than those estimated from the broad line luminosity by ∼ 0.2 dex.The gray y = x line demonstrates there is no strong discrepancy between virial and SED -based BH mass estimates for our sample using either method.We assume typical uncertainties of 0.4 dex for the virial BH masses(Shen 2013) and 0.3 dex for the SED-based BH masses(Paliya et al. 2020).

Figure 4 .
Figure 4. Left: Comparison between bolometric luminosities measured from the continuum of our NIR spectra using the 5100 Å luminosity (or 3000 Å luminosity for J1510+5702) (y-axis) and from the SED fitting approach (x-axis).The gray y = x line demonstrates that the spectral continuum-based bolometric luminosities are, on average, larger than the SED fitting -based luminosities by ∼ 0.2 dex.Right: Comparison between Eddington luminosity ratios derived from the bolometric luminosity measured from the continuum and broad-line BH masses from the M BH, vir, line method from our NIR spectra (y-axis) and from the SED fitting approach (x-axis).We assume typical uncertainties of 0.2 dex for spectral continuum -based bolometric luminosities (L bol, spec = 9.26 L 5100 or 5.15 L 3000 ;Richards et al. 2006;Runnoe et al. 2012) and 0.4 dex for the SED fitting -based bolometric luminosities (L bol, SED = 2 L disk ;Paliya et al. 2020).The gray y = x line and the large uncertainties in the Eddington ratio estimates (summed-in-quadrature from the BH mass and bolometric luminosity uncertainties) demonstrates that they are broadly consistent with each other and that the Eddington ratios are λ Edd ∼ 0.1-1.

Table 1 .
Paliya et al. (2020)r our blazar sample with GNIRS spectra.The SED-based masses (typical uncertainty of 0.3 dex) and accretion disk luminosities (typical uncertainty of 0.4 dex) are fromPaliya et al. (2020).The continuum luminosities are measured from our GNIRS spectra at 5100 Å or 1 3000 Å otherwise.The uncertainties on the continuum luminosities are measurement errors only.True uncertainties are dominated by systematic differences in data reduction and spectral fitting choices.Virial BH masses are given in Table2.

Table 2 .
Virial BH masses measured from broad-line detections (line S/N > 2) in our GNIRS (Mg II, Hβ) or SDSS (C IV) spectra.All uncertainties are measurement errors only.True uncertainties are dominated by systematic differences in data reduction and spectral fitting choices.Missing table entries either do not have sufficient line detections to estimate a BH mass or are not covered within the wavelength range of the spectrum.