First measurement of the Mg II forest correlation function in the Epoch of Reionization

One of the standard probes to study the intergalatic medium (IGM) is neutral hydrogen manifested as a blended collection of absoption lines known as the Ly$\alpha$ forest. However, due to the large scattering cross section of the Ly$\alpha$ transition, the absorption lines quickly saturate in the presence of residual neutral hydrogen as low as as 0.01%, thus rendering it challenging to utilize the Ly$\alpha$ forest as a probe of the IGM beyond $z\sim 6$. An alternative probe of the high-redshift IGM is required.

Heavy elements, or metals, are expected to be present in the IGM around the Epoch of Reionization due to their ejection by massive stars exploding as supernovae. Low-ionization metals lines (e.g. O i, Mg II, and Si ii) are good tracers of neutral hydrogen in the pre-reionized IGM due to their similar ionization energy, and their low number densities mean that they do not saturate as easily as the Ly$\alpha$ line. Before reionization, forests of these low-ionization lines are expected to arise from predominantly-neutral and overdense regions (Oh, 2002). During and after reionization, the hardening of the UV background causes the disappearance of the low-ionization forests (e.g. O i, Si ii, and Mg ii) and gives rise instead to forests of high-ionization lines (e.g. C iv, N v, and O vi). This phase transition signalling reionization has been tentatively observed in the circumgalactic medium (CGM) of galaxies (e.g. Becker et al., 2009, 2019; Cooper et al., 2019; D’Odorico et al., 2022). Absorption lines by metal ions therefore provide us with an additional tool with which to study the IGM and the reionization history.

Despite the expected presence of metals at early times, their abundance and distribution are unknown, and models do not agree on the typical metallicity of the IGM (e.g. Oppenheimer & Davé, 2006; Oppenheimer et al., 2009; Jaacks et al., 2018; Kirihara et al., 2020; Liu & Bromm, 2020). However, assuming reionization is powered by ionizing photons from massive stars, cosmic reionization and enrichment are then intimately linked, since the same massive stars that reionize the Universe would enrich the IGM to $Z\sim 10^{-3}-10^{-4}Z_{\odot}$ when they explode as supernovae (e.g., Miralda-Escudé & Rees, 1998; Ferrara, 2016). Existing observations of metal absorbers concentrate around $z\sim 2-4$, and they suggest a typical IGM metallicity of $Z\sim 10^{-3}-10^{-2}Z_{\odot}$ as traced by C iv (e.g. Songaila & Cowie, 1996; Songaila, 1997; Schaye et al., 2003; Simcoe, 2011) and exhibit a lack of evolution in the metallicity from $z\sim 4.3-2$ (Schaye et al., 2003). Higher-redshift observations of metal absorbers (e.g. Simcoe et al., 2011; Ryan-Weber et al., 2006; Bosman et al., 2017) are generally probing the CGM of galaxies rather than the diffuse IGM. This is because current observations are limited by the use of standard methods that rely on detecting discrete absorbing lines, which necessitates high signal-to-noise ratio (SNR) and high resolution measurements, without which one is likely probing overdense gas around galaxies. Additionally, observing in the near-IR, where metal absorbers redshift into at higher redshifts, faces myriad challenges arising from the higher sky background and worsening detector sensivity and spectral resolution.

Hennawi et al. (2021) (hereafter H2021) introduced using the auto-correlation of the metal-line forest, which is based on clustering analyses of the Ly$\alpha$ forest, to detect IGM absorbers. The proposed method treats the metal-line forest as a continuous random field and utilizes all spectral pixels to statistically average down the noise, resulting in far higher sensitivity to the expected weak absorption signals. By combining a large hydrodynamical simulation with a semi-numerical reionization topology, H2021 simulated the Mg ii forest in the pre-reionized IGM at $z=7.5$ and showed that its auto-correlation is sensitive to the neutral hydrogen fraction, $x_{\rm{\text{H\,{i}}}}$, and the IGM metallicity, [Mg/H]. The characteristic signal is a peak in the auto-correlation function at the velocity separation of the Mg ii doublet (768 km/s), where the peak amplitude increases with [Mg/H] and the large-scale power decreases with increasing $x_{\rm{\text{H\,{i}}}}$. They showed that $x_{\rm{\text{H\,{i}}}}$ and [Mg/H] can be constrained to within 5% and 0.02 dex, respectively, assuming a mock JWST NIRSpec dataset of 10 quasars. Tie et al. (2022) later extended this method to the C iv forest at $z=4.5$ focusing on the enrichment topology and using a more realistic model for the metal distribution, and showed that the model parameters can also be well-constrained (see Xin et al. in prep for a cross-correlation analysis between metal-line and Ly$\alpha$ forests). Most recently, Karaçaylı et al. (2023) introduced an analytical framework for metal absorption based on the halo model and attempted to measure their 1D power spectra in high- and low-resolution data.

In this work, we aim to jointly constrain the neutral fraction and metallicity of the IGM at $z\geq 6.8$ for the first time with the auto-correlation of the Mg ii forest in a sample of ten high-redshift quasars. We describe our quasar dataset in §2 and the reduction of their spectra in §3. In §4, we address contamination from CGM absorbers that can bias our signal and describe steps to mask them out. We present our results in §5, where we measure the correlation function of the Mg ii forest and describe our procedure to constrain the metallicity and neutral fraction of the IGM using forward models and Bayesian inference. We conclude in §6. Throughout this work, we adopt a $\Lambda$CDM cosmology with the following parameters: $\Omega_{m}$ = 0.3192, $\Omega_{\Lambda}$ = 0.6808, $\Omega_{b}$ = 0.04964, and $h$ = 0.67038, which agree with the cosmological constrains from the CMB (Planck Collaboration et al., 2020) within one sigma. All distances in this work are comoving, denoted as cMpc or ckpc, unless explicitly indicated otherwise. We define metallicity as $[\mathrm{Mg/H}]$ $\equiv\mathrm{log}_{10}(Z/Z_{\odot})$, where $Z$ is the ratio of the number of magnesium atoms to the number of hydrogen atoms and $Z_{\odot}$ is the corresponding Solar ratio, where we used $\mathrm{log}_{10}(Z_{\odot})=-4.47$ (Asplund et al., 2009).

Our dataset consists of ten quasars at $z\geq$ 6.80: five $z\geq$ 7 quasars observed with Keck/MOSFIRE (McLean et al., 2010, 2012), four $z\geq$ 6.8 quasars observed with Keck/NIRES (Wilson et al., 2004), and one ultra-deep VLT/XSHOOTER (Vernet et al., 2011) observation of a $z=7.1$ quasar (Bosman et al., 2017). Our sample of quasars and their properties are listed in table 1. J0313$-$1806 (Wang et al., 2021) is the current record holder for the highest-redshift quasar discovered, followed by J1342$+$0928 (Bañados et al., 2018), while J0038$-$1527 (Yang et al., 2019) is a broad-absorption line quasar (as is J0313$-$1806). We also include two new unpublished high-$z$ quasars, denoted as ‘newqso1’ (Bañados et al. in prep) and ‘newqso2’ (Yang et al. in prep) for the rest of this paper. We note that there are other ground-based observations of $z\geq$ 6.8 quasars, e.g. observed with VLT/XSHOOTER and GEMINI/GNIRS, but these are either too low in signal-to-noise ratio (SNR $\leq$ 3) or have too low spectral resolution in the case of GNIRS.

Our MOSFIRE observations were executed following an ABAB dither sequence in the $K$-band covering from $1.9-2.4$ $\mu$m and through a 0.7” slit, giving a nominal resolution of $R=3610$ or FWHM = 83 km/s, where the spectral sampling is 2.78 pixels per resolution element¹¹1https://www2.keck.hawaii.edu/inst/mosfire/grating.html, resulting in a pixel size of $dv=30$ km/s. The NIRES observations were acquired with an ABBA dither sequence and provide full wavelength coverage of the YJHK bands ($0.95-2.45$ $\mu$m) through a fixed 0.55” slit, where the mean resolution is $R=2700$ or FWHM = 111 km/s and the sampling is 2.7 pixels per resolution element²²2https://www2.keck.hawaii.edu/inst/nires/spectrometer.html, resulting in a pixel size of $dv=40$ km/s. The XSHOOTER observation of J1120$+$0641 was taken by Bosman et al. (2017) over 30 hours, observed through a 0.9” slit in the VIS and NIR arms. Bosman et al. (2017) measured the resolution to be $R=7000$ or FWHM = 43 km/s, which we adopt here, and we will use a sampling of 3.7 pixels³³3https://www.eso.org/sci/facilities/paranal/instruments/xshooter/inst.html, which gives a pixel size of $dv=12$ km/s.

We reduce our spectra with PypeIt¹⁹¹⁹19https://pypeit.readthedocs.io/en/latest/, a Python package for semi-automated reduction of astronomical slit-based spectroscopy (Prochaska et al., 2020b, a). After performing standard calibration steps and flux calibration on each exposure, we coadd all exposures of the object and apply telluric correction on the combined spectrum. Since the spectra have different resolutions, wavelength grids, and pixel sizes, we coadd them all onto a common grid with a pixel size of 40 km/s. Finally, we restrict our analyses to a minimum wavelength of 1.95 $\mu$m, corresponding to $z_{\text{Mg\,{ii}}}=5.96$.

We fit for the quasar continuum using a B-spline with a breakpoint at every-60^th pixel (amounting to a breakpoint at every 2400 km/s spacing) and we mask strong absorbers identified by-eye before fitting the continuum. We use the same breakpoint spacing for all the quasars. Figures 1 to 3 show the spectra before and after continuum normalization for representative members of our dataset. For each quasar, the gray shaded regions in the top panel indicate regions that are masked before continuum fitting due to visually-detected strong absorbers and masks provided by PypeIt from the data reduction. In the bottom panels, we mask regions due to proximity to quasars and where $z>z_{\rm{QSO}}$ (gray), low SNR due to telluric absorption (blue), and detected CGM absorbers (red; see §4); these are excluded from the correlation function analyses.

We also mask pixels that fall within a quasar proximity zone, as these pixels will be affected by the quasar own ionizing radiation and produce biased results on the IGM neutral fraction. For the proximity zone (PZ) mask, we remove pixels within rest-frame $\Delta v=7643$ km/s blueward from the quasar Mg ii emission line, where we define the Mg ii emission line to be the midpoint of the Mg ii $\lambda 2796$Å and $\lambda 2804$Å doublet at 2800 $\rm{\AA}$. Our choice of $\Delta v=7643$ km/s follows Bosman et al. (2021), who consider pixels at $<$ 1185 $\rm{\AA}$ for their Ly$\alpha$ optical depth study, which results in a proximity zone size of 7643 km/s, or 8.64 proper Mpc at the mean redshift of our quasars, which is $z$ = 7.16. We use a fixed PZ size for all quasars as there are indications that quasar PZ sizes do not evolve strongly with redshift (Eilers et al., 2017; Satyavolu et al., 2023), and even in the presence of evolution, our adopted value of 8.64 proper Mpc is rather conservative given that the typical PZ size for $z\sim 6$ quasars ranges from $2-7$ proper Mpc (Eilers et al., 2017; Eilers et al., 2020; Satyavolu et al., 2023). The unshaded regions are used for the correlation function measurements, and they range from $z_{\rm{\text{Mg\,{ii}}}}=5.96$ to $z_{\rm{\text{Mg\,{ii}}}}=7.42$, with $z_{\rm{\text{Mg\,{ii}},med}}=6.469$, where $z_{\rm\text{Mg\,{ii}}}\equiv(\lambda_{\rm{obs}}\,/\,2800\,\rm{\AA})-1$ is the pixel redshift and $z_{\rm{\text{Mg\,{ii}},med}}$ is the median value. Table 1 includes the pathlength of each quasar that contributes towards the correlation function measurement, where the total pathlength of our dataset is $\Delta z=9.75$ or $\Delta X=47.1$. $\Delta X$ is the absorption pathlength, where $X(z)=\int_{0}^{z}(1+z^{\prime})^{2}H_{0}/H(z^{\prime})\,dz^{\prime}$ (Bahcall & Peebles, 1969). We compute $\Delta X_{\rm{qso}}=X(z_{i})-X(z_{f})$ for each quasar and take the sum as the total absorption pathlength quoted above.

We discover that the noise in the quasar spectra produced by PypeIt are slightly over- and under-estimated and not normally distributed, and we apply correction factors to correct for this such that the noise distribution is Gaussian. To determine the correction factors, we compute the residual $\frac{1-F_{\rm{norm}}}{\sigma_{\rm{norm}}}$ after applying all masks (including CGM masks described in §4) and obtain the median absolute deviation (MAD) of the residual distribution for each quasar as the correction factor to the noise. We obtain correction factors of 0.972, 0.956, 0.908, 1.059, 1.086, 1.055, 0.93, 0.898, 0.88, 1.051 to the noise vectors of quasars no. 1 to no. 10 (see Table 1), respectively. After applying these corrections, the median SNR of our dataset is 13.7. Figure 4 zooms in on certain regions of the Mg ii forest from the normalized spectra, where the typical flux fluctuations given our data quality is $\sim$ 10%, or higher towards the redder part of the spectrum. Typical IGM absorbers are expected to produce percent-level fluctuations (see Figure 2 of H2021), so attempting to detect them at this SNR level will be challenging.

As we are interested in measuring properties of the IGM, Mg ii absorbers within the virial radii of galaxies, namely those in their circumgalactic media (CGM), can bias our measurement (Figure 15 of H2021). While the strongest absorbers can be identified by eye, weaker but more abundant absorbers are difficult to detect visually. We briefly describe our procedures to mask out the CGM absorbers here, but see H2021 or Tie et al. (2022) for more details.

The first masking procedure uses the flux decrement ($1-F$) probability distribution function to filter out CGM absorbers. The rationale behind this procedure is that the $1-F$ PDF of the Mg ii forest will exhibit a tail towards large values (strong absorptions), which is sourced by rare strong CGM absorbers, and a peak at smaller value due to IGM absorbers (Figure 14 of H2021). Figure 5 (left panel) shows the flux PDF for our data, where we mask all pixels with $1-F>0.3$ (shaded region). The top axis shows the the pixel equivalent width, $W_{\lambda,\mathrm{pix}}\equiv(1-F)\,\Delta\lambda$, where $\Delta\lambda$ is the spectral pixel width at rest-frame Mg ii 2796 $\mathrm{\AA}$, which depends on the spectral sampling. Using our fixed common grid with $dv=40$ km/s, $\Delta\lambda\sim 0.37\,\rm{\AA}$.

While a flux cut helps remove strongly-absorbed pixels, it is usually insufficient for removing all affected pixels due to the broad absorption wings of CGM absorbers. As such, we supplement the flux cut with a second masking procedure, which identifies absorbers using the transmission field convolved with the profile of the Mg ii doublet. The convolved transmission field is known as the significance or $\chi$ field,

where $\sigma_{F}^{2}$ is the flux variance (noise vector) of the data and $W(v)=1-e^{-\tau(v)}$ is a matched filter with $\tau(v)$ being the optical depth of a Mg ii doublet with a Gaussian velocity distribution, assuming $N_{\mathrm{\text{Mg\,{II}}}}=10^{13.5}$ cm^-2 and Doppler parameter $b=\sqrt{2}\times\mathrm{FWHM}/2.35$, in which we use the FWHM of the respective quasar spectrum, except for the XSHOOTER spectrum of J1120$+$0641. For J1120$+$0641, we use a larger FWHM of 150 km/s, considering that we rebin all spectra to a fixed 40 km/s grid. Note that $F$ is the continuum-normalized flux and $\sigma_{F}^{2}$ is the continuum-normalized variance. Given the $\chi$ spectrum, absorbers are identified by their extreme $\chi$ value. We mask out pixels above $\chi>3$ as well as pixels $\pm 300$ km/s around each masked pixel to account for their absorption wings. Figure 5 (right panel) shows the $\chi$ distribution of the data. In both panels, in addition to the distribution of each quasar (colored lines), we also plot the distribution for all quasars (black) and pure Gaussian noise sampled from the noise vectors of all quasars (blue). While the $1-F$ PDF from all quasars is similar to that of noise, the $\chi$ PDF shows evidence for fluctuations that are not pure noise at $\chi>3$, since the noise PDF drops off at $\chi\sim 3$. Both PDFs also peak at some characteristic values as expected (see Figure 13 of Hennawi et al., 2021), although the $1-F$ PDF does not flatten off at large values and the plateau at large $\chi$ is weak in the $\chi$ PDF. In summary, we mask 1) pixels with $1-\rm{F}>0.3$ and 2) pixels with $\chi>3$ and those $\pm 300$ km/s around each masked pixel.

Figure 6 shows the masked absorbers for the three highest-$z$ quasars in our dataset (we show the rest in Figure 18 and 19 in Appendix B for brevity). The green and magenta horizontal lines indicate our flux and $\chi$ cutoffs, respectively. We discover new strong absorber candidates in J0313$-$1806, J0038$-$1527, and newqso1, and recover known absorbers at $z=6.84$ (Simcoe et al., 2020) and $z=6.27$ in the spectrum of J1342$+$0928. We recover two of the three weak absorbers identified by Bosman et al. (2017) in J1120$+$0641, and for the one that we do not detect, its chi values are $\sim 2$ and it is masked manually. These absorbers are indicated by the blue ticks in Figure 18 in Appendix B.

The auto-correlation of the Mg ii forest, when treated as a continuous random field, has been shown to probe the neutral fraction and metallicity of the IGM (H2021). We measure this statistic for the first time using a sample of ten $z\geq 6.80$ quasars observed with the Keck/MOSFIRE, Keck/NIRES, and VLT/X-SHOOTER spectrographs, where the median redshift of our measurement is $z=6.469$. We also measure the correlation function over three redshift bins: “low-$z$” ($z_{\rm\text{Mg\,{ii}}}<6.469$, $z_{\rm{\text{Mg\,{ii}},med}}=6.235$), “high-$z$” ($z_{\rm\text{Mg\,{ii}}}\geq 6.469$, $z_{\rm{\text{Mg\,{ii}},med}}=6.715$), and “all-$z$” ($z_{\rm{\text{Mg\,{ii}},med}}=6.469$). Our masking procedures to remove CGM absorbers successfully recover known strong absorbers as well as weaker absorbers.

The correlation function of the Mg ii forest without the CGM absorbers masked validates our measurement technique. We detect strong peaks in the correlation functions at the characteristic velocity separation of the Mg ii doublet (768 km/s) in all redshift bins as well as a rise towards small velocity lags owing to the velocity structure of CGM absorbers. After masking out the identified CGM absorbers, our correlation function measurements yield a null result that is consistent with noise in all redshift bins. We performed Bayesian inference where forward-modeled spectra were used to syntehsize a model-dependent covariance matrix of the data. Our MCMC analysis yields an upper limit of [Mg/H] $<-3.73$ at 95% confidence at a median redshift $z_{\rm\text{Mg\,{ii}}}=6.47$, while we are unable to place any strong constraints on the hydrogen neutral fraction. We place our upper limits in context with other literature measurements of the IGM metallicity in Figure 15. Simcoe (2011) (red points in figure) obtained an estimate of [C/H] via Voigt profile fitting of a sample of discrete H i absorbers with log$(\rho/\bar{\rho})=0.47$, finding a median abundance [C/H] = $-3.55$ at $z\sim 4.3$, which is $2-3$ times lower than similar measurements at $z\sim 2.4$ using both C iv and O vi. On the other hand, Schaye et al. (2003) measured the distribution of C iv pixel optical depth as a function of the corresponding H i optical depth and converted this relation into an estimate of [C/H] as a function of density (the so-called pixel optical method). They measured a median [C/H] $=-3.47$ at $z=3$ in the low density IGM gas with log$(\rho/\bar{\rho})=0.5$ and found little evolution in [C/H] from $z=4$ to $z=2$, as indicated by their best-fit model in the gray shaded region. To establish the existence or lack thereof of an evolution in the IGM metallicity, it is crucial to perform more measurements at $z\sim 3$ and fill the redshift gap at $z>4.5$ and our results provide an important step in this direction.

Despite a null detection in this first study, a statistical method like the correlation function remains one of the best ways to measure the IGM metallicity at high redshift, which has never been attempted before due to observational difficulties and challenges posed by the traditional method of discrete line identification. Looking ahead to the future with spectra from JWST (James Webb Space Telescope; Gardner et al., 2006), we expect to place stronger constraints on primordial cosmic enrichment and reionization. For instance, using a sample of eight $z\geq 7$ quasar spectra with a SNR/pixel of $\sim$76 from JWST/NIRSpec²⁰²⁰20https://www.stsci.edu/jwst/science-execution/program-information?id=3526, we expect to jointly constrain [Mg/H] within 1$\sigma$ precision of 0.04 dex and $x_{\rm{\text{H\,{i}}}}$ within 1$\sigma$ precision of 10%, as well as place an upper limit of [Mg/H] $<-3.94$ on the IGM enrichment at 95% credibility. Finally, by combining very high-$z$ measurements from JWST that is probing the neutral IGM with ground-based measurements of, e.g. the C iv forest, at $z\sim 3-4$ that is probing the ionized IGM, we may be able to detect the transition from forests of low-ionization metal absorbers to forests of high-ionization lines, which is the smoking gun of reionization as seen through metal absorption lines.

We acknowledge helpful conversations with the ENIGMA group at UC Santa Barbara and Leiden University. This project has received funding from the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation program (grant agreement No 885301). JFH acknowledges support from the National Science Foundation under Grant No. 1816006. JO acknowledges support from grants PID2022-138855NB-C32 and CNS2022-135878 from the Spanish Ministerio de Ciencia y Tecnologia. This research made use of Astropy²¹²¹21http://www.astropy.org, a community-developed core Python package for Astronomy (Astropy Collaboration et al., 2013, 2018) and SciPy²²²²22https://scipy.org (Virtanen et al., 2020).

The data underlying this article will be shared on reasonable request to the corresponding author.

During our investigations, we discovered that our likelihood yields MCMC samples and best-fit models that appear visually inconsistent with the actual measurements and error bars. The source of this discrepancy is correlated velocity bins with negative values in the data that are not supposed to occur given the correlation structure of the covariance matrix. To rectify this, we mask the outlying velocity bins before performing the MCMC. To determine which bins should be masked, we compare values of the correlation function measured from the data with values measured from the forward models in which the IGM model has no detectable signal. This is because which velocity bins are outliers depend on the structure of the covariance matrix, which is model-dependent. Therefore, we choose a model with no signal to remove sensitivity towards the signal.

Figures 16 and 17 show the correlation functions in one velocity lag bin plotted against the values of all other velocity bins for the all-$z$ bin. The black points are measured from forward models with $x_{\rm{\text{H\,{i}}}}$ = 0.50 and [Mg/H] = $-4.50$, where the ellipses denote their 1, 2, and 3$\sigma$ contours. Measurements from the real data that fall outside of the 3$\sigma$ ellipse are marked as red while measurements that fall within the 3$\sigma$ ellipse are marked as green. We mask the velocity bin where 20 or more of the real measurements fall outside the 3$\sigma$ ellipse; Figure 16 shows two examples of bad bins that are masked, while Figure 17 shows two examples of good bins that we retain. In conclusion, we mask $dv$ = 370, 450, 1170, and 1490 km/s in the all-$z$ bin, $dv$ = 50 and 1090 km/s in the high-$z$ bin, and $dv$ = 50, 130, 290, 370, 450, 770, 930, and 1490 km/s in the low-$z$ bin (see Figures 12 - 14).

We generate the forward-modeled datasets and compute their covariance matrices for all bins and only apply the bin masking before the MCMC sampling, in which we mask the bad bins in both the correlation function and the covariance matrix. The best-fit models and upper limits are therefore obtained from the masked data and masked covariance matrices.

Here we show the masked absorbers for the rest of the quasars in our dataset that are not shown in §4, which are Figure 18 and 19.

Here we show the forward-modeled spectra for the rest of the quasars in our dataset that are not shown in §5, assuming an IGM model with (${x_{\rm{\text{H\,{i}}}},\rm{[Mg/H]}})=(0.50,-4.50)$. These are shown in Figure 20 and 21.