The Effect of Interstellar Absorption on Measurements of the Baryon Acoustic Peak in the Lyman-{\alpha} Forest

TL;DR

This paper investigates how interstellar absorption lines in the Milky Way might bias measurements of the baryon acoustic peak in the Lyman-alpha forest, finding current effects negligible but potentially significant for future surveys like DESI.

Contribution

It introduces a method to assess and mitigate interstellar absorption effects on baryon acoustic peak measurements in upcoming large-scale surveys.

Findings

Current interstellar absorption effects are too small to impact existing measurements.

Future surveys like DESI may need to account for these effects.

A stacking method was used to generate absorption line maps from over a million spectra.

Abstract

In recent years, the autocorrelation of the hydrogen Lyman-{\alpha} forest has been used to observe the baryon acoustic peak at redshift 2 < z < 3.5 using tens of thousands of QSO spectra from the BOSS survey. However, the interstellar medium of the Milky-Way introduces absorption lines into the spectrum of any extragalactic source. These lines, while weak and undetectable in a single BOSS spectrum, could potentially bias the cosmological signal. In order to examine this, we generate absorption line maps by stacking over a million spectra of galaxies and QSOs. We find that the systematics introduced are too small to affect the current accuracy of the baryon acoustic peak, but might be relevant to future surveys such as the Dark Energy Spectroscopic Instrument (DESI). We outline a method to account for this with future datasets.

Figures

Equations

$$T(z)=\frac{f(z)}{C_q(z)}$$

$${\delta }_F(z)=\frac{T(z)}{\bar{T}(z)}-1$$

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

$$w_{ij}\propto \frac{{\left[(1+z_i)(1+z_j)\right]}^{\gamma /2}}{{\xi }_{ii}^2{\xi }_{jj}^2}$$

$${\xi }_{ii}^2=\frac{{\sigma }_{\text{pipeline},i}^2}{\eta (z_i)}+{\sigma }_{\text{LSS}}^2(z_i);{\sigma }_{\text{pipeline},i}^2=\frac{1}{\mathrm{ivar}\cdot (C_q\bar{T}{)}^2}$$

$$w_i=(1+z_i{)}^{\gamma /2}{\left[\frac{1}{\mathrm{ivar}\cdot {\left(C_q{\bar{T}}_i\right)}^2\eta (z_i)}+{\sigma }_{\text{LSS}}^2(z_i)\right]}^{-1}$$

$$f_1{f}_2\frac{\Delta {\lambda }_1\Delta {\lambda }_2}{(\Delta {\lambda }_{\text{survey}}{)}^2}.$$

$$0.05^{2}\frac{($5\text{\,}\mathrm{\SIUnitSymbolAngstrom}$)^{2}}{($500\text{\,}\mathrm{\SIUnitSymbolAngstrom}$)^{2}}=$2.5\text{\times}{10}^{-6}\text{\,}$$

$

Full text

††pubyear: 2017††pagerange: The Effect of Interstellar Absorption on Measurements of the Baryon Acoustic Peak in the Lyman\Hyphdash*$\boldsymbol{\upalpha}$ Forest–The Effect of Interstellar Absorption on Measurements of the Baryon Acoustic Peak in the Lyman\Hyphdash*$\boldsymbol{\upalpha}$ Forest

1 Introduction

The baryon acoustic peak is a cosmological feature of a predictable physical size, which can be used to constrain cosmological parameters that describe the expansion of the universe. It was first discovered in the power spectrum of temperature differences in the Cosmic Microwave Background (CMB) radiation, where it is known as Baryon Acoustic Oscillations (BAO; de Bernardis et al. 2000, Keisler et al. 2011; Sievers et al. 2013; Hinshaw et al. 2013; Planck Collaboration et al. 2014). The same phenomenon can be seen as a single peak in the autocorrelation of the baryon density, and it has first been observed in the autocorrelation of galaxy number counts (Eisenstein et al., 2005).

The Sloan Digital Sky Survey (SDSS; York et al. 2000) has so far obtained spectra for about $2.5\text{\times}{10}^{6}\text{\,}$ galaxies, spread over about a third of the sky. As part of SDSS, the Baryon Oscillation Spectroscopic Survey (BOSS; Dawson et al. 2013), has observed almost $3\text{\times}{10}^{5}\text{\,}$ QSOs (Pâris et al., 2017). In recent years, there have been a number of works which used BOSS data to measure the baryon peak using absorption in the Lyman\Hyphdash*$\upalpha$ forest of QSOs at a redshift range of about $2<z<3.5$ (Busca et al., 2013; Slosar et al., 2013; Delubac et al., 2015).

The baryon acoustic peak is a weak feature, and the typical signal-to-noise ratio (SNR) of the Lyman\Hyphdash*$\upalpha$ forest is low in an individual BOSS spectrum. The aggregate signal from over ${\mathord{\sim}\mkern 1.0mu}{}10^{4}$ QSOs is required in order to achieve the required SNR in the autocorrelation function. This makes this process potentially sensitive to systematics that could be in the data, even if they cannot be detected in a single QSO spectrum. In this paper, we focus on one such effect, which is intervening absorption lines formed in the interstellar medium (ISM) of the Milky Way (MW). Interstellar emission (see e.g., Brandt & Draine 2012) will not be considered here.

ISM absorption at some wavelengths, if not accounted for, could mimic Lyman\Hyphdash*$\upalpha$ at particular redshifts. If there are correlations between absorption lines, or between sight lines, they could introduce a spurious signal to the 2D autocorrelation function of the Lyman\Hyphdash*$\upalpha$ forest, and bias the shape or position of the baryon acoustic peak. This possible effect was first noted in Lee et al. (2013), who provided a global correction for the spectra.

Poznanski et al. (2012) followed by Baron et al. (2015a, b), have used the integrated MW absorption spectrum in order to study correlations of absorption features with dust extinction, and the properties of the Diffuse Interstellar Bands (DIBs; see also Murga et al. 2015 and Lan et al. 2015). They achieved this by stacking hundreds of thousands of SDSS spectra, in which the lines are not detected individually. The SDSS coverage allowed them to probe different extinction regimes, as well as different directions on the sky. In this work we use the same methodology in order to measure the ISM features. In order to calculate the autocorrelation of the Lyman\Hyphdash*$\upalpha$ forest, we re-created a pipeline which is similar to previous works (Busca et al., 2013; Slosar et al., 2013; Delubac et al., 2015). Combining these two approaches, we study the ISM contribution in isolation from other effects.

Near the final stages of this work, Bautista et al. (2017) published an updated measurement of the Lyman\Hyphdash*$\upalpha$ BAO autocorrelation from 12th data release (DR12) of SDSS BOSS (Alam et al., 2015). There, they briefly discuss the possible systematic contribution of the MW ISM to the BAO signal, which we study here in more detail.

2 Baryon Peak Pipeline

In order to measure the ISM contribution to the autocorrelation of the Lyman\Hyphdash*$\upalpha$ forest, we need a pipeline to calculate the autocorrelation estimator. We follow the method used in Busca et al. (2013), Slosar et al. (2013), and Delubac et al. (2015), with some modifications as discussed below.

2.1 Data selection

We use $176\,944\text{\,}$ QSOs from BOSS DR12 with redshift $2.1<z<3.5$ as selected by Ross et al. (2012), and observed as part of the BOSS program with the upgraded SDSS spectrographs (Smee et al., 2013). This redshift range translates to a wavelength range of about $3600\text{\,}\mathrm{\SIUnitSymbolAngstrom}5500\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. The lower limit is set by atmospheric absorption of UV light, whereas the higher limit is set by diminishing observed QSO numbers at higher redshifts. Note that this wavelength range falls within the blue cameras of the BOSS spectrographs. In the rest frame, we restrict ourselves to the range $1040\text{\,}\mathrm{\SIUnitSymbolAngstrom}1200\text{\,}\mathrm{\SIUnitSymbolAngstrom}$.

Within the spectra we discard pixels as follows. First, SDSS provides a variety of flags that are derived from the subexposures, and then propagated to the reduced stacked spectrum with both AND and OR operators. The OR operator corresponds to pixels that are marked in a single sub exposure, while the AND operator corresponds to pixels that are marked in all the sub exposures of a given source.We discard all pixels that are flagged in all sub-exposures, but we use only specific OR mask flags. Out of those, the BRIGHTSKY flag is the most common, and causes the removal of $18\text{\,}$ percent of the data. A few other flags (BADTRACE, BADFLAT, MANYBADCOLUMNS) discard up to 2 additional percent. We also mask the telluric Hg i line with a $4356\text{\,}\mathrm{\SIUnitSymbolAngstrom}4364\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ mask. This line is flagged by the SDSS pipeline in most of the spectra, but we find some residuals in other spectra as well.

Broad Absorption Lines (BALs) in the rest frame of the QSOs overlap with the Lyman\Hyphdash*$\upalpha$ forest and introduce noise that we wish to avoid. We mask BALs in our spectra using the BOSS catalogue (Pâris et al., 2017, DR12Q_BAL). We mask out the following lines: C iv, O vi, NV , Si iv+O iv as well as Lyman\Hyphdash*$\upalpha$, based on the supplied BAL reference frames (minimum and maximum velocities). We add a fractional safety margin of $2\text{\times}{10}^{-3}\text{\,}$ of the wavelength to the line widths. For the same reason, we mask Damped Lyman\Hyphdash*$\upalpha$ Absorption systems (DLAs), using a catalogue provided by Garnett et al. (2016). Similarly to Delubac et al. (2015), we remove regions with more than $20\text{\,}$ percent absorption and correct the damping wings using a Lorentzian profile for Lyman\Hyphdash*$\upalpha$ and Lyman\Hyphdash*$\upbeta$, ignoring possible metal lines from these systems.

2.2 Continuum fit

We implement a continuum fit process based on MF-PCA (see Suzuki et al. 2005; Pâris et al. 2011; Lee et al. 2012), with some minor differences in the fit process. First, we do not perform a visual inspection of all the fitted continua (except for a small fraction of objects). Lee et al. (2012) fitted redshift correction, power-law correction and flux normalization parameters simultaneously with the PCA coefficients, using the Levenberg-Marquardt (LM) algorithm. We perform LM on power-law correction and flux normalization, with a linear least squares for the PCA coefficients in every iteration. We find that allowing the redshift to vary sometimes produces unreliable results. For every QSO, we calculate a ‘goodness of fit’ value, defined as the mean of the absolute fit errors, over a wavelength range of $1216\text{\,}\mathrm{\SIUnitSymbolAngstrom}1600\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. We bin the QSOs according to their SNR, and discard those with the lowest goodness of fit in every bin. To account for the lack of visual inspection, we discard the worst $10\text{\,}$ percent of the continuum fits, instead of the $5\text{\,}$ percent discarded in Lee et al. (2012). Assuming the number of pairs is quadratic with QSO density, we thus remove about 20 percent of the pixel pairs. At the end of this stage we discard 29662 QSOs due to bad fit, 5967 QSOs in which the continuum level is comparable to the noise, and 546 QSOs in which the relevant part of the spectrum contains less than 150 valid pixels. The Lyman\Hyphdash*$\upalpha$ forest autocorrelation is computed from the remaining 140769 QSOs. This is a similar number to the total used by Delubac et al. (2015) who used DR11 while we use DR12, but due to our more conservative cuts, and lack of visual inspection, we reject more QSOs.

2.3 The autocorrelation estimator

We obtain the transmittance, $T(z)$ for every QSO, where $z$ is the Lyman\Hyphdash*$\upalpha$ equivalent redshift of the wavelength, $f(z)$ is the QSO intrinsic flux, and $C_{q}(z)$ is the estimated continuum.

[TABLE]

$\bar{T}(z)$ is the mean transmittance, where the average is done over all QSOs, as a function of redshift. The fractional change in transmittance, $\delta_{F}$, is defined as:

[TABLE]

The overall slope predicted by the continuum fit is not exact, especially for higher redshift QSOs, where the Lyman\Hyphdash*$\upalpha$ forest spans a wider wavelength range. We define $\delta_{F}(D_{C}(z))$ to be $\delta_{F}(z)$ as a function of the line-of-sight comoving distance, $D_{C}(z)$. We detrend $\delta_{F}(D_{C}(z))$ of Lyman\Hyphdash*$\upalpha$ forests that span more than $500\text{\,}\mathrm{Mpc}$ ($h_{0}=0.7$), using a weighted nonuniform boxcar filter with a window size of $\pm 300\text{\,}\mathrm{Mpc}$. We find, by trial and error, that the window size is large enough to avoid affecting the 2D autocorrelation in a way that distorts the baryon peak.

After calculating $\delta_{F}$, for all QSOs, following Busca et al. (2013), we remove a residual bias, $\left<\delta_{F}(z)\right>$, in every redshift bin. We then calculate the autocorrelation estimator:

[TABLE]

where A is a bin of comoving parallel and transverse distances between pairs of pixels. The sum over $i,j$ represents all pixel pairs with the required distances.

2.4 Pixel weights

We use a model for the weights given by Busca et al. (2013), based on McDonald et al. (2006). The expression for the weight of a pair of pixels is:

[TABLE]

Where:

[TABLE]

This model takes into consideration:

• •

The measured redshift dependence of the correlation function $(1+z)^{\gamma}$ where $\gamma=3.8$.

• •

The flux error for each point in the Lyman\Hyphdash*$\upalpha$ forest $\sigma_{\text{pipeline},i}^{2}$, which in turn is based on flux error ($\operatorname{ivar}$).

• •

Variance introduced by the large scale structure $\sigma_{\text{LSS}}$.

• •

A scaling factor for the contribution of pipeline variance to the correlation variance $\eta$.

We use the values obtained by Busca et al. (2013) for $\sigma_{LSS}(z)$ and $\eta(z)$.

The contributions of $i$ and $j$ can be separated to the product $w_{ij}=w_{i}w_{j}$, so that weights can be calculated per pixel rather than for every pair:

[TABLE]

2.5 Implementation

Our pipeline is written in Python+NumPy, with the critical path in the autocorrelation calculation written in Cython. We use the reference Python implementation, which is a bytecode interpreter, and as such it can be 2–3 orders of magnitude slower than compiled languages such as C. Each instruction has to be translated at runtime, which incurs an overhead. Numpy offers improved performance by including pre-compiled versions of common mathematical calculations, like matrix operations. However, some problems do not map efficiently to the ‘recipes’ provided by Numpy. Cython uses a simple syntax similar to python, while achieving performance that is comparable to C. Using Cython, we can write loops explicitly and control the iteration boundaries, which allows us to ignore pixel pairs that are too far away in redshift.

Given two QSOs, $A$ and $B$, with Lyman\Hyphdash*$\upalpha$ forest sizes $N$ and $M$, respectively, calculating their correlation contribution using matrix operations requires all $N\times M$ combinations of pixel pairs to be computed. A typical sightline can have a range of $200\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}600\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$, while the total range $1.95<z<3.5$ corresponds to $1700\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$. Given a point in QSO A, the matching points within $200\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$ in QSO B all lie in a limited distance range. Therefore, the overlap between sightlines can be small, and most of the matrix elements represent points that are too far apart. Because distances are monotonically increasing, it is possible to choose the iteration boundaries of the inner loop (over QSO B) to avoid pairs which cannot be found within the required range of $200\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$ (see fig. 1). Suppose we know the maximal range of pixels in QSO B, $\left[B_{i-n},B_{i+k}\right]$ that are in range of $A_{i}$. It is possible to find the maximal range $\left[B_{i+1-n^{\prime}},B_{i+1+k^{\prime}}\right]$ corresponding to the next point in $A$, $A_{i+1}$, without iterating over all points in $B$ for every $A_{i}$, because the series of range start (or end) points is monotonically increasing as well.

In order to find the set of all QSO pairs, we choose QSOs that are separated by less than [math] which corresponds to about $200\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$ at $z=1.9$. The $1.5\text{\times}{10}^{5}\text{\,}$ QSOs form about $2\text{\times}{10}^{10}\text{\,}$ potential pairs, out of which only $5\text{\times}{10}^{7}\text{\,}$ unique pairs are within the required angular range. We use a KD-Tree implementation in the astropy package, which can perform a k-Nearest Neighbor Search in $O(k\log n)$. To parallelize the search, each core is assigned a different subset of the QSOs, to check against the full list of QSOs.

The pipeline supports parallelization through MPI, using the Mpi4Py package. Performance is about $3\text{\times}{10}^{7}\text{\,}$ pixel pairs per second per core. The BOSS Lyman\Hyphdash*$\upalpha$ forest yields ${\mathord{\sim}\mkern 1.0mu}$3\text{\times}{10}^{12}\text{,}$$ pixel pairs and can be calculated within an hour on a 32-core machine.

Every pixel pair in our pipeline is binned according to sky region (HealPix bin), parallel comoving separation, transverse comoving separation, and parallel comoving distance (the average distance of the pixel pair). The range and resolution on each axis can be changed to obtain different cuts of the results (a tradeoff is required to prevent the size of the 4-dimensional array from becoming too large).

We make the pipeline source code available at:

https://github.com/yishayv/lyacorr.

2.6 Lyman\Hyphdash*$\boldsymbol{\upalpha}$ autocorrelation results

Figure 2 shows the autocorrelation estimator, $\hat{\xi}$, as a function of parallel and transverse comoving separation, $({r_{\parallel}},{r_{\perp}})$. We include the autocorrelation estimator from Delubac et al. (2015) for comparison (received by private communication). Figure 3 shows the autocorrelation estimator multiplied by distance squared in three angular wedges where $\mu$ is the cosine of the angle of the coordinate $({r_{\parallel}},{r_{\perp}})$. There is faint trace of the BAO signal, at small angles (blue line), with an amplitude of about ${\mathord{\sim}\mkern 1.0mu}$5\text{\times}{10}^{-5}\text{,}$$ near $100\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$.

Since our purpose is not to perform a full cosmological measurement, we did not derive the covariance matrix of the autocorrelation estimator, so we do not have a robust error estimation, in contrast to previous works. However, the autocorrelation estimator values are similar to previous results with the same data, with slightly more noise along ${r_{\perp}}$, appearing as faint vertical stripes in fig. 2. We can therefore proceed to study the ISM imprint on the cosmological signal.

3 Stacked ISM Spectra

Poznanski et al. (2012) have used SDSS spectra stacked in the observer frame to measure the correlation of dust extinction with the Na i D absorption doublet in the ISM of the MW. Baron et al. (2015a, b) have expanded that process to the entire wavelength range of the SDSS spectrograph to detect multiple DIBs, study their properties, and measure the correlations between DIBs and with other species.

Baron et al. (2015a) used only data redder than $3800\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, corresponding to a Lyman\Hyphdash*$\upalpha$ redshift of 2.1, due to the limits of the original SDSS spectrographs (Smee et al., 2013). We create ISM spectra using a similar process, using only spectra from the BOSS spectrograph, with a limit of $3600\text{\,}\mathrm{\SIUnitSymbolAngstrom}$, or a Lyman\Hyphdash*$\upalpha$ redshift below 2. We use $1.5\text{\times}{10}^{6}\text{\,}$ galaxies and $3\text{\times}{10}^{5}\text{\,}\mathrm{Q}$SOs from the 12th data release (DR12). For every spectrum, we estimate the continuum using a Savitzky-Golay filter with a $80\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ window (Savitzky & Golay, 1964)111Savitzky-Golay is a generalization of ‘running mean’ where instead of fitting a constant in the window (the mean) one fits a polynomial.. Filters have edge effects, they are ill defined near the edges of a sequence (spectrum in our case), over a region proportional to the window size. We use a narrower window than in the previous works, in order to discard less blue pixels, our region of interest. For every object (galaxy or QSO), we divide the original spectrum by the continuum to obtain a (noisy) detrended spectrum of relative absorption. Pixels where the continuum was nonpositive or comparable to noise (below $0.5\text{\times}{10}^{-17}\text{\,}\mathrm{erg}\mathrm{/}\mathrm{s}\mathrm{/}\mathrm{c}\mathrm{m}^{2}\mathrm{/}\mathrm{\SIUnitSymbolAngstrom}$) were discarded.

We then calculate the median of the detrended spectra, at every wavelength, in 20 bins of extinction, using the maps of Schlegel et al. (1998). The bin boundaries are chosen so that there is a similar number of spectra in each bin. This stacking allows us to reach the SNR necessary to detect the weak absorption features we seek. As discussed below, these stacks allow us to address only some of aspects of the ISM bias on the BAO measurement.

Since, generally speaking, QSOs are blue and galaxies are red, the blue end of the median spectra is dominated by QSOs (fig. 4), and the red by galaxies, where their much larger number contributes in addition to their colour. The split point occurs at roughly $4400\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ which corresponds to a redshift of ${\mathord{\sim}\mkern 1.0mu}2.6$. The Lyman\Hyphdash*$\upalpha$ autocorrelation data spans $2<z<3.5$, but most of the statistical weight is near $z\sim 2.2$. This means that for the most part, our ISM estimation relies on the QSO spectra, which include the Lyman\Hyphdash*$\upalpha$ QSOs and lower redshift QSOs, and we cannot leverage significantly the power in the galaxy spectra. Lyman\Hyphdash*$\upalpha$ QSOs account for about half of the flux from QSOs in the $u$ and $g$ bands.

The stacked ISM spectra are shown in fig. 5. Each row represents a spectrum in a single extinction bin. The mean absorption has been removed to show only features that are correlated with extinction. One can see that the spectra are very noisy at short wavelengths, where there is a penury of photons, and that only very strong lines appear. The Ca ii H&K lines ($3968\text{\,}\mathrm{\SIUnitSymbolAngstrom}3933\text{\,}\mathrm{\SIUnitSymbolAngstrom}$) are clearly the most notable feature in our range. The artifact at ${\mathord{\sim}\mkern 1.0mu}$5580\text{,}\mathrm{\SIUnitSymbolAngstrom}$$ is caused by the transition between the red and blue spectrographs. The Na i doublet can be seen near $5890\text{\,}\mathrm{\SIUnitSymbolAngstrom}$. Also visible are the DIBs at $5780\text{\,}\mathrm{\SIUnitSymbolAngstrom}$ and, to a lesser extent, $5797\text{\,}\mathrm{\SIUnitSymbolAngstrom}$.

4 Estimating the ISM influence on the BAO signal

We proceed to estimate the influence of the ISM absorption lines on the estimated baryonic peak. We identify two potential ways for the ISM to produce a BAO-like feature in the autocorrelation function: The first is correlation between different lines at different wavelengths. Such correlations would mimic intrinsic Lyman\Hyphdash*$\upalpha$ correlations in the parallel direction, along the line of sight, and have no angular dependence. Additionally, there could be angular correlations in all (or some) of the line strengths, due to preferred scales in the MW gas distribution.

4.1 Wavelength/parallel component

For correlations between lines to bias the BAO measurement, lines would need to have wavelength separations that correspond to a line-of-sight comoving distance of the same order as the baryon peak. The spurious feature created by such lines would be at a comoving distance that is the mean of the apparent comoving distance of the lines (when wrongly interpreted as Lyman\Hyphdash*$\upalpha$).

To create a BAO-like signal of the same strength (${\mathord{\sim}\mkern 1.0mu}$5\text{\times}{10}^{-5}\text{,}) using absorption line pairs, we would need relative absorption which is the square root of the signal strength, ${\mathord{\sim}\mkern 1.0mu}$7\text{\times}{10}^{-3}\text{\,}, for a particular redshift value. However, because we average over a large redshift range, multiple pairs of different line transitions would have to conspire to create a single feature that would not be smoothed.

Consider two absorption lines with effective profile widths $\Delta\lambda_{1}$ and $\Delta\lambda_{2}$ and average relative absorptions $f_{1},f_{2}$, over a redshift range $\Delta\lambda_{\text{survey}}$. For simplicity, we treat the measured absorption profiles as boxcar functions. In addition, for simplicity, we disregard the nonlinear relation between wavelength/redshift and comoving distance.

Assuming the correlation function resolves the lines (in comoving space), their contribution (without taking noise into account) would be:

[TABLE]

The effect is attenuated by the fact that we measure absorption relative to the mean transmittance over the whole survey. On the other hand, the statistical weight contribution varies with wavelength. Therefore, it would be wise to take a smaller, more conservative estimate for $\Delta\lambda_{\text{survey}}$. For BOSS, most of the data is around redshift $2.1\text{\,}2.5\text{\,}$, so the wavelength range is about $500\text{\,}\mathrm{\SIUnitSymbolAngstrom}$.

For the strongest lines in our range of interest, the Ca ii H&K, we find in the stacked ISM spectra a relative absorption of $10\text{\,}$ percent, but less than $3\text{\,}$ percent change in absorption between the highest and lowest extinction bins. Using a rough value of $5\text{\,}$ percent, we obtain:

[TABLE]

This number is 20 times smaller than the BAO amplitude of ${\mathord{\sim}\mkern 1.0mu}$5\text{\times}{10}^{-5}\text{,}$$, however it assumes that our stacks based on the dust maps can recover the full potential correlation between the Ca ii H&K lines.

We estimate the bias that could be introduced to the BAO peak estimation by ISM lines. As we show in section 4.3, the current data, especially in the bluest part of the spectrum, are insufficient to clearly detect most lines, so we cannot rule out hypothetical correlations from lines that we know exist in this wavelength range (e.g., Morton 2003). As a proxy for these lines, we assume the amplitude of the Ca ii H&K-induced correlation (section 4.4), but allow it to shift in apparent comoving distance, in order to simulate the potential effect of ISM lines which may remain hidden in our noise. We fit a Gaussian curve to the sum of the real and spurious signals to find the displacement in the primary peak position. We find that the maximum bias is obtained when the lines are positioned ${\mathord{\sim}\mkern 1.0mu}$0.7\text{,}\sigma$$ from the primary peak, where $\sigma$ is the standard deviation of the BAO signal. We obtain a value of about $0.3\text{\,}$ percent to the BAO peak displacement, if indeed the ISM correlation is about 20 times smaller than the cosmological signal. Specifically for the Ca ii H&K lines, the maximum displacement, of order $0.05\text{\,}$ percent, occurs in the $0.1<\mu<0.5$ angular wedge. The baryon peak has a constant $r$, whereas ISM lines have a constant ${r_{\parallel}}$, which means their relative position depends on $\mu$. The exact bias depends on the details, mainly the location of the spurious peak, and to a certain extent the fit method, but it is likely insignificant for BOSS QSOs.

4.2 Angular component

Most of the statistical weight in BOSS data occurs around $z=2.2$. This can be seen in figs. 6 and 7, where we calculate the relative contribution of pixel-pairs at different redshifts to the autocorrelation function. The baryon peak scale of $100\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$ corresponds to about $1.5\text{\,}\mathrm{\SIUnitSymbolDegree}$ at this redshift. As a consequence of the rather narrow peak, having pairs at different redshifts might not be enough to completely erase transverse correlations that may exist at the \sim$$1\text{\,}\mathrm{\SIUnitSymbolDegree} scale in the ISM. The sky density of galaxies and QSOs in SDSS is insufficient to create ISM stacks at every position, and quantify this. Instead, we can use extinction maps as a proxy to ISM absorption lines, assuming the absorption lines are correlated with dust column (which we know is only true to first order, e.g., Baron et al. 2015a).

We calculate the angular autocorrelation of extinction in two ways: In Method 1, we use the foreground extinction through every QSO sightline using the maps of Schlegel et al. (1998). We compute the autocorrelation in angular bins for all possible pairing within a maximum separation of ${\mathord{\sim}\mkern 1.0mu}$3\text{,}\mathrm{\SIUnitSymbolDegree}$$. The resulting angular autocorrelation can be seen in fig. 8. Clearly, there is no preferred scale, as previously shown in Figure 9 of Schlegel et al. (1998).

In Method 2, we use the more recent extinction map from the Planck collaboration (Planck Collaboration et al., 2016). We calculate the angular autocorrelation for different regions of the sky, within the SDSS field. The result, as seen in fig. 9, behaves similarly to the one obtained with Method 1. While extinction varies greatly in magnitude at different galactic latitudes, the autocorrelation function in all fields has a smooth slope which decreases with angular separation. It appears there is no preferred scale in the dust column density. We therefore do not expect any significant transverse spurious correlation from the ISM. However, this is only a first order estimation, since ISM lines do not correlate perfectly with dust column density (Baron et al., 2015a).

4.3 Autocorrelation of stacked ISM spectra

Since the correlation computation is multi-staged, complex, and very non-linear, we proceed with a data driven approach to quantifying the ISM absorption effect on the autocorrelation estimator. Optimally we would want to derive an ISM spectrum along every sightline, however, since most of our signal in the blue is from QSOs, the source density in BOSS is insufficient. Instead, since to first order all the absorption lines correlate with extinction, we can make a reasonable estimate of the effect of ISM absorption, by using extinction as a proxy for line strength.

For every QSO sightline we use the dust extinction (from the maps of Schlegel et al. 1998) to substitute the actual spectrum and Lyman\Hyphdash*$\upalpha$ forest with an ISM spectrum stacked from sightlines with similar extinctions. The baryon peak pipeline then proceeds normally, removes the mean transmittance as before, and only features that are correlated with extinction remain.

We then compute the autocorrelation function in bins of comoving separation. We discard pixels pairs with $z<2.1$ because the ISM spectra are dominated by noise in this range. We plot the ISM autocorrelation estimator in a similar way to the Lyman\Hyphdash*$\upalpha$ forest, in figs. 10 and 11.

From this calculation we find the following. First, the amplitude of the signal from the ISM is about an order of magnitude weaker than the cosmological component. Secondly, the amplitude of the correlation drops by a factor of about 1.5 with ${r_{\perp}}$, with no preferred scale, as obtained before from dust maps (fig. 8 and fig. 9). Because we use stacks based on extinction, we cannot see changes along ${r_{\perp}}$ that would occur due to potential clustering of lines, as measured by Baron et al. (2015a), hence the straight lines. Instead we are mostly sensitive to correlations between different absorption lines that appear in ${r_{\parallel}}$. The Ca ii H&K doublet mimics a $27\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$ separation, is probably responsible for the positive correlation around this value of ${r_{\parallel}}$, and we confirm this in the following section. These results are consistent with our indirect analysis in section 4.1 and section 4.2.

4.4 Autocorrelation in distance bins

In order to study the effect of specific lines, most notably Ca ii H&K, and trace the source of the stripes seen in fig. 10, we divided the data into bins of parallel comoving distance (effectively redshift or wavelength bins).

Figure 12 shows how the Ca ii H&K lines create spurious correlations with each other as well as with other weaker features in the ISM spectrum, including a feature at ${\mathord{\sim}\mkern 1.0mu}$3720\text{,}\mathrm{\SIUnitSymbolAngstrom}$$, which we barely detect, nor can we constrain its carrier (it could perhaps be atomic iron, see Morton 2003). These Ca ii H&K correlations contribute to the strongest stripes in fig. 10, near $30\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$, $140\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$, and $170\text{\,}\mathrm{{}^{-1}}\text{\,}\mathrm{Mpc}$. The total correlation contribution of the Ca ii H&K peaks, weighted according to their corresponding redshift bin, is about $6\text{\times}{10}^{-6}\text{\,}$, comparable to our estimate in section 4.2, and an order of magnitude smaller than the BAO signal. Since the largest impact is from the Ca ii H&K lines, our previous use of the contribution from these lines to the autocorrelation as proxies for spurious signal seems justified.

5 Impact on future Surveys

While we find that the ISM currently has a negligible impact on the cosmological signal recovered from BOSS QSOs, we would like to quantify whether future surveys, such as the Dark Energy Spectroscopic Instrument (DESI), need to take this contaminant into consideration. DESI will obtain many more spectra, with a resolving power similar to that of BOSS, as well as a comparable SNR, with the larger aperture telescope compensating for the fainter QSOs (DESI Collaboration et al., 2016).

We use the QSO density estimate from the DESI final design report (DESI Collaboration, 2016, Table 2.7) to simulate a random Lyman\Hyphdash*$\upalpha$ forest which resembles what we expect from DESI, but without any large scale correlations. We choose Lyman\Hyphdash*$\upalpha$ forests from random SDSS QSOs with a survey area of about $10\,000\text{\,}\mathrm{{}^{2}}\deg$. This value lies between the proposed descoped DESI survey ($9000\text{\,}\mathrm{{}^{2}}\deg$) and the full DESI survey ($14\,000\text{\,}\mathrm{{}^{2}}\deg$).

We begin by taking the expected DESI QSO density as a function of redshift and using it as a probability distribution function (PDF). We smooth the PDF using linear interpolation. We draw $672\,000\text{\,}$ redshift values using the smoothed PDF, and match every value with an SDSS QSO at a similar redshift. We take the pre-calculated $\delta_{F}(z)$ values of every SDSS QSO, and rescale them to the chosen redshift value. Since most SDSS Lyman\Hyphdash*$\upalpha$ forests will appear more than once in our sample, we multiply every selected QSO by a random sign, to avoid contaminating the autocorrelation function with positive contribution from same-object pairs. We then draw another random SDSS QSO, and use its sky position, with a small random displacement as the position of our new artificial DESI Lyman\Hyphdash*$\upalpha$ forest. We then compute the autocorrelation function of the new collection of Lyman\Hyphdash*$\upalpha$ forests.

While the large scale structure is removed by the randomization process, the resulting autocorrelation bins contain statistical noise which is similar to what we expect from the real DESI survey, under the assumption that the distributions of SNR and instrumental uncertainty values in the Lyman\Hyphdash*$\upalpha$ forest are not too different. We also obtain the accumulated statistical weight as a function of separation and redshift.

The total number of pixels summed and the total statistical weight are both about 10 times larger than in BOSS. The standard deviation of individual autocorrelation bins is ${\mathord{\sim}\mkern 1.0mu}$5\text{\times}{10}^{-6}\text{,}$$, which is the same order of magnitude as the correlation arising from the Ca ii H&K lines, which are representative of the maximal spurious signal. In the near-parallel angular wedge, where a few bins are summed together, the Ca ii H&K lines might have a noticeable contribution, but this happens at a separation a few times smaller than the baryon acoustic peak. Overall using our method we find no indication the ISM will significantly influence DESI results. However, we caution that this is under the first-order assumption we had to make that the ISM follows the dust column density and has no additional covariance.

While Ca ii H&K lines could simply be masked by future surveys, this does not resolve correlations arising from other lines or possible structures in the ISM which we cannot probe due to our limited signal in the bluest part of the spectra. A better method to treat these biases would be to empirically derive the ISM spectrum around every QSO sightline, by stacking nearby lower-$z$ spectra in the observer frame, and removing the ISM contribution from every spectrum. As a consequence, Lyman\Hyphdash*$\upalpha$ cosmology would indirectly benefit from a large density of other unrelated sources.

6 Conclusions

By building an independent (though similar) pipeline we reproduce the detection of the BAO signal in the Lyman\Hyphdash*$\upalpha$ forest of BOSS QSOs. We analyze in different ways the possible biases that this measurement may suffer from, due to interloping lines from the MW ISM. The spurious signal could be produced by correlations between lines, i.e., in the parallel direction, or by angular correlations on the sky, affecting the transverse direction.

We find that the imprint of ISM absorption is too small to significantly affect current BAO results from the 3-D Lyman\Hyphdash*$\upalpha$ forest. As future surveys such as DESI contain more data with more precision in BAO, it might be preferable to explicitly correct for these effects.

While our current analysis suffers from a lack of signal in the blue wavelength range of interest, these future surveys should have a better handle on the ISM, allowing them, via the methods presented here, to better measure and correct for biases due to the MW ISM.

Acknowledgements

We would like to thank the 2nd referee A. Font-Ribera, and K.G. Lee for useful comments on this manuscript, and Roman Garnett for providing us early access to the DLA catalogue from Garnett et al. (2016). We further thank Julien Guy for his advice and for making numerical autocorrelation results available to us for comparison. This research was supported by Grant No. 2014413 from the United States-Israel Binational Science Foundation (BSF). 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 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 Astrofísica 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.