The Completed SDSS-IV Extended Baryon Oscillation Spectroscopic Survey: Baryon Acoustic Oscillations with Lyα Forests

This analysis benefits from more than 10 years of cosmological observations from SDSS (York et al. 2000) on the 2.5 m Sloan Foundation telescope (Gunn et al. 2006) at the Apache Point Observatory. Most of the tracer quasar spectra and the entirety of the background quasar spectra were gathered during SDSS-III (Eisenstein et al. 2011) in the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013) and in the eBOSS component of SDSS-IV. A small fraction of tracer quasars were observed during SDSS-I and SDSS-II and are found in the seventh data release (Schneider et al. 2010).

The DR16 data used in this paper (Ahumada et al. 2020) contain the complete 5 yr BOSS sample (Fall 2009–Spring 2014) and the 5 yr eBOSS sample (2014 July–2019 March), including its pilot program SEQUELS (Myers et al. 2015), which began in the final year of SDSS-III and concluded in the first year of SDSS-IV. The quasar target selection algorithms for BOSS are summarized in Ross et al. (2012), and those for eBOSS are summarized in Myers et al. (2015). An additional sample of targets was observed as part of the Time-Domain Spectroscopic Survey (TDSS; Morganson et al. 2015; Ruan et al. 2016) and the SPectroscopic IDentification of ERosita Sources survey (SPIDERS; Dwelly et al. 2017). An example high-S/N quasar spectrum is shown in Figure 1. The DR16 footprint is presented in Figure 2, where the color map reflects the statistical improvement between DR12 and DR16 for the autocorrelation of Lyα pixels.

Zoom In Zoom Out Reset image size

The tracer and background quasars are taken from the DR16 quasar catalog (DR16Q; Lyke et al. 2020), which includes the DR7 quasar catalog (Schneider et al. 2010). Most of the DR7 quasars at z > 2.1 were reobserved in BOSS or eBOSS; they are all included in DR16Q. For the few quasars that were not reobserved, we use them as tracer quasars, but not as background quasars, since their spectra were processed with another pipeline.

Each spectrum from BOSS or eBOSS was acquired using one of the two spectrographs, with wavelength coverage ranging from 360 to 1000 nm (Smee et al. 2013). Each spectrograph views 500 fibers, roughly 450 of which are dedicated to science targets and the remaining to standard F stars for flux calibration and to empty sky locations for sky background subtraction.

The data were processed by version v5_13_0 of the eBOSS pipeline. The reduction is organized in two steps. The pipeline initially extracts the 2D raw data into a 1D flux-calibrated spectrum. During this procedure, the spectra are wavelength and flux calibrated and the individual exposures of one object are co-added into a rebinned spectrum with ${\rm{\Delta }}\mathrm{log}(\lambda )={10}^{-4}$. The spectra are then classified as STAR, GALAXY, or QSO, and their redshift is estimated.

DLA systems and broad emission lines (BALs) were identified in the quasar spectra. The details of these searches are presented in Lyke et al. (2020). The BAL search used a procedure similar to that used by Guo & Martini (2019). DLAs were identified using a neural network as described in Parks et al. (2018).

The current DR16 pipeline slightly differs from its last public release (DR14; Abolfathi et al. 2018). Details of the changes can be found in Ahumada et al. (2020). We focus on two relevant changes for the Lyα analysis here. The first concerned improvements in the background estimates on the spectral CCD images. This allowed for more accurate determination of spectral densities f(λ). The second improves modeling of the spectra of calibration stars by using a new set of stellar templates. The set was produced for the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al. 2016) pipeline and was provided to eBOSS. These templates are able to reduce residuals in flux calibration by improving the modeling of spectral lines in F stars.

For the purposes of measuring Lyα correlations, we rebin three original pipeline ${\rm{\Delta }}{\mathrm{log}}_{10}(\lambda )\sim {10}^{-4}$ spectral pixels into one ${\rm{\Delta }}{\mathrm{log}}_{10}(\lambda )\sim 3\times {10}^{-4}$ "analysis pixel." Throughout this paper, the use of the word "pixel" refers to these rebinned analysis pixels, unless otherwise stated.

The DR16Q catalog provides up to four estimates of quasar redshift based on the position of the broad emission lines. The absence of visible narrow spectral lines that would be associated with the quasar host galaxies means that individual quasar redshifts have statistical uncertainties of order $100\,\mathrm{km}\,{{\rm{s}}}^{-1}$ and the different estimators have somewhat different systematic differences. As described in Appendix B, we studied in detail these differences using synthetic quasar spectra. Based on these studies, we choose to employ an estimator that does not use any spectral information in the vicinity of the Lyα emission line or at shorter wavelengths. This choice of rest-frame wavelength coverage mitigates redshift-dependent systematic errors due, for example, to evolution in the mean Lyα opacity.

The redshift range of useful quasars is set by the goal of measuring the quasar–Lyα cross-correlations up to separations of $200\,{h}^{-1}\,\mathrm{Mpc}$. The lowest-redshift Lyα pixel has z = 1.96, meaning that quasars with redshifts z > 1.77 can be used to sample the cross-correlations up to the maximum separation. We impose a maximum redshift of z = 4, beyond which the number of tracers is insufficient for a useful correlation measurement. Our final sample is thus composed of ${\rm{341,468}}$ tracer quasars, as summarized in Table 1. The redshift distribution is shown in the left panel of Figure 3.

Zoom In Zoom Out Reset image size

Table 1. Definitions and Statistics for the Matter Density Tracers

Tracer	${\lambda }_{\mathrm{RF}}^{\min }$	${\lambda }_{\mathrm{RF}}^{\max }$	${z}_{{\rm{q}},\min }$	${z}_{{\rm{q}},\max }$	Nq	Npix	$\overline{S/N}$
	$(\mathrm{nm})$	$(\mathrm{nm})$				$({10}^{6})$
Quasar			1.77	4	341,468
$\mathrm{Ly}\alpha \left(\mathrm{Ly}\beta \right)$	92	102	2.65	4	69,656	8.4	1.88
$\mathrm{Ly}\alpha \left(\mathrm{Ly}\alpha \right)$	104	120	2.10	4	210,005	34.3	2.56

Note. The rest-frame wavelength intervals are shown in the second and third columns, the associated redshift interval in the fourth and fifth columns, the number of background or tracer quasars in the sixth column, the total number of spectral pixels in the seventh column, and, in the last column, the mean S/N of the pixels, ${\rm{S}}/{\rm{N}}=\langle (\delta +1)* \sqrt{w}\rangle$, where the weight, w, is given by the inverse of Equation (4).

Download table as:  ASCIITypeset image

As was done in DR14 (dSA19, B19), we use two different spectral regions shown in Figure 1: the "Lyα region," ${\lambda }_{\mathrm{RF}}\,\in [104,120]\,\mathrm{nm}$, and the "Lyβ region," ${\lambda }_{\mathrm{RF}}\in [92,102]\,\mathrm{nm}$. The Lyα region is defined between the Lyβ+O vi quasar emission line and the Lyα emission line. Defined as such, we ensure that this region only has Lyα absorption and weak metal absorption (see Pieri et al. 2014). This is the same definition as in previous analyses. The Lyβ+O vi and Lyα emission-line wings are avoided because of their higher quasar-to-quasar diversity.

The Lyβ spectral region covers most of the range between the Lyman limit, ${\lambda }_{\mathrm{Ly}-\mathrm{limit}}=91.18\,\mathrm{nm}$, and the Lyβ+O vi quasar emission line. This increase in range over the $[97.4,102]\,\mathrm{nm}$ in DR14 to [92,102] nm increases the number of spectral pixels in the region by a factor of two.

This analysis uses pixels in the observed wavelength range λ ∈ [360,600] nm. The upper limit comes from the requirement z_q  < 4, which is motivated by the decreasing number of quasars and the increasing contamination from sky emission lines. The lower limit is motivated by high atmospheric absorption in the UV and the low CCD throughput at the lowest wavelengths. From the previously defined rest-frame wavelength range of the two spectral regions and the observed wavelength range, the minimum background quasar redshift is z_q  > 2 for the Lyα region and z_q  > 2.53 for the Lyβ region. This selection gives us ${\rm{265,328}}$ background quasars for the Lyα region and ${\rm{103,080}}$ background quasars for the Lyβ region. We removed BALs with BALPROB > 0.9 in the DR16Q Catalog (Lyke et al. 2020). This reduced the numbers to ${\rm{249,814}}$ and ${\rm{97,124}}$.

As in dSA19 and B19, all good spectral observations of a given background quasar are stacked to give a higher-S/N spectrum. The use of only good BOSS or eBOSS observations reduces the number of forests to ${\rm{244,514}}$ and ${\rm{96,110}}$ for the Lyα and Lyβ regions, respectively.

A spectral region is used only if it contains at least 50 pixels, since short lines of sight may be overfitted, leading to erasing of structure in a correlated manner. This requirement shifts the minimum quasar redshift to z_q  > 2.10 for the Lyα region and z_q  > 2.65 for the Lyβ region, leaving ${\rm{214,542}}$ and ${\rm{71,602}}$ background quasars, respectively. The fit of the continuum (Equation (2)) fails for a few low-S/N or outlier spectral regions. The final sample has ${\rm{210,005}}$ and ${\rm{69,656}}$ background quasars, respectively. We summarize the sample in Table 1.

Spectra and noise variance estimates are corrected for Galactic extinction using the dust map of Schlegel et al. (1998). This had a very small effect in the final results, since extinction corrections are very smooth and can be absorbed in the continuum fitting.

Following B17 and du Mas des Bourboux et al. (2019), we correct for small residual flux calibration errors from the eBOSS pipeline. To do so, we designate a "calibration spectral region" on the red side of the Mg ii emission line, ${\lambda }_{\mathrm{RF}}\,\in [290,312]\,\mathrm{nm}$. We compute the mean flux as a function of observed wavelength $\overline{{f}_{\mathrm{calib}.}}({\lambda }_{i})$, correcting for the shape of the quasar continuum (see below). The new flux is then defined to be $f({\lambda }_{i})\to f({\lambda }_{i})/\overline{{f}_{\mathrm{calib}.}}({\lambda }_{i})$, and its inverse variance ${ivar}({\lambda }_{i})\to {ivar}({\lambda }_{i})\times {\overline{{f}_{\mathrm{calib}.}}}^{2}({\lambda }_{i})$. We find that the correction varies by at most 3%. This process captures residual errors from incomplete modeling of standard F stars and sky emission lines, as well as the Ca ii H and K absorption of the Milky Way.

The next step is to mask out spectral intervals, in observed wavelength, where the variance increases sharply owing to unmodeled emission lines from the sky. However, a balance has to be found between pixel quantity and pixel quality. We use the previously defined calibration spectral region to compute the variance as a function of observed wavelength. Comparing to the intrinsic variance of LSS in the Lyα and Lyβ regions, we mask the following six intervals: λ ∈ [404.30, 405.13], λ ∈ [435.31, 436.51], λ ∈ [545.68, 546.73], λ ∈ [557.05,559.00], λ ∈ [588.33, 590.33], and λ ∈ [629.48, 631.68] nm. This mask produces less data loss compared to previous studies, in large part due to the improved sky calibration in the latest eBOSS pipeline. In addition, we mask both Ca ii H and K absorption of the Milky Way, λ ∈ [393.01, 393.82] and λ ∈ [396.55, 397.38] nm. Even though we do not observe a significant excess of variance in these regions, masking them avoids spurious large-scale angular correlated absorption between lines of sight.

For each spectral region in each line of sight, q, the flux transmission field, δ_q (λ), at the observed wavelength, λ, is obtained from the ratio of the observed flux, f_q (λ), to the mean expected flux, $\overline{F}(z){C}_{q}(\lambda )$:

$$\begin{eqnarray}&&{\delta }_{q}(\lambda )=\displaystyle \frac{{f}_{q}(\lambda )}{\overline{F}(\lambda ){C}_{q}(\lambda )}-1.\end{eqnarray} \tag{ 1 }$$

Here C_q (λ) is the unabsorbed quasar continuum and $\overline{F}(\lambda )$ is the mean transmission. Their product is taken to be a universal function of the rest-frame wavelength, $\overline{C}({\lambda }_{\mathrm{RF}})$, corrected by a first-degree polynomial in $\mathrm{log}\lambda$:

$$\begin{eqnarray}&&\overline{F}(\lambda ){C}_{q}(\lambda )=\overline{C}({\lambda }_{\mathrm{RF}})\left({a}_{q}+{b}_{q}{\rm{\Lambda }}\right).\quad {\rm{\Lambda }}\equiv \mathrm{log}\lambda .\end{eqnarray} \tag{ 2 }$$

The coefficients (a_q , b_q ) are fitted separately for the Lyα and Lyβ spectral regions of each line of sight by maximizing the likelihood function

$$\begin{eqnarray}&&2\mathrm{ln}L=-\displaystyle \sum _{i}\displaystyle \frac{{[{f}_{i}-\overline{F}{C}_{q}({\lambda }_{i},{a}_{q},{b}_{q})]}^{2}}{{\sigma }_{q}^{2}({\lambda }_{i})}-\mathrm{ln}[{\sigma }_{q}^{2}({\lambda }_{i})],\end{eqnarray} \tag{ 3 }$$

where the sum is over all forest pixels for the quasar q and ${\sigma }_{q}^{2}(\lambda )$ is the estimated variance of the flux, f_i . Our estimate of the ${\sigma }_{q}(\lambda )$ depends on (a_q , b_q ), so we include this dependence in the likelihood function of Equation (3).

The variance σ_q ² receives a contribution from both instrumental noise (readout and photostatistics) and large-scale structure (LSS), and we model it as the sum of three terms:

$$\begin{eqnarray}&&\displaystyle \frac{{\sigma }_{q}^{2}(\lambda )}{{\left(\overline{F}{C}_{q}(\lambda )\right)}^{2}}=\eta (\lambda ){\tilde{\sigma }}_{\mathrm{pip},q}^{2}(\lambda )+{\sigma }_{\mathrm{LSS}}^{2}(\lambda )+\displaystyle \frac{\epsilon (\lambda )}{{\tilde{\sigma }}_{\mathrm{pip},q}^{2}(\lambda )}.\end{eqnarray} \tag{ 4 }$$

The first term is the instrumental noise, ${\tilde{\sigma }}_{\mathrm{pip},q}(\lambda )\,={\sigma }_{\mathrm{pip},q}(\lambda )/\overline{F}{C}_{q}(\lambda )$, where ${\sigma }_{\mathrm{pip},q}^{2}$ is the pipeline estimate of the flux variance. We multiply this by a wavelength-dependent correction, η(λ), that is found to range from 1.04 at λ = 360 nm to 1.20 at λ = 580 nm. The LSS variance, ${\sigma }_{\mathrm{LSS}}^{2}$, gives a minimum value of the variance equal to the intrinsic variance of the flux transmission field. Finally, we add an ad hoc factor proportional to $1/{\tilde{\sigma }}_{\mathrm{pip},q}^{2}$. It describes the observed increase of variance, at high S/N, most probably due to quasar-to-quasar spectral diversity. The three terms in Equation (4) dominate for different ranges of flux and wavelength: for our sample of quasars, the instrumental and LSS terms dominate, respectively, for wavelengths less than or greater than λ ≈ 450 nm. The term is significant only for the brightest quasars, corresponding to ≈10⁻⁴ of the pixels.

The functions $\overline{C}$, η, ${\sigma }_{\mathrm{LSS}}^{2}$, and are computed iteratively. One starts with an initial estimate of $\overline{C}({\lambda }_{\mathrm{RF}})$ as the mean spectrum using as weights an initial estimate of the $1/{\sigma }_{q}^{2}(\lambda )$. The first estimates of the quasar parameters a_q and b_q are then calculated. The resulting δ_q (λ) are calculated and their variance determined in bins of ${\tilde{\sigma }}_{\mathrm{pip},q}$ and λ. The functions η(λ), (λ), and σ_LSS(λ) are then determined by fitting the variance of δ_q (λ) as a function of ${\tilde{\sigma }}_{\mathrm{pip},q}$ and λ. The mean spectrum, $\overline{C}({\lambda }_{\mathrm{RF}})$, is then recalculated with the new weights. This process is repeated until stable values are obtained after about five iterations. This process is executed independently for the Lyα and Lyβ regions. Figure 4 shows histograms of ${\sigma }_{q}^{-2}$ after correction for bias evolution (Equation (7)).

Zoom In Zoom Out Reset image size

This analysis does not determine separately the mean transmission, $\overline{F}(z)$, and the quasar continuum, C_q (λ). In Figure 1, the latter is calculated assuming the $\overline{F}(z)$ measured by Calura et al. (2012) and shown purely for illustration purposes.

The δ_q (λ) for forests with identified DLAs are calculated masking pixels where a DLA reduces the transmission by more than 20%. The masked pixels are not used for the calculation for correlation functions. The absorption in the wings is corrected using a Voigt profile following the procedure of Noterdaeme et al. (2012).

The catalog of δ_q (λ) is publicly available with a data format described in Appendix G. The catalog for data used in this analysis consists of 34.3 and 8.4 million (analysis) pixels in the Lyα and Lyβ regions, respectively. The summary of the Lyα and Lyβ regions is given in Table 1. The redshift distribution of the pixels, assuming Lyα absorption, is presented in the left panel of Figure 3.

As demonstrated in the DR12 Lyα analysis (B17, dMdB17), the quasar fitting of Equation (2) biases the mean and the spectral slope of each individual δ_q toward zero for each line of sight, resulting in a distortion of the correlation function. This distortion can be modeled if we make the biases exact by redefining the δ_q , per spectral region:

$$\begin{eqnarray}&&{\delta }_{q}({\lambda }_{i})\to \displaystyle \sum _{j}{\eta }_{{ij}}^{q}{\delta }_{q}({\lambda }_{j}),\end{eqnarray} \tag{ 5 }$$

where

$$\begin{eqnarray}&&{\eta }_{{ij}}^{q}={\delta }_{{ij}}^{K}-\displaystyle \frac{{w}_{j}}{{\displaystyle \sum }_{k}{w}_{k}}-\displaystyle \frac{{w}_{j}\left({{\rm{\Lambda }}}_{i}-\overline{{{\rm{\Lambda }}}_{q}}\right)\left({{\rm{\Lambda }}}_{j}-\overline{{{\rm{\Lambda }}}_{q}}\right)}{{\displaystyle \sum }_{k}{w}_{k}{\left({{\rm{\Lambda }}}_{k}-\overline{{{\rm{\Lambda }}}_{q}}\right)}^{2}},\end{eqnarray} \tag{ 6 }$$

and $\overline{{{\rm{\Lambda }}}_{q}}$ is the mean of ${\rm{\Lambda }}=\mathrm{log}\lambda$ for spectrum q. In the definition of the projection (Equation (6)), the weights are the ones used in the measurement of the correlation functions:

$$\begin{eqnarray}&&{w}_{i}={\sigma }_{q}^{-2}({\lambda }_{i}){\left(\displaystyle \frac{1+{z}_{i}}{1+2.25}\right)}^{{\gamma }_{\mathrm{Ly}\alpha }-1},\end{eqnarray} \tag{ 7 }$$

where the redshift evolution of the Lyα bias is taken into account (${\gamma }_{\mathrm{Ly}\alpha }=2.9$; McDonald et al. 2006). In Section 3.5, we will use Equations (5) and (6) to establish the relation between the measured correlations and the true correlations.

The mapping applied by Equation (5) changes the evolution of mean δ_q per wavelength bin slightly, allowing it to deviate from zero. We reintroduce this property for the cross-correlation (Section 3.3) by redefining the δ_q per observed wavelength bin:

$$\begin{eqnarray}&&{\delta }_{q}({\lambda }_{i})\to {\delta }_{q}({\lambda }_{i})-\overline{\delta (\lambda )},\end{eqnarray} \tag{ 8 }$$

where the weighted average, $\overline{\delta (\lambda )}$, is computed using the same weights as the ones used when computing the correlation function (Equation (7)). This modification ensures that the cross-correlation approaches zero at large scales, whatever the redshift distribution of the quasars.

In principle, the correlation functions can be computed as a function of redshift and angular separation, ξ(Δz, Δθ), as those are directly observed quantities. However, because of the redshift dependence of the comoving angular diameter distance ${D}_{M}$(z) = (1 + z)D_A (z) and of the Hubble distance ${D}_{H}$(z) = c/H(z), this would widen the BAO peak unless the data were sorted into multiple redshift bins. To avoid degradation of the BAO feature, we convert angular separations to ${r}_{\perp }$ and redshifts separations to ${r}_{\parallel }$ by adopting a "fiducial" cosmology. This acts as an optimal data compression designed to maximize the BAO signal if the fiducial cosmology is correct. As demonstrated in Appendix A, the conversion from observed to comoving coordinates does not bias the measurements of the BAO distance scale.

The fiducial cosmology that we adopt is the ΛCDM cosmology of Planck Collaboration et al. (2016, hereafter Planck 2016). If the fiducial cosmology approximates the true cosmology, then the BAO peak will be at a constant separation, ${r}_{\mathrm{BAO}}\sim 100\,{h}^{-1}\,\mathrm{Mpc}$, for all redshifts. The cosmological parameters for this model are shown in the first part of Table 2 and the derived parameters in the second part; they are computed using the "Code for Anisotropies in the Microwave Background" (CAMB; Lewis et al. 2000). The same assumed cosmology is used to produce and analyze the mock data of Section 5. This cosmology is the same as the one used in DR12 and DR14 studies of Lyα BAOs.

Table 2. Parameters of the Flat-ΛCDM Cosmological Model, from Planck Collaboration et al. (2016), Used for the Production and Analysis of the Mock Data and the Analysis of the Data

Parameter	Planck (2016) Cosmology
	$\left(\mathrm{TT}+\mathrm{lowP}\right)$
${{\rm{\Omega }}}_{m}{h}^{2}$	0.14252
$={{\rm{\Omega }}}_{c}{h}^{2}$	0.1197
$\,+{{\rm{\Omega }}}_{b}{h}^{2}$	0.02222
$\,+{{\rm{\Omega }}}_{\nu }{h}^{2}$	0.0006
h	0.6731
${N}_{\nu }$	3
${\sigma }_{8}$	0.8299
ns	0.9655
${{\rm{\Omega }}}_{m}$	0.31457
${{\rm{\Omega }}}_{r}$	$7.975\,{10}^{-5}$
${r}_{d}\,[\mathrm{Mpc}]$	147.33
${r}_{d}\,[{h}^{-1}\,\mathrm{Mpc}]$	99.17
${D}_{H}(z=2.334)/{r}_{d}$	8.6011
${D}_{M}(z=2.334)/{r}_{d}$	39.2035
$f(z=2.334)$	0.9704

Note. The first part of the table gives the cosmological parameters; the second part gives derived quantities used in this paper. They are computed using CAMB (Lewis et al. 2000).

Download table as:  ASCIITypeset image

The separation between two tracers is determined along and across the line of sight, ${\boldsymbol{r}}=({r}_{\parallel },{r}_{\perp })$. For tracers i and j, of redshift z_i and z_j and offset by an observed angle Δθ, the separation is computed as follows:

$$\begin{eqnarray}&&{r}_{\parallel }=\left[{D}_{c}({z}_{i})-{D}_{c}({z}_{j})\right]\cos \left(\displaystyle \frac{{\rm{\Delta }}\theta }{2}\right)\end{eqnarray} \tag{ 9 }$$

and

$$\begin{eqnarray}&&{r}_{\perp }=\left[{D}_{M}({z}_{i})+{D}_{M}({z}_{j})\right]\sin \left(\displaystyle \frac{{\rm{\Delta }}\theta }{2}\right),\end{eqnarray} \tag{ 10 }$$

where ${D}_{c}(z)={\int }_{0}^{z}{dz}/H(z)$ is the comoving distance. We will also refer to (r, μ) in this study; they are defined as ${r}^{2}={r}_{\parallel }^{2}+{r}_{\perp }^{2}$ and $\mu ={r}_{\parallel }/r$. The quantity μ is the cosine of the angle formed by the median line of sight of both tracers and the vector ${\boldsymbol{r}}$.

For quasars, the redshift is defined to be Z_LYAWG  (see Section 2.2). We test the effects of other estimators of redshifts on the measurement of the BAO in Appendix D. For spectral pixels, the main absorber, and the one used to measure BAO, is Lyα. We thus define the redshift from its rest-frame wavelength. For each spectral pixel i, of observed wavelength ${\lambda }_{i}$, the redshift is then given by ${z}_{i}\,={\lambda }_{i}/{\lambda }_{\mathrm{Ly}\alpha }-1={\lambda }_{i}/121.567-1$. The consequences of the presence of other absorption, like Si ii(126) or the C iv doublet, are discussed in Section 4.

As given in Table 2, the BAO peak is expected at the separation ${r}_{d}\sim 100\,{h}^{-1}\,\mathrm{Mpc};$ furthermore, we know from linear theory that the peak has a width (FWHM) of $\approx 20\,{h}^{-1}\,\mathrm{Mpc}$. For both these reasons, we compute the different correlations up to $\pm 200\,{h}^{-1}\,\mathrm{Mpc}$ along and across the line of sight, i.e., twice the expected BAO scale, with a bin size of 4 h⁻¹ Mpc. At the effective redshift of this study, ${z}_{\mathrm{eff}.}=2.33$, the BAO scale is ${\rm{\Delta }}{\theta }_{\mathrm{BAO}}\sim 1.5\,\deg$ across the line of sight, and ${\rm{\Delta }}{z}_{\mathrm{BAO}}\sim 0.12$ (corresponding to ≈50 analysis pixels) along the line of sight.

For the estimator of the 3D autocorrelation of Lyα absorption in spectral pixels, we use

$$\begin{eqnarray}&&{\xi }_{A}=\displaystyle \frac{{\displaystyle \sum }_{(i,j)\in A}{w}_{i}{w}_{j}\,{\delta }_{i}{\delta }_{j}}{{\displaystyle \sum }_{(i,j)\in A}{w}_{i}{w}_{j}}.\end{eqnarray} \tag{ 11 }$$

This is the standard "covariance" estimator when the mean has been subtracted: $\left\langle {\delta }_{i}\right\rangle =0$, by definition of the projection of Equation (5).

In this equation, i and j refer to two spectral pixels of the flux transmission field, δ_i and δ_j , from Equation (1). The weights w_i and w_j are defined by Equation (7), where the factor ${(1+z)}^{\gamma }$ ${\gamma }_{\mathrm{Ly}\alpha }=2.9$ is chosen to favor high-redshift pixels where the amplitude of the correlation function is greatest (McDonald et al. 2006). We have tested that small changes in the assumed redshift evolution of the weights do not translate into noticeable differences in the BAO results presented in the next sections.

As explained in Section 3.1, the computation of the correlation is done for all possible pairs of pixels (i, j), of separation $({r}_{\parallel },{r}_{\perp })$, within $[0,200]\,{h}^{-1}\,\mathrm{Mpc}$ in both directions. Each bin A is $4\,{h}^{-1}\,\mathrm{Mpc}$ wide in both directions. As a result, each correlation function has N_bin = 50 × 50 = 2500 bins. The sum runs over all possible pairs of pixels from different lines of sight. We exclude pairs of pixels from the same line of sight because of correlated continuum errors that could bias our measurement of the correlation function.

In previous measurements of the Lyα autocorrelation, e.g., B17 and dSA19, pairs of pixels involving the same spectrograph and ${r}_{\parallel }\lt 4\,{h}^{-1}\,\mathrm{Mpc}$ were avoided to minimize spurious correlations due to the sky subtraction procedure. In this analysis, we keep these pairs, allowing us to model their direct effect on the lowest ${r}_{\parallel }$ bins and their indirect effect on other bins through distortion due to continuum fitting. This procedure is detailed in Section 4.5.

The resulting correlations have 6.8 × 10¹¹ pairs of pixels for the Lyα(Lyα)$\,\times \,$Lyα(Lyα) autocorrelation and 3.8 × 10¹¹ for the Lyα(Lyα)$\,\times \,$Lyα(Lyβ) autocorrelation. The lower number of pixels in the Lyβ region, combined with the increased absorption fluctuations and Poisson noise, translates to a variance in the Lyα(Lyα)$\,\times \,$Lyα(Lyβ) correlation function that is approximately three times larger than that in the Lyα(Lyα)$\,\times \,$Lyα(Lyα) function.

The right panel of Figure 3 shows the normalized redshift distribution of the pairs within the BAO region. The Lyα(Lyα) × Lyα(Lyα) distribution is presented in orange, Lyα(Lyα) × Lyα(Lyβ) in green, and the sum of the two in blue. Figure 2 presents the comparison of the weighted number of pairs between this analysis and the DR12 autocorrelation (B17) reproduced using the different improvements of this analysis. The ratio is shown over the BOSS+eBOSS survey footprint, sampled by HEALPix pixels (Górski et al. 2005).

Figure 5 presents the Lyα(Lyα) × Lyα(Lyα) autocorrelation for the data on the left and for the best-fit model on the right (Section 6). The correlation is multiplied by the separation $| r|$ for visual purposes, and the color bar is saturated and symmetric around zero. Even though each individual bin of the correlation is noisy when presented in 2D, the BAO scale is seen in the data at large μ ($r\approx {r}_{\parallel }\approx 100\,{h}^{-1}\,\mathrm{Mpc}$) by the transition from blue (negative values) to white (zeros).

Zoom In Zoom Out Reset image size

The measured autocorrelation, Lyα(Lyα) × Lyα(Lyα), is also presented in the top four panels of Figure 6. These panels show the 2D correlation of Figure 5 reduced to a weighted 1D correlation for four different wedges of $| \mu | =| {r}_{\parallel }/r|$. In the same figure, the best-fit model is shown in red and is discussed in Section 6. The BAO scale peak at $r\sim 100\,{h}^{-1}\,\mathrm{Mpc}$ is visible, especially for $\mu \gt 0.8$. Four similar panels in Appendix F present the autocorrelation Lyα(Lyα) × Lyα(Lyβ) (Figure F1).

Zoom In Zoom Out Reset image size

The autocorrelation, Lyα(Lyα)$\,\times \,$Lyα(Lyα), is also measured in two redshift bins. This allows us to look at the evolution of the different parameters of the best-fit models, e.g., the bias of Lyα, along with allowing systematic tests. In their DR14 study, dSA19 showed that for the Lyα autocorrelation the optimal way to assign pixel pairs to the high- or low-redshift measurement is to cut on the mean of the maximum redshift of the two forests, rather than simply on the redshift of the two pixels. This minimizes the cross-covariance between the low-z and high-z autocorrelations (Section 3.4). As in dSA19, we chose z_cut = 2.5, in order to approximately get the same weighted number of pairs in both redshift splits. The resulting autocorrelation function is presented in the four top panels of Figure F2 in Appendix F.

The covariance matrix of each autocorrelation function has ${N}_{\mathrm{bin}}^{2}=2500\times 2500={\rm{6,250,000}}$ elements. For two bins A and B of the correlation function ξ, the covariance is defined by

$$\begin{eqnarray}&&{C}_{{AB}}=\langle {\xi }_{A}{\xi }_{B}\rangle -\langle {\xi }_{A}\rangle \langle {\xi }_{B}\rangle .\end{eqnarray} \tag{ 12 }$$

In this study, the covariance matrix is estimated by dividing the sky into subsamples defined by HEALPix pixels and computing their weighted covariance. Using nside = 16 over the eBOSS footprint, we get around 880 subsamples, each covering 3.7 × 3.7 = 13.4 deg² on the sky. This solid angle is equivalent to a $250\times 250\,{({h}^{-1}\mathrm{Mpc})}^{2}$ patch at ${z}_{\mathrm{eff}}=2.33$. Lyα pixel pairs are assigned to HEALPix pixels according to the Lyα pixel with the smallest right ascension. The covariance is then given by the following, neglecting the small correlations between subsamples:

$$\begin{eqnarray}&&{C}_{{AB}}=\displaystyle \frac{1}{{W}_{A}{W}_{B}}\displaystyle \sum _{s}{W}_{A}^{s}{W}_{B}^{s}\left[{\xi }_{A}^{s}{\xi }_{B}^{s}-{\xi }_{A}{\xi }_{B}\right].\end{eqnarray} \tag{ 13 }$$

Here s is a subsample with summed weight W_A ^s and measured correlation ${\xi }^{s}$, and ${W}_{A}={\sum }_{s}{W}_{A}^{s}$.

The covariance is dominated by the diagonal elements that are of order ${C}_{{AA}}\approx \langle {\delta }^{2}{\rangle }^{2}/{N}_{A}$, where $\langle {\delta }^{2}\rangle$ is the pixel variance and N_A is the number of pixel pairs in the bin A. Deviations from this simple expression are due to intraforest correlations that make the effective number of independent pairs less than N_A . We find

$$\begin{eqnarray}&&{\mathrm{Var}}_{A}\approx \displaystyle \frac{\langle {\delta }^{2}{\rangle }^{2}}{{{fN}}_{A}^{\mathrm{pair}}},\end{eqnarray} \tag{ 14 }$$

where $\langle {\delta }^{2}\rangle \approx 0.13\,(0.24)$ and f ≈ 0.4 (0.8) for the Lyα and Lyβ regions. For the Lyα(Lyα)$\,\times \,$Lyα(Lyα) correlation this is equivalent to ${\mathrm{Var}}_{A}\approx 1.5\times {10}^{-10}(100\,{h}^{-1}\,\mathrm{Mpc}/{r}_{\perp })$, where the ${r}_{\perp }$ dependence reflects the approximate proportionality between N_A and ${r}_{\perp }$.

Off-diagonal elements of the covariance matrix are due to intraforest correlations that break the independence of pixel pairs in different bins, $A\ne B$. ²⁹ Intraforest correlations reflect both physical correlations of absorption and the effect of continuum fitting. Imperfections in the continuum fitting increase both $\langle {\delta }^{2}\rangle$ and the intraforest correlations and are therefore reflected in the covariance matrix, as long as the imperfections themselves are not correlated between different forests.

The size of the off-diagonal elements is made clearer through the normalized covariance matrix, i.e., the correlation matrix, with elements in [−1, 1]:

$$\begin{eqnarray}&&{\mathrm{Corr}}_{{AB}}=\displaystyle \frac{{C}_{{AB}}}{\sqrt{{\mathrm{Var}}_{A}{\mathrm{Var}}_{B}}},\end{eqnarray} \tag{ 15 }$$

where ${\mathrm{Var}}_{A}={C}_{{AA}}$ is the variance. The largest values are Corr_AB  ≈ 0.4 for ${r}_{\perp }^{A}={r}_{\perp }^{B}$ and $| {r}_{\parallel }^{A}-{r}_{\parallel }^{| }=4{h}^{-1}\,\mathrm{Mpc}$. Elements with ${r}_{\perp }^{A}\ne {r}_{\perp }^{B}$ are very small, <0.03.

The estimates of the off-diagonal elements of the correlation matrix from Equation (15) are noisy, and we smooth them by modeling them as a function of the difference of separation along and across the line of sight: ${\mathrm{Corr}}_{{AB}}={Corr}({r}_{\parallel }^{A},{r}_{\perp }^{A},{r}_{\parallel }^{B},{r}_{\perp }^{B})\,={Corr}({\rm{\Delta }}{r}_{\parallel },{\rm{\Delta }}{r}_{\perp })$, with ${\rm{\Delta }}{r}_{\parallel }=| {r}_{\parallel }^{A}-{r}_{\parallel }^{B}|$ and ${\rm{\Delta }}{r}_{\perp }=| {r}_{\perp }^{A}-{r}_{\perp }^{B}|$. They are presented in the left panel of Figure 7, where they are shown to decrease rapidly as a function of ${\rm{\Delta }}{r}_{\parallel }$ and ${\rm{\Delta }}{r}_{\perp }$.

Zoom In Zoom Out Reset image size

We have also calculated the covariance matrix using other techniques presented in Appendix C. We find no significant change from the estimates of the covariance matrix that is used in the main analysis.

For the Lyα–quasar cross-correlation we use the same estimator as the one from previous studies (Font-Ribera et al. 2012b, 2013, dMdB17, B19). It is defined to be the weighted mean of the flux transmission field at a given separation from a quasar:

$$\begin{eqnarray}&&{\xi }_{A}=\displaystyle \frac{{\displaystyle \sum }_{(i,j)\in A}{w}_{i}{w}_{j}\,{\delta }_{i}}{{\displaystyle \sum }_{(i,j)\in A}{w}_{i}{w}_{j}}.\end{eqnarray} \tag{ 16 }$$

In this equation, i indexes a flux pixel and j a quasar. The weights w_i are as defined in Equation (7) for the Lyα absorption fluctuations. In their study, du Mas des Bourboux et al. (2019) fitted the redshift evolution of the quasar bias and found a best-fit power law of index ${\gamma }_{\mathrm{quasar}}=1.44\pm 0.08$. The quasar weights are then defined to be

$$\begin{eqnarray}&&{w}_{j}={\left(\displaystyle \frac{1+{z}_{j}}{1+2.25}\right)}^{{\gamma }_{\mathrm{quasar}}-1}.\end{eqnarray} \tag{ 17 }$$

In a similar way as for the autocorrelation, the cross-correlation is computed for all pixel–quasar pairs, though omitting pairings of pixels with their own background quasar whose mean correlation vanishes owing to the continuum fitting procedure. The correlation is computed for all pairs within ${r}_{\perp }\in [0,200]\,{h}^{-1}\,\mathrm{Mpc}$. Unlike the autocorrelation, the cross-correlation is not symmetric by permutation of the two tracers. We thus have the opportunity to define positive values of the separation along the line of sight, r_∥, when the tracer quasar is in front of the Lyα pixel tracer, i.e., ${z}_{\mathrm{Ly}\alpha }\gt {z}_{\mathrm{quasar}}$. The line-of-sight separation then ranges over ${r}_{\parallel }\in [-200,200]\,{h}^{-1}\,\mathrm{Mpc}$. With a width of $4\,{h}^{-1}\,\mathrm{Mpc}$, the correlation is computed on ${N}_{\mathrm{bin}}=100\times 50=5000$ bins.

The weighted distribution of the pair redshifts in the BAO region is similar for the cross-correlations to what is presented for the autocorrelation in the right panel of Figure 3. The comparison on the sky of the weighted number of pairs between DR12 (dMdB17) and this study is also similar to what is shown for the autocorrelation in Figure 2.

Figure 8 presents in 2D the measured cross-correlation, Lyα(Lyα) × quasar, and its best-fit model (Section 6). In such a display, the BAO scale would be seen as a half ring of radius $r\sim 100\,{h}^{-1}\,\mathrm{Mpc}$. However, the signal is difficult to see owing to the noise in the data. We show in the bottom four panels of Figure 6 the same cross-correlation but averaged over four different wedges of $| \mu |$. There the BAO scale can be clearly observed as a dip at $r\sim 100\,{h}^{-1}\,\mathrm{Mpc}$. We note that the correlation is reversed, i.e., negative, with respect to the matter correlation function, because the bias of Lyα absorption is negative and the bias of quasars is positive, giving a negative product of biases. The cross-correlation Lyα(Lyβ) × quasar is presented in Figure F1 of Appendix F.

Zoom In Zoom Out Reset image size

As in B19, we alternatively compute the cross-correlation, Lyα(Lyα) × quasar, in two redshift bins. For this measurement, it is important that $\langle \delta \rangle =0$ per bin of observed wavelength. To ensure this, we recompute the δ of Section 2.3, independently for the two redshift bins. We split lines of sight by their background quasar redshift, z_q , and select the splitting redshift, z_cut = 2.57, to have approximately the same weighted number of pairs in each binned cross-correlation. This definition of redshift bins allows us to minimize the cross-covariance between the two samples (Section 3.4). Figure F2 of Appendix F presents, in its four bottom panels, the two redshift bins of the cross-correlation, in four wedges of μ.

The covariance matrix is computed using the same estimators as in Section 3.2 for the autocorrelation. The covariance matrix is dominated by its diagonal elements, i.e., the variance. The latter is approximately inversely proportional to the number of pairs:

$$\begin{eqnarray}&&{\mathrm{Var}}_{A}\approx \displaystyle \frac{\langle {\delta }^{2}\rangle }{0.7{N}_{A}^{\mathrm{pair}}},\end{eqnarray} \tag{ 18 }$$

where the factor 0.7 gives the effective loss of number of pairs, due to correlations between neighboring pixels. For the DR16 data set, the resulting variance for the Lyα(Lyα)$\,\times \,$quasar correlation is ${\mathrm{Var}}_{A}\approx 8.6\times {10}^{-8}(100\,{h}^{-1}\,\mathrm{Mpc}/{r}_{\perp })$. There are a total of ≈1.2 × 10⁹ quasar–pixel pairs in the BAO region, $80\lt r\lt 120\,{h}^{-1}\,\mathrm{Mpc}$.

The off-diagonal elements of the covariance matrix are presented in the right panel of Figure 7. As with those for the autocorrelation, they decrease rapidly as a function of ${\rm{\Delta }}{r}_{\parallel }$ and ${\rm{\Delta }}{r}_{\perp }$.

We modify the estimator of the covariance matrix via subsampling defined by Equation (13) to compute the cross-covariance between the different measured correlation functions:

$$\begin{eqnarray}&&{C}_{{AB}}^{12}=\displaystyle \frac{1}{{W}_{A}{W}_{B}}\displaystyle \sum _{s}{W}_{A}^{s}{W}_{B}^{s}\left[{\xi }_{A}^{1,s}{\xi }_{B}^{2,s}-{\xi }_{A}^{1}{\xi }_{B}^{2}\right].\end{eqnarray} \tag{ 19 }$$

In this equation, ξ¹ and ξ² are two different measured correlation functions, e.g., ξ¹ could be the autocorrelation Lyα(Lyα) × Lyα(Lyα) and ξ² the cross-correlation Lyα(Lyα) × quasar.

The different correlation functions are found to be marginally correlated. For example, the correlations between the Lyα(Lyα)$\,\times \,$Lyα(Lyα) and Lyα(Lyα)$\,\times \,$quasar functions are less than 1%, as shown in Figure 11.

The continuum fitting (Equation (2)) and the projection (Equation (5)) mix the ${\delta }_{q}$ within a given forest and thereby modify significantly the correlation function. The transformation (Equation (8)) also has a small effect on the correlations. Under the assumption that the δ transformations are linear, the true and distorted correlation functions are related by a "distortion matrix," D_AB :

$$\begin{eqnarray}&&{\hat{\xi }}_{\mathrm{distorted}}(A)=\displaystyle \sum _{B}{D}_{{AB}}{\xi }_{\mathrm{true}}(B).\end{eqnarray} \tag{ 20 }$$

To the extent that the projection given by Equations (5) and (6) captures the full distortion, the distortion matrices for the auto- and cross-correlations can be read off from the coefficients in Equation (6):

$$\begin{eqnarray}&&{D}_{{AB}}^{\mathrm{auto}}={W}_{A}^{-1}\displaystyle \sum _{{ij}\in A}{w}_{i}{w}_{j}\,\left(\displaystyle \sum _{i^{\prime} j^{\prime} \in B}{\eta }_{{ii}^{\prime} }{\eta }_{{jj}^{\prime} }\right),\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{D}_{{{AA}}^{{\prime} }}^{\mathrm{cross}}={W}_{A}^{-1}\displaystyle \sum _{(i,j)\in A}{w}_{i}{w}_{j}\left(\displaystyle \sum _{({i}^{{\prime} },j)\in {A}^{{\prime} }}{\eta }_{{{ii}}^{{\prime} }}\right).\end{eqnarray} \tag{ 22 }$$

In the fits of the data described in Section 6, the physical model of the correlation function (Section 4) is multiplied by the distortion matrix before comparison with the data. This procedure is validated with the analysis of the mock data sets (Section 5), where it is shown that the BAO parameters are accurately recovered in the presence of distortion.

The components of the correlation function that reflect the underlying matter correlations are given by the Fourier transform of the tracer biased power spectrum:

$$\begin{eqnarray}&&\hat{P}({\boldsymbol{k}})={b}_{i}{b}_{j}\left(1+{\beta }_{i}{\mu }_{k}^{2}\right)\left(1+{\beta }_{j}{\mu }_{k}^{2}\right){P}_{\mathrm{QL}}({\boldsymbol{k}}){F}_{\mathrm{NL}}({\boldsymbol{k}})G({\boldsymbol{k}}),\end{eqnarray} \tag{ 27 }$$

where the vector ${\boldsymbol{k}}=({k}_{\parallel },{k}_{\perp })=(k,{\mu }_{k})$ of modulus k has components along and across the line of sight, $({k}_{\parallel },{k}_{\perp })$, with ${\mu }_{k}={k}_{\parallel }/k$. The bias and RSD parameters, (b, β), are for the tracer i or j. P_QL is the quasi-linear power spectrum defined below, F_NL corrects for nonlinear effects at large ${\boldsymbol{k}}$, and $G({\boldsymbol{k}})$ is a damping term that accounts for averaging of the correlation function in individual $({r}_{\parallel },{r}_{\perp })$ bins. For the autocorrelation, the main tracer is the Lyα absorption, i = j, and for the cross-correlation two main tracers play a role: Lyα absorption and quasars, $i\ne j$.

For the cross-correlation power spectrum we do not include relativistic effects that lead to an asymmetry in ${r}_{\parallel }$ (quasar more distant or less distant than the forest). These effects have been calculated to be sufficiently small to be negligible for this study (Lepori et al. 2020). They were included in the study of B19 and were shown to be partially degenerate with the parameter (${\rm{\Delta }}{r}_{\parallel ,{QSO}}$) describing systematic errors in quasar redshift.

The quasi-linear power spectrum provides for the aforementioned decoupling of the peak component:

$$\begin{eqnarray}&&{P}_{\mathrm{QL}}({\boldsymbol{k}},z)={P}_{\mathrm{sm}}(k,z)+\exp \left[-\displaystyle \frac{{k}_{\parallel }^{2}{{\rm{\Sigma }}}_{\parallel }^{2}+{k}_{\perp }^{2}{{\rm{\Sigma }}}_{\perp }^{2}}{2}\right]{P}_{\mathrm{peak}}(k,z).\end{eqnarray} \tag{ 28 }$$

The smooth component, P_sm, is derived from the linear power spectrum, ${P}_{{\rm{L}}}(k,z)$, via the sideband technique (Kirkby et al. 2013), which is implemented in picca, with the help of nbodykit ³¹ (Hand et al. 2018). The CAMB ${P}_{{\rm{L}}}({z}_{\mathrm{eff}})$ is Fourier transformed into the correlation function where two sidebands are defined on both side of the BAO. A smooth function is then fitted to connect the two sidebands, allowing us to produce a correlation function, without the BAO peak. The latter correlation is Fourier transformed back to a smooth power spectrum, P_sm, that does not have the BAO wiggle features. The peak power spectrum is thus defined as ${P}_{\mathrm{peak}}\,={P}_{{\rm{L}}}-{P}_{\mathrm{sm}}$.

The correction for nonlinear broadening of the BAO peak (Eisenstein et al. 2007) is parameterized by $({{\rm{\Sigma }}}_{\parallel },{{\rm{\Sigma }}}_{\perp })$, with ${{\rm{\Sigma }}}_{\perp }=3.26\,{h}^{-1}\,\mathrm{Mpc}$ and

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{\Sigma }}}_{\parallel }}{{{\rm{\Sigma }}}_{\perp }}=1+f,\end{eqnarray} \tag{ 29 }$$

where $f=d(\mathrm{ln}g)/d(\mathrm{ln}a)\approx {{\rm{\Omega }}}_{{\rm{m}}}^{0.55}(z)$ is the linear growth rate of structure, resulting in ${{\rm{\Sigma }}}_{\parallel }=6.42\,{h}^{-1}\,\mathrm{Mpc}$.

The function ${F}_{\mathrm{NL}}({\boldsymbol{k}})$ in Equation (27) accounts for nonlinear effects at small scales. For the autocorrelation, the most important effects are thermal broadening, peculiar velocities, and nonlinear structure growth. We use Equation (3.6) of Arinyo-i-Prats et al. (2015) with parameter values from their Table 7 interpolated to our effective redshift z = 2.334. For the cross-correlation, the most important effect is that of quasar nonlinear velocities, and, following Percival & White (2009), we adopt

$$\begin{eqnarray}&&{F}_{\mathrm{NL}}^{\mathrm{cross}}({k}_{\parallel })=\displaystyle \frac{1}{1+{\left({k}_{\parallel }{\sigma }_{v}\right)}^{2}},\end{eqnarray} \tag{ 30 }$$

where ${\sigma }_{v}$ is a free parameter. This function also takes into account statistical quasar redshift errors.

The final term in Equation (27) accounts for the effect of the binning of the correlation function on the separation grid. We assume the distribution to be homogeneous on each bin ³² and compute the function $G({\boldsymbol{k}})$ as the product of the Fourier transforms of the rectangle functions that model a uniform square bin:

$$\begin{eqnarray}&&G({\boldsymbol{k}})=\mathrm{sinc}\left(\displaystyle \frac{{k}_{\parallel }{R}_{\parallel }}{2}\right)\mathrm{sinc}\left(\displaystyle \frac{{k}_{\perp }{R}_{\perp }}{2}\right),\end{eqnarray} \tag{ 31 }$$

where R_∥ and R_⊥ are the radial and transverse widths of the bins, respectively.

The dominant Lyα absorption can be viewed as the sum of contributions from the diffuse IGM and from HCD ³³ systems. HCD absorbers are expected to trace the underlying density field, and their effect on the flux transmission field depends on whether they are identified and given the special treatment described in Section 2. If they are correctly identified with the total absorption region masked and the wings correctly modeled, they can be expected to have no significant effect on the field. Conversely, if they are not identified, the measured correlation function will be modified because their absorption is spread along the radial direction. This broadening effect introduces a ${k}_{\parallel }$ dependence of the effective bias (Font-Ribera & Miralda-Escudé 2012):

$$\begin{eqnarray}&&\left\{\begin{array}{l}{b}_{\mathrm{Ly}\alpha }^{{\prime} }={b}_{\mathrm{Ly}\alpha }+{b}_{\mathrm{HCD}}{F}_{\mathrm{HCD}}({k}_{\parallel })\\ {b}_{\mathrm{Ly}\alpha }^{{\prime} }{\beta }_{\mathrm{Ly}\alpha }^{{\prime} }={b}_{\mathrm{Ly}\alpha }{\beta }_{\mathrm{Ly}\alpha }+{b}_{\mathrm{HCD}}{\beta }_{\mathrm{HCD}}{F}_{\mathrm{HCD}}({k}_{\parallel }),\end{array}\right.\end{eqnarray} \tag{ 32 }$$

where $({b}_{\mathrm{Ly}\alpha },{\beta }_{\mathrm{Ly}\alpha })$ and $({b}_{\mathrm{HCD}},{\beta }_{\mathrm{HCD}})$ are the bias parameters associated with the IGM and HCD systems, respectively, and ${F}_{\mathrm{HCD}}$ is a function that depends on the number and column density distribution of HCDs. Rogers et al. (2018) numerically calculated ${F}_{\mathrm{HCD}}$ using hydrodynamical simulations, and we find that a simple exponential form, ${F}_{\mathrm{HCD}}=\exp (-{L}_{\mathrm{HCD}}{k}_{\parallel })$, well approximates their results. Here ${L}_{\mathrm{HCD}}$ is a typical length scale for unmasked HCDs. For eBOSS spectral resolution, DLA identification is possible for DLA widths (wavelength interval for absorption greater than 20%) greater than ∼2.0 nm, corresponding to $\sim 14\,{h}^{-1}\,\mathrm{Mpc}$ in our sample. Degeneracies between ${L}_{\mathrm{HCD}}$ and other anisotropic parameters make the fitting somewhat unstable. We therefore impose ${L}_{\mathrm{HCD}}=10\,{h}^{-1}\,\mathrm{Mpc}$ while fitting for the bias parameters ${b}_{\mathrm{HCD}}$ and ${\beta }_{\mathrm{HCD}}$. We have verified that setting ${L}_{\mathrm{HCD}}$ in the range $7\lt {L}_{\mathrm{HCD}}\lt 13\,{h}^{-1}\,\mathrm{Mpc}$ does not significantly change the inferred BAO peak position. Cuceu et al. (2020) also showed that fixing ${L}_{\mathrm{HCD}}$ has a minimal impact on BAO results when doing a Bayesian analysis of Lyα BAOs.

In our measured correlations functions, different $({r}_{\parallel },{r}_{\perp })$ bins have mean redshifts that vary over the range $2.32\lt \bar{z}\lt 2.39$. Therefore, in order to fit a unique function to the data, we need to assume a redshift dependence of the bias parameters. Following McDonald et al. (2006), we assume that the product of ${b}_{\mathrm{Ly}\alpha }$ and the growth factor of structures varies with redshift as ${(1+z)}^{{\gamma }_{\alpha }-1}$, with ${\gamma }_{\alpha }=2.9$, while we make use of the approximation that ${\beta }_{\mathrm{Ly}\alpha }$ does not depend on redshift. Because the fit of the cross-correlation is only sensitive to the product of the quasar and Lyα biases, we choose to adopt a value for the quasar bias. Following the analysis of du Mas des Bourboux et al. (2019) based on a compilation of quasar bias measurements (Croom et al. 2005; Shen et al. 2013; Laurent et al. 2016, 2017), we set ${b}_{{\rm{q}}}\equiv {b}_{{\rm{q}}}({z}_{\mathrm{eff}})=3.77$ with a redshift dependence given by

$$\begin{eqnarray}&&{b}_{{\rm{q}}}(z)=3.60{\left(\displaystyle \frac{1+z}{3.334}\right)}^{1.44}.\end{eqnarray} \tag{ 33 }$$

The quasar RSD, assumed to be redshift independent in the matter-dominated epochs explored here, is

$$\begin{eqnarray}&&{\beta }_{{\rm{q}}}=\displaystyle \frac{f}{{b}_{{\rm{q}}}}.\end{eqnarray} \tag{ 34 }$$

Setting f = 0.9704 for our fiducial cosmology yields ${\beta }_{{\rm{q}}}\,=0.269$, which we take to be redshift independent.

We model the power spectrum for correlations involving metals with the same form we use for Lyα–Lyα and Lyα–quasar correlations (Equation (27)) except that HCD effects are neglected. Each metal species, n, has bias parameters (b_n , β_n ). The Fourier transform of the power spectrum ${P}_{{mn}}({\boldsymbol{k}},z)$ for the absorber pair (m, n) is then the model correlation function of the pair (m, n): ${\xi }_{\mathrm{mod}}^{m-n}({\tilde{r}}_{\parallel },{\tilde{r}}_{\perp })$, where $({\tilde{r}}_{\parallel },{\tilde{r}}_{\perp })$ are the separations calculated using the correct rest-frame wavelengths, $({\lambda }_{m},{\lambda }_{n})$. For the cross-correlation, the same form is used except that the bias parameters for the metal n are replaced with the quasar bias parameters and the separations are calculated using the quasar redshift.

Since we assign redshifts to pixels assuming Lyα absorption, the rest-frame wavelength we ascribe to a metal transition is not equal to the true rest-frame wavelength. For $m\ne n$ and n = Lyα, this misidentification results in a shift of the model correlation function with vanishing separation of absorbers corresponding to reconstructed separations ${r}_{\perp }=0$ and ${r}_{\parallel }\,\approx (1+z){D}_{H}(z)({\lambda }_{m}-{\lambda }_{\mathrm{Ly}\alpha })/{\lambda }_{\mathrm{Ly}\alpha }$. The values of ${r}_{\parallel }$ for the major metal transition are given in Table 3. For m = n, the reconstructed separations are scaled from the true separations by a factor ${D}_{H}(z)/{D}_{H}({z}_{\alpha })$ for ${r}_{\parallel }$ and by a factor ${D}_{M}(z)/{D}_{M}({z}_{\alpha })$ for ${r}_{\perp }$, where z and z_α are the true and reconstructed redshifts, respectively.

Table 3. Metal Lines Included in the Synthetic Spectra

Metal Line	${\lambda }_{m}$ (nm)	Relative	${r}_{\parallel }$
		Strength ($\times {10}^{3}$)	$({h}^{-1}\,\mathrm{Mpc})$
Si iii	120.7	1.892	−21
Si iia	119.0	0.642	−64
Si iib	119.3	0.908	−56
Si iic	126.0	0.354	+111

Note. For each metal line their transition wavelength and their relative strength with respect to the Lyα optical depth are given. The last column gives the position of maximum apparent correlation ${r}_{\perp }=0$ and ${r}_{\parallel }\approx (1+z){D}_{H}(z)({\lambda }_{m}-{\lambda }_{\mathrm{Ly}\alpha })/\lambda$, corresponding to metal and Lyα absorption at the same physical position but reconstructed assuming only Lyα absorption at ${z}_{\mathrm{eff}}=2.334$.

Download table as:  ASCIITypeset image

For each pair (m, n) of contaminants, we compute the shifted–model correlation function with respect to the unshifted–model correlation function, ${\xi }_{\mathrm{mod}}^{m-n}$, by introducing a metal matrix M_AB (Blomqvist et al. 2018), such that

$$\begin{eqnarray}&&{\xi }_{\mathrm{mod}}^{m-n}(A)\to \displaystyle \sum _{B}{M}_{{AB}}{\xi }_{\mathrm{mod}}^{m-n}({\tilde{r}}_{\parallel }(B),{\tilde{r}}_{\perp }(B)),\end{eqnarray} \tag{ 35 }$$

where

$$\begin{eqnarray}&&{M}_{{AB}}=\displaystyle \frac{1}{{W}_{A}}\displaystyle \sum _{(m,n)\in A,(m,n)\in B}{w}_{m}{w}_{n},\end{eqnarray} \tag{ 36 }$$

where ${W}_{A}={\sum }_{(m,n)\in A}{w}_{m}{w}_{n}$, $(m,n)\in A$ refers to pixel separations computed assuming z_α , and $(m,n)\in B$ to pixel separations computed using the redshifts of the m and n absorbers, z_m and z_n . We take into account the redshift dependence of the weights in the computation of w_m and w_n .

In the fits of the data, we include terms corresponding to metal–Lyα correlations for the four silicon transitions listed in Table 3. Since these correlations are only visible in $({r}_{\parallel },{r}_{\perp })$ bins corresponding to small physical separations, b_m and β_m cannot be determined separately. We therefore fix β_m  = 0.50 as done in B17 based on measurements of the cross-correlation between DLAs and the Lyα forest (Font-Ribera et al. 2012b). Correlations between C iv and Lyα give a contribution outside the ${r}_{\parallel }$ range studied here, but we include effects of the Civ–Civ autocorrelation and fix ${\beta }_{{\rm{C}}{\rm\small{IV}}({eff})}=0.27$ (Blomqvist et al. 2018).

The term in Equation (26) representing the transverse proximity effect takes the form (Font-Ribera et al. 2013)

$$\begin{eqnarray}&&{\xi }^{\mathrm{TP}}={\xi }_{0}^{\mathrm{TP}}{\left(\displaystyle \frac{1{h}^{-1}\mathrm{Mpc}}{r}\right)}^{2}\exp (-r/{\lambda }_{\mathrm{UV}}).\end{eqnarray} \tag{ 37 }$$

This form supposes isotropic emission from the quasars. We fix ${\lambda }_{\mathrm{UV}}=300\,{h}^{-1}\,\mathrm{Mpc}$ (Rudie et al. 2013) and fit for the amplitude ${\xi }_{0}^{\mathrm{TP}}$.

The correlation ${\xi }^{\mathrm{sky}}$, of Equation (25), models the correlations induced by the sky subtraction procedure. The sky subtraction for each spectrum is done independently for each spectrograph, i.e., per half-plate, for the 500 fibers (450 science fibers and ∼40 sky fibers). The Poisson fluctuations in the sky spectra that are subtracted induce correlations in spectra obtained with the same spectrograph at the same observed wavelength, leading to an excess correlation in ${r}_{\parallel }=0$ bins. This was observed in the DR12 autocorrelation measurement of B17, where they decided to reject such same-spectrograph pairs in their measurement. This effectively removes the excess correlation at ${r}_{\parallel }=0$. However, because of continuum fitting, the excess correlation at ${r}_{\parallel }=0$ generates a smooth distortion of the correlation function for all ${r}_{\parallel }$, and this was not removed by the procedure of B17. Here we do not remove the same-spectrograph pairs, which allows us to fit for its amplitude and thereby take into account the smooth component induced by continuum fitting. Apart from improving the fit to the data, this procedure is necessary because the DR16 analysis combines all measurements of each quasar, so spectra no longer correspond to a unique spectrograph.

The sky-fiber-induced correlations are easily seen at rest-frame wavelength longer than the Lyα quasar emission line, where correlations due to absorption are small. An example is shown in Figure 9, showing correlations in the "Mg ii(11)" spectral region, ${\lambda }_{\mathrm{RF}}\in [260,276]\,\mathrm{nm}$ (du Mas des Bourboux et al. 2019), where the pixels are immediately blueward of the Mg ii quasar emission line, ${\lambda }_{\mathrm{RF}}=279.6\,\mathrm{nm}$. The correlation is consistent with zero for pairs of pixels taken with different spectrographs but is significant for pairs from the same spectrograph. This spurious correlation decreases rapidly with increasing angular separation. This decrease is due to the eBOSS pipeline that adds broadband functions to the sky calibration, modeling its variation across the size of the plate and thus decreasing the correlations with increasing angular separation.

Zoom In Zoom Out Reset image size

Since, as the angular separation increases, the number of same-spectrograph pairs decreases and the number of different-spectrograph pairs increases, the weighted sum of the two, shown as the red curve in the figure, can be modeled as a Gaussian function of ${r}_{\perp }$:

$$\begin{eqnarray}{\xi }^{\mathrm{sky}}({r}_{\parallel },{r}_{\perp })=\left\{\begin{array}{ll}\tfrac{{A}_{\mathrm{sky}}}{{\sigma }_{\mathrm{sky}}\sqrt{2\pi }}\exp \left(-\tfrac{1}{2}{\left(\tfrac{{r}_{\perp }}{{\sigma }_{\mathrm{sky}}}\right)}^{2}\right), & \mathrm{if}\,{r}_{\parallel }=0\\ 0, & \mathrm{if}\,{r}_{\parallel }\ne 0\end{array}\right..\end{eqnarray} \tag{ 38 }$$

The two free parameters ${(A,\sigma )}_{\mathrm{sky}}$ give the scale and the width of the correlation.

We fit the measured autocorrelation, Mg ii(Mg ii(11)) × Mg ii(Mg ii(11)), via the distorted model of Equation (40) and the model of the sky correlation of Equation (38). The fit is performed in the same condition as for the Lyα autocorrelation, hence $r\in [10,180]\,{h}^{-1}\,\mathrm{Mpc}$ (Section 6), resulting in 1590 bins. We get a best fit with ${\chi }^{2}/(\mathrm{DOF})=1551.94/(1590-2)$, corresponding to a probability p = 0.76 when using only the best observation. We find ${\chi }^{2}/(\mathrm{DOF})=1531.49/(1590-2)$ and p = 0.86 when using all observations. In both cases, we find ${(A,\sigma )}_{\mathrm{sky}}\sim (1\times {10}^{-3},20\,{h}^{-1}\,\mathrm{Mpc})$. A null test on the autocorrelation, using all observations, proves the very high significance of these two parameters: ${\chi }^{2}=2824.85$ and Δχ² = 1293.36, for a difference of two parameters. Interestingly, a similar null test, removing the ${r}_{\parallel }=0$ bins, still yields a significant measurement of the two parameters, Δχ² = 223.50. This test shows that the effect of the sky subtraction residuals is mainly in the ${r}_{\parallel }=0$ bins but that a nonnegligible fraction is brought to other bins by the effects of continuum fitting.

For the Lyα(Lyα)$\,\times \,$Lyα(Lyα) correlation function, the spurious correlation from the sky calibration in the ${r}_{\parallel }=0$ bins reaches ≈20% of the physical correlation at ${r}_{\perp }\approx 50\,{h}^{-1}\,\mathrm{Mpc}$. For ${r}_{\parallel }\ne 0$, the distortion-induced correlation is smooth, so it does not affect our BAO measurement. However, it has a nonnegligible effect on the measurement of the Lyα bias parameters (Table D1), because of its μ dependence.

An important test of systematic effects in the position of the BAO peak is performed by adding polynomial "broadband" terms to the correlation function. We follow the procedure and choice of broadband forms used by B17 and adopt the form

$$\begin{eqnarray}&&B(r,\mu )=\displaystyle \sum _{j=0}^{{j}_{\max }}\displaystyle \sum _{i={i}_{\min }}^{{i}_{\max }}{a}_{{ij}}\displaystyle \frac{{L}_{j}(\mu )}{{r}^{i}}\ \,(j\,\mathrm{even}),\end{eqnarray} \tag{ 39 }$$

where the L_j are Legendre polynomials. Following B17, we fit with $({i}_{\min },{i}_{\max })=(0,2)$ corresponding to a parabola in ${r}^{2}{\xi }_{\mathrm{smooth}}$ underneath the BAO peak. We set ${j}_{\max }\,=\,6$, giving four values of j corresponding to approximately independent broadbands in each of the four angular ranges.

The expected measured correlation function, $\hat{\xi }$, is related to the true ${\xi }^{t}$, by the distortion matrix:

$$\begin{eqnarray}&&{\hat{\xi }}_{A}=\displaystyle \sum _{{A}^{{\prime} }}{D}_{{{AA}}^{{\prime} }}{\xi }_{{A}^{{\prime} }}^{t}.\end{eqnarray} \tag{ 40 }$$

The distortion matrix, ${D}_{{{AA}}^{{\prime} }}$, is computed in Section 3.5 for the auto- and cross-correlations. The matrix accounts for the correlations introduced by the continuum fitting procedure. In this equation, A' is a bin of the model and A is a bin of the measurement. ³⁴

Font-Ribera et al. (2012a) presented a method to efficiently simulate Lyα forest spectra for a given survey configuration. This method was used in multiple analyses of the Lyα autocorrelation in BOSS (Slosar et al. 2011, 2013; Busca et al. 2013; Delubac et al. 2015; Bautista et al. 2017), and it enabled multiple tests of potential systematics and the validation of the BAO analysis pipeline.

However, because the quasar positions behind the simulated forests did not correspond to density peaks, these simulated data sets did not have the correct cross-correlation between the quasars and the Lyα forest. This made them unfit for BAO studies using this cross-correlation. For this reason, the DR12 analysis of the cross-correlation in dMdB17 presented a new set of simulations using a different algorithm described in Le Goff et al. (2011), and that included a realistic cross-correlation. This allowed the first validation of BAOs in cross-correlation using mock data, and it showed that the auto- and cross-correlation measurements of the BAO scale have a very small correlation and can therefore be combined. These simulations, however, did not have the correct survey geometry and assumed instead that the lines of sight were parallel.

We have developed two different sets of simulations that have been generated in collaboration with the Lyα working group of DESI. The results presented in this section use simulations described in Farr et al. (2020), developed primarily at University College London and referred to here as the London mocks. A second set (T. Etourneau et al. 2020, in preparation) was developed primarily at CEA Saclay and is referred to here as the Saclay mocks. These mocks were completed after the London mocks and benefit from refinements inspired by experience with the London mocks. They will be used for future studies of the non-BAO component of the correlation function, including the effect of HCDs.

5.1.1. London Mocks

5.1.2. Saclay Mocks

Once we have the simulated transmitted flux fraction along each line of sight, we use these to simulate synthetic quasar spectra with the relevant astrophysical and instrumental artifacts. The methodology here is similar to the one used in the BOSS Lyα analyses, described in Bautista et al. (2015), and it can be summarized as follows:

• 1.  
  Add contaminant absorption, including higher-order hydrogen lines, metal absorbers, and HCD systems.
• 2.  
  Multiply each transmitted flux fraction by a simulated quasar continuum.
• 3.  
  Convolve the simulated spectrum with the resolution of the spectrograph, pixelize it, and add Gaussian noise simulating eBOSS observations.

Figure 10 shows a simulated spectrum with different levels of complexity.

Zoom In Zoom Out Reset image size

5.2.1. Adding Contaminants

5.2.2. Simulating Quasar Continua

5.2.3. Simulating Instrumental Effects

As described in Farr et al. (2020), the synthetic Lyα spectra were tuned in order to reproduce large-scale measurements from the BOSS DR12 results, as well as reproducing the line-of-sight power spectrum over the relevant scales. In Appendix F we present figures comparing the measured correlations in the mocks and in the data, including for the first time comparisons of measurements of Lyα absorption in the Lyβ region.

An important role of the mocks is to verify the subsampling procedure for calculating the covariance matrix of the correlation functions. This can be done by comparing the subsampling calculation with the mock-to-mock variations of the correlation functions. We see no difference between the two at more than the 1% level. For example, in the left panel of Figure 11 we illustrate off-diagonal elements of the correlation (normalized covariance) matrix, which has been smoothed as described in Section 3.2 for the Lyα(Lyα)$\,\times \,$Lyα(Lyα) and Lyα(Lyα)$\,\times \,$quasar correlation functions. The correlations measured by subsampling and by mock-to-mock variations are in excellent agreement, verifying our procedure.

Zoom In Zoom Out Reset image size

Figure 11 shows that the covariances for mocks and data have a very similar structure, the differences being due in part to the smaller metal component in the mocks than in the data, leading to less correlation between nearby wavelengths. The mock variances are a factor of two smaller for the mocks for the autocorrelation and a factor of 1.5 for the cross-correlations. This is due to fewer low-flux spectra in the mocks compared to the data.

In Table 4 we present the results of our BAO analysis on multiple synthetic realizations of our data set. In order to study the effect of the different systematics included in the mock spectra, we run our analysis on different versions of each synthetic data set:

• 1.  
  Lyα: analysis run directly on the simulated transmitted flux fractions, including only Lyα absorption. This analysis does not need to do a continuum fitting and is therefore unaffected by the continuum distortions discussed in Section 3.5.
• 2.  
  + continuum + noise: analysis run on simulated spectra, with quasar continua and instrumental noise.
• 3.  
  + metals: analysis run on the same simulated spectra as above, but including also absorption from metal lines.
• 4.  
  + HCDs + ${\sigma }_{v}$: analysis run on the same simulated spectra as above, but including also absorption from HCDs and including random redshift errors for the quasars. When analyzing the real data, we mask HCDs identified with an algorithm that is able to find the largest systems, with $\mathrm{log}{N}_{{\rm{H}}{\rm\small{I}}}\gt 20.3$. Similarly, when computing the delta field in the mock spectra, HCDs with $\mathrm{log}{N}_{{\rm{H}}{\rm\small{I}}}\gt 20.3$ are corrected using the input N_{H i} and pixels with absorption larger than 20% are masked. ³⁹

Table 4. Fit Results for Mock Data Sets

Mock Set	$\overline{{\alpha }_{\parallel }}\,\,\overline{\sigma }$	$\overline{{\alpha }_{\perp }}\,\,\overline{\sigma }$	$\overline{{b}_{\eta ,\mathrm{Ly}\alpha }}\,\,\overline{\sigma }$	$\overline{{\beta }_{\mathrm{Ly}\alpha }}\,\,\overline{\sigma }$	$\overline{{\chi }_{\min }^{2}}/\mathrm{DOF},\overline{\mathrm{proba}}$
Lyα(Lyα) × Lyα(Lyα)
raw	1.012 0.021	0.985 0.028	−0.2 0.001	1.568 0.021	1629.67/(1590-4) 0.27
+cont+noise	1.003 0.027	0.995 0.04	−0.201 0.002	1.486 0.028	1603.57/(1590-4) 0.41
+metals	1.012 0.029	0.987 0.05	−0.202 0.002	1.485 0.03	1600.32/(1590-8) 0.43
+HCD+${\sigma }_{v}$	1.004 0.029	1.001 0.041	−0.205 0.003	1.48 0.051	1597.99/(1590-10) 0.42
Lyα(Lyα) × quasar
raw	1.008 0.025	0.999 0.024	−0.189 0.003	1.568 0.041	3227.31/(3180-6) 0.29
+cont+noise	1.008 0.029	0.992 0.033	−0.192 0.004	1.491 0.061	3203.5/(3180-6) 0.39
+metals	1.006 0.029	0.994 0.033	−0.193 0.004	1.51 0.063	3194.72/(3180-10) 0.4
+HCD+${\sigma }_{v}$	1.003 0.033	0.998 0.033	−0.199 0.007	1.48 0.081	3212.8/(3180-10) 0.35
Lyα(Lyα) × Lyα(Lyβ)
raw	1.005 0.025	0.996 0.034	−0.2 0.002	1.588 0.026	1609.71/(1590-4) 0.41
+cont+noise	1.014 0.049	0.983 0.069	−0.202 0.003	1.509 0.05	1592.39/(1590-4) 0.47
+metals	1.02 0.049	0.994 0.065	−0.203 0.004	1.528 0.054	1589.8/(1590-8) 0.46
+HCD+${\sigma }_{v}$	1.009 0.054	1.019 0.087	−0.206 0.004	1.502 0.085	1609.88/(1590-10) 0.35
Lyα(Lyβ) × quasar
raw	1.028 0.042	1.009 0.044	−0.189 0.005	1.595 0.073	3189.37/(3180-6) 0.45
+cont+noise	1.008 0.07	1.015 0.082	−0.193 0.01	1.527 0.146	3183.49/(3180-6) 0.47
+metals	0.994 0.071	1.002 0.093	−0.19 0.01	1.495 0.149	3216.12/(3180-10) 0.36
+HCD+${\sigma }_{v}$	1.011 0.08	1.013 0.099	−0.192 0.015	1.447 0.186	3195.57/(3180-10) 0.4
all combined
raw	1.009 0.012	0.995 0.014	−0.203 0.001	1.628 0.015	9948.08/(9540-7) 0.02
+cont+noise	1.005 0.017	0.992 0.022	−0.206 0.002	1.553 0.023	9788.16/(9540-7) 0.1
+metals	1.01 0.018	0.989 0.023	−0.206 0.002	1.558 0.025	9822.24/(9540-11) 0.08
+HCD+${\sigma }_{v}$	1.005 0.019	0.998 0.023	−0.205 0.002	1.464 0.036	9625.56/(9540-13) 0.31

Note. Shown are the mean values and mean standard deviations (${\rm{\Delta }}{\chi }^{2}=1$) for ${\alpha }_{\parallel }$, ${\alpha }_{\perp }$, ${b}_{\eta ,\mathrm{Ly}\alpha }$, and ${\beta }_{\mathrm{Ly}\alpha }$. Results are shown for "raw" mocks (no quasar continuum or noise added to the transmission field) and for three progressively more realistic mocks (adding continua, noise, metals, HCDs, and quasar velocity dispersion). Means are based on 100 mocks for (+cont+noise) and (HCD+${\sigma }_{v}$) and 10 otherwise.

Download table as:  ASCIITypeset image

We ran all four analyses (Lyα autocorrelations and quasar cross-correlations using the Lyα and the Lyβ regions) on 10 realizations for each of the synthetic data sets described above. We analyze 90 extra realizations for two of the settings above, (+continuum+noise) and (+HCDs+${\sigma }_{v}$), allowing us to test possible systematic biases on the BAO measurements at a level 10 times smaller than our statistical uncertainty.

The main conclusion of this analysis is that we are able to recover the right BAO scale even in the presence of contaminants, and that our model is able to describe the correlations measured on the synthetic data sets.

5.4.1. Correlations in a BAO Measurement

Anisotropic BAO studies provide a measurement of the line-of-sight BAO scale, ${\alpha }_{\parallel }$, and of the transverse scale, ${\alpha }_{\perp }$. In BAO analysis using galaxy clustering these parameters are anticorrelated at the 40% level, i.e., they have a correlation coefficient around $\rho \left({\alpha }_{\parallel },{\alpha }_{\perp }\right)=-0.4$.

In Table 5 we use the BAO analyses on the 100 synthetic data sets to measure the correlation coefficient $\rho \left({\alpha }_{\parallel },{\alpha }_{\perp }\right)$ in the four BAO measurements discussed in this publication. One can see that the results are also anticorrelated, with values ranging from −0.347 to −0.615. The correlations involving $({\alpha }_{\parallel },{\alpha }_{\perp })$ from distinct pairs of correlation functions are consistent with zero.

Table 5. Correlations between the BAO Parameters Measured in the Different Correlation Functions as Measured in Mocks

Correlation	$\rho \left({\alpha }_{\parallel },{\alpha }_{\perp }\right)$
Functions
auto × auto	−0.572 ± 0.093
cross × cross	−0.470 ± 0.093
auto Lyβ × auto Lyβ	−0.615 ± 0.064
cross Lyβ × cross Lyβ	−0.347 ± 0.052

Download table as:  ASCIITypeset image

5.4.2. Covariance between Different Large-scale Correlations

The correlation functions are fit for separations $r\,\in [10,180]\,{h}^{-1}\,\mathrm{Mpc}$, and for all directions, $\mu \in [0,1]$ for the autocorrelations and μ ∈ [−1, 1] for the cross-correlations. This gives 1590 bins for the two autocorrelations and twice as many, 3180, for both cross-correlations. The different best-fit models differ in their number and the nature of the free parameters. However, they all share the two BAO parameters, $({\alpha }_{\parallel },{\alpha }_{\perp })$, that are the focus of this study.

The autocorrelation baseline model is composed of the Kaiser model for biased matter density tracers, the effect of correlated sky residuals, the Kaiser model for metal-induced correlations, and finally the effect of HCDs. The model has 13 free parameters when fitting separately the Lyα(Lyα)$\,\times \,$Lyα(Lyα) or Lyα(Lyα)$\,\times \,$Lyα(Lyβ) functions. The simultaneous fit for the two functions has 15 free parameters because the two sky subtraction parameters are independent for the two functions, while the other parameters are common to them.

The cross-correlation baseline model is composed of the Kaiser model for quasars, Lyα, and the contamination by metals, along with the modeling of the quasar systemic and statistical redshift errors. Because the RSD parameter of quasars, ${\beta }_{\mathrm{QSO}}$, is highly correlated to the bias and beta parameters of Lyα, we fix it, following Equations (33) and (34), at the effective redshift of the fit. We fix the HCD parameters at the values determined by the combined (auto + cross) fits since they are highly correlated with other parameters, especially the quasar velocity offset. We also fix the quasar proximity effect parameters to those determined by the combined fit. This results in a best-fit model with 10 parameters, for both cross-correlations. Since all these parameters are shared between Lyα(Lyα) × quasar and Lyα(Lyβ) × quasar, the combined fit is composed also of 10 free parameters.

The combined fit to the four measured correlation functions, two auto- and two cross-correlations, is composed of the different models cited in the two previous paragraphs. Running a combined fit to the auto- and cross-correlations allows us to break degeneracies and thus allows us to free ${\beta }_{\mathrm{QSO}}$, to fit for HCDs in the cross-correlation and take into account the proximity effect of quasars onto Lyα absorption. The resulting best-fit model has 19 parameters.

When fitting simultaneously more than one of the four correlation functions, we neglect the cross-covariance between them. The subsampling covariance of Section 3.4 was at the level of 1%, and this was confirmed by the mock-to-mock variations shown in Figure 11. The expected low correlation between BAO parameters measured with different correlation functions was confirmed by the mock-to-mock variations discussed in Section 5.

We follow the prescriptions of dSA19 (their Appendix A) and compute the effective redshift of the measurement to be ${z}_{\mathrm{eff}.}=2.334$, when running the combined fit to all four correlations. This effective redshift is defined to be the pivot redshift for the two BAO parameters ${\alpha }_{\parallel }$ and ${\alpha }_{\perp }$. All parameters are assumed to be redshift independent except for the bias parameters (Section 4.2), which the fitter returns at the effective redshift.

The best-fit model of the individual fit to the autocorrelation Lyα(Lyα) × Lyα(Lyα) is given in 2D in the right panel of Figure 5 and in four wedges of $| \mu |$ in the top four panels of Figure 6. In a similar way the best-fit model of the cross-correlation Lyα(Lyα) × quasar is given in 2D in the right panel of Figure 8 and for four wedges of $| \mu |$ at the bottom of Figure 6. Finally, both auto- and cross-correlations, using the Lyβ region, Lyα(Lyα) × Lyα(Lyβ), and Lyα(Lyβ) × quasar, are given in Figure F1 of Appendix F.

The results of the seven different fits are given in Table 6. The first four columns give the results for each individual measured correlation. The last three columns give combined fits: "full auto," "full cross," and "all combined." ⁴¹ In the first section, the first two lines give the BAO parameters ${\alpha }_{\parallel }$ and ${\alpha }_{\perp }$ then follows the different parameters that describe the overall shape of the correlation functions. The second section gives the attributes of the fit. The last two lines give properties of two different models. The second-to-last line is for a model where the BAO parameters are fixed to the cosmology of Planck (2016). The first number is the difference of χ², and the second number in parentheses gives the conversion to the significance in number of trials, as computed using Monte Carlo realizations (Table 7). The last line gives the significance of the BAO peak.

Table 7. Values of ${\rm{\Delta }}{\chi }^{2}$ Corresponding to $\mathrm{CL}=(68.27,95.45 \% )$

Parameter	${\rm{\Delta }}{\chi }^{2}$ $(68.27 \% )$	${\rm{\Delta }}{\chi }^{2}$ $(95.45 \% )$
Lyα(Lyα)$\,\times \,$Lyα(Lyα)
${\alpha }_{\parallel }$	1.06 ± 0.06	5.06 ± 0.35
${\alpha }_{\perp }$	1.24 ± 0.08	4.51 ± 0.25
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.64 ± 0.09	7.07 ± 0.33
Lyα(Lyα)$\,\times \,$Lyα(Lyβ)
${\alpha }_{\parallel }$	1.25 ± 0.09	5.1 ± 0.41
${\alpha }_{\perp }$	1.39 ± 0.07	5.03 ± 0.26
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.91 ± 0.11	7.9 ± 0.4
Lyα(Lyα)$\,\times \,$quasar
${\alpha }_{\parallel }$	1.21 ± 0.06	4.4 ± 0.22
${\alpha }_{\perp }$	1.22 ± 0.08	4.98 ± 0.26
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.59 ± 0.12	6.79 ± 0.31
Lyα(Lyβ)$\,\times \,$quasar
${\alpha }_{\parallel }$	1.43 ± 0.07	4.5 ± 0.18
${\alpha }_{\perp }$	1.41 ± 0.07	4.74 ± 0.2
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.93 ± 0.13	6.69 ± 0.22
auto all
${\alpha }_{\parallel }$	1.06 ± 0.07	4.83 ± 0.46
${\alpha }_{\perp }$	1.14 ± 0.08	4.86 ± 0.29
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.42 ± 0.1	7.16 ± 0.33
cross all
${\alpha }_{\parallel }$	1.23 ± 0.06	4.3 ± 0.19
${\alpha }_{\perp }$	1.11 ± 0.09	5.02 ± 0.41
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.64 ± 0.11	6.97 ± 0.36
combined
${\alpha }_{\parallel }$	1.05 ± 0.05	4.07 ± 0.37
${\alpha }_{\perp }$	1.05 ± 0.08	4.37 ± 0.27
(${\alpha }_{\parallel }$, ${\alpha }_{\perp }$)	2.4 ± 0.07	6.4 ± 0.27

Note. Values are derived from 1000 Monte Carlo simulations of the correlation function that are fit using the model containing only Lyα absorption. Confidence levels are the fractions of the generated data sets that have best fits below the ${\rm{\Delta }}{\chi }^{2}$ limit. The uncertainties are statistical and estimated using the bootstrap technique.

Download table as:  ASCIITypeset image

In Table 6, parameters without error bars are fixed in the fit. This is the case for the BAO parameters using only the Lyβ region, where the BAO peak is not detected with high significance and we choose to fix $({\alpha }_{\parallel },{\alpha }_{\perp })=(1,1)$. In fits for the cross-correlation we fix the values of some parameters to the values found in the combined fit. For the HCD bias parameters, this is motivated by their high correlation with the quasar velocity distribution parameter σ_v . The parameters ${b}_{\eta ,{\rm{C}}{\rm\small{IV}}(\mathrm{eff})}$, ${\beta }_{\mathrm{QSO}}$, and ${\xi }_{0}^{\mathrm{TP}}$ are fixed because they are not significantly constrained by the cross-only fits.

We can notice that both widths, ${\sigma }_{\mathrm{sky}}$, of the model of the correlated calibration residuals are of the same order of magnitude as that found for the Mg ii(Mg ii(11)) × Mg ii(Mg ii(11)) autocorrelation in Section 4.5. The strength, ${A}_{\mathrm{sky}}$, is, however, 1 mag stronger, due to the difference in effective redshift of the correlation functions between Lyα and Mg ii.

In another aspect we can note that the parameter describing the systematic quasar redshift error, ${\rm{\Delta }}{r}_{\parallel ,\mathrm{QSO}}=0.31\,\pm 0.11\,{h}^{-1}\,\mathrm{Mpc}$, is compatible with zero at the 3σ level. Using DR12 data, dMdB17 measured this parameter to be $-0.79\pm 0.13\,{h}^{-1}\,\mathrm{Mpc}$, which was 6σ from zero. This difference is linked to the change in the quasar redshift estimator (Section 2.2). Other estimators give even larger values, as shown in Table B1. Though this suggests that our new quasar redshift estimator is less biased, it is not possible to draw any definitive conclusion. Indeed, B19 showed that this parameter is highly correlated with models for relativistic effects (Lepori et al. 2020).

Combining the two measurements of the Lyα autocorrelation yields

$$\begin{eqnarray}&&\left\{\begin{array}{l}{D}_{H}(z=2.334)/{r}_{d}=8.93{\,}_{-0.27}^{+0.28}{\,}_{-0.63}^{+0.64}\\ {D}_{M}(z=2.334)/{r}_{d}=37.6{\,}_{-1.9}^{+1.9}{\,}_{-4.0}^{+4.2}\\ \rho \left({D}_{H}(z)/{r}_{d},{D}_{M}(z)/{r}_{d}\right)=-0.49\end{array}\right..\end{eqnarray} \tag{ 41 }$$

Combining the two measurements of the Lyα–quasar cross-correlation yields

$$\begin{eqnarray}&&\left\{\begin{array}{l}{D}_{H}(z=2.334)/{r}_{d}=9.08{\,}_{-0.34}^{+0.34}{\,}_{-0.60}^{+0.61}\\ {D}_{M}(z=2.334)/{r}_{d}=37.3{\,}_{-1.6}^{+1.7}{\,}_{-3.6}^{+4.2}\\ \rho \left({D}_{H}(z)/{r}_{d},{D}_{M}(z)/{r}_{d}\right)=-0.43\end{array}\right..\end{eqnarray} \tag{ 42 }$$

Finally, combining the two autocorrelations along with the two cross-correlations yields

$$\begin{eqnarray}&&\left\{\begin{array}{l}{D}_{H}(z=2.334)/{r}_{d}=8.99{\,}_{-0.19}^{+0.20}{\,}_{-0.38}^{+0.38}\\ {D}_{M}(z=2.334)/{r}_{d}=37.5{\,}_{-1.1}^{+1.2}{\,}_{-2.3}^{+2.5}\\ \rho \left({D}_{H}(z)/{r}_{d},{D}_{M}(z)/{r}_{d}\right)=-0.45\end{array}\right..\end{eqnarray} \tag{ 43 }$$

These three results are presented in Figure 12. The black point gives the Planck (2016) cosmology, while the contours give the 1 and 2 confidence levels, corresponding to (68.27%, 95.45%).

Zoom In Zoom Out Reset image size

As first shown in dMdB17, the BAO parameters $({\alpha }_{\parallel },{\alpha }_{\perp })$ have non-Gaussian uncertainties. This means that, e.g., Δχ² = 1 does not correspond to 68% of trials for one degree of freedom. To estimate the mapping between (68.27%, 95.45%) confidence levels to Δχ² values, we generate 1000 fast Monte Carlo (fastMC) realizations of our seven fits. Each fit is reduced to the simple Kaiser model, in order to build enough statistics given the available computation time. A complete model would give similar results, as shown by dMdB17. In each realization, we generate a random noisy measurement of the best-fit model, given the covariance matrix, with BAO parameters set to 1. Then, the fit is performed in four different combinations: for the BAO parameters free, another for both of them fixed, then two more for one fixed and the other free. This allows for each of the seven results to give the error bars for α_∥ and α_⊥ independently. Table 7 gives the results of the estimation of the error bars. The BAO error bars given in Table 6 are the results of this fastMC correction.

The distribution of ${\rm{\Delta }}{\chi }^{2}({\alpha }_{\parallel }={\alpha }_{\perp }=1)$ in the 1000 fastMC allows us to compute the confidence levels in two dimensions. Furthermore, it allows us to estimate how significant the shift of our best-fit measurement is against the Planck (2016) cosmology. This result is given in Table 6 in the second-to-last line. Our final measurement, "all combined," is 1.5σ from the cosmology of Planck (2016).

The significance of the BAO detection is estimated by leaving free the parameter A_BAO that describes the size of the BAO peak relative to what is expected according to the Planck (2016) cosmology. We give the results for this extension to the baseline model in the last line of Table 6. Using fastMC realizations, we verify that this parameter is linear in its mapping from Δχ² to confidence levels. Using 1000 fastMC, we verify this property up to 3σ and then assume that it holds true for higher confidence levels.

Our measurements of the autocorrelation, Lyα(Lyα)$\,\times \,$Lyα(Lyα), and cross-correlation, Lyα(Lyα)$\,\times \,$quasar, both have a BAO peak that has a significance of higher than 4σ. It is barely 3σ for the autocorrelation Lyα(Lyα)$\,\times \,$Lyα(Lyβ), and barely 2σ for the Lyα(Lyβ)$\,\times \,$quasar. Because of this low significance of the BAO measurement in these two extra correlations, the 2σ error bars are very nonlinear, with Lyα(Lyβ)$\,\times \,$quasar having none. We thus choose not to give the BAO measurement for these two correlations because it should not be used alone, but only in combination with other measurements. The final BAO measurement that combines all measured correlation functions has a BAO significance of more than 7σ.

Our measured $({\alpha }_{\parallel },{\alpha }_{\perp })$ are marginally correlated to other parameters but significantly correlated, ρ ∼ −50%, with one another. The only other parameter that is correlated with $({\alpha }_{\parallel },{\alpha }_{\perp })$ is the bias parameter of the Si ii(126) metal line. Indeed, the correlations of absorption by this metal and Lyα absorption at the same physical position generate an apparent peak in the reconstructed correlation function at $({r}_{\parallel },{r}_{\perp })\,\sim (+111,0)\,{h}^{-1}\,\mathrm{Mpc}$ (Table 3). This correlation is ≈ − 22% with ${\alpha }_{\parallel }$ and +9% with α_⊥ in the case of the "full auto." Since in the cross-correlations we have access to positive and negative distances along the line of sight, this correlation is reduced to −14% with α_∥ and +5% with α_⊥ for the "full cross." Combining all results gives a correlation of −19% with α_∥ and +8% with α_⊥ for the "all combined" fit. All other parameters are less than ±3% correlated with the BAO parameters.

In Appendix D, we test for systematic errors in our measurement of the BAO parameters, along with changes in the estimation of their error bars. We either modify the best-fit model in Table D1 or change some aspects of the analysis or study data splits in Table D2. In Table D1 the strongest effect on both the recovered BAO best-fit value and error bars happens when taking the effect of metals into account, which results in a 0.5σ shift in ${\alpha }_{\parallel }$. As previously discussed, the Si ii(126) absorption line produces a scale along the line of sight of $+111\,{h}^{-1}\,\mathrm{Mpc}$. The change in the error bars is explained by the effect of preventing the best-fit model to build false BAO significance from this line.

Of the fits in Table D1, those that add polynomial "broadband" curves to the model are of special importance because they test the sensitivity of the BAO parameters to unidentified systematic errors in the model. We performed two fits, allowing the broadband curves different amounts of freedom in each. The first placed "physical priors" on $({b}_{\mathrm{Ly}\alpha },{\beta }_{\mathrm{Ly}\alpha },{b}_{\mathrm{HCD}})$ in the form of a Gaussian of mean and width of the fit without broadband terms. Such priors ensured that the broadband terms were relatively small perturbations to the physical model. The second type of fit placed no priors on $({b}_{\mathrm{Ly}\alpha },{\beta }_{\mathrm{Ly}\alpha },{b}_{\mathrm{HCD}})$. The results of these fits are given in Table D1. We see that the addition of such terms does not change significantly the values of $({\alpha }_{\parallel },{\alpha }_{\perp })$ in any of the fits.