Measurement of BAO correlations at $z=2.3$ with SDSS DR12 \lya-Forests

TL;DR

This study measures baryon acoustic oscillations at redshift 2.3 using SDSS quasar spectra, improving analysis methods to refine cosmological distance measurements and test cosmological models.

Contribution

It introduces enhanced analysis techniques for BAO detection in Lyman-alpha forests, providing more precise measurements and systematic error assessments at high redshift.

Findings

Measured Hubble distance at z=2.33 as 9.07 ± 0.31 times the sound horizon.

Found the combined distance measure to be 13.94 ± 0.35, consistent with flat- extLambda CDM.

Constrained cosmological parameters: Ω_m=0.296 ± 0.029, Ω_Λ=0.699 ± 0.100, Ω_k≈0.

Abstract

We use flux-transmission correlations in \Lya forests to measure the imprint of baryon acoustic oscillations (BAO). The study uses spectra of 157,783 quasars in the redshift range $2.1\le z \le 3.5$ from the Sloan Digital Sky Survey (SDSS) Data Release 12 (DR12). Besides the statistical improvements on our previous studies using SDSS DR9 and DR11, we have implemented numerous improvements in the analysis procedure, allowing us to construct a physical model of the correlation function and to investigate potential systematic errors in the determination of the BAO peak position. The Hubble distance, $\DHub=c/H(z)$, relative to the sound horizon is $\DHub(z=2.33)/r_d=9.07 \pm 0.31$. The best-determined combination of comoving angular-diameter distance, $\DM$, and the Hubble distance is found to be $\DHub^{0.7}\DM^{0.3}/r_d=13.94\pm0.35$. This value is $1.028\pm0.026$ times the prediction of the flat-\lcdm model consistent with the cosmic microwave background (CMB) anisotropy spectrum. The errors include marginalization over the effects of unidentified high-density absorption systems and fluctuations in ultraviolet ionizing radiation. Independently of the CMB measurements, the combination of our results and other BAO observations determine the open-\lcdm density parameters to be $\om=0.296 \pm 0.029$, $\ol=0.699 \pm 0.100$ and $\Omega_k = -0.002 \pm 0.119$.

Tables

Table 1. Table 1: Parameters of the flat Λ Λ \Lambda CDM cosmological model used for the production and analysis of the mock spectra and the CMB-based flat Λ Λ \Lambda CDM model from Planck Collaboration et al. ( 2016 ) used for the analysis of the data. The models are defined by the cold dark matter, baryon, and massive neutrinos densities, the Hubble constant, the number of light neutrino species, and the inhomogeneity parameters, σ 8 subscript 𝜎 8 \sigma_{8} and n s subscript 𝑛 𝑠 n_{s} . The sound horizon, r d subscript 𝑟 𝑑 r_{d} , is calculated using CAMB (Lewis et al., 2000 ) .

	Mocks	Planck
		(TT+lowP)
${\Omega }_M{h}^2$	0.1323	0.1426
$={\Omega }_C{h}^2$	0.1096	0.1197
$+{\Omega }_B{h}^2$	0.0227	0.02222
$+{\Omega }_{\nu }{h}^2$	0	0.0006
$h$	0.7	0.6731
$N_{\nu }$	3	3
${\sigma }_8$	0.795	0.829
$n_s$	0.97	0.9655
${\Omega }_M$	0.27	0.3147
$r_d$ [Mpc]	149.7	147.33
	(104.80 $h^{-1}$)	(99.17 $h^{-1}$)
$D_M\left(2.33\right)/r_d$	38.63	39.15
$D_H\left(2.33\right)/r_d$	8.744	8.612

Table 2. Table 2: Transitions contributing to absorption by metals in the mock spectra. The second column shows r ∥ = ( 1 + z ¯ ) D H ( z ¯ ) ( λ m ​ e ​ t − λ Ly ​ α / ) λ r_{\parallel}=(1+\bar{z})D_{H}(\bar{z})(\lambda_{met}-\lambda_{\rm Ly\alpha}/)\lambda , the reconstructed separation corresponding to correlations of Ly α 𝛼 \alpha  and metal absorption at the same physical position.

transition	$r_{\parallel }$ [$h^{-1}\mathrm{Mpc}$]
SiIII(120.7 nm)	21
SiIIa(119.0 nm)	64
SiIIb(119.3 nm)	56
SiIIc(126.0 nm)	111

Table 3. Table 3: Parameters of the physical model for the correlation function and their best-fit values for the data fit over the range 10 < r < 180 ​ h − 1 ​ Mpc 10 𝑟 180 superscript ℎ 1 Mpc 10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc} . The parameters include the BAO peak-position parameters, ( α ⟂ , α ∥ ) subscript 𝛼 perpendicular-to subscript 𝛼 parallel-to (\alpha_{\perp},\alpha_{\parallel}) , the bias parameters ( b , β ) 𝑏 𝛽 (b,\beta) at redshift z = 2.3 𝑧 2.3 z=2.3 for the various components, the length scale, L HCD subscript 𝐿 HCD L_{\rm HCD} , for unidentified HCDs, and the bias parameter, b Γ subscript 𝑏 Γ b_{\Gamma} for UV fluctuations. The biases b HCD subscript 𝑏 HCD b_{\rm HCD} and b Γ subscript 𝑏 Γ b_{\Gamma} are relative to b Ly ​ α subscript 𝑏 Ly 𝛼 b_{\rm Ly\alpha} . Parameters without uncertainties are fixed in the fit. These include the UV absorber response bias, b a ′ superscript subscript 𝑏 𝑎 ′ b_{a}^{\prime} , the UV attenuation length, λ UV subscript 𝜆 UV \lambda_{\rm UV} , the BAO amplitude parameter, A peak subscript 𝐴 peak A_{\rm peak} , the bias redshift-evolution parameter, γ 𝛾 \gamma , and the peak-broadening parameters Σ @vec ⁡ Σ \@vec{\Sigma} .

parameters	best fit	best fit
${\alpha }_{\parallel },{\alpha }_{⟂}$	$1.053\pm 0.036$	$0.965\pm 0.055$
$b\left(1+\beta \right),\beta$ Ly$\alpha$	$-0.325\pm 0.004$	$1.663\pm 0.085$
$b,\beta$ Si3	$-0.0033\pm 0.0013$	$0.500$
$b,\beta$ Si2a	$-0.0044\pm 0.0009$	$0.500$
$b,\beta$ Si2b	$-0.0035\pm 0.0009$	$0.500$
$b,\beta$ Si2c	$-0.0015\pm 0.0012$	$0.500$
$b,\beta$ CIV	$-0.0233\pm 0.0114$	$0.500$
$b,\beta$ HCDs	$-0.0288\pm 0.0043$	$0.681\pm 0.175$
$L_{\mathrm{HCD}}$	$24.3410\pm 1.1308$
$b_{\Gamma },b_a^{\prime }$	$0.1125\pm 0.0512$	$-2/3$
${\lambda }_{\mathrm{UV}}\left[h^{-1}\mathrm{Mpc}\right]$	$300.00$
$A_{\mathrm{peak}}$	$1.00$
$\gamma$	$2.90$
${\Sigma }_{⟂},{\Sigma }_{\parallel }$	$3.26$	$6.41$

Table 4. Table 4: Weighted mean of fit parameters and mean uncertainties for six sets of 100 mocks of increasing realism: Ly α 𝛼 \alpha  only; including a quasar continuum; including metal absorption (for Met1 or Met2); including HCDs (and Met1) and masking those with N HI > 10 20 ​ cm − 2 subscript 𝑁 HI superscript 10 20 superscript cm 2 N_{\rm HI}>10^{20}{\rm cm^{-2}} or > 10 21 ​ cm − 2 absent superscript 10 21 superscript cm 2 >10^{21}{\rm cm^{-2}} . Fits are over the range 10 < r < 180 ​ h − 1 ​ Mpc 10 𝑟 180 superscript ℎ 1 Mpc 10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc} . The bias parameter, b Ly ​ α subscript 𝑏 Ly 𝛼 b_{\rm Ly\alpha} , refers to the reference redshift z = 2.25 𝑧 2.25 z=2.25 . The input values for the mock generation were b ​ ( 1 + β ) Ly ​ α = − 0.336 𝑏 subscript 1 𝛽 Ly 𝛼 0.336 b(1+\beta)_{\rm Ly\alpha}=-0.336 and β Ly ​ α = 1.4 subscript 𝛽 Ly 𝛼 1.4 \beta_{\rm Ly\alpha}=1.4 .

mock set	$\overline{{\alpha }_{\parallel }}\left(\overline{{\sigma }_{{\alpha }_{\parallel }}}\right)$	$\overline{{\alpha }_{⟂}}\left(\overline{{\sigma }_{{\alpha }_{⟂}}}\right)$	$\overline{b{\left(1+\beta \right)}_{\mathrm{Ly}\alpha }}\left(\overline{{\sigma }_{b\left(1+\beta \right)}}\right)$	$\overline{{\beta }_{\mathrm{Ly}\alpha }}\left(\overline{{\sigma }_{{\beta }_{\mathrm{Ly}\alpha }}}\right)$	$\overline{{\chi }_{\min }^2}/DOF$    $\overline{Prob.}$
Ly$\alpha$ only	$0.998\left(0.014\right)$	$1.002\left(0.020\right)$	$-0.341\left(0.001\right)$	$1.360\left(0.017\right)$	$1632.5/\left(1590-4\right)p=0.292$
+continuum	$1.000\left(0.023\right)$	$0.993\left(0.040\right)$	$-0.338\left(0.002\right)$	$1.339\left(0.034\right)$	$1589.1/\left(1590-4\right)p=0.487$
+metals (Met1)	$0.998\left(0.025\right)$	$0.993\left(0.040\right)$	$-0.337\left(0.002\right)$	$1.334\left(0.038\right)$	$1583.8/\left(1590-12\right)p=0.473$
(or) +metals (Met2)	$1.002\left(0.023\right)$	$0.992\left(0.039\right)$	$-0.339\left(0.002\right)$	$1.334\left(0.037\right)$	$1582.8/\left(1590-12\right)p=0.480$
+HCDs (20)	$0.995\left(0.032\right)$	$0.999\left(0.058\right)$	$-0.351\left(0.035\right)$	$1.338\left(0.172\right)$	$1590.0/\left(1590-15\right)p=0.429$
(or) +HCDs (21)	$0.995\left(0.032\right)$	$0.995\left(0.056\right)$	$-0.359\left(0.026\right)$	$1.321\left(0.145\right)$	$1594.1/\left(1590-15\right)p=0.417$

Table 5. Table 5: Results of fits of the data assuming increasingly complicated forms for ξ ph ​ ( r ⟂ , r ∥ ) superscript 𝜉 ph subscript 𝑟 perpendicular-to subscript 𝑟 parallel-to \xi^{\rm ph}(r_{\perp},r_{\parallel}) : Ly α 𝛼 \alpha  absorption only; including metal absorption; including unidentified HCDs; and including fluctuations in the ionizing UV flux (i.e. the complete physical model whose complete set of parameters is given in Table 3 ). Fits are over the range 10 < r < 180 ​ h − 1 ​ Mpc 10 𝑟 180 superscript ℎ 1 Mpc 10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc} . The bias parameter, b Ly ​ α subscript 𝑏 Ly 𝛼 b_{\rm Ly\alpha} , refers to the reference redshift z = 2.3 𝑧 2.3 z=2.3 .

analysis	${\alpha }_{\parallel }$	${\alpha }_{⟂}$	$b_{\mathrm{Ly}\alpha }\left(1+{\beta }_{\mathrm{Ly}\alpha }\right)$	${\beta }_{\mathrm{Ly}\alpha }$	${\chi }_{\min }^2/DOF$, prob
Ly$\alpha$	$1.040\pm 0.033$	$0.975\pm 0.056$	$-0.326\pm 0.002$	$1.246\pm 0.044$	$1763.8/\left(1590-4\right)p=0.001$
+metals	$1.050\pm 0.035$	$0.967\pm 0.054$	$-0.330\pm 0.002$	$1.275\pm 0.045$	$1644.5/\left(1590-9\right)p=0.130$
+HCD	$1.053\pm 0.036$	$0.962\pm 0.054$	$-0.321\pm 0.003$	$1.656\pm 0.086$	$1561.4/\left(1590-12\right)p=0.612$
+UV	$1.053\pm 0.036$	$0.965\pm 0.055$	$-0.326\pm 0.003$	$1.666\pm 0.085$	$1556.5/\left(1590-13\right)p=0.639$

Table 6. Table 6: Best-fit values of ( α ⟂ , α ∥ ) subscript 𝛼 perpendicular-to subscript 𝛼 parallel-to (\alpha_{\perp},\alpha_{\parallel}) for fits including metals, HCDs and UV fluctuations, with and without power-law broadbands (BBs) of the form ( i min , i max , j max ) = ( 0 , 2 , 6 ) subscript 𝑖 min subscript 𝑖 max subscript 𝑗 max 0 2 6 (i_{\rm min},i_{\rm max},j_{\rm max})=(0,2,6) . The fit without broadbands is over 10 < r < 180 ​ h − 1 ​ Mpc 10 𝑟 180 superscript ℎ 1 Mpc 10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc} and those with broadbands are over 40 < r < 180 ​ h − 1 ​ Mpc 40 𝑟 180 superscript ℎ 1 Mpc 40<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc} . The broadband fit with “physical priors” requires that the best-fit values for parameters be within two standard deviations of those found without the broadband and imposes A peak = 1 subscript 𝐴 peak 1 A_{\rm peak}=1 . The broadband fit with “no additional priors” adds no constraints beyond fixing the parameters that are fixed in the fit without broadband (see Table 3 ).

analysis	${\alpha }_{\parallel }$	${\alpha }_{⟂}$	${\alpha }_{\parallel }^{0.7}{\alpha }_{⟂}^{0.3}$	${\chi }_{\min }^2/DOF$, prob.
no BB	${1.053}_{-0.035-0.071}^{+0.037+0.077}$	${0.965}_{-0.051-0.102}^{+0.061+0.146}$	${1.026}_{-0.024-0.047}^{+0.025+0.052}$	$1556.5/\left(1590-13\right)\mathrm{\hspace{0.33em}0.639}$
BB, physical priors	${1.057}_{-0.037-0.073}^{+0.039+0.081}$	${0.962}_{-0.053-0.106}^{+0.063+0.155}$	${1.028}_{-0.025-0.049}^{+0.026+0.054}$	$1475.8/\left(1515-25\right)\mathrm{\hspace{0.33em}0.598}$
BB, no additional priors	${1.043}_{-0.033-0.068}^{+0.035+0.075}$	${1.000}_{-0.060-0.114}^{+0.095+0.346}$	${1.030}_{-0.021-0.043}^{+0.023+0.063}$	$1472.6/\left(1515-26\right)\mathrm{\hspace{0.33em}0.614}$

Table 7. Table 7: Systematic shifts of the Ly α 𝛼 \alpha  bias parameters and BAO peak-position parameters due to various spurious correlations discussed in Sect. 5 . Individual shifts in ( α ⟂ , α ∥ ) subscript 𝛼 perpendicular-to subscript 𝛼 parallel-to (\alpha_{\perp},\alpha_{\parallel}) are less than σ / 10 𝜎 10 \sigma/10 and the total shifts are consistent with zero. Their uncertainty is dominated by the statistical uncertainty on the sky model.

	${\beta }_{\mathrm{Ly}\alpha }$	$b\left(1+\beta \right)$	${\alpha }_{\parallel }$	${\alpha }_{⟂}$
Sky model noise	$-0.026$	$-0.002$
Calibration noise	$+0.047$	$+0.002$		$+0.001$
Fiber cross-talk	$+0.003$
ISM absorption	$+0.003$
Sum	$+0.027$		$+0.001$
Quadratic sum	$+0.054$	$+0.002$		$+0.001$

Equations

$$104.0<{\lambda }_{\mathrm{RF}}<120.0\text{leavevmode}\mathrm{nm}.$$

$$F({\delta }_G)=\exp \left[-a(z)e^{b(z){\delta }_G}\right].$$

$${\delta }_q(\lambda )=\frac{f_q(\lambda )}{C_q(\lambda )\overline{F}(z)}-1.$$

$$C_q(\lambda )=C({\lambda }_{\mathrm{RF}})(a_q+b_q\log (\lambda ))$$

$$L=\prod _{q,\lambda }P(f_q(\lambda )|C_q(\lambda )).$$

$$P(f_q(\lambda )|C_q(\lambda ))\propto {\int }_0^{1}dFD(F,z)\exp \left[\frac{-(C_{q}F-f_q(\lambda ){)}^2}{2{\sigma }_q(\lambda {)}^2}\right].$$

$${\tilde{\delta }}_q(\lambda )={\delta }_q(\lambda )-\overline{\delta (\lambda )}$$

$${\hat{\delta }}_q(\lambda )={\tilde{\delta }}_q(\lambda )-\overline{{\tilde{\delta }}_q}-(\Lambda -\overline{{\Lambda }_q})\frac{\overline{(\Lambda -\overline{{\Lambda }_q})\tilde{\delta }}}{\overline{(\Lambda -\overline{{\Lambda }_q}{)}^2}}\Lambda \equiv \log \lambda$$

$${\hat{\xi }}_A=\frac{{\sum }_{ij\in A}{w}_i{w}_j{\hat{\delta }}_i{\hat{\delta }}_j}{{\sum }_{ij\in A}{w}_i{w}_j},$$

$${\hat{\xi }}_A=\sum _B{D}_{AB}{\xi }_B$$

$$D_{AB}=w_A^{-1}\sum _{i,j\in A}{w}_i{w}_j\left[\sum _{i^{\prime },j^{\prime }\in B}{\eta }_{i{i}^{\prime }}{\eta }_{j{j}^{\prime }}\right]$$

$${\eta }_{i{i}^{\prime }}={\delta }_{i{i}^{\prime }}-\frac{w_{i^{\prime }}}{{\sum }_k{w}_k}-\frac{({\Lambda }_i-\overline{\Lambda })w_{i^{\prime }}({\Lambda }_{i^{\prime }}-\overline{\Lambda })}{{\sum }_k{w}_k({\Lambda }_k-\overline{\Lambda }{)}^2}\Lambda =\log \lambda$$

$$C_{AB}=\frac{1}{W_A{W}_B}\sum _{ij\in A}\sum _{kl\in B}{w}_i{w}_j{w}_k{w}_l[⟨{\hat{\delta }}_i{\hat{\delta }}_j{\hat{\delta }}_k{\hat{\delta }}_l⟩-⟨{\hat{\delta }}_i{\hat{\delta }}_j⟩⟨{\hat{\delta }}_k{\hat{\delta }}_l⟩]$$

$$C_{AB}=\frac{1}{W_A{W}_B}\sum _s{W}_A^s{W}_B^s[{\hat{\xi }}_A^s{\hat{\xi }}_B^s-{\hat{\xi }}_A{\hat{\xi }}_B],$$

$$C_{AB}=\frac{1}{W_A{W}_B}\sum _{ij\in A}\sum _{kl\in B}{w}_i{w}_j{w}_k{w}_l{\xi }_{1d}({\lambda }_i/{\lambda }_k){\xi }_{1d}({\lambda }_j/{\lambda }_l)$$

$$C_{AB}=\overline{{\hat{\xi }}_A{\hat{\xi }}_B}-\overline{{\hat{\xi }}_A}\overline{{\hat{\xi }}_B},$$

$$C(r_{\perp },r_{\parallel },r_{\perp },r_{\parallel })\sim \frac{0.043}{N_{pairs}}\sim \frac{2.5\times 10^{-10}}{r_{\perp }/100\text{leavevmode}{h}^{-1}\mathrm{Mpc}}.$$

$$f_q(\lambda )=A_q(\lambda )[f_q^{\mathrm{m}}(\lambda )-B_q(\lambda )].$$

$${\xi }_A=\sum _{A^{\prime }}{D}_{A{A}^{\prime }}\left[{\xi }_{A^{\prime }}^{\mathrm{ph}}+B_{A^{\prime }}\right]$$

$${\xi }_{\mathrm{Ly}\alpha }(r_{\perp },r_{\parallel },{\alpha }_{\perp },{\alpha }_{\parallel })={\xi }_{\mathrm{smooth}}(r_{\perp },r_{\parallel })+{\xi }_{\mathrm{peak}}({\alpha }_{\perp }{r}_{\perp },{\alpha }_{\parallel }{r}_{\parallel })$$

$${\alpha }_{\parallel }=\frac{\left[D_H(\bar{z})/r_d\right]}{{\left[D_H(\bar{z})/r_d\right]}_{\mathrm{fid}}}\mathrm{and}{\alpha }_{\perp }=\frac{\left[D_M(\bar{z})/r_d\right]}{{\left[D_M(\bar{z})/r_d\right]}_{\mathrm{fid}}}\text{leavevmode},$$

$$P_{\mathrm{Ly}\alpha }(\text{@vec}k,z)=P_{\mathrm{QL}}(\text{@vec}k,z)b_{\mathrm{Ly}\alpha }^2(1+{\beta }_{\mathrm{Ly}\alpha }{\mu }_{\text{@vec}k}^2{)}^2{F}_{\mathrm{NL}}(\text{@vec}k)G(\text{@vec}k)$$

$$P_{\mathrm{QL}}(\text{@vec}k,z)=P_{\mathrm{smooth}}(k,z)+e^{-k^2{\Sigma }^2({\mu }_k)/2}{A}_{\mathrm{peak}}{P}_{\mathrm{peak}}(k,z)$$

$${\xi }_A^{\mathrm{ph}}=\sum _{m,n}{\xi }_A^{mn}$$

$$\frac{{\xi }^{\mathrm{q}-\mathrm{CIV}}(z=1.65)}{{\xi }^{\mathrm{q}-\mathrm{Ly}\alpha }(z=2.3)}=0.08\pm 0.02r=10h^{-1}\mathrm{Mpc}$$

$$\frac{{\xi }_{z=1.65}^{\mathrm{q}-\mathrm{CIV}}}{{\xi }_{z=2.3}^{\mathrm{q}-\mathrm{Ly}\alpha }}=\frac{[b_{\mathrm{CIV}}{b}_{\mathrm{q}}{\xi }_{\mathrm{L}}{]}_{z=1.65}}{[b_{\mathrm{Ly}\alpha }{b}_{\mathrm{q}}{\xi }_{\mathrm{L}}{]}_{z=2.3}}\frac{1+({\beta }_{\mathrm{CIV}}+{\beta }_{\mathrm{q}})/3+{\beta }_{\mathrm{CIV}}{\beta }_{\mathrm{q}}/5}{1+({\beta }_{\mathrm{Ly}\alpha }+{\beta }_{\mathrm{q}})/3+{\beta }_{\mathrm{Ly}\alpha }{\beta }_{\mathrm{q}}/5}$$

$$B(r,\mu )=\sum _{j=0}^{j_{\max }}\sum _{i=i_{\min }}^{i_{\max }}{a}_{ij}\frac{L_j(\mu )}{r^i}(j\mathrm{even}),$$

$${\alpha }_{\parallel }=1.053\pm 0.036\Rightarrow \frac{D_H(z=2.33)}{r_d}=9.07\pm 0.31$$

$${\alpha }_{\perp }=0.965\pm 0.055\Rightarrow \frac{D_M(z=2.33)}{r_d}=37.77\pm 2.13.$$

$${\alpha }_{\parallel }^{0.7}{\alpha }_{\perp }^{0.3}=1.028\pm 0.026,$$

Full text

1 Introduction

The sound waves that propagated in the pre-recombination universe produced a pronounced peak in the two-point correlation function of the cosmological-density field (Peebles & Yu, 1970; Sunyaev & Zeldovich, 1970). This “baryon-acoustic oscillation” (BAO) peak, first observed by Eisenstein et al. (2005) and Cole et al. (2005), is centered on a comoving distance equal to the sound horizon at the drag epoch, $r_{d}$. Observed at a redshift $z$, the BAO peak position in the transverse (angular) direction determines the ratio $D_{M}(z)/r_{d}$, where $D_{M}(z)=(1+z)D_{A}(z)$ is the “comoving” angular-diameter distance ($D_{A}$ is the “traditional” angular-diameter distance). In the radial (redshift) direction, the peak position determines $D_{H}(z)/r_{d}$, where $D_{H}(z)=c/H(z)$ is the Hubble distance. Since $D_{M}$, $D_{H}$ and $r_{d}$ all have simple dependencies on the cosmological parameters, observation of the BAO feature constrains those parameters, especially when combined with cosmic microwave background (CMB) data (Planck Collaboration et al., 2014, 2016; Aubourg et al., 2015). In particular, one can constrain models beyond the flat-$\Lambda$CDM model that describes CMB data, deriving constraints on cosmological curvature and the dark-energy equation of state. Furthermore, the shape of the spectrum of CMB anisotropies can be used to calculate the value of $r_{d}$ to percent-level precision, $r_{d}\sim 147$ Mpc and the use of this value allows one to derive $D_{M}(z)$ and $D_{H}(z)$ from BAO measurements. These absolute distances can be combined with relative distances determined with type Ia supernovae (Betoule et al., 2014), to extrapolate to $z=0$, yielding a “top-down” measurement of $H_{0}$ (Aubourg et al., 2015).

Most studies of the BAO peak have used galaxies at redshifts $z<0.8$ as tracers of the density. The first observations used $z\sim 0.35$ data from the Sloan Digital Sky Survey (SDSS) (Eisenstein et al., 2005) and $z\sim 0.2$ data from the Two-Degree Field Redshift Survey (Cole et al., 2005), and the combination of these sets (Percival et al., 2007, 2010). Since then, results of increasingly higher precision have been obtained, most significantly in the redshift range $0.35<z<0.65$ from the Baryon Oscillation Spectroscopy Survey (BOSS) of SDSS-III (Anderson et al., 2012, 2014b, 2014a) with the results of the complete survey being summarized in Alam et al. (2016). Results at other redshifts have been obtained by 6dFGRS at $z\sim 0.11$ (Beutler et al., 2011), WiggleZ at $0.4<z<0.8$ (Blake et al., 2011), SDSS-I at $z\sim 0.35$ (Padmanabhan et al., 2012; Mehta et al., 2012; Chuang & Wang, 2012; Xu et al., 2013) and SDSS-I at $z\sim 0.15$ (Ross et al., 2015). There is an impressive agreement of the results of these studies with the expectations of flat-$\Lambda$CDM models based on CMB data, as emphasized by Planck Collaboration et al. (2016).

At higher redshifts, BAO correlations can be seen using quasars and their Lyman-$\alpha$ (Ly$\alpha$) forests as mass tracers (McDonald & Eisenstein, 2007). The correlations in the Ly$\alpha$-forest flux-transmission field of BOSS quasars were first studied in Slosar et al. (2011) and the BAO peak was seen in SDSS data release DR9 (Busca et al., 2013; Slosar et al., 2013; Kirkby et al., 2013) and DR11 (Delubac et al., 2015). Cross-correlations of the Ly$\alpha$ absorption with the distribution of quasars were detected in DR9 (Font-Ribera et al., 2013), and the first BAO detection was presented in Font-Ribera et al. (2014).

In this paper, we study the auto-correlation function, $\xi$, of the Ly$\alpha$ flux-transmission field using the SDSS data release DR12 and we update the cosmological constraints reported in Delubac et al. (2015). Our study of the quasar-forest cross-correlation function will be presented in a future publication (du Mas des Bourboux et al., in preparation). In addition to using the $15\%$ increase of survey area of DR12 over DR11, the following improvements over the analysis of Delubac et al. (2015) have been implemented:

• •

Spectroscopic pipeline improvements (Sect. 2) that includes a new algorithm for the extraction of spectra from the CCD images that results in a more linear flux response. This modification allows us to correct for the mean distortion of the flux-transmission field due to imperfect spectral modeling of standard stars. We also correct for the differential positioning of quasar and stellar fibers in the focal plane due to optimization at different wavelengths (Margala et al., 2016).

• •

The use of mock spectra with improved modeling of metal absorbers (Sect. 3), including both Ly$\alpha$-metal and metal-metal correlations.

• •

Modeling of the distortions of the correlation function due to quasar-continuum fitting (Sect. 4). This allows us to construct a physical model of the correlation function over the range $10<r<200\leavevmode\nobreak\ h^{-1}{\rm Mpc}$.

• •

Modeling of spurious correlations introduced by the pipeline (Sect. 5), calculation of their effect on the correlation function, and searches for unidentified spurious correlations using the Carbon-IV (CIV) forest, $142.0<\lambda_{\rm RF}<152.0$ nm.

• •

Fits of the data (Sect. 6) that marginalize over the contributions to the correlation function of metals, unidentified high-column-density systems, and fluctuations of ionizing UV flux.

None of these improvements induce significant changes in the derived values of $D_{H}/r_{d}$ and $D_{M}/r_{d}$; the differences (Sect. 8) with those of Delubac et al. (2015) are consistent with statistical fluctuations induced by the increased sample size.

This paper is organized as follows. Section 2 describes the DR12 data used in this analysis. Section 3 gives a brief description of the mock spectra used to test the analysis procedure, with a more detailed description being found in Bautista et al. (2015). Section 4 presents our method of estimating the flux-transmission field, its correlation function, and the associated covariance matrix. Section 5 studies spurious correlations induced by the pipeline. In Section 6 we fit the mocks and data to derive the BAO peak position parameters, $D_{M}(z=2.33)/r_{d}$ and $D_{H}(z=2.33)/r_{d}$. Section 7 investigates possible systematic errors in the measurement. In Section 8 we compare our measured peak position with the predictions of $\Lambda$CDM models and derive constraints on the associated cosmological parameters. Section 9 presents a brief summary.

2 The BOSS quasar sample and data reduction

The quasar sample was gathered over a five-year period by the SDSS-III Collaboration (Eisenstein et al., 2011; Gunn et al., 1998, 2006; Smee et al., 2013). We use the data from the twelfth Data Release (DR12) of SDSS as presented in Alam et al. (2015). The associated quasar catalog is described in Pâris et al. (2016). Most of the quasar spectra were obtained by the Baryon Oscillation Spectroscopic Survey, BOSS (Dawson et al., 2013), but DR12 also includes six months of data from the SEQUELS 111http://www.sdss.org/dr12/algorithms/ancillary/boss/sequels/ program. The DR12 celestial footprint covering $\sim\pi\leavevmode\nobreak\ \rm{sr}\sim 10^{4}\leavevmode\nobreak\ {\rm deg^{2}}$ is displayed in Fig. 1. An example of a quasar spectrum in the Ly$\alpha$-forest region is shown in Fig. 2. The redshift distribution of measurement pairs in the forest is shown in Fig. 3.

The quasar target selection used in BOSS, summarized in Ross et al. (2012), combines different targeting methods described in Yèche et al. (2010), Kirkpatrick et al. (2011), and Bovy et al. (2011). The selection algorithms use SDSS photometry and, when available, data from the GALEX survey (Martin et al., 2005) in the UV; the UKIDSS survey (Lawrence et al., 2007) in the NIR, and the FIRST survey (Becker et al., 1995) in the radio.

The DR12 data were processed using a new software package that differs from the standard DR12 SDSS-III pipeline (Bolton et al., 2012) and which has become the standard pipeline for SDSS DR13 (SDSS Collaboration et al., 2016). For each object, both pipelines provide a flux calibrated spectrum, $f(\lambda)$, errors, and an object classification (galaxy, quasar, star). A model spectrum is fit to $f(\lambda)$ providing a redshift estimate. For this study, we use the “coadded” spectra constructed from typically four exposures of 15 minutes resampled at wavelength pixels of width $\Delta\log_{10}\lambda=10^{-4}$ ($c\Delta\lambda/\lambda\sim{\rm 69\leavevmode\nobreak\ km\,s^{-1}}$). For the small number of quasars with repeated observations, the coadded spectra can include exposures widely separated in time.

An important difference with respect to the DR12 pipeline is that pixels on the CCD image are combined to give a flux $f(\lambda)$ with pixel-weights determined only by the CCD readout noise. While this method is sub-optimal because it ignores photo-electron Poisson noise, compared to the DR12 method it yields an unbiased estimate of $f(\lambda)$ since the weights do not depend on the observed CCD counts which are needed to estimate Poisson noise. A more detailed description of the changes to the extraction pipeline is given in Appendix A.

The ratio of observed flux to model flux, averaged over all spectra, is an important diagnostic of pipeline systematic errors. Figure 4 shows the average ratio, $R(\lambda)$, as a function of observed wavelength for quasar spectra on the red side of Ly$\alpha$ emission. Since model imperfections for individual quasar spectra are averaged over in this figure (because of the range of quasar redshifts), the deviations from unity reflect imperfections in the spectrograph flux calibration. The figure reveals percent-level deviations that are mostly due to imperfect modeling of photo-spectroscopic standard stars. The calcium H and K lines from interstellar absorption are also visible. In the analysis to be presented here, the $f(\lambda)$ are divided by $R(\lambda)$ to correct on average for these artifacts. This procedure is effective only with the fluxes from the DR13 pipeline since the non-linearities present in the DR12 pipeline make the correction flux-dependent. We show in Section 5 that after this global correction, remaining calibration artifacts (due to their time-dependence) are sufficiently small to have a negligible effect on the measurement of the correlation function.

The spectra of all quasar targets were visually inspected (Pâris et al., 2012, 2014, 2016) to correct for misidentifications, to flag broad absorption lines (BALs), and to establish the definitive quasar redshift. Damped Ly$\alpha$ troughs (DLAs) were visually flagged, but also identified and characterized automatically (Noterdaeme et al., 2012). The visual inspection of DR12 confirmed 297,301 quasars, of which 181,719 are in the redshift range appropriate for this study, 2.1 $\leq$ $z_{q}$ $\leq$ 3.5. We discarded quasars with visually identified BALs (to avoid the necessity of modeling their profiles in the forest) leaving 160,868 quasars. A further cut requiring a minimum number of unmasked forest pixels (50 “analysis pixels”; see below) yielded a sample of 157,922 quasars. Finally, 139 spectra failed the continuum-fitting procedure (Sect. 4.1), leaving 157,783 spectra compared to 137,562 in the Delubac et al. (2015) investigation.

For the measurement of the flux transmission, we use the rest-frame wavelength interval

[TABLE]

As illustrated in Fig. 2, this range is bracketed by the emission lines $\lambda_{\rm Ly\beta}=102.572$ nm and $\lambda_{\rm Ly\alpha}=121.567$ nm. This region was chosen as the maximum range that avoids the large pixel variances on the slopes of the two lines due to quasar-to-quasar diversity of line-emission strength and profile. The absorber redshift, $z=\lambda/\lambda_{\rm Ly\alpha}-1$, is required to lie in the range $1.96<z<3.44$. The lower limit is set by the requirement that the observed wavelength be greater than 360 nm, below which the system throughput is less than 10% of its peak value. The upper limit is produced by the maximum quasar redshift of 3.5, beyond which the BOSS surface density of quasars is not sufficient to be useful for this study. The weighted distribution of redshifts of absorber pairs near the BAO peak position is shown in Fig. 3. The distribution has a mean of $\overline{z}=2.33$.

For the determination of the correlation function, we use analysis pixels that are the inverse-variance-weighted flux average over three adjacent pipeline pixels. Throughout the rest of this paper, “pixel” refers to analysis pixels unless otherwise stated. The width of these pixels is $207\leavevmode\nobreak\ {\rm km\leavevmode\nobreak\ s^{-1}}$ corresponding at $z\sim 2.33$ to an observed-wavelength width $\sim 0.28\leavevmode\nobreak\ {\rm nm}$ and a comoving radial distance of $\sim 2.0\leavevmode\nobreak\ h^{-1}{\rm Mpc}$. The total sample of 157,783 spectra thus provides $\sim 3\times 10^{7}$ measurements of Ly$\alpha$ absorption over an effective volume of $\sim 50\leavevmode\nobreak\ {\rm Gpc}^{3}$.

3 Mock quasar spectra

In order to test the analysis procedure and investigate statistical and possible systematic errors, we created 100 sets of mock spectra that reproduce the essential physical and instrumental characteristics of the BOSS spectra. The basic method for the production of the mocks with Ly$\alpha$ absorption is described in Font-Ribera et al. (2012). Except for the implementation of absorption due to metals and high column-density systems (HCDs, i.e., damped Ly$\alpha$-systems and Lyman-limit systems), the mocks used in this study are identical to the DR11 mocks (Bautista et al., 2015) that were used in Delubac et al. (2015). The DR11 mocks therefore cover $\sim 15\%$ less solid angle than the DR12 data.

For each set of spectra, the background quasars were assigned the angular positions and redshifts of the DR11 quasars. The unabsorbed spectra (continua) of the quasars were generated using the Principal Component Analysis eigenspectra of Suzuki et al. (2005). The amplitudes for each eigenspectrum were randomly drawn from Gaussian distributions with a dispersion equal to that of the corresponding eigenvalues in Table 1 of Suzuki (2006). The overall normalization was chosen by fitting the mock spectrum to the corresponding observed spectrum.

The Ly$\alpha$ absorption field was generated using the technique described in Font-Ribera et al. (2012) in which a Gaussian random field, $\delta_{G}$, is first defined at the positions of the $\sim 10^{7}$ observed forest pixels. Flux transmissions were defined by the non-linear transformation

[TABLE]

The two functions, a(z) and b(z) are chosen so that the resulting mean transmission and variance are near those observed in the data. The correlations of the field $\delta_{G}$ are chosen so that the correlations of $F(\delta_{G})$ follow the linear correlations of the “mock” cosmology of Table 1 modified by non-linear effects (McDonald, 2003). For pixels randomly distributed in space, this procedure would involve inverting a $10^{7}\times 10^{7}$ matrix. To reduce the problem to a manageable size, use was made of the fact that the forest pixels are nearly parallel, allowing a separate treatment of radial and transverse coordinates.

In the final step, the spectra were modified to include the effects of the BOSS spectrograph point spread function (PSF), readout noise, photon noise, and flux systematic errors.

For each of the 100 mock data sets, four types of spectra were produced and analyzed. The first type consists simply of Ly$\alpha$ flux-transmission fraction, $F_{\rm Ly\alpha}(\lambda)$, modified for the wavelength resolution but without multiplication by a quasar continuum spectrum, $C(\lambda)$. Analysis of this mock type allowed us to study the recovery of the BAO peak position under the most favorable conditions. With the introduction of the quasar continuum, the second type consists of more realistic spectra, $F_{\rm Ly\alpha}(\lambda)C(\lambda)$. Analysis of this type tests our ability to fit the quasar continuum and to model the resulting distortion of the correlation function.

The final two types of spectra add to $F_{\rm Ly\alpha}(\lambda)$ absorption that is due to HCDs and to metals. Following Font-Ribera & Miralda-Escudé (2012), HCDs are placed at randomly chosen pixels where the optical depth was above a chosen threshold. The neutral-hydrogen column densities were drawn randomly with $N_{\rm HI}>10^{17.2}\leavevmode\nobreak\ {\rm cm^{-2}}$, assuming an intrinsic power-law distribution corrected for self-shielding and normalized to match the observations of Prochaska et al. (2005). Because the HCDs are placed in redshift space, the resulting HCDs have correlations that trace the underlying density field but with a redshift-space distortion parameter that is not necessarily equal to that of physical HCDs.

Absorption due to metal transitions was also simulated by adding absorption in proportion to the Ly$\alpha$ absorption at the same redshift. The important transitions can be seen in the “one dimensional” correlation function, $\xi_{\rm 1d}(\lambda_{1}/\lambda_{2})$ for fluxes at two wavelengths within the same forest, as shown in Fig. 5. In this figure, the peaks correspond to absorption by two different transitions by material at the same physical position but a different wavelength. The most prominent features are due to pairs comprising Ly$\alpha$ and one metal transition, but there are also features due to metal-metal pairs. More detailed information can be obtained by stacking spectra around strong absorbers in the forest (Pieri et al., 2010, 2014). Table 2 lists the transitions included in the mock spectra.

We have tested two procedures for adding metal absorption to the primary Ly$\alpha$ absorption. The methods are defined by the assumed relation between Ly$\alpha$ absorption and metal absorption at the same physical position. For each wavelength, $\lambda$, with Ly$\alpha$ absorption $A_{\rm Ly\alpha}(\lambda)=1-F_{\rm Ly\alpha}(\lambda)$, an appropriate metal absorption $A_{m}(\lambda+\Delta\lambda_{m})$ for each transition, $m$, is chosen. A simple procedure would be to choose $A_{m}\propto A_{\rm Ly\alpha}$ with a proportionality constant chosen to reproduce the Ly$\alpha$-metal features in the observed $\xi_{\rm 1d}$. However, this procedure does not produce significant features corresponding to metal-metal correlations. To do this it is necessary to add, for a randomly chosen small fraction of wavelengths, a much larger metal absorption.

The first procedure (hereafter Met1) aims to match the observed line-of-sight correlation function $\xi_{\rm 1d}$, in particular the amplitudes of the peaks associated to Ly$\alpha$-metal and metal-metal correlations. We first transform the pure mock Ly$\alpha$ flux into “optical depth” $\tau_{\rm Ly\alpha}$ and define the metal absorption of transition $m$ as $\tau_{m}=a_{m}\tau_{\rm Ly\alpha}$. A quadratic term in $\tau_{\rm Ly\alpha}$ is added for a fraction of the strong Ly$\alpha$ absorbers to simulate strong metal absorption. The parameters $a_{m}$, the quadratic terms and rates of strong absorption are set to match the observed amplitudes in the data $\xi_{\rm 1d}$.

The second procedure (hereafter Met2) aims to match the observed stack of high signal-to-noise ratio Ly$\alpha$ absorbers, while maintaining consistency with the $\xi_{\rm 1d}$ in the resulting mock $\xi_{\rm 1d}$. Following Pieri et al. (2014), stacks were produced using the DR12 Ly$\alpha$ forest sample. As a function of the Ly$\alpha$ transmission, we measured amplitudes for the metal absorption features. The obtained “flux decrements” caused by metals were treated as target average flux decrements. Metals are implemented in Met2 mocks as a mix of weak and strong metal absorption, generating a mock $\xi_{\rm 1d}$ consistent with the observed $\xi_{\rm 1d}$ and consistent with the absorption frame measurements of Pieri et al. (2014).

The metal absorption added by these methods is due to the presence of metals at the same physical position as HI. Because our mocks provide the HI density only for redshifts $z>1.9$, we cannot generate absorption at lower redshifts that nevertheless appears in the Ly$\alpha$ forest because of transition wavelengths much greater than that of Ly$\alpha$. An important example is CIV doublet ($\overline{\lambda}=154.9$ nm) where the absorption at an observed wavelength of 400 nm is due to material at $z\sim 1.6$. Fortunately this absorption has little effect on the correlation function, as we will see in Sect. 6.

The statistical properties of metal absorption in the mocks are determined by the underlying density field. However, the analysis procedure interprets absorption at a given wavelength as absorption due to the Ly$\alpha$ transition. Because of this behavior, the metal contribution to the measured correlation function is shifted and deformed in $(r_{\perp},r_{\parallel})$ space. In particular, the large correlation due to HI and metals at the same physical position is seen at $(r_{\perp}\sim 0,r_{\parallel})$ with $r_{\parallel}=(1+z)D_{H}(z)\Delta\lambda/\lambda$ where $\Delta\lambda/\lambda$ is the relative wavelength separation of the metal feature with respect to Ly$\alpha$. This leads to a large correlation at this $r_{\parallel}$ (and $r_{\perp}\sim 0$), so if the amplitude is significant, it can be measured. Table 2 lists the apparent $r_{\parallel}$ corresponding to vanishing Ly$\alpha$-metal separation.

4 Measurement of flux-transmission field and its correlation function

In this Section we describe the measurement of the correlation function of the transmitted flux fraction:

[TABLE]

Here, $f_{q}(\lambda)$ is the observed flux density for quasar $q$ at observed wavelength $\lambda$, $C_{q}(\lambda)$ is the unabsorbed flux density (the so-called “continuum”) and $\overline{F}(z)$ is the mean transmitted fraction at the absorber redshift, $z(\lambda)=\lambda/\lambda_{\rm Ly\alpha}-1$. Measurement of the flux-transmission field $\delta_{q}(\lambda)$ requires estimates of the product $C_{q}(\lambda)\overline{F}(z)$ for each quasar. An example is shown in Fig. 2. The estimation procedures, described in this Section, differ slightly from those of our previous studies (Busca et al., 2013; Delubac et al., 2015). One important modification is that we now calculate the distortion of the correlation function due to continuum-fitting (Sect. 4.3).

4.1 Estimation of $\delta_{q}(\lambda)$

As in our previous investigations, we assume the quasar continuum, $C_{q}(\lambda)$, is the product of a universal function of the rest-frame wavelength, $\lambda_{\rm RF}=\lambda/(1+z_{q})$, and a linear function of $\log\lambda$, included to account for quasar spectral diversity:

[TABLE]

with $C(\lambda_{\rm RF})$ being normalized so that its integral over the forest is equal to unity. The $(a_{q},b_{q})$ and $C(\lambda_{\rm RF})$ are determined by maximizing the likelihood function given by

[TABLE]

Here $P(f_{q}(\lambda)\,|\,C_{q}(\lambda))$ is the probability to observe a flux $f_{q}(\lambda)$ for a given continuum found by convolving the intrinsic probability, $D(F=f_{q}(\lambda)/C_{q}(\lambda),z)$, with the observational resolution assumed to be Gaussian:

[TABLE]

Here, $\sigma_{q}(\lambda)^{2}$ is the variance due to readout noise and photon statistics.

A continuum-determination method is defined by the assumed form of $D(F,z)$. For this work, it is taken to be the log-normal model of absorption used to generate the mock data, corresponding to “method 2” of our previous studies (Delubac et al., 2015). As a check, we also use “method 1”, which is equivalent to using a narrow Gaussian for $D(F,z)$, thereby producing only the product $C_{q}(\lambda)\overline{F}(\lambda)$ for each forest.

In practice, we maximize the likelihood iteratively by assuming a $C(\lambda_{\rm RF})$ to determine the $(a_{q},b_{q})$. The mean absorption $\overline{F}(\lambda)$ is then calculated by an appropriately weighted average of $f_{q}(\lambda)/C_{q}(\lambda)$ (for fixed $\lambda$) after which $C(\lambda_{\rm RF})$ is recalculated as an average of $f_{q}(\lambda)/\overline{F}(z)$ (for fixed $\lambda_{\rm RF}$). The procedure is stopped after ten iterations, at which point a stability of $\sim 10^{-4}$ is reached for $C(\lambda_{\rm RF})$ and $\overline{F}(\lambda)$. Figure 2 shows a spectrum with its $C_{q}(\lambda)$ and $C_{q}\overline{F}$.

The function $\overline{F}(z)$ to be used in Eq. 2 to calculate $\delta_{q}(\lambda)$ is the mean of $f_{q}(\lambda)/C_{q}(\lambda)$ and can thus differ slightly from $\overline{D(F,z)}$ used in the model to estimate the $C_{q}(\lambda)$. The use of the mean of $f_{q}(\lambda)/C_{q}(\lambda)$ ensures that the mean of $\delta$ is zero at each redshift (or wavelength).

Figure 6 displays the calculated $\overline{F}(z)$, both for the DR13 pipeline used here and the DR12 pipeline used previously. The use of the new pipeline removes most of the artifacts in the old analysis. We emphasize, however, that the derived $\overline{F}(z)$ is dependent on the assumed form of $D(F,z)$ and should therefore not be considered as a measurement of the mean absorption. For example, the flattening of $\overline{F}(z)$ for $z<2.1$ ($\lambda<377\leavevmode\nobreak\ {\rm nm}$) suggests that we have slightly underestimated $\overline{F}(z)$ at these redshifts. Since, by construction, the mean $\delta$ vanishes at each redshift, this implies that our procedure makes a compensating overestimate of the $C_{q}(\lambda)$. Since it is the product $\overline{F}(z)C_{q}(\lambda)$ that determines $\delta_{q}(\lambda)$, the measured correlation function is therefore not strongly affected.

Those forests with identified DLAs are given a special treatment. All pixels where the absorption due to the DLA is higher than 20% are not used. The absorption in the wings is corrected using a Voigt profile following the procedure of Noterdaeme et al. (2012).

We denote as $\tilde{\delta}_{q}(\lambda)$ the estimate of $\delta_{q}(\lambda)$ using the relation (2). Because forest data is used to fit for $C_{q}(\lambda)$ and $\overline{F}(z)$ the measured $\tilde{\delta}_{q}(\lambda)$ is not equal to the original $\delta_{q}(\lambda)$. We can identify two effects. First, the use of $\overline{F}(\lambda)$ is equivalent to the transformation

[TABLE]

where the over-bar refers to the average over forests at fixed $\lambda$. Second, the fitting of $(a_{q},b_{q})$ with the data biases toward zero the mean $\tilde{\delta}_{q}(\lambda)$ and mean $(\log\lambda-\overline{\log\lambda})\tilde{\delta}_{q}(\lambda)$. To simplify this effect and facilitate its correction we explicitly subtract the mean and first moment for each forest

[TABLE]

where the over-bars refer to averages within a given forest using the weights defined in the next Section. The accompanying distortion of the correlation function is non-negligible, as we will demonstrate in Sect. 6.2 with the mock spectra.

The transformation (7) has two interesting effects on the measured flux-transmission field. Most importantly, it makes it simple to calculate the distortion of the correlation function (Sect. 4.3) and thus simplifies the relation between the underlying physical model and the measured correlation function. Second, it nearly eliminates the difference between the correlations functions calculated with the two continuum fitting methods used in Delubac et al. (2015) with the r.m.s. difference in the two $\xi(r_{\perp},r_{\parallel})$ being 0.056 of the r.m.s. uncertainty per $(r_{\perp},r_{\parallel})$ bin.

4.2 Estimation of the correlation functions

For the estimator of the flux auto-correlation function, we adopt a simple weighted sum of products of the $\hat{\delta}$:

[TABLE]

where the $w_{i}$ are weights (see below) and each $i$ (or $j$) indexes a measurement on a quasar $q_{i}$ at wavelength $\lambda_{i}$. The sum over $(i,j)$ is understood to run over all pairs of pixels within a bin $A$ in the space of pixel separations, $(A)\rightarrow(r_{\perp},r_{\parallel})$. We exclude pairs of pixels from the same quasar to avoid the correlated errors in $\hat{\delta}_{i}$ and $\hat{\delta}_{j}$ arising from the estimate of $C_{q}(\lambda)$ for the spectrum of the quasar. The bins $A$ are defined by a range of width $4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ of the components perpendicular and parallel to the line of sight, $r_{\perp}$ and $r_{\parallel}$. We use 50 bins in each component, spanning the range from [math] to $200\leavevmode\nobreak\ h^{-1}{\rm Mpc}$; the total number of bins used for evaluating the correlation function is therefore 2500. Separations in observational pixel coordinates (RA,Dec,$\lambda$) are transformed to $(r_{\perp},r_{\parallel})$ in units of $h^{-1}{\rm Mpc}$ by assuming that absorption is due to the Ly$\alpha$ transition and using the cosmological parameters from Table 1 (Planck cosmology for the data and the mock cosmology for the mocks).

As described in Delubac et al. (2015), the weights, $w_{i}$, are chosen so as to account for both Poisson noise in the flux measurement and for the intrinsic fluctuations in $\hat{\delta}_{i}$ due to cosmological large-scale structure. The weights are set to zero for pixels flagged by the pipeline as having problems due, for example to sky emission lines or cosmic rays. To reduce the pipeline systematics discussed in Sect. 5, we also do not use pairs of pixels that have nearly the same wavelength ($r_{\parallel}<4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$) and that were taken on the same exposures.

4.3 The distortion matrix

The transformations (6) and (7) mix pixels so that the correlation $\langle\hat{\delta}_{q}(\lambda)\hat{\delta}_{q^{\prime}}(\lambda^{\prime})\rangle$ is equal to the original $\langle\delta_{q}(\lambda)\delta_{q^{\prime}}(\lambda^{\prime})\rangle$ plus a linear combination of the correlations of other pixel-pairs $\langle\delta_{q^{\prime\prime}}(\lambda^{\prime\prime})\delta_{q^{\prime\prime\prime}}(\lambda^{\prime\prime\prime})\rangle$. This approach means that the measured correlation function $\hat{\xi}$ is a “distorted” version of the true correlation function $\xi$. Since the transformations (6) and (7) are linear, the relation between measured and true correlation functions is given by a distortion matrix $D_{AB}$:

[TABLE]

where $A$ and $B$ refer to bins in pixel separation space. Writing $\hat{\delta}_{i}=\sum_{i^{\prime}}\eta_{ii^{\prime}}\delta_{i^{\prime}}$ produces

[TABLE]

where $w_{A}=\sum_{i,j\in A}w_{i}w_{j}$. We ignore the small effect of transformation (6), in which case $\eta_{ii^{\prime}}=0$ unless the pixels $i$ and $i^{\prime}$ are in the same forest:

[TABLE]

where the two sums over $k$ include only pixels in the same forest as that of $i$ and $i^{\prime}$ and the over-bars refer to averages in that forest. The matrix $D_{AB}$ thus depends only on the geometry and weights of the survey. We will see its effect on the mock correlation function in Sect. 6.2 (Fig. 11).

Previous analyses (Busca et al., 2013; Slosar et al., 2013; Delubac et al., 2015; Blomqvist et al., 2015) have dealt with the distortions introduced by continuum-fitting in different ways. Busca et al. (2013) and Delubac et al. (2015) model it as an additive power law in $r$ and $\mu$ with 12 free parameters. Blomqvist et al. (2015) assumed that continuum-fitting reduces the observed amplitude of long-wavelength modes in the direction parallel to the line of sight. They then model it as a multiplicative function of $k_{\parallel}$ that tends to zero at large scales (with $k_{\parallel}\lesssim 2\pi/L$, where $L\sim 380\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ is the typical length of a forest) and to one at small scales ($k_{\parallel}\gg 2\pi/L$). They tune the shape of the function using simulated data to ultimately reduce the number of free parameters to one.

In our approach we do not introduce free parameters to account for the effects of continuum-fitting. Instead, we follow the assumption, first proposed by Slosar et al. (2013), that at each line of sight the continuum-fit delta field differs from the true delta field by a linear function in $\Lambda$. Slosar et al. (2013) then de-weight these “linear modes” in their covariance matrix. Alternatively, we use the transformation from $\delta$ to $\hat{\delta}$ to remove the “linear modes” from the model via the distortion matrix.

4.4 The covariance matrix

The covariance matrix associated with $\hat{\xi}$ is:

[TABLE]

where $W_{A}=\sum_{ij\in A}w_{i}w_{j}$. Following Delubac et al. (2015), we evaluate $C_{AB}$ by dividing the BOSS footprint into sub-samples and measuring $\xi^{s}_{A}$ and $\xi^{s}_{B}$ in each sub-sample $s$. Neglecting the small correlations between sub-samples, and replacing the four-point function by the measured product of correlations in subsamples, the covariance is given by

[TABLE]

where the $W_{A}^{s}$ are the sums of weights in the sub-sample $s$. As in Delubac et al. (2015), the SDSS plates define the sub-samples.

As a check of the subsampling method, the sum (12) can also be estimated by neglecting inter-forest correlations, in which case the four-point function vanishes unless the four pixels are drawn from just two spectra:

[TABLE]

The sum can then be estimated from a random sample of forest pairs. Because neighboring forests are nearly parallel, the sum necessarily gives $C_{AB}=0$ unless $r_{\perp}^{A}\sim r_{\perp}^{B}$.

Finally, for the analysis of the mock data, the covariance matrix can also by calculated from the mock-to-mock variations of $\xi$:

[TABLE]

where the over-bar refers to averages over the set of 100 mocks. Given the approximations used in the calculation of $C_{AB}$ via eqns. 13 or 14, it is of great importance that the mock-to-mock variations confirm the accuracy of the other two methods when applied to mock data. This comparison is shown in Fig. 7.

The $2500\times 2500$ element matrix, $C_{AB}$, has a relatively simple structure. By far the most important elements are the diagonal elements which are, to good approximation, inversely proportional to the number of pixel pairs used in the calculation of the correlation function; the number of pairs is roughly proportional to $r_{\perp}$:

[TABLE]

The variance is about twice as large as what one would calculate assuming all (analysis) pixels used to calculate $\xi(r_{\perp},r_{\parallel})$ are independent. This decrease in the effective number of pixels is due to the physical and instrumental correlations between neighboring pixels in a given forest (eqn. 14).

The off-diagonal elements of the covariance matrix also have a simple structure. As previously noted, the covariance is mostly due to the two-forest part of the four-point function which, because neighboring forests are nearly parallel, only contribute to the covariance matrix elements with $r_{\perp}=r_{\perp}^{\prime}$. This behavior is illustrated in Fig. 7, which displays the mean values of the correlation matrix elements as a function of $r_{\parallel}-r_{\parallel}^{\prime}$ for the smallest values of $r_{\perp}-r_{\perp}^{\prime}$.

The statistical precision of the sub-sampling calculation is $\sim 0.02$ for individual elements of the correlation matrix. Figure 7 reveals that only correlations with $\Delta r_{\perp}\sim 0$ and $\Delta r_{\parallel}<20\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ are greater than the statistical precision and therefore sufficiently large for individual matrix elements to be measured accurately by sub-sampling. As in Delubac et al. (2015), we therefore average the sub-sampling correlation matrix over $(r_{\perp}^{\prime},r_{\parallel}^{\prime})$ and use the resulting covariance matrix that is a function of $(r_{\perp}-r_{\perp}^{\prime},r_{\parallel}-r_{\parallel}^{\prime})$.

5 Correlations introduced by the optics and data pipeline

Spurious correlations in $\delta_{q}(\lambda)$ are introduced by the telescope and spectrometer optics, and by the pipeline reductions of the data. These correlations are superimposed on the physical $\xi(r_{\perp},r_{\parallel})$ that we wish to measure; in this Section we estimate the various contributions.

5.1 Optical cross-talk

At the optical level, correlations are introduced by signals from one object “scattering” into the spectra of other objects, that is through optical cross-talk. There is negligible cross-talk introduced in the telescope focal plane by photons from one object entering into the fiber of another object. However, there can be measurable cross-talk between neighboring fibers downstream of the wavelength dispersion where they are focused on the CCD. This contamination arises through imperfect modeling of the point-spread function when transforming the two-dimensional CCD image into a series of one-dimensional spectra, $f(\lambda)$. We have measured the cross-talk as a function of fiber separation by fitting the signal in sky fibers to a sky model and a weighted sum of the spectra of the objects (quasars or galaxies) in neighboring fibers. Consistent with the results of Croft et al. (2016), we find the cross-talk for the neighboring fibers is $\sim 0.2\%$ and $\sim 0.05\%$ for fibers separated by two or three rows respectively, and is consistent with zero for larger separations.

The cross-talk directly introduces correlations between pixels at the same wavelength but these pixel pairs are, at any rate, not used in the analysis. The cross-talk introduced between pixels $\lambda_{1}$ and $\lambda_{2}$ of two quasars is proportional to the product of the cross-talk amplitude and to $\xi_{\rm 1d}(\lambda_{1}/\lambda_{2})$. We have verified that these correlations are insignificant compared to the measured $\xi$. This fact is in part due to the fiber-assignment strategy which avoided placing two quasar candidates on neighboring fibers.

5.2 Pipeline-induced correlations

The pipeline treatment of the spectrometer data transforms flat-field corrected CCD counts, $f_{q}^{\rm m}(\lambda)$, to fluxes, $f_{q}(\lambda)$. This process requires subtracting a sky contribution, $B_{q}(\lambda)$, and multiplying by a calibration vector, $A_{q}(\lambda)$, that corrects for the wavelength-dependent throughput of the system:

[TABLE]

Both $A_{q}(\lambda)$ and $B_{q}(\lambda)$ are determined from spectra taken simultaneously with the data: spectra of spectro-photometric standard stars for $A_{q}$, and spectra of “sky fibers” (fibers pointing to empty sky regions) for $B_{q}$. The spectra from the 1000 fibers of a given exposure are treated independently for the two 500-fiber spectrographs because of possible differential variations of instrumental properties such as throughput, PSF, and scattered light. (See Smee et al. (2013) for the relevant details of the spectrometer construction.) Hereafter, the collection of fibers of a given plate assigned to one of the two spectrographs will be referred to as a “half-plate.” Inaccuracies in the determinations of $A_{q}$ and $B_{q}$ will lead directly to correlated inaccuracies in the $f_{q}(\lambda)$ of a given half-plate.

Errors in the sky subtraction, $B_{q}(\lambda)$, are simple to model because they are mostly due to well-understood Poisson fluctuations of photo-electron counts in the $\sim 50$ sky fibers per half-plate. The pipeline interpolates these measurements in the focal plane to produce the sky background to be subtracted from a given quasar. The limited number of sky-fibers results in correlations between $f_{q}(\lambda)$ and $f_{q^{\prime}}(\lambda^{\prime})$ for $q\neq q^{\prime}$ and $\lambda=\lambda^{\prime}$ since the sky subtractions where calculated using the same noisy sky-fibers. These correlations for $\lambda=\lambda^{\prime}$ generate small correlations for $\lambda\neq\lambda^{\prime}$ through continuum fitting, as discussed in Sect. 4.2. Because the correlations are primarily for $\lambda=\lambda^{\prime}$ they contribute a spurious $\xi(r_{\perp},r_{\parallel}=0)$ for pairs of pixels on the same half-plate.

Errors on the calibration vector $A_{q}$ are also Poissonian in that they are due to the fluctuations in the imperfect modeling of $\sim 10$ calibration stars per half-plate. The mean over all spectra of the imperfections are visible in the mean “transmission” of the unabsorbed parts of quasar spectra shown in Fig. 4. As mentioned in Section 2, the mean imperfections are removed by dividing forest spectra by $R(\lambda)$ from Fig. 4. However, fluctuations from exposure to exposure remain because of the spectral variations of the small number of calibration stars used. The statistical properties of these variations have been estimated by comparing the $A_{q}$ measured with random sub-samples of stars and, independently, by comparing the $A_{q}$ on different half-plates for the same exposure. Unlike the correlations due to sky-subtraction, which are confined to $\lambda=\lambda^{\prime}$, the errors in the calibration vectors are strongly correlated over the characteristic wavelength range of stellar spectral features, $\Delta\lambda\sim 1$ nm.

We constructed a model for pipeline-induced correlations that uses the measured statistical properties of the $A_{q}(\lambda)$ and $B_{q}(\lambda)$. The contribution of the sky noise was estimated using a fit of random realizations of the signal in the sky fibers (using a degree 2 polynomial fit as a function of the fiber number). The realizations were based on the statistical uncertainty estimation from the pipeline (which we have found to be accurate at the 5% level, see Fig. 19, bottom panel). A normalization factor was applied to account for the average difference of calibration between the sky and target fibers (a $\sim 10\%$ correction). Performing a fit allowed us to capture accurately the resulting correlations as a function of fiber separation.

The contribution of the calibration was estimated using the fact that because the flux calibration of the two spectrographs is determined independently, the difference of those two calibration solutions for each exposure provides us with an estimate of their statistical fluctuation (up to a normalization factor $\sqrt{2}$) while retaining the correlation of those fluctuations as a function of wavelength. Those differences for all DR12 plates was computed. Their average per observation run was subtracted in order to account for the actual difference of throughput of the two spectrographs and their evolution during the course of the survey. This data set provided us with a library of calibration uncertainties that was used to calculate calibration-induced correlations between spectra.

This model was used to calculate the expected correlation between pixels in the Ly$\alpha$ forest of one quasar and pixels in the CIV forest ($142<\lambda_{\rm RF}<152$ nm) of another quasar. The physical correlation of the two forests is due mostly to the auto-correlation of the weak CIV absorption in the two forests, but any pipeline-induced correlations should be present at full strength. Since this cross-correlation is designed to isolate pipeline-induced effects, we denote this correlation function as $\xi^{\rm pl}(\lambda_{2}-\lambda_{1},\theta)$ where $(\lambda_{1},\lambda_{2})$ are the wavelengths in the two forests and $\theta$ is the angular separation. We also distinguish between $\xi^{\rm pl}$ and correlations between flux pairs on the same half-plate, $\xi^{\rm pl}({hp=hp^{\prime}})$, and correlations between flux pairs on different half-plates, $\xi^{\rm pl}({hp\neq hp^{\prime}})$. While these correlations are naturally given as a function of wavelength and angular separations, for convenience, the $\xi^{\rm pl}$ will be given as functions of pseudo-separations $(\tilde{r}_{\perp},\tilde{r}_{\parallel})$ calculated using the Ly$\alpha$ rest-frame wavelength to define redshifts in both spectra. With this approach, $\tilde{r}_{\parallel}=0$ corresponds to absorption at the same observed wavelength.

This correlation for same-half-plate pairs (measured using the techniques of Sect. 4.2) is shown with the red points in Fig. 8 for the first $\tilde{r}_{\parallel}$ bin ($\tilde{r}_{\parallel}<4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$). Superimposed on the data is the prediction of the model of the pipeline-induced correlations. For the first $\tilde{r}_{\parallel}$ bin, the correlation is dominated by those induced by to the sky-subtraction model. There is good agreement between this simple model and the observed correlations.

The pipeline-induced correlations that we have considered do not contribute to correlations between pixels observed on different half-plates. The different-half-plates correlation functions for the $\tilde{r}_{\parallel}<4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ bin is shown by the blue points in Fig. 8. The correlations for different half-plates (blue points) are clearly far less than those for same half-plates (red points). For $\tilde{r}_{\parallel}<4\leavevmode\nobreak\ h^{-1}{\rm Mpc},\tilde{r}_{\perp}<\leavevmode\nobreak\ 60\leavevmode\nobreak\ h^{-1}{\rm Mpc}$, $\chi^{2}=28.6$ ($N_{\rm dof}=15$) for the no-correlation hypothesis.

For the $\tilde{r}_{\parallel}>4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ bins, the model predicts much smaller correlations. In particular, the sky-subtraction model noise induces non-zero correlations only because of the continuum fit which distorts the original correlation function as described in Section 4.2. This behavior is illustrated in Fig. 9 showing $\xi^{\rm pl}$ in four ranges of $\mu$. Because of the low level of absorption fluctuations redward of Ly$\alpha$ emission, the variance of the points is a factor $\sim 7$ smaller than the corresponding variance for the forest-forest correlation function.

The model for pipeline-induced correlations assumes that the sky-subtraction errors are entirely due to Poisson statistics of photo-electron counts, neglecting possible wavelength- and position-dependent systematic mis-estimates of the sky flux. While the changes applied to the pipeline (Appendix A) considerably reduced the systematic sky residuals, significant residuals remain on bright sky lines due to an imperfect PSF model, PSF variations and displacement of spectral traces (due to changes of temperature and variations of the gravity load of the spectrographs which are at the Cassegrain focus of the telescope). We tested the effect of those sky residuals by computing the correlation function of fake Ly$\alpha$ forests consisting of the residual signal in the nearest sky fiber in the CCD divided by the quasar continuum model (we made sure to use only once each sky fiber in this process to avoid introducing an artificial correlation). The measured correlation signal could be entirely explained by the Poissonian sky model noise described above. We hence conclude that the systematic sky residuals induce a negligible contamination to the Ly$\alpha$ correlation function.

This analysis presented above evaluates the contamination due to any additive signal to the Ly$\alpha$ forest that leads to a correlated signal in the CCD. It includes the systematic sky residuals, their fluctuations from plate to plate, but also the potential effect of scattered light in the spectrograph.

We have used our model of pipeline-induced correlations to calculate their effect on the determination of the BAO peak position. Since the model contains no scale near the BAO scale, it is not surprising that it predicts no measurable influence, as reported in Table 7 of Sect. 7. Furthermore, to facilitate the comparison of the measured correlation function with the physical model developed in the next section, we do not use pairs of pixels on the same exposure that would contribute to the $r_{\perp}<4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ bins of $\xi(r_{\perp},r_{\parallel})$.

6 Fits for the BAO peak position

To determine the position of the BAO peak in the transverse and radial directions, we fit the measured $\xi(r_{\perp},r_{\parallel})$ to functions that describe the underlying large-scale-structure correlations. These correlations are primarily due to Ly$\alpha$ absorption in the intergalactic medium (IGM), but we also include absorption by metals in the IGM and neutral hydrogen in high-column-density systems (HCDs). The physical correlations are corrected for distortions that are introduced by the procedure for determining the quasar continuum (Sect. 4.3). The fitting routine we use gives results consistent with those found using the publicly available baofit222http://darkmatter.ps.uci.edu/baofit/ package (Kirkby et al., 2013; Blomqvist et al., 2015) modified to include the effects described in this Section. The parameters of the fits are described below and in Table 3. The best-fit model correlation function is shown in Fig. 10.

6.1 The model of the correlation function

The general form for the model correlation function $\xi_{A}$ of the $(r_{\perp},r_{\parallel})$ bin $A$ is a distorted sum of the cosmological, or “physical” correlation function, $\xi^{\rm ph}_{A}$, and a slowly varying function, $B_{A}$, used to test for systematics:

[TABLE]

Here, $D_{AA^{\prime}}$ is the distortion matrix (eqn. 10) that models the effects of continuum fitting.

The physical component of the model is dominated by the auto-correlation due to Ly$\alpha$ absorption. It is assumed to be a biased version of the total matter auto-correlation of the appropriate flat-$\Lambda$CDM model (Table 1) modified to free the position of the BAO peak (Kirkby et al., 2013):

[TABLE]

where the BAO peak-position parameters to be fit are

[TABLE]

and where the subscript “fid” refers to the fiducial cosmological model from Table 1 used to transform angle differences and redshift differences to $(r_{\perp},r_{\parallel})$. The nominal correlation function, $\xi_{\rm Ly\alpha}(r_{\perp},r_{\parallel},\alpha_{\perp}=\alpha_{\parallel}=1)$, is derived from its Fourier transform

[TABLE]

where $P_{\rm QL}$ is the (quasi) linear power spectrum decomposed into a smooth component and a peak component corrected for non-linear broadening of the BAO peak:

[TABLE]

The smooth component is derived from the CAMB-calculated linear power spectrum, $P_{\rm L}$, via the side-band technique (Kirkby et al., 2013) and $P_{\rm peak}=P_{\rm L}-P_{\rm smooth}$. The correction for non-linear broadening of the BAO peak is parameterized by $\Sigma^{2}=\mu_{k}^{2}\Sigma_{\parallel}^{2}+(1-\mu_{k}^{2})\Sigma_{\perp}^{2}$. The nominal values used are $\Sigma_{\parallel}=6.41\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ and $\Sigma_{\perp}=3.26\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ (Kirkby et al., 2013).

In equation 21, the bias, $b_{\rm Ly\alpha}$, is assumed to have a redshift dependence $b_{n}\propto(1+z)^{\gamma}$ with $\gamma=2.9$ (so that $P(k,z)b_{\rm Ly\alpha}(z)^{2}\propto(1+z)^{3.8}$), and $\beta_{\rm Ly\alpha}$ is assumed redshift independent. The function $F_{\rm NL}$ corrects for non-linear effects at large $k$ due to the isotropic enhancement of power due to non-linear growth, the isotropic suppression of power due to gas pressure, and the suppression of power due to line-of-sight non-linear peculiar velocity and thermal broadening. We use the form given by equation 21 and Table 1 of McDonald (2003). The forms proposed by Arinyo-i-Prats et al. (2015) produce nearly identical results for the range of $(r_{\perp},r_{\parallel})$ used in this study.

The last term in equation 21, $G(\@vec{k})=G_{\parallel}(k_{\parallel})G_{\perp}(k_{\perp})$, accounts for the binning in $(r_{\perp},r_{\parallel})$ space which effectively averages the correlation function over a bin. (The large width of these bins renders unnecessary a term that accounts for the spectrometer resolution.) If the distribution of observed $(r_{\perp},r_{\parallel})$ were uniform in a bin, $G(\@vec{k})$ would be the Fourier transform of the bin. The distribution is approximately uniform in the radial direction, which implies $G_{\parallel}={\rm sinc}^{2}(R_{\parallel}k_{\parallel}/2)$ where $R_{\parallel}=4\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ (the bin width). In the perpendicular direction, the distribution is approximately proportional to $r_{\perp}$ so $G_{\perp}=[(J_{1}(k_{\perp}r_{+})-J_{1}(k_{\perp}r_{-}))/k_{\perp}\Delta r]^{2}$ where $r_{+}$ and $r_{-}$ are the extrema of a given $r_{\perp}$ bin and $\Delta r=r_{+}-r_{-}$ is the bin width. The expression is bin-dependent and performing a Fourier transform on each bin would be too time-consuming. We therefore replace $G_{\perp}$ by ${\rm sinc^{2}}(k_{\perp}\Delta r)$ and evaluate the correlation function at the average $r_{\perp}$ of the pairs in a bin. This procedure produces a sufficiently accurate correlation function.

We allow for fluctuations of ionizing UV radiation (Pontzen, 2014; Gontcho A Gontcho et al., 2014) which lead to a scale-dependent bias of Ly$\alpha$ absorption given by eqn. 12 of Gontcho A Gontcho et al. (2014). The effect of the fluctuations is to increase $b_{\rm Ly\alpha}$ from its nominal value at small scale to a different value at large scale. The transition scale is determined by the UV photon mean-free-path which we set to a comoving value of $300\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ (Rudie et al., 2013). We then fit for one parameter, $b_{\Gamma}$ corresponding to the $b_{\Gamma}(b_{s}-b_{a})$ of Gontcho A Gontcho et al. (2014); it determines the change in $b_{\rm Ly\alpha}$ between large and small scales. A second bias, $b_{a}^{\prime}$, that determines the precise dependence of the bias on scale, is set to the nominal value of $-2/3$ used by Gontcho A Gontcho et al. (2014).

In our previous investigations, $\xi_{\rm Ly\alpha}$ was assumed to be a sufficiently accurate approximation for the correlation function. In this study, we include absorption by unidentified High-Column-Density Systems (HCDs) and by metals. The total physical correlation function is then the sum of the auto- and cross-correlations of the various components:

[TABLE]

Each component, $n$ or $m$, has its own bias and redshift-space distortion parameter.

The correlation function for HCDs and their cross-correlation with Ly$\alpha$ absorption is derived from a power spectrum as in eqn. 21 but with different bias parameters (Font-Ribera & Miralda-Escudé, 2012). While these absorbers are expected to trace the underlying density field their effect on the flux-transmission field depends on whether they are identified and given the special treatment described in Sect. 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. We model this by multiplying the bias for unidentified HCDs, $b_{\rm HCD}$, by a function of $k_{\parallel}$, $F_{\rm HCD}(k_{\parallel})$, that reflects the typical column densities of the unidentified HCDs. Following the studies with the mock spectra described in Sect. 3, we adopted the form $F_{\rm HCD}(k_{\parallel})={\rm sinc}(L_{\rm HCD}k_{\parallel})$, where $L_{\rm HCD}$ is a free parameter. The parameter $\beta_{\rm HCD}$ is poorly determined by the fit and we impose a Gaussian prior $\beta_{\rm HCD}=0.5\pm 0.2$

Because there is little absorption by metals the treatment of metal components is simplified without the separation into peak and smooth components. The fiducial correlation function is directly used to calculate in real space both metal-metal and metal-Ly$\alpha$ correlations. The parameters $(b,\beta)$ for each species (Table 3) are included by using a multipole expansion.

Absorption by metals is complicated by the fact that in the data analysis the Ly$\alpha$ transition is naturally assumed when transforming wavelength differences to position differences. Because of this procedure, for the cross-correlation between absorption by transition $n$ and transition $m$, the nominal $(r_{\perp},r_{\parallel})$ of the bin $A$ is not equal to the true separation coordinates $(r_{\perp}^{mn},r_{\parallel}^{mn})$. The model correlation, $\xi^{mn}_{A}$ for the bin $A$ is thus evaluated at the true, rather than nominal, separation.

Because amplitudes for SiII and SiIII are mostly determined by the excess correlation at $(r_{\perp}\sim 0,r_{\parallel}\neq 0)$, the $\beta$ for each metal is poorly determined. We therefore fix their value to $\beta=0.5$ corresponding to host halos with bias of 2, the value found for DLAs (Font-Ribera & Miralda-Escudé, 2012) which is also typical of star-forming galaxies.

The bias, $b_{\rm CIV}$, of foreground ($z\sim 1.65$) CIV absorption is poorly determined by the fit of the Ly$\alpha$-forest auto correlation. This is because the large wavelength difference between the Ly$\alpha$ and CIV transitions places the peak corresponding to vanishing metal-Ly$\alpha$ separation outside the forest. We have therefore measured $b_{\rm CIV}$ by correlating the flux in the CIV forest with quasars at $z\sim 1.65$. We compared these correlations with the analogous ones between the $\delta_{q}(\lambda)$ in the Ly$\alpha$ forest with $z\sim 2.3$ quasars. The resulting $\mu$-averaged correlation ratio is

[TABLE]

Because the quasar-CIV flux correlation is linear in $b_{\rm CIV}$ while the contribution to the flux auto-correlation is proportional to $b_{\rm CIV}^{2}$, this result implies a negligible ($<1\%$) contamination of the Ly$\alpha$-forest auto-correlation by CIV auto-correlation.

To transform the ratio (24) into a constraint on $b_{\rm CIV}$ that is needed for the fits, we equate it to the appropriate function of bias parameters:

[TABLE]

where $\xi_{\rm L}$ is the linear correlation function at $10\leavevmode\nobreak\ h^{-1}{\rm Mpc}$. To determine $b_{\rm CIV}(1.65)$ we use the values $(\beta_{\rm Ly\alpha},\beta_{\rm q},\beta_{\rm CIV})=(1.4,0.27,0.5)$, a growth factor $\xi_{\rm L}(1.65)/\xi_{\rm L}(2.3)=1.5$, and the redshift evolution of the quasar bias from Croom et al. (2005), $b_{\rm q}(1.65)]/b_{\rm q}(2.3)=0.70$. The resulting value of $b_{\rm CIV}(1.65)/b_{\rm Ly\alpha}(2.3)$ is multiplied by the measured value of $b_{\rm Ly\alpha}$ (Table LABEL:fitautotable) and then used as a Gaussian prior in the fits: $b_{\rm CIV}(2.3)=-0.0224\pm 0.0112$ (and $b_{\rm CIV}<0$), where we have conservatively doubled the uncertainty in (24).

To test for systematic errors, the model in equation (18) includes the broadband term $B_{A}$, a smoothly varying function of $(r_{\perp},r_{\parallel})$:

[TABLE]

where the $L_{j}$ are Legendre polynomials. In our previous studies, $B_{A}$ had a central role because we did not attempt to model distortions due to continuum fitting as we now do with $D_{AA^{\prime}}$. We include $B_{A}$ now as an option to study how the fit position of the BAO peak depends on the assumed form of the smooth component of the correlation function. By including this function, we can assure that the measurement of the BAO parameters is due only to the peak position, and is not significantly influenced by the assumed form of the smooth component.

6.2 Fits with the mock data sets

The 100 mock sets were used to test the fitting procedure and to search for possible biases in the fit BAO parameters. Studies were made on individual mock sets and for the stack of the 100 sets. The former test our ability to reliably extract $(\alpha_{\perp},\alpha_{\parallel})$ in data sets similar to the BOSS observations. The latter test the accuracy with which we model the analysis procedure including the effects of weighting, $(r_{\perp},r_{\parallel})$ binning, and continuum fitting.

Four types of mocks were produced and analyzed. The results are summarized in table LABEL:fitmocktable and Figs. 11 and 12. The four types and the analysis procedure are as follows:

Only Ly$\alpha$ absorption: in the fits the distortion matrix is set to the identity matrix and only the parameters $(\alpha_{\perp},\alpha_{\parallel},b_{\rm Ly\alpha},\beta_{\rm Ly\alpha})$ are fit. 2. 2.

A quasar continuum is added: the fits use the distortion matrix and the parameters $(\alpha_{\perp},\alpha_{\parallel},b_{\rm Ly\alpha},\beta_{\rm Ly\alpha})$ are fit. 3. 3.

Metal absorption is added: the additional parameters are the $(b_{i},\beta_{i})$ for each metal transition, $i$. 4. 4.

High column density systems are added: three additional parameters $(b_{\rm HCD},\beta_{\rm HCD},L_{\rm HCD})$ are fit.

Figure 11 shows the correlation function in four ranges of $\mu$. The black points are the stack of the measured correlation function for type 1 mocks (only Ly$\alpha$ absorption). The red points are the stacked measurements for type 2 mocks (quasar continuum added). The difference between the black and red points shows the effect of continuum-fitting. The heavy curves are the best fit model for the cosmology used to generate the mocks. The good agreement between the red curve and points indicates the accuracy of the modeling of the continuum-fitting distortion.

Fits with type 3 mocks (including metals) were used to search for possible influences of metal absorption on derived BAO parameters. Metal absorption is mostly confined to the small $r_{\perp}$ region and this is seen in the two top panels of Fig. 12. The most radial sector, $\mu>0.95$, displays the effects of metal absorption, most prominently the peak due to SiIIa,b absorption around $r\sim 60\leavevmode\nobreak\ h^{-1}{\rm Mpc}$. The peak is not visible in the next sector, $0.8<\mu<0.95$.

Two sets of fits for type 4 mocks (HCDs) were also performed. In one, HCDs with $N_{\rm HI}<10^{20}{\rm cm^{-2}}$ were not given the DLA treatment described in Section 4.1. This cut approximates the threshold for DLA detection either automatically or visually. In the second set, only HCDs with $N_{\rm HI}>10^{21}\leavevmode\nobreak\ {\rm cm^{-2}}$ were treated specially, testing our sensitivity to an underestimation of the efficiency to flag DLAs.

Table LABEL:fitmocktable lists the mean fit values for the 100 mocks and for the 4 types. Fig. 13 presents the corresponding distributions for type 4 mocks. For all types, the mean of the best-fit values of $(\alpha_{\perp},\alpha_{\parallel})$ is consistent with unity at the sub-percent level. Since the precision of our measurement with the data is of order three percent, this means that the mocks do not indicate any significant bias in the measurement of the BAO peak position.

6.3 Fits of the data with physical correlation function

The data were fit over the range $10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ under various hypotheses of increasing complexity concerning the form of $\xi^{\rm ph}(r_{\perp},r_{\parallel})$. The results are summarized in Table LABEL:fitautotable. The first fit assumes only Ly$\alpha$ absorption, the second includes absorption by metals, and the third absorption by unidentified HCDs. The fourth fit also assumes fluctuation in the flux of ionizing UV radiation. Inclusion of each effect yields a significantly better $\chi^{2}$; however, the values of $(\alpha_{\perp},\alpha_{\parallel})$ do not change significantly with each step. The data and fit correlation functions are presented for the four angular ranges in Fig. 14.

In what follows, we will focus our attention on the fourth fit that will be referred to as the “complete physical model”. The BAO peak-position parameters for this model are

[TABLE]

[TABLE]

where we have used the Planck values from Table 1 that were used in the (angle,redshift) to $(r_{\perp},r_{\parallel})$ transformation. The one, two and three standard deviation contours (corresponding to $\Delta\chi^{2}=2.3,6.18,11.83$) are shown as the red contours in Fig. 15. Our result is within one standard deviation of the prediction of flat-$\Lambda$CDM model of Table 1 that is in agreement with the Planck data.

In the next Section we will study the robustness of the determination of the BAO-peak position parameters by adding power-law broadband terms to the complete physical model. We can also estimate the significance of the BAO peak by modifying this model by setting $A_{\rm peak}=0$. This degrades the quality of the fit by $\Delta\chi^{2}=27.7$ corresponding to a BAO peak detection of $5.2\sigma$. Adding broadband terms reduces this to $\Delta\chi^{2}=21.8$ ($4.7\sigma$).

Other than $(\alpha_{\perp},\alpha_{\parallel})$, the parameters of the model listed in Table 3 have no direct cosmological implications. However, the fact that the best-fit values are physically reasonable reinforces our confidence in the model and, therefore, on the cosmological conclusions. Most importantly, the best-fit values of $(b_{\rm Ly\alpha},\beta_{\rm Ly\alpha})$ in table LABEL:fitautotable are both relatively stable and near the expected values. For example, the model of Arinyo-i-Prats et al. (2015) predicts $\beta_{\rm Ly\alpha}\sim 1.4$ and $|b_{\rm Ly\alpha}|(1+\beta_{\rm Ly\alpha})$ in the range 0.25 to 0.33. In Appendix B we discuss in more detail the comparison of our measurement of $(b_{\rm Ly\alpha},\beta_{\rm Ly\alpha})$ with theory and with previous measurements.

The best-fit values of most of the other parameters in Table 3 are near those found in the mock fits. For the metal biases, this only means that we have placed a reasonable amount of metal absorption in the mock spectra. Of more significance is the value of $L_{\rm HCD}=(24.3\pm 1.13)\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ which is, as expected, roughly the width of a DLA of $N_{\rm HI}=10^{20}\leavevmode\nobreak\ {\rm cm^{2}}$, corresponding to our estimated threshold for good efficiency in identifying DLAs. In the mocks, best-fit values of $L_{\rm HCD}$ when masking only those DLAs with $N_{\rm HI}>10^{20}\leavevmode\nobreak\ {\rm cm^{2}}$ cluster near $L_{\rm HCD}\sim 25$ in about half the fits. For the other half, the fits find $L_{\rm HCD}\sim 0$, corresponding to the addition of a second Ly$\alpha$ component but with a different $\beta$. This indicates that more than one way exists to model sufficiently well the correlation function with unidentified HCDs. The bias parameter found in the data, $b_{\rm HCD}=-0.0288\pm 0.0043$, is about a factor of 7 higher than that found in those mocks that find $L_{\rm HCD}\sim 25$. This indicates that there are unidentified HCDs in the data with $N_{\rm HI}>10^{20}{\rm cm^{2}}$ and or that the HCDs that we have placed in the mocks are less biased than in reality.

We have not generated mocks with fluctuations in the flux of ionizing radiation. However, we note that the bias parameter $b_{\Gamma}=0.1125\pm 0.0511$ is near the value of $0.13$ suggested by Gontcho A Gontcho et al. (2014).

6.4 Fits of the data with power-law broadbands added

In addition to statistical uncertainties, it is important to determine to what extent the derived values of $(\alpha_{\perp},\alpha_{\parallel})$ are influenced by the assumed form of $\xi^{\rm ph}$. The amplitude and form of the metal components are well understood and constrained by the low-$r_{\perp}$ bins. However, the HCD and UV parameters have a more diffuse effect on the correlation function and are not strongly constrained by independent data. For example, the UV background fluctuation model we have used includes only fluctuations that follow the underlying density field. Additional correlations (Gontcho A Gontcho et al., 2014) due to the discrete nature of sources of UV radiation would modify the absorption in neighboring forests since they are influenced by the same discrete sources. This “shot-noise” term would give an additive, isotropic and positive contribution to the Ly$\alpha$ auto-correlation function. Its shape and amplitude as a function of comoving separation are however poorly known as they depend on the quasar luminosity distribution and variability.

To study the effects of possible uncertainties in the form of the non-BAO component, $\xi^{\rm ph}_{\rm smooth}$, we have performed fits adding to the complete physical model broadband functions in the form of power laws (26) defined by $(i_{\rm min},i_{\rm max},j_{\rm max})$. With the addition of new components of $\xi^{\rm ph}_{\rm smooth}$, the connection between the BAO peak amplitude and the smooth components is lost so we also performed fits releasing the constraint $A_{\rm peak}=1$.

The parameters of slowly varying broadband forms are necessarily somewhat degenerate with the parameters of the slowly varying part of the physical model: $(b,\beta)$ and the parameters describing HCDs and UV fluctuations. As such, fits with power-law broadband terms can produce unphysical best-fit values of these parameters. Conversely, we expect that a small number of power-law terms will not significantly affect the BAO-peak parameters. Fits with the mock spectra indicate that this is the case if the power-law broadbands are restricted to the ranges $-2\leq i_{\rm max}\leq 3$, $(i_{\rm max}-i_{\rm min})\leq 3$ and $j_{\rm max}\leq 6$. Of course in the case of the mocks we know precisely the form of the physical correlation function while for the data we do not. Shifts of ($\alpha_{\perp},\alpha_{\parallel}$) in the data when using power-law terms could indicate that the underlying physical correlation function has not been modeled with sufficient precision.

We want to ensure that the power-law terms model variations of the slowly-varying part of the correlation function under the BAO peak. We therefore perform these fits only over the restricted range $40<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc}$, avoiding undue influence of the $10<r<40\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ range on the amplitudes of the power laws. Using this range, the fit with the complete physical model without power-law broadbands yields BAO parameters indistinguishable from those found over the full range $10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ (Table LABEL:fitautotable).

To choose appropriate values of $(i_{\rm min},i_{\rm max},j_{\rm max})$, we note that apart from the BAO peak, $r^{2}\xi$ does not vary wildly over the range $40<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc}$. A reasonable choice is therefore $(i_{\rm min},i_{\rm max})=(0,2)$ corresponding to a parabola and which amounts to freeing the value and first two derivatives of $r^{2}\xi_{\rm smooth}$ underneath the BAO peak. A reasonable choice for $j_{\rm max}$ is $j_{\rm max}=6$ corresponding to approximately independent broadbands in each of the four angular ranges in Fig. 14. The best-fit correlation function for $(i_{\rm min},i_{\rm max},j_{\rm max})=(0,2,6)$ and $A_{\rm peak}$ free is shown in the figure and the BAO-peak parameters on the third line of Table LABEL:fitbroadbandtable. The best-fit value of $\alpha_{\parallel}$ is $\sim 1\%$ (i.e. $\sim 0.3\sigma$) lower than the value for the complete physical model and uncertainty is unchanged. The value of $\alpha_{\perp}$ is about 3% higher than the standard value with the error increasing slightly from $\sim 0.055$ to $\sim 0.075$. The one- and two-standard deviation contours are shown in green in Fig. 15. We emphasize that these results are obtained without placing any priors on the physical parameters, other than those already used in the standard fits.

We have experimented with other combinations of $(i_{\rm min},i_{\rm max},j_{\rm max})$ in the range $j_{\rm max}\leq 6$, $-2<i_{\rm min}<0$, $-1<i_{\rm max}<3$. While most combinations yield values of $(\alpha_{\perp},\alpha_{\parallel})$ that are close to the nominal values, there are cases where there are significant changes. These cases always correspond to fits yielding unphysical values for $(b,\beta)$ of one or more components. In all cases, the large changes in $(\alpha_{\perp},\alpha_{\parallel})$ can be avoided by placing a “physical priors” on the $(b,\beta)$ that requires them to take values within two standard deviations of those found in the standard fits over the range $10<r<180\leavevmode\nobreak\ h^{-1}{\rm Mpc}$ and imposing $A_{\rm peak}=1$. The results for the $(0,2,6)$ model are listed oneline two of Table LABEL:fitbroadbandtable and the contours shown in blue in Fig. 15.

6.5 Summary of data fits

Table LABEL:fitbroadbandtable and Fig. 15 summarizes the BAO parameters found for fits with our most complete physical model of the correlation function, with and without power-law broadbands. The results of the fits with physical priors imply that the uncertainty in the form of $\xi^{\rm ph}_{\rm smooth}$ does not translate into a significant uncertainty in $(\alpha_{\perp},\alpha_{\parallel})$. Because of the priors placed on the physical parameters, this conclusion follows only if our basic physical model is not too far removed from reality. Without such priors, the uncertainty of $\alpha_{\perp}$ can be significantly increased (line 3 of Table LABEL:fitbroadbandtable).

The limits on $\alpha_{\parallel}$ are relatively stable with the inclusion of broadbands. A quantity that is more robust than either $\alpha_{\parallel}$ or $\alpha_{\perp}$ is $\alpha_{\parallel}^{\gamma}\alpha_{\perp}^{1-\gamma}$ with the exponent $\gamma$ chosen to minimize the variance. For the observed correlation and variances of $(\alpha_{\perp},\alpha_{\parallel})$, one calculates $\gamma\sim 0.7$ with

[TABLE]

where the uncertainty is chosen to cover the one-standard-deviations of the three fits in Table LABEL:fitbroadbandtable. Using the Planck model values from Table 1, we deduce

[TABLE]

This result is largely independent of the addition of power-law broadbands. It contains most of the cosmological constraints of our measurement if used in combination with constraints from the forest-quasar cross-correlation, since the latter has more constraining power on $D_{M}/r_{d}$ than the auto-correlation result.

6.6 Comparison with previous results

We can compare the results (27) and (28) with those of Delubac et al. (2015): $D_{H}(2.34)/r_{d}=9.18\pm 0.28$, $D_{M}/r_{d}=37.67\pm 2.17$ and $D_{H}^{0.7}D_{M}^{0.3}/r_{d}=14.02\pm 0.3$. The change in $D_{H}/r_{d}$ is $\sim 0.5\sigma$. To determine whether this change is due to the change in analysis procedure or in the data set we performed fits using only data used in the DR11 analysis. The results for a HCDS,UV analysis are $D_{H}(2.33)/r_{d}=9.14\pm 0.28$ and $D_{M}(2.34)/r_{d}=37.54\pm 2.17$, within $0.15\sigma$ and $0.06\sigma$ of the results of Delubac et al. (2015). This result implies that most of the change is due to the change in the data set.

The $0.5\sigma$ change in $D_{H}/r_{d}$ is typical of the variations we can expect when increasing the sample size from DR11 to DR12. We confirmed this expectations by using jackknife samples of DR12 to simulate DR11 data sets and observing changes in $D_{H}/r_{d}$. Roughly 30% of the jackknife samples had a change as large as that observed for DR11 and we therefore conclude that the DR12-DR11 difference is primarily statistical.

7 Systematic Errors in the BAO-peak position

Systematic errors in $(\alpha_{\perp},\alpha_{\parallel})$ can result if the measured correlations function is not sufficiently well modeled by the fitting function (18). Such effects can be introduced by the pipeline and analysis procedures, by improper modeling for the flux-transmission field, and by astrophysical absorption not related to cosmological large-scale structure.

Spurious correlations introduced by the pipeline were discussed in Section 5. The most important correlations were confined to the lowest $r_{\parallel}$ bin with negligible effects on the fitting of the correlation function. No evidence was found for any additional instrumental correlations that could have an influence of the measurement of $(\alpha_{\perp},\alpha_{\parallel})$. Since our model of pipeline-induced correlations introduces no scales near the BAO scale, it is not surprising that it predicts negligible effects on the BAO scale. Upper limits on the pipeline-induced shifts in $(\alpha_{\perp},\alpha_{\parallel})$ are given in Table 7.

Biases could also be introduced at the analysis level, for example through the small correlations between weights and the measured $\delta_{q}(\lambda)$. These correlations are estimated to have a negligible impact on the BAO parameters. If they did not, they would have produced biased estimates of $(\alpha_{\perp},\alpha_{\parallel})$ in the mock data sets, while no such biases are seen.

The second type of systematic error would result from incorrect modeling of the various contributions to the physical correlations. The physical correlations are parameterized by the $(b_{i},\beta_{i})$ for each transition. These parameters are marginalized over and their uncertainties therefore contribute to the reported statistical errors on $(\alpha_{\perp},\alpha_{\parallel})$.

Modeling of the contribution of unidentified HCDs and of UV fluctuations also contain free-parameters that are marginalized over. We have verified that the best-fit BAO parameters are insensitive to assumption about non-linearities and redshift evolution of the bias parameters. We have also investigated possible influence of foreground absorbers in addition to CIV. Finally, we have included power-law broadband parameters to place limits on our sensitivity to unmodeled slowly-varying contributions.

In our physical model of the correlation function we have ignored two effects that could potentially bias our results at some level. First, we have assumed that the lines of sight sample random positions of the Universe, while in a real survey we only have information in pixels in front of a quasar. This makes our measured correlation sensitive to the three-point function of quasar-lya-lya, but given that this is averaged over pixel pairs at different separations from the quasars we expect this contamination to be very smooth on BAO scales. The second effect is that the transmitted flux fraction could be affected by the relative velocity between baryons and dark matter in the early Universe (Tseliakhovich & Hirata, 2010). These relative velocities are quite small in the redshift range of interest, but it is possible that the different velocities at early times left an imprint in the distribution of neutral hydrogen even at low redshift, similar to what has been suggested in the case of galaxy clustering (Dalal et al., 2010; Yoo et al., 2011; Blazek et al., 2016). In order to study this effect in the Lyman alpha forest we would require detailed hydrodynamical simulations, and this is beyond the scope of this paper. Given that the density of quasars should have a different dependency on the relative velocity, differences in the BAO scale measured from the auto-correlation and from the cross-correlation with quasars (Font-Ribera et al., 2014) could be used to constraint the amplitude of this effect.

The third type of systematic error would be due to correlations of non-cosmological origin. One possibility is atomic or molecular absorption in the Galactic interstellar medium (ISM). The average ISM absorption is included in the fit mean flux transmission, or the subsequent subtraction of the average $\delta_{q}(\lambda)$. However, residuals correlated across the sky might still be present. Current studies show that the ISM absorption is highly correlated with dust column density. We hence use the estimated Milky-Way dust extinction as a proxy for potential ISM absorption and fit for the correlation of $\delta_{q}(\lambda)$ with this dust extinction parameter. The resulting correlation is shown in Fig. 16. The Ca H and K lines are the only known ISM lines to appear significantly. Therefore ISM absorption has a negligible impact on the Ly$\alpha$ correlation function. We have verified this conclusion by computing the auto-correlation function of the quantity $E(B-V)\times d\delta_{q}(\lambda)/dE(B-V)$.

A final source of correlations that we have considered would originate in our use of a unique quasar-continuum template which we adapt to individual quasars by the $(a_{q},b_{q})$ parameters of equation 3. Templates could be improved by making them depend on a non-forest variable, like quasar redshift or CIV emission-line strength. These improved templates could change the derived flux correlation function if the non-forest variables of the quasars are themselves correlated at large distances. This is not expected to be the case of CIV line strength since this strength is believed to be due to local conditions and therefore not correlated at large distances. We have tested this by using two templates, one for strong CIV emitters (relative to side-bands) and one for weak emitters. This slightly improves the pixel variance (by 3%) but makes no significant change in the measured correlation function and less than $0.1\%$ difference in the derived $(\alpha_{\perp},\alpha_{\parallel})$.

A continuum that evolves with redshift could be more problematic because the correlation function measurement necessarily requires quasar pairs that are near in redshift. We searched for the effect of such evolution by using a quasar continuum template that evolves linearly in redshift. The results do not differ significantly from those for the redshift-independent template.

We have searched for unsuspected systematics by dividing the data into two redshift bins, $z<2.295$ and $z>2.295$. The 95%CL ranges on $\alpha_{\parallel}$ in the high and low redshift ranges are $0.90<\alpha_{\parallel}<1.094$ and $1.02<\alpha_{\parallel}<1.15$. For $\alpha_{\perp}$ the ranges are $0.74<\alpha_{\perp}<1.5$ and $0.86<\alpha_{\perp}<1.01$. As with the DR11 results of Delubac et al. (2015), the results are consistent with the expected absence of evolution at the two-standard-deviation level.

8 Cosmological implications

Our measurements of $D_{M}/r_{d}$ and $D_{H}/r_{d}$ at $z=2.33$ are in good agreement with the values predicted by the flat-$\Lambda$CDM model consistent with the CMB data. The precision is, however, less than that of the BOSS galaxy BAO data (Alam et al., 2016), and for simple cosmological models, parameter constraints based on combinations of CMB and BAO data are not strongly changed by our results. Our constraints however, can provide tests of more general models with complicated expansion dynamics (Aubourg et al., 2015).

When combined with other BAO data, our measurement provides interesting results using only low-redshift data. For example, the low-redshift expansion history can be directly measured through the measurements of $D_{H}(z)/r_{d}$. The Hubble distance (27) can be written as

[TABLE]

Figure 17 displays this result along with those of Alam et al. (2016); Beutler et al. (2011); Ross et al. (2015) and Font-Ribera et al. (2014). The expected phases of acceleration at $z<0.7$ and deceleration at $z>0.7$ are striking.

Combining $D_{H}/r_{d}$ and $D_{M}/r_{d}$ measurements, the cosmological parameters can be measured independently of CMB data. Non-flat $\Lambda$CDM models (o$\Lambda$CDM ) have three parameters $(\Omega_{M},\Omega_{\Lambda},H_{0}r_{d})$ that determine $D_{M}/r_{d}$ and $D_{H}/r_{d}$ and constraints on $(\Omega_{M},\Omega_{\Lambda})$ can be derived by marginalizing over $H_{0}r_{d}$. This procedure is equivalent to the use of SNIa to measure $(\Omega_{M},\Omega_{\Lambda})$ by marginalizing over the parameters that determine supernova luminosities. Fig. 18 shows the constraints derived using various combinations of BAO results. By themselves, the low-redshift BAO data provide no strong evidence for $\Omega_{\Lambda}>0$. The inclusion of the Ly$\alpha$-forest data results in $\sim 10\%$ precision on $(\Omega_{M},\Omega_{\Lambda})$:

[TABLE]

consistent with spatial flatness:

[TABLE]

The fit has $\chi^{2}=12.27$ for 12 data points and three parameters. The values of $(\Omega_{M},\Omega_{\Lambda})$ agree well with the flat-$\Lambda$CDM values found using CMB anisotropies $(\Omega_{M},\Omega_{\Lambda})=(0.315,0.685)$ (Planck Collaboration et al., 2016).

9 Conclusions

In this paper, we present an analysis of the complete SDSS-3 (DR12) set of Ly$\alpha$ forests. The position of the BAO peak in the flux correlation function is in good agreement with that expected in flat-$\Lambda$CDM models favored by CMB anisotropies.

Our analysis represents a significant improvement over our previous work, both in the understanding of instrumental effects on the flux transmission and in the data analysis. The modeling of the distortion of the correlation function due to continuum fitting allows fitting the entire correlation function using a physical model without resorting to arbitrary “broadband terms”. In the future, this approach may allow alternative ways to study the cosmological parameters, for instance by eliminating the need to separate the correlation function into “smooth” and “peak” components.

Now that we have direct access to the full correlation function, further improvements in the analysis will profit from better modeling of absorption mechanisms. In particular, it will be important to fully use external constraints on HCD and metal absorption.

The improvements in modeling and data analysis that we have implemented here will clearly be important for future projects with reduced statistical errors. For example, the DESI project (DESI Collaboration et al., 2016) will have up to 15 times the number of forest-pixel pairs resulting in a significant decrease of statistical errors of the correlation function. The understanding the physics of the broadband correlation function may be the limiting factor in the extraction of BAO parameters.

The cosmological information from the Ly$\alpha$-forest auto-correlation function studied here is complemented by the study of the cross-correlation between the forest and neighboring quasars. Results and implications of a new measurement of this function will be presented in a forthcoming publication (du Mas des Bourboux et al., in preparation).

Acknowledgements.

Funding for SDSS-III has been provided by the Alfred P. Sloan Foundation, the Participating Institutions, the National Science Foundation, and the U.S. Department of Energy Office of Science. The SDSS-III web site is http://www.sdss3.org/. SDSS-III is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS-III Collaboration including the University of Arizona, the Brazilian Participation Group, Brookhaven National Laboratory, Carnegie Mellon University, University of Florida, the French Participation Group, the German Participation Group, Harvard University, the Instituto de Astrofisica de Canarias, the Michigan State/Notre Dame/JINA Participation Group, Johns Hopkins University, Lawrence Berkeley National Laboratory, Max Planck Institute for Astrophysics, Max Planck Institute for Extraterrestrial Physics, New Mexico State University, New York University, Ohio State University, Pennsylvania State University, University of Portsmouth, Princeton University, the Spanish Participation Group, University of Tokyo, University of Utah, Vanderbilt University, University of Virginia, University of Washington, and Yale University. The French Participation Group of SDSS-III was supported by the Agence Nationale de la Recherche under contracts ANR-08-BLAN-0222 and ANR-12-BS05-0015-01. NB and JG acknowledge support from the French Programme National Cosmologie et Galaxies. MB, MMP and IP were supported by the A*MIDEX project (ANR- 11-IDEX-0001-02) funded by the “Investissements d’Avenir” French Government program, managed by the French National Research Agency (ANR), and by ANR under contract ANR-14-ACHN-0021. The work of JB and KD was supported in part by the U.S. Department of Energy, Office of Science, Office of High Energy Physics, under Award Number DE-SC0009959. N.P.R acknowledges support from the S.T.F.C. and the Ernest Rutherford Fellowship scheme.

Appendix A Modifications applied to the BOSS spectral extraction pipeline

Each of the two BOSS spectrograph slitheads consists of 500 fibers organized in 25 aligned blocks (or bundles) of 20 fibers (see Smee et al. 2013 for a description of the BOSS spectrographs). Fibers in a bundle are separated by 260 $\mu$m whereas the separation between fibers of adjacent bundles is 624 $\mu$m. This design allows one to treat independently each fiber bundle in the data reduction: whereas the overlap of light coming from fibers of the same bundle must be accounted for, one can safely ignore interbundle contamination.

The spectral processing consists in a row-by-row extraction from each CCD, where the signal of a row of pixels illuminated by a fiber bundle is fit independently from the rest of the data (wavelength are dispersed approximately along CCD columns). The extraction consists in minimizing the following quantity

[TABLE]

where $D$ are the pixel values (after preprocessing, i.e. after the subtraction of the bias, conversion from digital units to electrons, and flat-fielding), $w$ the weights assigned to each pixel, $P$ the true cross-dispersion PSF profile of a fiber and $\hat{P}$ the profile used (a Gaussian). The quantity $\hat{N}$ is the number of collected electrons from a fiber, to be estimated. The last term is a low-order polynomial contribution used to fit any residual light in the row; it compensates both for scattered light and the imperfect modeling of the cross-dispersion profile tails.

For the extraction to be optimal, one must have $w_{i}=var_{i}^{-1}$, where $var_{i}$ is the variance of $D_{i}$. This pixel variance is, however, not known a priori. In the BOSS DR12 pipeline, the pixel data were directly used as an estimate of the expected signal: $var_{i}=var_{read}+D_{i}$ if $D_{i}>0$ and $var=var_{read}$ if not, where $var_{read}$ is the readout noise variance, and $D_{i}$ is here an estimate of the Poisson noise.

This use of the data to determine the variance results in a biased estimate of the flux. As an illustration, if one considers a single fiber extraction, without the background polynomial contribution, the $\chi^{2}$ minimization yields

[TABLE]

To first order, assuming $D_{i}=NP_{i}+\delta_{i}>0$, and noting $\bar{w}_{i}=1/var(\delta_{i})$ the expectation value of $w_{i}$, one has $w_{i}\simeq\bar{w}_{i}(1-\bar{w}_{i}\delta_{i})$ and

[TABLE]

The first term produces an unbiased estimate of $N$ if $\hat{P}=P$, which is not strictly the case. This bias is in principle corrected for by the calibration, but since the weights $\bar{w}$ depend on the amount of signal, the calibration obtained with bright spectra does not perfectly correct the bias of faint spectra. The second term is a systematic underestimation of the signal; it can be interpreted as an effective number of pixels contributing to the measurement (if the profile is normalized, with $\sum_{i}P_{i}=1$).

The fractional bias is more important at low fluxes, and as a consequence i) this bias has a non trivial wavelength dependence correlated with bright sky emission lines and ii) its amplitude relative to the calibrated flux depends on the throughput and hence the observation conditions. Those two effects could lead to correlated biases in the Ly$\alpha$ forest.

We have rerun the data processing using $w_{i}=1/var_{read}$ (this weight is in practice not a constant quantity per CCD amplifier because of the flat-field correction). We also improved the outlier rejection procedure that discards unmasked pixels affected by cosmic rays.

This new approach to extraction is unbiased but less optimal. Expressing this procedure using a matrix representation, if we define $X$ as the parameters of the model (fluxes and polynomial corrections) such that we can interpret the data as $D=HX+\delta$ (where $\delta$ is the noise), the optimal extraction results in a covariance of the parameters $C_{X}=(\sum_{i}w_{i}H_{i}H_{i}^{T})^{-1}$, whereas the sub-optimal extraction produces a larger covariance $C_{X}^{\prime}=C_{X}(\sum w_{i}^{2}var_{i}H_{i}H_{i}^{T})C_{X}$. We use this latter expression in the subsequent steps of the pipeline, using the estimated flux (and not directly the pixel values) to evaluate the Poisson noise contribution to $var_{i}$.

The ratio of uncertainties does not exceed 5% over the dynamic range of BOSS observations of quasars. It reaches a higher value of 10% for the brightest measurements (see Fig. 19).

Subsequent steps of the data processing (fiber flat-field, sky spectrum model, photometric calibration) also involve the weighted average of the data. Similarly to the extraction, using directly the flux inverse variances as weights results in a bias because of the flux-inverse-variance correlation. We hence took care of replacing the inverse variances of spectra by their average over many wavelength (using a spline fit) to mitigate this bias and to simultaneously minimize the loss of optimality. Those changes have been incorporated in the pipeline used for the SDSS data release DR13 (SDSS Collaboration et al., 2016).

Appendix B The bias parameters of the Ly$\alpha$ forest

The physical model of the Ly$\alpha$ forest autocorrelation function involves two independent parameters for the linear biasing: the density bias factor $b_{\rm Ly\alpha}$, and the redshift distortion factor $\beta_{\rm Ly\alpha}$. The values for our best fit with errors are given in Table 5. As usual, we give the value of $b_{\rm Ly\alpha}(1+\beta_{\rm Ly\alpha})$ because its error is smaller than the error for $b_{\rm Ly\alpha}$ and less correlated with $\beta_{\rm Ly\alpha}$ (see Slosar et al. (2011)). Some variations of these parameters occur as the model is improved to include the effects of metal lines, and possible effects of HCDs and ionizing intensity fluctuations. As expected, there is a degeneracy in the values of various parameters of the model that affect broadband terms. In particular, the value of $\beta_{\rm Ly\alpha}$ shows a substantial change from 1.2 to 1.7 depending on the model that is used. The first model, assuming that only Ly$\alpha$ forest correlation is present with a constant ionizing intensity, does not yield a good fit. Taking into account the effect of metal lines is the main improvement of the model that allows for an acceptable fit, and does not change appreciably the value of $\beta_{\rm Ly\alpha}$ although errors are substantially increased. The addition of the effects of HCDs and UV fluctuations further improves the model and increases substantially the best fit value of $\beta_{\rm Ly\alpha}$.

To compare our values of $b$ and $\beta_{\rm Ly\alpha}$ with the previous determination from the DR11 data from Blomqvist et al. (2015), we first translate their result to the cosmological model we use in this paper, taking into account that the measurement of the Ly$\alpha$ autocorrelation depends only on $b\sigma_{8}(\bar{z})$, where $\bar{z}=2.33$ is the mean redshift of our Ly$\alpha$ forest correlation measurement. Blomqvist et al. (2015) measured $\beta_{\rm Ly\alpha}=1.39\pm 0.11$, and $b_{\rm Ly\alpha}(1+\beta_{\rm Ly\alpha})=-0.374\pm 0.007$, using a model with $\Omega_{m}=0.27$ and $\sigma_{8}=0.7877$, implying $\sigma_{8}(\bar{z})=0.3072$. The model used in this paper has $\Omega_{m}=0.3147$ and $\sigma_{8}=0.8298$, implying $\sigma_{8}(\bar{z})=0.3131$. Therefore the Ly$\alpha$ bias of Blomqvist et al. (2015) translated to our model is $b_{\rm Ly\alpha}(1+\beta_{\rm Ly\alpha})=-0.367\pm 0.007$. This should be compared with our result on the first row of Table 5, because Blomqvist et al. (2015) did not consider the effects from metals, HCDs or UV fluctuations

The main reason for the difference with the values of $b_{\rm Ly\alpha}(1+\beta_{\rm Ly\alpha})$ and $\beta_{\rm Ly\alpha}$ we obtain here for our Ly$\alpha$ forest model applied to DR12 given in Table 5 is the different radial range for the fit: Blomqvist et al. (2015) restricted their fit to $r>40\leavevmode\nobreak\ h^{-1}{\rm Mpc}$, whereas we use $r>10\leavevmode\nobreak\ h^{-1}{\rm Mpc}$. There are other details that contribute to the difference, e.g., the continuum distortion and bin smoothing modeling, but the radial range of the fit is most important. As mentioned earlier, the $\chi^{2}$ value shows that our Ly$\alpha$ fit is not good, reflecting the fact that the shape of the correlation cannot be properly fit over our broad radial range. This is why the values of $\|b_{\rm Ly\alpha}\|\,(1+\beta_{\rm Ly\alpha})$ and $\beta_{\rm Ly\alpha}$ decrease when extending the fit to smaller separations. The better fits obtained for the more complete models in Table 5 show that $b_{\rm Ly\alpha}(1+\beta_{\rm Ly\alpha})$ is relatively stable, but the value of $\beta_{\rm Ly\alpha}$ can be more model dependent.

These results can also be compared to theoretical predictions from hydrodynamic simulations of the intergalactic medium, which were compared to the observational measurement of Blomqvist et al. (2015) in Section 7 of Arinyo-i-Prats et al. (2015). There, it was noticed that while the predicted value of $\beta_{\rm Ly\alpha}$ is roughly in agreement with the observational determination, the value of $\|b_{\rm Ly\alpha}\|\,(1+\beta_{\rm Ly\alpha})$ was predicted to be lower than observed. Now, the lower value of $\|b_{\rm Ly\alpha}\|\,(1+\beta_{\rm Ly\alpha})$ we measure from DR12, correcting for the effect of metals and extending the fit to smaller separations, is in much better agreement with the range of the expected theoretical values. Despite the considerable uncertainties that are still present both in the predictions from hydrodynamic simulations of the Ly$\alpha$ forest (as discussed in Arinyo-i-Prats et al. (2015)), and in the observational determinations that are affected by the modeling of metals, HCDs and ionizing intensity fluctuations, we think that the rough agreement we now find of two independently predicted linear bias factors of the Ly$\alpha$ forest is fairly remarkable. In the future, we expect to be able to better measure these linear factors and their redshift evolution.