Insights into physical conditions and magnetic fields from high redshift quasars

As shown in Figure 1, the redshifting of strong emission lines, such as H$\alpha$, [NII], H$\beta$, and [OIII], in $z>5$ galaxies, can introduce a photometric excess in certain bands. These lines induce a significant difference between the observed photometry ($f_{obs}$) and the underlying continuum flux density ($f_{temp}$), especially if the equivalent width of the lines is high. Thus, if we can fit for the underlying continuum flux density using templates, we can obtain an estimate of the emission line fluxes using the excess emission in the different bands. We note that the IRAC 5.8 and 8 $\mu$m bands and the WISE 3.4$\mu m$ and 4.6$\mu m$ bands, are much noisier than the IRAC 3.6 and 4.5 $\mu$m bands. So an accurate estimate of the photometric uncertainty in each band is crucial for fitting the underlying continuum robustly and to obtain robust signal to noise estimates of the derived line fluxes.

In addition, we make the underlying assumption that variability is not important. In reality, AGN may vary since the multiwavelength data were taken from different surveys at different times by different groups and even in a single band, the data is a combination over multiple epochs of observations. However, variability should result in either an increase or decrease in the flux density in a band - if the dominant trend in the flux density residuals after fitting has a positive excess, it would imply that there is an additional contribution to the emission on top of the intrinsic variability of the source. Secondly, since the sources are unresolved, the flux densities we are measuring are a combination of the emission from the AGN broad line region, the narrow line region and the host galaxy. Although the sources are spectroscopically confirmed to be AGN, through rest-frame ultraviolet spectroscopy, it is unclear how much the integrated rest-frame optical emission would actually vary. Nevertheless, we assess the role of variability on the measurements later in this section.

We first compile 23 AGN-dominated continuum templates obtained from local (Brown et al., 2019) and composite AGN (Polletta et al., 2007; Assef et al., 2010). We exclude strong emission lines from the original templates and smooth the templates to remove any artificial patterns. We redshift the templates to the redshift of the quasars and then apply a correction for absorption by the intergalactic medium (IGM) using a recipe from Inoue et al. (2014). Before fitting templates to the photometry, systematic uncertainties (3% for ground-based near-IR, and 5% for the Spitzer and WISE bands) are added in quadrature to the flux density uncertainties. We also limit the continuum fits to wavelengths where the photometric $SNR>5$.

As an example, for the z=5.81 quasar shown in Figure 2, H$\alpha$ and [NII] line fluxes can be estimated from the photometry in the IRAC 2 and WISE 2 bands, while H$\beta$ and [OIII] can be estimated with IRAC 1 and WISE 1 bands, as shown in Figure 1. However, with typical quasar line widths, because H$\alpha$ is strongly contaminated with the [NII] emission line, it is difficult to separately constrain H$\alpha$ fluxes from [NII] with the broad IRAC/WISE bandpasses. Recent studies (Schindler et al., 2020; de Rosa et al., 2011) have shown through rest-frame ultraviolet spectroscopy that the metallicity of high redshift quasars seems to be very similar to those at lower redshifts, indicating that the stars in their host galaxies have already built up most of their metal abundance. Based on this result, we adopt [NII]/ H$\alpha$ =1 as estimated from local AGN (e.g., SDSS DR12 in Figure 4) and z$\sim$2.5 quasars (Kewley et al., 2013). Since we cannot directly fit the spectra to measure line widths, we fix FWHM = 3400km/s for all lines. FWHM = 3400 km/s is chosen as a median FWHM of the MgII line in 82 quasars at $4.5<z<6.4$ (Kurk et al., 2007; Willott et al., 2010; Trakhtenbro et al., 2011; de Rosa et al., 2011; Chehade et al., 2018) and assuming a correlation of FWHM between MgII and H$\alpha$ lines (McLure & Dunlop, 2004; Jun et al., 2015; Zuo et al., 2015).

With these assumptions, we then evaluate the line fluxes as follows:

1. 1.

   We find the best-fit continuum template of each source by fitting to the high SNR near-IR bands ($y$, J, H, and K) and IRAC 3 and 4 channels, if available, by minimizing $\chi^{2}$ between observed photometry and the continuum templates convolved through the corresponding bandpass.

2. 2.

   We use the excess emission at IRAC 2 and/or WISE 2, to compute H$\alpha$. We assume a Gaussian profile for the line with FWHM=3400km/s and [NII]/H$\alpha$=1. Due to the different shapes of the IRAC and WISE band transmission (as shown in Figure 2), the H$\alpha$ and [NII] contribution in each band should be different, and depends on the redshift. We, therefore, calculate the fractional contribution of H$\alpha$ and [NII] in the IRAC 2 and WISE 2 bands for each object separately.

   Note that the uncertainty weighted line flux values from both IRAC 2 and WISE 2 are derived, when available, for quasars at $z>5.16$ (see Figure 1).

3. 3.

   H$\beta$ is derived by assuming an H$\alpha$/H$\beta$$=2.86$ (i.e. no dust). This is consistent with the H$\alpha$/H$\beta$ ratio of 2.6 that we find for a $z=5.37$ quasar in our sample. This object is unique since only the [OIII] line is in IRAC 1; we can therefore independently measure H$\beta$ after subtracting the [OIII] contribution from the WISE 1 excess. Although we only have this one constraint on the H$\alpha$/H$\beta$ ratio, we note that any dust in the quasar host would only increase the [OIII]/H$\beta$ ratio as we discuss later.

4. 4.

   We derive [OIII] after subtracting the H$\beta$ contribution from the excess emission in IRAC 1 and/or WISE 1. Depending on the redshift, the final line flux is the weighted line flux from both the IRAC 1 and WISE 1 excess, when available, as shown in Figure 1. At $z<5.3$, for instance, we only rely on the WISE 1 photometry to estimate [OIII] after subtracting the H$\beta$ contribution.

To produce robust uncertainty estimates on the fitted emission lines, we use a Monte Carlo simulation. We modify the true photometry in each of the bands by adding a random uncertainty value to the fluxes. The uncertainties are drawn from a normal Gaussian distribution with $\sigma$ corresponding to the photometric uncertainty in that band. We then repeat the fits with the synthetic photometry. The above calculation is repeated with the number of trials, $N_{tri}=600$, to obtain the $1\sigma$ uncertainty and the median of the emission line fluxes for each quasar in the sample. Note that the result does not vary significantly with larger $N_{tri}$. About 38% of 53 sources have negative $(f_{obs}-f_{continum})$ at IRAC and/or WISE bands, suggesting that the emission lines are weak in these systems or that variability is dominating the fitting residuals. As a result, we are able to constrain H$\beta$ and [OIII] line fluxes for only a subset of 33 quasars at $5.03<z<6.3$.

As mentioned earlier, since the multi-wavelength photometry of the quasars is taken at different epochs, it is possible that the variability of the sources between the different observing epochs can cause a significant difference between observed photometry and derived continuum flux. In general, since the quasars were first detected in short wavelength photometry, those bands would then be skewed brighter relative to the continuum in the other bands. To assess this, we investigate the distributions of $(f_{obs}-f_{temp})/\sigma_{obs}$ in all 10 filters, $yJHK$, IRAC 1, 2, 3, 4 channels and WISE bands in Figure 3. We find that the distributions of residuals in yJHK are nearly centered on zero indicating that the residuals are mostly due to photometric noise. In contrast, the bands in which the expected strength of emission lines is strong (i.e. IRAC 1, 2 and WISE bands), are mostly skewed positive. This further indicates that the contamination of the photometry due to emission lines is real and the methodology we have adopted to derive emission line constraints is robust.

In order to assess if the photometry itself may be biased due to systematics, ideally one would compare the line flux derived from photometry with spectroscopy. Unfortunately, measuring spectral lines at these redshifts is only becoming possible now with the James Webb Space Telescope. Such measurements do not exist at the present time and our work is intended to provide a prediction as to the expected line fluxes. However, we can characterize the robustness of the methodology by measuring the photometry of quasars at lower redshifts ($z\sim 3$), where the line flux does not contaminate the Spitzer and WISE bands. This is discussed in the Appendix where we find that lower redshift quasars show no noticeable excess in the IRAC bands, implying that systematics do not affect the measured photometry.

In order to investigate the physical characteristics of z $>5$ quasars and compare them with the local sample, we plot the derived line ratios of our sample on the [NII]/H$\alpha$ versus [OIII]/H$\beta$ diagnostics diagram (Baldwin, Phillips, & Terlevich (1981); BPT diagram). As discussed earlier, we only have line flux constraints on 33 of the quasars in the sample due to large photometric uncertainties in most of the cases, and inadequacy of the templates in a small subset of cases. In order to derive line ratios, we further need to restrict the sample to those which have high SNR line detections from our Monte-Carlo analysis. By restricting the sample to those which have $SNR>3$ in $H\alpha$ and $[OIII]$, we have a final sample of 8 quasars at $5.03<z<6.23$ with a median redshift of $z=5.94$ with which we can study the line ratios. In the Appendix, the catalog of 8 quasars is available in Table A3, A3, and A3. The median bolometric luminosity of this sample of quasars is $2\times 10^{46}$ erg/s with a range of $3.9\times 10^{45}<L_{AGN}~{}[erg/s]<4.6\times 10^{46}$. The median derived rest-frame equivalent widths of H$\alpha$, H$\beta$ and [OIII] are 400Å, 100Å and 440Å. The median derived [OIII] luminosity is $2\times 10^{45}$ erg/s.

Although we have to assume an [NII]/H$\alpha$ ratio as discussed earlier based on published UV spectroscopy, we can for the first time, constrain the location of the eight $z\sim 6$ quasars (black dots) on the BPT diagram (Figure 4). The median [OIII]/H$\beta$ ratio of the eight quasars at $5.03<z<6.23$ is 4.46 (red star with an uncertainty) with an assumption of solar metallicity ([NII]/H$\alpha$$=1$) and no dust. If dust in the host galaxy were present, it would extinct H$\beta$ slightly more than [OIII], and would increase the H$\alpha$/H$\beta$ ratio. This in turn would only serve to increase the derived contribution of [OIII] to the excess in the photometric bandpass and thereby the [OIII]/H$\beta$ ratio. We investigate this by assuming a H$\alpha$/H$\beta$ ratio of 4.5. This has been measured by Jun et al. (2015) in 155 luminous quasars at $3.3<z<6.4$. We find that in this case, we derive a factor of 1.7 vertical increase in the [OIII]/H$\beta$ ratio. We note that simulations suggest that the stellar masses of these objects are likely to be $>10^{10}$ M_☉ (Marshall et al., 2020), so the assumption of solar metallicity and high [NII]/H$\alpha$ is well motivated.

If in contrast, the quasar hosts have built up their black holes faster than their underlying stellar population as is suggested by some models and observations (e.g DeGraf et al., 2015; Vayner et al., 2021), they would be expected to have lower metallicities, and low [NII]/H$\alpha$ ratios, as well. We repeat the calculations described earlier for this unlikely scenario, which we note is inconsistent with the UV derived metallicities of these quasars. We find that if we adopt an [NII]/H$\alpha$=0.25, the high H$\alpha$ line flux translates to a high H$\beta$ line flux as well and therefore reduces the [OIII] contribution to the photometry. This is shown as the purple point in Figure 4. We can therefore say with high confidence, that the [OIII]/H$\beta$ ratios of $z\sim 6$ quasars are likely to be between $2.4-4.5$ with a strong preference towards the high end of these values.

For comparison, we select 14 of the most luminous local AGN (orange circles) having L${}_{AGN}>10^{45}erg/s$ with SNR$>$10 in H$\alpha$, [NII], H$\beta$, and [OIII] (hereafter, the luminous Sloan AGNs). These AGN selected from the Sloan Digital Sky Survey have been classified with the BPT diagnostic (data from SDSS DR12³³3http://skyserver.sdss.org/dr12/en/home.aspx). We find that our sample of $z\sim 6$ quasars, if they have solar metallicity, have unusually strong [OIII]/H$\beta$ which is above the median of local quasars (about 1.4$\sigma$). Their line ratio is comparable to the most extreme local systems as shown by the three AGNs having log([OIII]/H$\beta$) $>0.5$. Even if their metallicities are substantially sub-solar (purple dot), the [OIII]/H$\beta$ ratio is still higher than most luminous AGN in the local Universe. Our result indicates that high-z quasars are located at the high end of the [OIII]/H$\beta$ distribution found among local AGNs. It is possible that this is due to the selection biases introduced from the high SNR requirements for detecting and characterizing lines with broadband photometry which only upcoming spectroscopic observations with JWST will be able to mitigate.

It is worth comparing the derived equivalent widths of the high-z quasars with other samples of quasars. We find rest-frame H$\alpha$ equivalent widths of $\sim$300-800Å at luminosities of $\sim$10⁴⁶ erg s^-1. This is in good agreement with the values derived for $z>5$ quasars, using a similar methodology by Leipski et al. (2014), and similar to the H$\alpha$ equivalent widths of $\sim$5000 Sloan quasars (Shen et al., 2011). The difference is we extend the work of Leipski et al. (2014) to [OIII] and H$\beta$ as well. A comparison between our derived H$\beta$ equivalent widths which assumes a constant H$\beta$/H$\alpha$ ratio and those of Shen et al. (2011), is also worthwhile. The median H$\beta$ rest-frame equivalent width of the local SDSS quasars is 100Å, similar to our sample, implying that much of the excess in the bluer bandpass is because of the [OIII] doublet as discussed earlier and shown in Table A3. We note that this corresponds to a substantially larger [OIII] equivalent width than that observed in $z\sim 2$ quasars of similar luminosities (Vietri et al., 2020).

Nebular line ratios in narrow-line regions of AGN are typically explained by strong radiative shocks resulting from outflows into the surrounding ISM (e.g Allen et al., 2008). It is generally argued that the scatter in the BPT diagram is driven by the strength and hardness of the ionizing radiation field, metallicity, gas density and shock parameters. Stern & Laor (2013) show that in Type I AGN which our systems are, the [OIII]/H$\beta$ remain constant even if one separates the NLR region emission from the BLR region, while the [NII]/H$\alpha$ typically decreases. Since we do not have the kinematic information to distinguish the relative intensities of the two components, we characterize the physical properties derived from the integrated emission. This provides a lower limit to the emission line ratios of the narrow-line region.

Due to the high [OIII] to bolometric luminosity ratio of $\sim$0.1 seen in our sample, it appears that the mechanism for powering the line emission is not the same as the broad line-width quasar-powered outflows seen at $z\sim 2$. Those objects, showed [OIII] to bolometric luminosity ratios of $\sim$10^-3 and had [OIII] luminosities of a few times 10⁴³ erg s^-1 . They showed [OIII] narrow to broad line ratios of 1:2 or 1:3 at comparable bolometric luminosities. If the bulk of our observed [OIII] emission is instead from star-formation driven winds, the ionizing photons from star-formation would also contribute to the H$\alpha$ and H$\beta$ emission. Since high redshift starbursts are thought to have low [NII]/H$\alpha$ ratios and high [OIII]/H$\beta$ ratios (e.g. Shim & Chary, 2013), subtracting a component arising from the quasar itself would shift our measured ratios by at most 0.2 dex in either axis, up and to the right (the direction as the arrow in Figure 5-a), but this is challenging to constrain robustly - it is instead more illustrative to compare the derived properties with rare objects that have similar [OIII]/H$\beta$ line ratios in the local Universe, as we do below.

We use a library of fully radiative shock models calculated with the MAPPINGS III shock and photoionization code from Allen et al. (2008) to interpret our measures of emission line flux ratios. The library consists of grids of models with shock velocities, magnetic fields, and pre-shock density for the solar abundance model. The typical volume average gas density in the NLR is expected to be $n>100cm^{-3}$ which we adopt (Baskin & Laor, 2005). Here, we use the precursor + shock set of models with $n=100cm^{-3}$ and solar metallicity since they most closely match realistic expectations for the physical conditions. Since the grids of models are sparsely distributed in a discretized grid, we do a linear interpolation to fill the gap between grids. We then find the best-matched shock velocity (V_s) and magnetic field (B) by minimizing the $\chi^{2}$ statistics between the shock models and our observed line ratios as shown in Figure 5-a). As a result, we find B = 8 $\mu$G and $V_{s}$ = 415 $km/s$ for the median line ratios of the $z\sim 6$ quasars. We note that these are galaxy-wide field strengths averaged over the narrow-line region of the halo and do not correspond to the Gauss-strength fields in the immediate vicinity of the supermassive black hole which has been probed by the Event Horizon Telescope (Event Horizon Telescope Collaboration et al., 2021). Note that if the gas density is higher, we find B= 316 $\mu$G and $V_{s}$ = 255 $km/s$ with $n=1000\,cm^{-3}$.

For the median of the narrow-line regions of the luminous Sloan AGNs, we derive B = 25.1 $\mu$G and $V_{s}$ = 187.5 $km/s$, with the values for the individual luminous Sloan AGNs shown in Figure 5-b). Comparing values with similar assumptions on gas densities, we find the typical luminous Sloan AGNs in the local Universe appear to have larger B-fields and lower shock velocities, as derived from our fits to the integrated emission line ratios.

However, for the subsample of Sloan AGN (blue dots) having log([OIII]/H$\beta$) $>0.5$, we find shock velocities similar to those that we find in high-z quasars, while the best-fit B-field values are in the range 1.5 to 50 $\mu$G, straddling the best fit that we derive for the high-z quasars. Due to the fact that the [FeII]/[MgII] ratios in luminous quasars appear to have a narrow range across a wide range of redshifts, and simulations argue for high stellar masses and solar metallicities, we conclude that the physical parameters that we derive for this subsample of Sloan AGN is representative of the physical parameters for the high-z quasars and conclude a lower limit to the magnetic field strength to be 1.5 $\mu$G with the best estimate being 8 $\mu$G. Our results argue for high shock velocities and strong magnetic fields in the primordial Universe, similar to that seen in local luminous AGNs.

Our result places the first constraints on the amplitude of magnetic fields in massive halos that harbor $>10^{10}$ M_☉ galaxies and $\sim 10^{9}$ M_☉ super-massive black holes within 1 billion years of the Big Bang. We interpret these results in the context of the expected evolution of seed fields into galactic magnetic fields.

These halos are likely only about 1 kpc in radius, indicating that they have only undergone about $\sim$200 dynamical timescales over their lifetimes by $z\sim 6$. If primordial seed fields from the inflation era were present, they are thought to be of order 10${}^{-17}-10^{-20}$ G and would only be amplified by a small factor within 1 Gyr. Such an evolution would be unable to account for the microGauss fields that we find in these halos.

In contrast, if there were turbulence in the interstellar medium (ISM) during structure formation, this can result in the formation of current eddies which produce a protogalactic seed magnetic field with amplitudes of order 0.01-1 $\mu$G (Biermann & Schlüter, 1951; Ruzmaikin et al., 1988; Kulsrud et al., 1997). This field is large enough, that it can easily get amplified to our derived seed field strengths, as baryons collapse into dark matter halos as has been estimated by Howard & Kulsrud (1997) and is shown in Figure 6. As the seed field strengthens through structure formation, dissipative processes come into play which cause the amplified field strength to decline. Our derived field strength implies that the seed field had to be as high as 0.1 $\mu$G to be amplified to $>$1 $\mu$G by $z\sim 6$.

A weaker seed field ($\sim$0.01 $\mu$G) with a higher amplification factor by the alpha-$\Omega$ dynamo wildly overpredicts the derived fields for low-z quasars as shown in Figure 6. Such a seed evolves as (Ruzmaikin et al., 1988):

where B_seed is the seed field, $t$ is the age of the Universe, $\gamma$ is the number of dynamical timescales while $\tau_{\rm rot}$ is the typical winding timescale expected to be $\sim 10^{7}$ yrs. A stronger seed field ($\sim$6 $\mu$G) with a longer dynamical time would be consistent with both measurements. However, it is challenging to account for such a slow dynamical timescale since even a lower mass dark matter halo like our Galaxy would have undergone at least that many rotations over the central kpc. The dark matter halo of a quasar, which is at least an order of magnitude more massive would likely have far exceeded the 750 dynamical times that a higher seed field would require.