First Measurement of the Mg II Forest Correlation Function in the Epoch of Reionization

MNRAS 535, 223–246 (2024) https://doi.org/10.1093/mnras/stae2193 Advance Access publication 2024 September 24 First measurement of the Mg II forest correlation function in the Epoch of Reionization Suk Sien Tie , 1 ‹Joseph F. Hennawi , 1 , 2 ‹Feige Wang, 3 Silvia Onorato , 2 Jinyi Yang, 3 Eduardo Ba ˜ nados , 4 Frederick B. Davies 4 and Jose O ˜ norbe 5 1 Department of Physics, Broida Hall, University of California, Santa Barbara, Santa Barbara, CA 93106-9530, USA 2 Leiden Observatory, Leiden University, PO Box 9513, NL-2300 RA Leiden, the Netherlands 3 Steward Observatory, University of Arizona, 933 North Cherry Avenue, Tucson, AZ 85721, USA 4 Max Planck Institut f ¨ ur Astronomie, K ¨ onigstuhl 17, D-69117, Heidelberg, Germany 5 Facultad de F ´ ısicas, Universidad de Sevilla, Avda. Reina Mercedes s/n. Campus Reina Mercedes. E-41012, Seville, Spain Accepted 2024 September 5. Received 2024 August 18; in original form 2023 August 22 A B S T R A C T Strong low-ionization transitions like the Mg II λ2796 , 2804 Å doublet are believed to produce a detectable ‘metal-line forest’, if metals pollute the neutral intergalatic medium (IGM). We measure the autocorrelation of the Mg II forest transmission using 10 ground-based z ≥6 . 80 quasar spectra probing the redshift range 5 . 96 < z Mg II < 7 . 42 ( z Mg II , median = 6 . 47). The correlation function exhibits strong small-scale clustering and a pronounced peak at the doublet velocity ( v = 768 km s −1 ) arising from discrete absorbers in the circumgalactic medium of galaxies. After these strong absorbers are identified and masked the signal is consistent with noise. Our measurements are compared to a suite of models generated by combining a large hydrodynamical simulation with a seminumerical reionization topology, assuming a simple uniform enrichment model. We obtain a 95 per cent credibility upper limit of [ Mg / H ] < −3 . 73 at z Mg II , median = 6 . 47, assuming uninformative priors on [Mg/H] and the IGM neutral fraction x H I . Splitting the data into lowz (5 . 96 < z Mg II < 6 . 47; z Mg II , median = 6 . 235) and highz (6 . 47 < z Mg II < 7 . 42; z Mg II , median = 6 . 72) subsamples again yields null detections and 95 per cent upper limits of [ Mg / H ] < −3 . 75 and < −3 . 45, respectively. These first measurements set the stage for making the Mg II forest an emerging tool to precisely constrain the Universe reionization and enrichment history. K ey words: methods: observ ational – intergalactic medium –quasars: absorption lines –cosmology: observations. 1 INTRODUCTION One of the standard probes to study the intergalatic medium (IGM) is neutral hydrogen manifested as a blended collection of absorption lines known as the Ly αforest. However, due to the large scattering cross section of the Ly αtransition, the absorption lines quickly saturate in the presence of residual neutral hydrogen as low as as 0.01 per cent, thus rendering it challenging to utilize the Ly αforest as a probe of the IGM beyond z ∼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 supernov ae. Lo w-ionization metals lines (e.g. O I , Mg II , and Si II ) are good tracers of neutral hydrogen in the prereionized IGM due to their similar ionization energy, and their low number densities mean that they do not saturate as easily as the Ly αline. Before reionization, forests of these low-ionization lines are expected to arise from predominantly neutral and o v erdense regions (Oh 2002 ). During and after reionization, the hardening of the ultraviolet background causes the disappearance of the low-  E-mail: [email protected] (SST); [email protected].edu (JFH) 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, Rauch & Sargent 2009 ; Becker et al. 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 ´ e 2006 ; Oppenheimer, Dav ´ e & Finlator 2009 ; Jaacks et al. 2018 ; Kirihara et al. 2020 ; Liu & Bromm 2020 ). Ho we ver, 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 ∼10 −3 −10 −4 Z when they explode as supernovae (e.g. Miralda-Escud ´ e & Rees 1998 ; Ferrara 2016 ). Existing observations of metal absorbers as traced by C IV concentrate around z ∼2 −4 and they suggest a typical IGM metallicity of Z ∼10 −3 −10 −2 Z  (e.g. Songaila & Cowie 1996 ; Songaila 1997 ; Schaye et al. 2003 ; Simcoe 2011 ) and exhibit a lack of evolution in the metallicity from z ∼4 . 3 −2 (Schaye et al. 2003 ). Based on a catalogue of © 2024 The Author(s). Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https:// creativecommons.org/ licenses/ by/ 4.0/ ), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 224 S. Tie et al. MNRAS 535, 223–246 (2024) C IV absorbers that is 50 per cent complete at log( N C IV /cm −2 ) = 13.2, which is assembled from 42 high signal-to-noise ratio (SNR, 10 −30 pixel −1 ) and high-resolution spectra of z ∼6 quasars from XSHOOTER (Vernet et al. 2011 ), Davies et al. ( 2023a , b ) observed a rapid decline in C I v from 4 . 3  z  6 . 3, which they concluded to plausibly arise from chemical evolution alone. Intervening Mg II absorbers are detected as discrete systems from redshift 2 to 7 (Chen et al. 2017 ), where these absorbers are assumed to arise in the CGM due to their weakly ionized state. Using 42 lines of sight from XQR-30 (Davies et al. 2023a ), Sebastian et al. ( 2024 ) hav e unco v ered a significant population ( N= 131) of weak Mg II absorbers (0 . 03 < W < 0 . 1 Å) that display a flat evolution in d n/ d X from redshift 1 . 9 < z < 6 . 4. In addition, Christensen et al. ( 2023 ) used JWST NIRSpec to probe four lines of sight and showed that the flat d n/ d X of all discrete Mg II absorbers ( W > 0 . 3 Å) extends to redshift 7, including the highest redshift intervening Mg II absorber disco v ered to date at z = 7.443. Their work reports significant differences in both number and equivalent width of discrete JWST NIRSpec absorbers compared with XSHOOTER spectra, which can be attributed to the difference in spectral resolution of the two instruments (see also Vanni et al. 2024). High-redshift observations of metal absorbers (e.g. Ryan-Weber, Pettini & Madau 2006 ; Simcoe et al. 2011 ; Bosman et al. 2017 ; Codoreanu et al. 2018 ) 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 SNR and high-resolution measurements, without which one is likely probing o v erdense 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 sensitivity and spectral resolution. Hennawi et al. ( 2021 , hereafter H2021 ) introduced using the autocorrelation of the metal-line forest, which is based on clustering analyses of the Ly α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 seminumerical reionization topology, H2021 simulated the Mg II forest in the pre-reionized IGM at z = 7 . 5 and showed that its autocorrelation is sensitive to the neutral hydrogen fraction, x H I , and the IGM metallicity, [Mg/H]. The characteristic signal is a peak in the autocorrelation function at the velocity separation of the Mg II doublet (768 km s −1 ), where the peak amplitude increases with [Mg/H] and the large-scale power decreases with increasing x H I . They showed that x H I and [Mg/H] can be constrained to within 5 per cent and 0.02 de x, respectiv ely, assuming a mock JWST NIRSpec data set 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 (for a cross-correlation analysis between metal-line and Ly αforests, see Xin et al. (in preparation), for its use in constraining the IGM enrichment and Meyer et al. 2019 for its use in constraining the sources of reionization). Most recently, Kara c ¸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 ≥6 . 8 for the first time with the autocorrelation of the Mg II forest in a sample of 10 high-redshift quasars. We describe our quasar data set in Section 2 and the reduction of their spectra in Section 3 . In Section 4 , we address contamination from CGM absorbers that can bias our signal and describe steps to mask them out. We present our results in Section 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 Section 6 . Throughout this work, we adopt a lambda-cold dark matter cosmology with the following parameters: m = 0.3192,  = 0.6808, b = 0.04964, and h = 0.67038, which agree with the cosmological constraints from the cosmic microwave background (Planck Collaboration VI 2020 ) within 1 σ. All distances in this work are comoving, denoted as cMpc or ckpc, unless explicitly indicated otherwise. We define metallicity as [ Mg / H ] ≡log 10 ( Z/Z ), where Z is the ratio of the number of magnesium atoms to the number of hydrogen atoms and Z is the corresponding Solar ratio, where we used log 10 ( Z ) = −4 . 47 (Asplund et al. 2009 ). We would like to mention the caveat that magnesium can be α-enhanced as a result of Type II SNe and so should not be treated as being representative. 2 OBSERVATION Our data set consists of 10 quasars at z ≥6.80: five z ≥7 quasars observed with Keck/MOSFIRE (McLean et al. 2010 , 2012 ), four z ≥ 6.8 quasars observed with Keck/NIRES (Wilson et al. 2004 ), and one ultradeep 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 disco v ered, followed by J1342 + 0928 (Ba ˜ nados 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 highz quasars, J0410 −0139 (Ba ˜ nados et al. 2024 ) and J0038 −0653 (Yang et al. in preparation). We note that there are other groundbased observations of z ≥6.8 quasars, for example, observed with VLT/XSHOOTER and GEMINI/GNIRS, but these are either too low in SNR ( ≤3) or have too low spectral resolution in the case of GNIRS. Our MOSFIRE observations were e x ecuted following an ABAB dither sequence in the K-band co v ering from 1 . 9 −2 . 4 μm and through a 0.7 arcsec slit, giving a nominal resolution of R = 3610 or FWHM = 83 km s −1 , where the spectral sampling is 2.78 pixels per resolution element, 1 resulting in a pixel size of d v = 30 km s −1 . The NIRES observations were acquired with an ABBA dither sequence and provide full-wavelength coverage of the YJHK bands (0 . 95 −2 . 45 μm) through a fixed 0.55 arcsec slit, where the mean resolution is R = 2700 or FWHM = 111 km s −1 and the sampling is 2.7 pixels per resolution element, 2 resulting in a pixel size of d v = 40 km s −-1 . The XSHOOTER observation of J1120 + 0641 was taken by Bosman et al. ( 2017 ) o v er 30 h, observed through a 0.9 arcsec slit in the VIS and NIR arms. Bosman et al. ( 2017 ) measured the resolution to be R = 7000 or FWHM = 43 km s −1 , which we adopt here, and we will use a sampling of 3.7 pixels, 3 which gives a pixel size of d v = 12 km s −1 . 1 https://www2.k eck.haw aii.edu/inst/mosfire/grating.html 2 https://www2.k eck.haw aii.edu/inst/nires/spectrometer.html 3 https:// www.eso.org/ sci/ facilities/ paranal/ instruments/ xshooter/ inst.html Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 225 MNRAS 535, 223–246 (2024) Table 1. Highz quasars used in this work. No. Name z em J mag K mag t exp SNR a SNR corr b z Mg II Instr. u eff k c u eff k u eff k (AB) (AB) (s) (Allz) (Highz) (Lowz) 1 J0313 −1806 d 7.642 20.94 ±0 . 13 19.96 ±0 . 13 25 920 18.04 18.56 1.45 MOSFIRE 0.332 0.360 0.314 2 J1342 + 0928 e 7.541 20.30 ±0 . 02 20.03 ±0 . 12 12 960 14.38 15.04 1.35 MOSFIRE 0.139 0.187 0.101 3 J1007 + 2115 f 7.515 20.22 ±0 . 18 19.75 ±0 . 08 7920 9.82 10.82 1.33 NIRES 0.036 0.046 0.032 4 J1120 + 0641 g 7.085 20.35 ±0 . 15 —114 000 14.7 13.88 0.91 XSHOOTER 0.089 0.035 0.130 5 J0252 −0503 h 7.001 20.19 ±0 . 07 19.92 ±0 . 08 13 680 20.63 19.0 0.82 MOSFIRE 0.117 0.072 0.151 6 J0038 −1527 i 7.034 19.69 ±0 . 07 19.33 ±0 . 05 7200 22.86 21.67 0.83 MOSFIRE 0.224 0.242 0.207 7 J0411 −0907 j 6.826 20.02 ±0 . 14 19.56 ±0 . 21 5760 12.58 13.53 0.66 NIRES 0.028 0.017 0.036 8 J0319 −1008 k 6.827 20.98 ±0 . 24 —18 720 7.61 8.48 0.66 NIRES 0.004 0.003 0.005 9 J0410 −0139 l 7.0 —20.39 m 9360 5.36 6.09 0.83 NIRES 0.001 0.001 0.001 10 J0038 −0653 n 7.1 —19.61 o 5760 13.29 12.65 0.90 MOSFIRE 0.028 0.036 0.024 Notes. a Median values calculated from the unmasked regions on the bottom panels of Figs 1–6. b We disco v er that the noise vectors produced by the PYPEIT pipeline could be slightly o v er- and underestimated. We correct for this by applying correction factors to the noise vectors of the quasars as described in the text. c Ef fecti ve weights of quasars towards the correlation function in each redshift bin (see definition following equation 4 ). d Wang et al. ( 2021 ). e Ba ˜ nados et al. ( 2018 ). f Yang et al. ( 2020 ). g Mortlock et al. ( 2011 ). h Yang et al. ( 2019 ). i Wang et al. ( 2018 ). j Pons et al. ( 2019 ); Wang et al. ( 2019 ). k Yang et al. ( 2019 ). l Ba ˜ nados et al. ( 2024 ). m Pri v ate communication. n Yang et al. (in preparation). o Pri v ate communication. 3 DATA REDUCTION We reduce our spectra with PYPEIT , 4 a PYTHON package for semi-automated reduction of astronomical slit-based spectroscopy (Prochaska et al. 2020a , b ). After performing standard calibration steps and flux calibration on each exposure, we co-add 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 co-add them all onto a common grid with a pixel size of 40 km s −1 . Finally, we restrict our analyses to a minimum wavelength of 1.95 μm, corresponding to z Mg II = 5 . 96. We fit for the quasar continuum using a B-spline with a breakpoint at ev ery-60th pix el (amounting to a breakpoint at every 2400 km s −1 spacing) and we mask absorbers identified by-eye before fitting the continuum. We use the same breakpoint spacing for all the quasars. Figs 1 –4 show the spectra before and after continuum normalization for representative members of our data set. For each quasar, the grey shaded regions in the top panel indicate regions that are masked before continuum fitting due to visually detected 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 QSO (grey), low SNR due to telluric absorption (blue), and detected CGM absorbers (red; see Section 4); these are excluded from the correlation function analyses. We also mask pixels that fall within a quasar proximity zone (PZ), as these pixels will be affected by the quasar own ionizing radiation and produce biased results on the IGM neutral fraction. For the PZ mask, we remo v e pix els within rest-frame v = 7643 km s −1 blueward from the quasar Mg II emission line, where we define the Mg II emission line to be the mid-point of the Mg II λ2796 Åand 4 https:// pypeit.readthedocs.io/ en/ latest/ λ2804 Å doublet at 2800 Å. Our choice of v = 7643 km s −1 follows Bosman et al. ( 2022 ), who consider pixels at < 1185 Åfor their Ly α optical depth study, which results in a PZ size of 7643 km s −1 , 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 ; Satya v olu et al. 2023 ), and even in the presence of evolution, our adopted value of 8.64 proper Mpc is rather conserv ati ve gi ven that the typical PZ size for z ∼6 quasars ranges from 2 −7 proper Mpc (Eilers et al. 2017 , 2020 ; Satya v olu et al. 2023 ).The unshaded regions are used for the correlation function measurements, and they range from z Mg II = 5 . 96 to 7 . 42, with z Mg II , med = 6 . 469, where z Mg II ≡( λobs / 2800 Å) −1 is the pixel redshift and z Mg II , med is the median value. Table 1 includes the path length of each quasar that contributes towards the correlation function measurement, where the total path length of our data set is z = 9 . 75 or X = 47 . 1. X is the absorption path length, where X( z) =  z 0 (1 + z  ) 2 H 0 /H ( z  ) d z  (Bahcall & Peebles 1969 ). We compute X qso = X( z i ) −X( z f ) for each quasar and take the sum as the total absorption path length quoted abo v e. We disco v er that the noise in the quasar spectra produced by PYPEIT are slightly o v er- and underestimated 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 1 −F norm σnorm after applying all masks (including CGM masks described in Section 4 ) and obtain the median absolute deviation 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, and 1.051 to the noise vectors of quasars numbers 1–10 (see Table 1 ), respectively. After applying these corrections, the median SNR of our data set is 13.7. Fig. 5 zooms-in on certain regions of the Mg II forest from the Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 226 S. Tie et al. MNRAS 535, 223–246 (2024) Figure 1. Three example spectra from our data set, J0313 −1806 (top), J1342 + 0928 (middle), and J1007 + 2115 (bottom). In the top panels of each quasar, the red line is the fitted quasar continuum, the faint black line is the spectral noise, and the shaded regions indicate masks obtained from PYPEIT reduction and from visually detected absorbers, which are applied before continuum fitting. The bottom panel of each quasar shows the continuum-normalized spectrum. Also shown is the telluric spectrum in blue. The multicoloured shaded regions indicate masks due to PYPEIT reduction (if available, in green), proximity to quasars and where z > z QSO (grey), low SNR due to telluric absorption (blue), and detected CGM absorbers (red; see Section 4), where we use the unshaded regions for the correlation function measurement and inference. The vertical dashed lines indicate the quasar redshift. Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 227 MNRAS 535, 223–246 (2024) Figure 2. Same as Fig. 1 , but for J1120 + 0641 (top), J0252 −0503 (middle), and J0038 −1527 (bottom). Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 228 S. Tie et al. MNRAS 535, 223–246 (2024) Figure 3. Same as Fig. 1 , but for J0411 −0907 (top), J0319 −1008 (middle), and J0410 −0319 (bottom). Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 229 MNRAS 535, 223–246 (2024) Figure 4. Same as Fig. 1 , but only showing the normalized spectrum of J0038 −0653 (upon permission by Yang et al. in preparation). normalized spectra, where the typical flux fluctuations given our data quality is ∼10 per cent, or higher towards the redder part of the spectrum. Typical IGM absorbers are expected to produce percentlevel fluctuations (see fig. 2 of H2021 ), so attempting to detect them at this SNR level will be challenging. 4 MASKING MG II ABSORBERS OF THE CGM 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 (fig. 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 (fig. 14 of H2021 ). Fig. 6 (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 equi v alent width, W λ, pix ≡(1 −F ) λ, where λ is the spectral pixel width at rest-frame Mg II 2796 Å, which depends on the spectral sampling. Using our fixed common grid with d v = 40 km s −1 , λ ∼0 . 37 Å. While a flux cut helps remo v e 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 χfield, χ( v) =  [1 −F ( v  )] W ( | v −v  | )d v   σ2 F ( v  ) W 2 ( | v −v  | )d v  , (1) where σ2 F is the flux variance (noise vector) of the data and W ( v) = 1 −e −τ( v) is a matched filter with τ( v) being the optical depth of a Mg II doublet with a Gaussian velocity distribution, assuming N Mg II = 10 13 . 5 cm −2 and Doppler parameter b = √ 2 ×FWHM / 2 . 35, in which we use the FWHM (full width at half-maximum) 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 −1 , considering that we rebin all spectra to a fixed 40 km s −1 grid. Note that F is the continuum-normalized flux and σ2 F is the continuum-normalized v ariance. Gi ven the χspectrum, absorbers are identified by their extreme χvalue. We mask out pixels abo v e χ> 3 as well as pixels ±300 km s −1 around each masked pixel to account for their absorption wings. Fig. 6 (right panel) shows the χ distribution of the data. In both panels, in addition to the distribution of each quasar (coloured 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 χPDF sho ws e vidence for fluctuations that are not pure noise at χ> 3, since the noise PDF drops off at χ∼3. Both PDFs also peak at some characteristic values as expected (see fig. 13 of Hennawi et al. 2021 ), although the 1 −F PDF does not flatten off at large values and the plateau at large χis weak in the χPDF. In summary, we mask (1) pixels with 1 −F > 0 . 3 and (2) pixels with χ> 3 and those ±300 km s −1 around each masked pixel. Fig. 7 shows the masked absorbers for the three highestz quasars in our data set (we show the rest in Figs B1, B2 and B3 in Appendix B for brevity). The green and magenta horizontal lines indicate our flux and χcut-of fs, respecti v ely. We disco v er new absorber candidates in J0313 −1806, J0038 −1527, and newqso1, and reco v er known absorbers at z = 6 . 84 (Simcoe et al. 2020 ) and z = 6 . 27 in the spectrum of J1342 + 0928. We reco v er 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 ∼2 and it is masked manually. These absorbers are indicated by the blue ticks in Fig. B1 in Appendix B. 5 RESULTS 5.1 Autocorr elation measur ement of the Mg II for est We perform the correlation function measurement o v er three redshift bins: (i) 5 . 96 ≤z Mg II < 6 . 469 (‘lowz’ bin, where z Mg II , med = 6 . 235). (ii) 6 . 469 ≤z Mg II ≤7 . 421 (‘highz’ bin, where z Mg II , med = 6 . 715). (iii) 5 . 96 ≤z Mg II ≤7 . 421 (‘allz’ bin, with z Mg II , med = 6 . 469). Recall that z Mg II ≡( λobs / 2800 Å) −1 is the redshift of the pixels that contributed to the correlation function and z Mg II , med is the median value. We first compute the relative flux fluctuation for each Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 230 S. Tie et al. MNRAS 535, 223–246 (2024) Figure 5. Zooming in on regions of the Mg II forest, where the top three panels are for the z > 7.5 quasars and the rest for the z ? 7.5 quasars in our dataset. Red lines are the normalized noise, which are offset by some amount from zero for visualization purpose, while the shaded regions are masked (see Section 2 ), although here we do not color code them like we did for Figs 1 –3 . continuum-normalized quasar spectrum ( F ) as δf ≡F −¯ F ¯ F , (2) where ¯ F is the mean flux of the entire data set in each redshift bin. We compute ¯ F as the weighted average of the normalized flux for the entire data set using the inverse variance of the corrected measurement noise (see below). We use the suboptimal unbiased estimator for the correlation function as the weighted averaged of all pixel pairs i, j (Slosar et al. 2011 ) ξ( v) =  i,j w i w j δi δj  i,j w i w j , (3) where  i,j means summing o v er pix el pairs i, j separated by the velocity lag v. Rewriting the sum as a sum o v er all pix el pairs from each quasar k, the correlation function can also be computed as the Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 231 MNRAS 535, 223–246 (2024) Figure 6. Masking schemes to remo v e CGM absorbers based on the PDF of the flux decrement (left) and significance (right). The coloured lines are the distributions for each quasar and the black lines are the distributions from all quasars. The dark blue lines are the distributions of Gaussian noise drawn from the corrected noise array from all quasars (see Section 2 ). To remo v e CGM absorbers, we mask all pixels with 1 −F > 0 . 3 (shaded region in left panel), as well as those with χ> 3 (shaded region in right panel) and ±300 km s −1 around each masked pixel. weighted average over the individual quasar correlation functions ξ( v) =  k  i,j w ik w jk δik δjk  k  i,j w ik w jk =  k u k ξk ( v)  k u k =  k u eff k ξk ( v) (4) where u k ≡ i,j w ik w jk and u eff k ≡u k  k u k is the result of normalizing u k (i.e.  k u eff k = 1), and ξk ( v) =  i,j w ik w jk δik δjk /u k is the weighted estimate of the correlation function of each quasar k. Along similar lines, ¯ F is computed as a weighted average via ¯ F =  ki w ik F ik  ki w ik , (5) where F ik is the continuum normalized flux of the ith spectral pixel in the kth quasar. The values of ¯ F are 1.0015, 1.001, and 1.0018 for the allz , highz , and lo wz bins, respecti vely, after masking out CGM absorbers. For the weights, we use the inverse variance of the corrected measurement noise, w ik = 1 /σ 2 ik , where σ2 ik is the error of the continuum normalized flux. While one should technically also include the variance due to the intrinsic fluctuations in the Mg II forest, we neglect that term here as we are noise-dominated – based on H2021 , for an IGM model with [Mg/H] = −3 . 50 ( −4 . 50), the typical fluctuation in the Mg II forest is ∼5 per cent (2 per cent), which requires a SNR of at least 20 (50), whereas the median SNR of our data set is 13.7 (considering only the top five highest SNR quasars, the mean SNR is still ∼18). We measure the correlation function of the Mg II forest from v min = 10 to v max = 3500 km s −1 . Given that the correlation function is expected to rise at small scales and peak at the doublet separation (768 km s −1 ), and flatten at large velocity separations, a fine bin size is needed to capture the small-scale rise and the correlation function peak, whereas a coarse bin size is sufficient at large scales. Therefore, to save computation time during the forward modelling step (Section 5.2 ), we adopt a heterogenous binning where we use a bin size of 80 km s −1 at v < 1500 km s −1 and a bin size of 200 km s −1 at v ≥1500 km s −1 . Since we are using the weighted average over all quasars as the measurement for the data set, giv en the hetere geneous SNR and path length of the quasars, they do not contribute equally to the final measurement. Table 1 sho ws the ef fecti ve weights, u eff k , which quantifies the contribution of each quasar to the total correlation function for each redshift bin. For the allz bin, J0313 −1806 has the highest weight (33 per cent), followed by J0038 −1527 (22 per cent), J1342 + 0928 (14 per cent), and J0252 −0503 (12 per cent). Similarly, the highz bin measurement is dominated by J0313 −1806, J0038 −1527, and J1342 + 0928, while the lowz bin measurement is dominated by J0313 −1806, J0038 −1527, and J0252 −0503. Given the unequal contributions of individual quasars to the correlation function measurement, we can compute a so-called ef fecti ve sample size (Kish 1995 ), which is an estimate of the sample size required for a random process where all data points carry uniform weight, to achieve the same level of variance as a random process with non-uniform weights. The ef fecti ve sample size is given by n eff = (  i n u eff i ) 2  i n ( u eff ) 2 , where u eff i is the weight of the ith data point. We find values of n eff = 4.9, 4.3, and 5.1 for the allz , highz , and lo wz bins, respecti vely, which is significantly less than the actual number of quasars in the data set (10). Figs 8 –10 show the correlation functions before and after we filter out CGM absorbers (see Section 4 ) in the three redshift bins. The coloured points are the measurements for each quasar and the black points are the average from all quasars (equation 4 ). In the left panels of the figures, we detect a strong peak in the correlation function at the Mg II doublet separation (768 km s −1 ) that is due to strong CGM absorbers, as well as a rise at small velocity lags due to the small-scale velocity structure of the absorbers. The detection of the Mg II peak at the correct location provides a sanity check of our methodology. We roughly estimate the error bars on the correlation function in the left panels as the dispersion between the individual correlation functions divided by  N q , with N q = 10 being the total number of quasars in our data set. We note that this error estimation is technically incorrect as our final correlation function measurement is a weighted one, but since we are only using these errors for visual purposes and given that Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 238 S. Tie et al. MNRAS 535, 223–246 (2024) 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 (Gardner et al. 2006 ), we expect to place stronger constraints on primordial cosmic enrichment and reionization. For instance, using a sample of eight z ≥7 quasar spectra with a SNR/pixel of ∼76 from JWST /NIRSpec, 5 we expect to jointly constrain [Mg/H] within 1 σprecision of 0.04 dex and x H I within 1 σprecision of 10 per cent, as well as place an upper limit of [Mg/H] < −3 . 94 on the IGM enrichment at 95 per cent credibility. Finally, by combining very highz measurements from JWST that is probing the neutral IGM with ground-based measurements of, for example, the C IV forest, at z ∼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. ACKNOWLEDGEMENTS 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 , 6 a community-developed core PYTHON package for Astronomy (Astropy Collaboration 2013 , 2018 ) and SCIPY 7 (Virtanen et al. 2020 ). DATA AVAILABILITY The data underlying this article will be shared on reasonable request to the corresponding author. REFERENCES Almgren A. S. , Bell J. B., Lijewski M. J., Luki ´ c Z., Van Andel E., 2013, ApJ , 765, 39 Asplund M. , Grevesse N., Sauval A. J., Scott P., 2009, ARA&A , 47, 481 Astropy Collaboration , 2013, A&A , 558, A33 Astropy Collaboration , 2018, AJ , 156, 123 Bahcall J. N. , Peebles P. J. E., 1969, ApJ , 156, L7 Ba ˜ nados E. et al.. 2024, preprint ( arXiv:2407.07236 ) Ba ˜ nados E. et al., 2018, Nature , 553, 473 Becker G. D. et al., 2019, ApJ , 883, 163 Becker G. D. , Rauch M., Sargent W. L. W., 2009, ApJ , 698, 1010 Bosman S. E. I. et al., 2022, MNRAS , 514, 55 Bosman S. E. I. , Becker G. D., Haehnelt M. G., Hewett P. C., McMahon R. G., Mortlock D. J., Simpson C., Venemans B. P., 2017, MNRAS , 470, 1919 Chen S.-F. S. et al., 2017, ApJ , 850, 188 Christensen L. et al., 2023, A&A , 680, A82 5 https:// www.stsci.edu/ jwst/ science-e x ecution/program-information?id= 3526 6 http://www .astropy .org 7 https://scipy.org Codoreanu A. , Ryan-Weber E. V., Garc ´ ıa L. ´ A., Crighton N. H. M., Becker G., Pettini M., Madau P., Venemans B., 2018, MNRAS , 481, 4940 Cooper T. J. , Simcoe R. A., Cooksey K. L., Bordoloi R., Miller D. R., Furesz G., Turner M. L., Ba ˜ nados E., 2019, ApJ , 882, 77 D’Odorico V. et al., 2016, MNRAS , 463, 2690 D’Odorico V. et al., 2022, MNRAS , 512, 2389 Davies F. B. , Furlanetto S. R., 2022, MNRAS , 514, 1302 Davies R. L. et al., 2023a, MNRAS , 521, 289 Davies R. L. et al., 2023b, MNRAS , 521, 314 Eilers A.-C. et al., 2020, ApJ , 900, 37 Eilers A.-C. , Davies F. B., Hennawi J. F., Prochaska J. X., Luki ´ c Z., Mazzucchelli C., 2017, ApJ , 840, 24 Ferrara A. , 2016, in Mesinger A.ed., Astrophysics and Space Science Library Vol. 423, Understanding the Epoch of Cosmic Reionization: Challenges and Progress. Springer International Publishing, Switzerland, p. 163 F oreman-Macke y D. , Hogg D. W., Lang D., Goodman J., 2013, PASP , 125, 306 Fox A. J. , 2011, ApJ , 730, 58 Gardner J. P. et al., 2006, Space Sci. Rev. , 123, 485 Hennawi J. F. , Davies F. B., Wang F., O ˜ norbe J., 2021, MNRAS , 506, 2963 (H2021) Jaacks J. , Thompson R., Finkelstein S. L., Bromm V., 2018, MNRAS , 475, 4396 Kara c ¸aylı N. G. et al., 2023, MNRAS , 522, 5980 Kirihara T. , Hase ga wa K., Umemura M., Mori M., Ishiyama T., 2020, MNRAS , 491, 4387 Kish L. , 1995, Surv e y Sampling. Wiley, New York Liu B. , Bromm V., 2020, MNRAS , 497, 2839 Luki ´ c Z. , Stark C. W., Nugent P., White M., Meiksin A. A., Almgren A., 2015, MNRAS , 446, 3697 McLean I. S. et al., 2010, in McLean I. S., Ramsay S. K., Takami H.eds, Proc. SPIE Conf. Ser. Vol. 7735, Ground-based and Airborne Instrumentation for Astronomy III. p. 77351E McLean I. S. et al., 2012, in McLean I. S., Ramsay S. K., Takami H.eds, Proc. SPIE Conf. Ser. Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV. p. 84460J Mesinger A. , Furlanetto S., Cen R., 2011, MNRAS , 411, 955 Meyer R. A. , Bosman S. E. I., Kakiichi K., Ellis R. S., 2019, MNRAS , 483, 19 Miralda-Escud ´ e J. , Rees M. J., 1998, ApJ , 497, 21 Mortlock D. J. et al., 2011, Nature , 474, 616 Oh S. P. , 2002, MNRAS , 336, 1021 Oppenheimer B. D. , Dav ´ e R., 2006, MNRAS , 373, 1265 Oppenheimer B. D. , Dav ´ e R., Finlator K., 2009, MNRAS , 396, 729 Planck Collaboration VI , 2020, A&A , 641, A6 Pons E. , McMahon R. G., Simcoe R. A., Banerji M., Hewett P. C., Reed S. L., 2019, MNRAS , 484, 5142 Prochaska J. X. et al., 2020a, pypeit/PypeIt: Release 1.0.0 , Zenodo, doi:10.5281/zenodo.3743493 Prochaska J. X. et al., 2020b, J. Open Source Softw. , 5, 2308 Ryan-Weber E. V. , Pettini M., Madau P., 2006, MNRAS , 371, L78 Satya v olu S. et al., 2023, MNRAS , 522, 4918 Schaye J. , Aguirre A., Kim T.-S., Theuns T., Rauch M., Sargent W. L. W., 2003, ApJ , 596, 768 Sebastian A. M. et al., 2024, MNRAS , 530, 1829 Simcoe R. A. et al., 2011, ApJ , 743, 21 Simcoe R. A. , 2011, ApJ , 738, 159 Simcoe R. A. , Onoue M., Eilers A.-C., Banados E., Cooper T. J., Furesz G., Hennawi J. F., Venemans B., 2020, preprint ( arXiv:2011.10582 ) Slosar A. et al., 2011, J. Cosmol. Astropart. Phys. , 2011, 001 Songaila A. , 1997, ApJ , 490, L1 Songaila A. , Cowie L. L., 1996, AJ , 112, 335 Tie S. S. , Hennawi J. F., Kakiichi K., Bosman S. E. I., 2022, MNRAS , 515, 3656 Vernet J. et al., 2011, A&A , 536, A105 Virtanen P. et al., 2020, Nat. Methods , 17, 261 W ang F . et al., 2018, ApJ , 869, L9 W ang F . et al., 2019, ApJ , 884, 30 Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 239 MNRAS 535, 223–246 (2024) W ang F . et al., 2021, ApJ , 907, L1 Wilson J. C. et al., 2004, in Moorwood A. F. M., Iye M.eds, Proc. SPIE Conf. Ser. Vol. 5492, Ground-based Instrumentation for Astronomy. p. 1295 Yang J. et al., 2019, AJ , 157, 236 Yang J. et al., 2020, ApJ , 897, L14 APPENDIX A During our investigations,we discovered that our likelihood yields MCMC samples and best-fitting models that appear visually inconsistent with the actual measurements and error bars. The source of this discrepancy is correlated velocity bins with negati ve v alues 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 remo v e sensitivity towards the signal. Figs A1 and A2 show the correlation functions in one velocity lag bin plotted against the values of all other velocity bins for the allz bin. The black points are measured from forward models with x H I = 0.50 and [Mg/H] = −4 . 50, where the ellipses denote their 1 σ, 2 σ, and 3 σcontours. Measurements from the real data that fall outside of the 3 σellipse are marked as red, while measurements that fall within the 3 σellipse are marked as green. We mask the velocity bin where 20 or more of the real measurements fall outside the 3 σellipse; Fig. A1 shows two examples of bad bins that are masked, while Fig. A2 shows two examples of good bins that we retain. In conclusion, we mask d v = 370, 450, 1170, and 1490 km s −1 in the allz bin, d v = 50 and 1090 km s −1 in the highz bin, and d v = 50, 130, 290, 370, 450, 770, 930, and 1490 km s −1 in the lowz bin (see Figs 13 –15 ). We generate the forward-modelled data sets 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-fitting models and upper limits are therefore obtained from the masked data and masked covariance matrices. Figure A1. Values of the allz redshift bin correlation function at d v = 450 km s −1 (left) and d v = 1170 km s −1 (right) plotted against the values of all other bins. The black points are computed from forward models with x HI = 0.50 and [Mg/H] = −4 . 50 and the ellipses denote their 1 σ, 2 σ, and 3 σcontours. Red points are measurements from the real data that fall outside of the 3 σellipse, where these two bins are examples of catastrophic failure as the real data fall outside the 3 σellipse of all the other bins. Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 240 S. Tie et al. MNRAS 535, 223–246 (2024) Figure A2. Same as Fig. A1 , but for d v = 530 km s −1 (left) and d v = 1410 km s −1 (right). Red (green) points are measurements from the real data that fall outside of (within) the 3 σellipse. In these two bins, most the real data fall within the 3 σellipse of the other bins. APPENDIX B Here, we show the masked absorbers for the rest of the quasars in our data set that are not shown in Section 4 , which are Figs B1, B2 and B3 . APPENDIX C Here, we show the forward-modelled spectra for the rest of the quasars in our data set that are not shown in Section 5 , assuming an IGM model with ( x H I , [Mg / H] ) = (0 . 50 , −4 . 50). These are shown in Figs C1, C2 and C3 . Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 241 MNRAS 535, 223–246 (2024) Figure B1. Continued from Fig. 7 , for J1120 + 0641 (top), J0252 −0503 (middle), and J0038 −1527 (bottom). For J1120 + 0641, we also manually mask the locations of weak absorbers found by Bosman et al. ( 2017 ), as indicated by the blue ticks. Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 242 S. Tie et al. MNRAS 535, 223–246 (2024) Figure B2. Continued from Fig. 7 , for J0411 −0907 (top), J0319 −1008 (middle), and J0410 −0139 (bottom). Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 243 MNRAS 535, 223–246 (2024) Figure B3. Continued from Fig. 7 for J0038 −0653. Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 244 S. Tie et al. MNRAS 535, 223–246 (2024) Figure C1. Continued from Fig. 12 , for J1120 + 0641 (top), J0252 −0503 (middle), and J0038 −1527 (bottom). Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 Mg II forest in the EoR 245 MNRAS 535, 223–246 (2024) Figure C2. Continued from Fig. C1 , for J0411 −0907 (top), J0319 −1008 (middle), and J0410 −0139 (bottom). Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025 246 S. Tie et al. MNRAS 535, 223–246 (2024) Figure C3. Continued from Fig. C1 for J0038 −0653. This paper has been typeset from a T E X/L A T E X file prepared by the author. © 2024 The Author(s). Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License ( https://cr eativecommons.or g/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. Downloaded from https://academic.oup.com/mnras/article/535/1/223/7774407 by guest on 26 March 2025