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

The targets are the farthest known Fermi-detected $\gamma$-ray blazars. A $\gamma$-ray detection is a definitive signature for the presence of a closely aligned relativistic jet. The targets are the nine $\gamma$-ray detected blazars presented in (Paliya et al., 2020). They were identified in a systematic search of $\gamma$-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 $\gamma$-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 $\gamma$-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.

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 $\lambda$1549 BH mass measurements.

We conducted Gemini-N/GNIRS NIR spectroscopy for the nine $z>3$ Fermi blazars with the goal to measure H$\beta$ (Mg II $\lambda$2800) in the $K$ ($J$) band for the eight targets at redshift $3.0<z<3.7$ and Mg II $\lambda$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\approx 800$ and covers the full J, H, and K bands in the observed 0.8–2.5 $\mu$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.

We used the pypeit version 1.10.0 package (Prochaska et al., 2020a, b, c) 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 night-sky 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$\beta$ and Mg II $\lambda$2800 emission line complex regions locally. We define narrow Gaussians as having FWHM $<1800$ km s^-1. The narrow line widths are tied together and their wavelengths are locked together allowing for a small systematic shift in each line complex. Additionally, we updated the Paliya et al. (2020) spectroscopic redshifts by eye for our sources to better match the narrow lines, especially the well-detected [O III] $\lambda$5007 line. 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_{\rm 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_{\rm br}>2$, which we justify using simulations. We simulate a single $\sim 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_{\rm br}$ is then compared to the true input flux $f_{i}$ by measuring the $S/N_{\rm br}$ from each simulated broad line. The $S/N_{\rm br}>2$ threshold is determined as the approximate point where the error on the recovered flux is within $\sim 20$, which is comparable to the systematic uncertainties from data reduction, flux calibration, and spectral modeling choices, as shown in Figure 2.

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 $\lambda L_{\lambda}$ and FWHM_br are the continuum luminosity and broad-line full-width-at-half-maximum (FWHM) with an intrinsic scatter of $\sim 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) used in Shen et al. (2011):

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) from Shaw et al. (2012):

We found that the masses derived from the continuum luminosities are $\sim 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.

For typical (non-beamed) AGNs, one can estimate the Eddington ratios by diving the AGN bolometric luminosity by the Eddington luminosity ($\lambda_{\rm{Edd}}=L_{\rm{bol}}/L_{\rm{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 BC${}_{5100}=9.26$, BC${}_{3000}=5.15$. We ignore the small inclination factor of $i=\frac{\cos\ 0^{\circ}}{\cos\ 30^{\circ}}=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_{\rm{bol}}=2\ L_{\rm{disk}}$. The results shown in Figure 4 show that our spectral continuum-based bolometric luminosities are significantly larger than the bolometric luminosities inferred from the disk luminosity by SED fitting, thought to be due to contamination from boosted optical emission from the jet.

Alternatively, one can estimate the disk bolometric luminosity directly from the broad-line luminosity as $L_{\rm{bol}}=\kappa L_{\rm{br}}$. This method assumes the broad line-emitting gas is ionized by continuum emission from the accretion disk. Therefore, the disk luminosity can be related to the broad-line luminosity by the correction factors given by Calderone et al. (2013):

Then, the Eddington ratios can be self-consistently calculated using the BH masses estimated from the broad line luminosities. This approach again avoids the issue of continuum contamination due to jetted emission and has a scatter of about a factor of 2 (Calderone et al., 2013). Given the systematically larger continuum luminosity-based BH masses and bolometric luminosities found previously, we consider this approach the most robust for our sample blazars. We have performed these calculations and show our robust Eddington ratio estimates calculated from only the line luminosities and widths, denoted $\lambda_{\rm{Edd,line}}$, in Figure 5.

In Figure 4, we show that the spectral continuum-based bolometric luminosities are, on average, larger than the SED fitting-based luminosities by $\sim 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 ($\sim 0.2$ dex; Runnoe et al. 2012) is larger for blazars, or that the bolometric corrections for blazars are different. This could be explained by contamination of the quasar continuum emission by optical emission from a relativistic jet. Taken at face value, the discrepancy of about a factor of $\sim 3$ in luminosity is small given the Doppler factors of $\delta=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).