A New Precision Measurement of the Small-scale Line-of-sight Power Spectrum of the Lyα Forest

Our measurement of the power spectrum was performed using 38 high-resolution quasar spectra (see Table 1) from Dall'Aglio et al. (2008) observed with the Ultraviolet and Visual Echelle Spectrograph (UVES; Dekker et al. 2000) at the Very Large Telescope (VLT), and 36 spectra (see Table 2) from the Keck Observatory Database of Ionized Absorption toward Quasars (KODIAQ) project (Lehner et al. 2014) observed with the High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) at Keck. For the latter, we used the highest S/N (signal-to-noise ratio) part of DR1 (O'Meara et al. 2015) and additional data beyond DR1 (mostly early reductions of objects in DR2; O'Meara et al. 2017) reduced in the same way. Reduced and continuum fitted spectra of all UVES and KODIAQ DR1 data used here are available in the igmspec package (Prochaska 2017). KODIAQ DR2 data will be available in future igmspec releases and on the KODIAQ database web page.⁷

Table 1.  UVES Spectra from Dall'Aglio et al. (2008) Used for Our Analysis

Object	zQSO	Median S/N per $6\,\mathrm{km}\,{{\rm{s}}}^{-1}$
HE1341−1020	2.137	58.1
Q0122−380	2.192	56.4
PKS1448−232	2.222	57.4
PKS0237−23	2.224	102.4
HE0001−2340	2.278	65.9
Q0109−3518	2.406	70.0
HE1122−1648	2.407	171.6
HE2217−2818	2.414	93.6
Q0329−385	2.437	58.4
HE1158−1843	2.459	66.7
Q2206−1958	2.567	74.5
Q1232+0815	2.575	45.8
HE1347−2457	2.615	62.0
HS1140+2711	2.628	88.9
Q0453−423	2.663	77.6
PKS0329−255	2.705	48.0
Q1151+068	2.758	49.1
Q0002−422	2.768	75.0
HE0151−4326	2.787	98.1
Q0913+0715	2.788	54.4
Q1409+095	2.843	24.7
HE2347−4342	2.886	152.3
Q1223+178	2.955	33.4
Q0216+08	2.996	36.8
HE2243−6031	3.011	118.8
CTQ247	3.026	69.1
HE0940−1050	3.089	69.6
Q0420−388	3.120	116.2
CTQ460	3.141	40.9
Q2139−4434	3.208	31.2
Q0347−3819	3.229	83.9
PKS2126−158	3.285	63.6
Q1209+0919	3.291	30.2
Q0055−269	3.665	75.7
Q1249−0159	3.668	69.7
Q1621−0042	3.708	77.7
Q1317−0507	3.719	42.0
PKS2000−330	3.786	150.8

Download table as:  ASCIITypeset image

Table 2.  HIRES Spectra from KODIAQ (O'Meara et al. 2015) Used for Our Analysis

Object	zQSO	Median S/N per $6\,\mathrm{km}\,{{\rm{s}}}^{-1}$
J122824+312837	2.200	87.3
J110610+640009	2.203	58.5
J162645+642655	2.320	103.7
J141906+592312	2.321	36.7
J005814+011530c	2.495	36.2
J162548+264658c	2.518	43.9
J121117+042222	2.526	33.6
J101723−204658	2.545	70.3
J234628+124859	2.573	75.1
J101155+294141a	2.620	129.9
J082107+310751	2.625	64.0
J121930+494052	2.633	90.3
J143500+535953	2.635	65.0
J144453+291905	2.669	133.7
J081240+320808	2.712	48.8
J014516−094517A	2.730	76.8
J170100+641209c	2.735	81.8
J155152+191104	2.830	30.2
J012156+144820	2.870	54.5
Q0805+046b	2.877	26.8
J143316+313126	2.940	53.8
J134544+262506	2.941	34.7
J073621+651313	3.038	25.7
J194455+770552a	3.051	30.4
J120917+113830	3.105	31.4
J114308+345222a	3.146	31.9
J102009+104002	3.168	35.9
J1201+0116b	3.233	30.1
J080117+521034a	3.236	43.2
J095852+120245a	3.298	44.8
J025905+001126	3.365	26.3
Q2355+0108b	3.400	58.3
J173352+540030	3.425	57.3
J144516+095836a	3.530	24.6
J142438+225600a	3.630	29.3
J193957−100241	3.787	65.5

Notes.

^aObjects are part of DR2, but a pre-DR2 reduction has been used. ^bObjects are not part of DR1 or DR2, but reduced in the same way. ^cObjects are part of DR1, but a pre-DR1 reduction has been used.

Download table as:  ASCIITypeset image

Most of the UVES spectra have full coverage between the atmospheric cutoff at ${\lambda }_{\mathrm{obs}}\sim 3100\,\mathring{\rm A}$ and ${\lambda }_{\mathrm{obs}}\sim 1\,\mu {\rm{m}}$. This allows us to use a large range of spectrum redward of the Lyα forest to search for metal lines and enables us to search for Lyman limit systems (LLSs) using higher-order Lyman series transitions, in many cases even exploiting coverage of the Lyman limit.

For the KODIAQ data, the typical red spectral coverage ends at ${\lambda }_{\mathrm{obs}}\sim 6000\,\mathring{\rm A}$, while the blue spectral cutoff is comparable to UVES ${\lambda }_{\mathrm{obs}}\sim 3100\,\mathring{\rm A}$. For a few cases in both data sets, however, even the Lyα forest was not fully covered.

The objects used in our analysis were chosen to have a median S/N > 20 per $6\,\mathrm{km}\,{{\rm{s}}}^{-1}$ interval inside the Lyα forest region covered by the spectra. We also chose to omit spectra with known broad absorption lines (BALs). Finally, we omitted sight lines with 3.0 < z < 3.5 that were color selected to avoid potential biases due to their increased abundance of LLSs (Worseck & Prochaska 2011).

The distribution of S/N for the data set used in our analysis is shown in Figure 1. Many objects have a much larger S/N than our cut (up to about a median S/N of 150 per 6 km s⁻¹ inside the Lyα forest in a sight line). This very high S/N enables us to perform the measurement without strong systematics due to the limited accuracy of our noise model affecting our measurement at the scales of interest.

Zoom In Zoom Out Reset image size

The nominal FWHM spectral resolution of the data varies between $3.1$ and $6.3\,\mathrm{km}\,{{\rm{s}}}^{-1}$ with a typical value around $6\,\mathrm{km}\,{{\rm{s}}}^{-1}$. Therefore all Lyα absorption lines are resolved, and we expect thermal broadening and broadening due to pressure smoothing to be the effects determining the smoothness of lines, and not smoothing due to finite spectroscopic resolution (see, e.g., Bolton et al. 2014 for a measurement of line width for thermally broadened lines).

The data were already reduced and continuum fit, and we briefly summarize the details here. The UVES data were fit with both a global power law in nonabsorbed regions and local cubic splines fitted automatically with spline point separations depending on the continuum slope. Systematic biases in this technique were estimated and corrected using a Monte Carlo analysis on mock data by Dall'Aglio et al. (2008). The HIRES continua were hand-fitted by John O'Meara one echelle order at a time by placing Legendre polynomial anchor points at nonabsorbed positions (for estimates of continuum fitting errors, see Kirkman et al. 2005; Faucher-Giguère et al. 2008b). Afterward, a fourth- to 12th-order polynomial was fit through the anchors. Further details about the reduction and continuum fitting techniques can be found in the respective data papers (Dall'Aglio et al. 2008; O'Meara et al. 2015).

We will use this high-resolution data set to measure the small-scale power spectrum of the Lyα forest in redshift bins of size Δz = 0.2 with central redshifts between $\bar{z}\,=\,1.8$ and $\bar{z}=3.4$. Each of the bins contains at least eight quasar spectra, with the majority (all but the low- and high-redshift edges of our sample) containing more than 14 quasar spectra. The data sets also contain eight spectra that cover higher redshifts that we did not analyze because of the small amount of data available (two spectra cover z ≳ 4.0; the rest only cover the forest for z < 3.6, just half a bin further than our analysis).

In this work and our companion paper (M. Walther et al. 2017b, in preparation), we compare our power spectrum measurements to measurements from data with lower spectral resolution from BOSS (PD+13) and from the XQ-100 survey (Iršič et al. 2017a) based on X-SHOOTER data, and we also conduct joint model fits. To facilitate this comparison, we use the same binning in redshift from $\bar{z}=2.2$ to $\bar{z}=3.4$.

To prepare the data for the power spectrum computation, we restrict our attention to the rest-frame wavelength range of $1050\,\mathring{\rm A} \lt {\lambda }_{r}\lt 1180\,\mathring{\rm A}$. This was done to exclude the Lyα proximity zone, accounting for possibly large redshift errors, as well as to exclude the Lyβ and O vi λ 1035 emission lines and possible blueshifted absorption from these as well as increased continuum fitting errors close to emission lines. This is the same range used in PD+13 and is considered a conservative choice for the Lyα forest region.

We masked parts of the spectrum to reduce contaminations due to low-quality data, high-column-density absorbers such as DLAs, and metal lines. If a pixel is already masked during the reduction (due to, e.g., cosmic rays or gaps in spectral coverage), it stays masked. According to McDonald et al. (2005), excluding DLAs and LLSs from power spectrum calculations only changes the power by <2% on the scales of interest for our analysis. We nevertheless excluded absorbers with clearly visible damping wings, that is, DLAs and super-LLSs, from our spectra to make sure those do not influence the result. For this we masked the core and wings of these strong absorbers by eye until the wings are below the noise level of the spectra. In most cases, this leads to the exclusion of big continuous spectral regions at the boundary of a redshift bin or the whole data set of this object falling inside a redshift bin. Where a DLA mask removed only the bin center, we used the longer of the two remaining spectral regions at the bin boundaries to compute the power. Therefore, spectra with DLAs are only shorter, but with no additional gaps in the data. Finally, we removed spectra that were shorter than 10% of a redshift bin corresponding to a minimum length of Δz ∼ 0.02 (or $3500\,\mathrm{km}\,{{\rm{s}}}^{-1}$) to avoid noisy contributions from short spectra that consist of only a few absorption lines.

Metal lines in the Lyα forest are expected to increase the power primarily on small scales ($k\gtrsim 0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$, see McDonald et al. 2000; Lidz et al. 2010), due to their narrower widths compared to Lyα. But due to correlations induced by the rest-frame velocity separations of different transitions (e.g., between Si iii and H i-Lyα or between the C iv doublet lines), there is contamination on larger scales as well (Croft et al. 1999; McDonald et al. 2006).

To reduce the impact of metal absorption inside our sight lines, we masked metal lines in the forest region. We first identified metal absorption lines in the Lyα forest by having two of the authors (H. Hiss and M. Walther) visually inspect the spectra and mask metal contamination in several ways. First, we look for DLAs inside the forest and mask all strong metal absorption lines corresponding to the DLA redshifts. For this we used all metal transitions in Table 3 and masked a region of $60\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in each direction around each identified absorber. Then we search for absorption from common doublet transitions (Si iv, C iv, Mg ii, Al iii, Fe ii) redward of the forest where the spectrum is mostly clean and mask all associated metal lines analogous to our procedure for DLAs (using the same metal catalog and mask width), until there are no doublet features left redward of the forest. We also masked out a region of $200\,\mathrm{km}\,{{\rm{s}}}^{-1}$ in each direction around redshift zero to get rid of metal contamination from the Milky Way (in this case the only relevant transition is Ca ii for $2.33\lt {z}_{\mathrm{QSO}}\lt 2.78$).

Table 3.  List of Metal Transitions That Were Masked

Ion	${\lambda }_{\mathrm{rest}}/\mathring{\rm A}$	Metal Transition	${\lambda }_{\mathrm{rest}}/\mathring{\rm A}$
O via	1031.9261	Si iva	1402.770
C ii	1036.3367	Si ii	1526.7066
O vi	1037.6167	C iva	1548.195
N ii	1083.990	C iva	1550.770
Fe iii	1122.526	Fe ii	1608.4511
Fe ii	1144.9379	Al ii	1670.7874
Si ii	1190.4158	Al iii	1854.7164
Si ii	1193.2897	Al iii	1862.7895
N i	1200.7098	Fe ii	2344.214
Si iiia	1206.500	Fe ii	2374.4612
N v	1238.821	Fe ii	2382.765
N v	1242.804	Fe ii	2586.6500
Si iia	1260.4221	Fe ii	2600.1729
O i	1302.1685	Mg ii	2796.352
Si ii	1304.3702	Mg ii	2803.531
C ii	1334.5323	Mg i	2852.9642
C ii*	1335.7077	Ca i	3934.777
Si iva	1393.755	Ca i	3969.591

Note.

^aStrongest transitions, therefore used for all of the masking techniques.

Download table as:  ASCIITypeset image

Additionally, we used an automated partial LLS (pLLS) finder written by John O'Meara to identify strong absorption systems. This finder works by searching for pixels with zero flux (within some threshold) at the corresponding positions of Lyα, β, γ, and higher Lyman transitions if available and grouping them into systems of LLS candidates. The candidates were then visually inspected by one of the authors (John O'Meara) and compared to theoretical line profiles of absorbers with $\mathrm{log}({N}_{{\rm{H}}{\rm{I}}})=15,16,17$ in Lyα to Lyγ, and systems that appeared consistent were positively identified as pLLSs. For these systems, associated metal absorbers from a reduced line list (see Table 3 absorbers marked with "a") were masked with the same velocity window size as above. Note that the hydrogen absorption arising from the pLLSs identified in this way was not masked regardless of their N_{H i}.

Lastly, we perform a line-fitting analysis (see Hiss et al. 2017) using a semiautomatic wrapper around VPFIT (Carswell & Webb 2014) on the same set of Lyα spectra. The result of this is a distribution of line widths b and column density ${N}_{{\rm{H}}{\rm{I}}}$ for all fitted lines assuming that all absorption is due to hydrogen. For H i gas at a given column density, lines are broadened both thermally (due to finite pressure broadening and instantaneous temperature Doppler broadening) and hydrodynamically by local gas motions. There is a minimum broadening ${b}_{\mathrm{cut}}({N}_{{\rm{H}}{\rm{I}}})$ populated by absorbers with zero line-of-sight peculiar velocity that are purely thermally broadened (see, e.g., Schaye et al. 2000; Hiss et al. 2017, for more details). Therefore, $b({N}_{{\rm{H}}{\rm{I}}})$ should cut off for $b\lt {b}_{\mathrm{cut}}({N}_{{\rm{H}}{\rm{I}}})$, and all remaining lines narrower than this cutoff can be attributed to fitting artifacts at the edges of strong absorbers, noise fluctuations, or narrow metal absorption lines. This cutoff is fit using an iterative procedure similar to the one used in Rudie et al. (2012). We identified all of the lines in previously unmasked spectral regions fulfilling $b\lt 11\,\mathrm{km}\,{{\rm{s}}}^{-1}\,{({N}_{{\rm{H}}{\rm{I}}}/{10}^{12.95}{\mathrm{cm}}^{-2})}^{0.15}$ as these are narrower than the thermal cutoff (see Hiss et al. 2017 for details). For each of these lines (that are in the Lyα forest region), we checked if they could be identified using a second metal transition clearly lining up at the same redshift. In practice, these identifications were most easily made when both metal transitions form a doublet. Given a positive identification, masking was then performed as for metal absorbers that we identified redward of the Lyα forest, that is, by determining the redshift of the absorber and using the full list of metal lines in Table 3.

Neither of these techniques produces a fully metal-free Lyα forest, and we briefly elaborate on these limitations. The search for metals redward of the forest only finds transitions if both doublet counterparts are visible at redder wavelength than the Lyα emission line. If the doublet counterpart of a line falls into a spectral region that was not covered or is blended with a different line (from either a different absorber redshift or, e.g., a telluric absorption feature), the contaminating system will often go undetected. For example, for a C ivλ1548/1550 system at z > 2.5 that might contaminate the forest, the Mg ii λ2796/2803 doublet would land at ${\lambda }_{\mathrm{obs}}\gt 1\,\mu {\rm{m}}$ and hence outside of the spectral coverage of our spectra, so only if the absorber shows Fe ii or Al iii doublets would this C iv absorber be masked. On the other hand, especially for the higher redshift bins, a large fraction of the spectral range redward of the forest is contaminated by telluric absorption, making doublet identification very challenging. The largest problem for this method is therefore limited usable spectral coverage in the red.

For the automated pLLS finder, at least Lyα to Lyγ need to be detectable to identify a pLLS. This leads to a minimal redshift of z ≳ 2.1 for absorbers that can be identified as they need to be at observed wavelengths higher than the atmospheric cutoff ${\lambda }_{\mathrm{obs}}\approx 3000\,\mathring{\rm A}$. For reduced spectral coverage in the blue, this minimal redshift is correspondingly higher.

The last of our metal masking procedures only recognizes metal absorption lines significantly narrower than the cutoff in the $b({N}_{{\rm{H}}{\rm{I}}})$ distribution and requires a second metal transition for identification. Therefore, singlet lines are not masked by this technique unless another transition from a different metal species clearly lines up with them. Also, metal absorbers were not removed if all components are broadened above the $b({N}_{{\rm{H}}{\rm{I}}})$ threshold adopted, and they are therefore not recognized as metal candidates.

In Figure 2 we show an example to illustrate how much of a typical spectrum is masked. The black line shows the nonmasked part of the spectrum, while the red line shows parts that are masked by possible metal contamination. One can see that many narrow lines in the spectrum are masked, but also that many regions of the forest coincidentally overlap with positions where metal lines associated with our identified absorbers could lie, but that do not actually contain any visible metal absorption.

Zoom In Zoom Out Reset image size

In Figure 3 we present an overview of the data set. For each quasar spectrum, the emission redshift is shown as well as the full coverage of the spectrum after masking preexisting gaps in the data and DLAs (which therefore appear as white gaps). The remaining spectrum is divided into the data used (dark colors), masked data due to possible metal contamination (bright colors), and spectra not used due to failing our requirement of spectral extent (yellow). We also show the usable Lyα forest path length after applying the full masking procedure compared to the available path length when not masking metals in Figure 4. Up to ∼40% of the forest gets removed by applying our metal masking procedure, strongly reducing our data set size and therefore the precision of our measurement results. Although not perfect, our procedure significantly reduces the metal line contamination in our spectra, which decreases the amount of small-scale power compared with an unmasked data set, as we will see in Section 4.1. However, this masking also changes the power spectrum due to the application of a complex window function in configuration space. To correct for this effect, we forward model the masking (see Section 3.4).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In principle, the estimation of quasar continua in the data is subject to errors as well. We perform our power spectrum measurement on the flux contrast

$$\begin{eqnarray}&&{\delta }_{F}=\displaystyle \frac{F-\bar{F}}{\bar{F}},\end{eqnarray} \tag{ 2 }$$

with $\bar{F}$ being the mean transmission of the Lyα forest. Because of this, any global misplacement of the continuum will be divided out as long as the mean transmission measurement we divide by is measured on the same spectrum. In addition, incorrect placement of the continuum could lead to gradients or wiggles in the data that could source additional large-scale power (on scales of $k\lt 0.001\,{\rm{s}}\,{\mathrm{km}}^{-1}$). However, this will not strongly impact the small-scale power measurement we want to perform in this work⁸ (Lee et al. 2012).

Nevertheless, we perform a mean flux regulation on the data set using the technique of Lee et al. (2012), which enables us to easily divide out the mean transmission as well as possible gradients in the continuum fits. For this, the transmission of the Lyα forest region (excluding masked pixels as well as possible proximity regions, as discussed earlier) of each quasar sight line is first fit by a linear relation f(z) times the mean transmission function $\bar{F}=\exp (-{\tau }_{\mathrm{eff}})$ with

$$\begin{eqnarray}&&{\tau }_{\mathrm{eff}}=C+{\tau }_{0}{\left(\displaystyle \frac{1+z}{1+{z}_{0}}\right)}^{\beta }\end{eqnarray} \tag{ 3 }$$

following the functional form and parameters from Becker et al. (2013). Afterwards, f(z) is divided out so that the mean flux evolution of each spectrum follows the same relationship.

The flux contrast ${\delta }_{F}$ is now easily obtained by using Equation (2) on mean flux regulated spectra using Equation (3) for the mean flux evolution instead of dividing each spectrum by a mean transmission estimate for this spectrum. This allows us to divide out the mean flux at each pixel analytically based on the fit of the mean flux evolution across each spectrum. While this in principle leads to a reduced large-scale power (Lee et al. 2012), we are not measuring the affected scales in this work.

The resulting spectra are finally divided into our redshift bins of Δz_chunk = 0.2 (with the first bin starting at z = 1.7) to increase redshift resolution to a level where we could closely monitor thermal evolution. We will henceforth call these pieces of spectra "spectral chunks."

We calculate the power spectrum from the flux contrast δ_F defined in Equation (2). As our spectra are not periodic and are not regularly sampled because of masking, we Fourier transform δ_F using a Lomb–Scargle periodogram (Lomb 1976; Scargle 1982), allowing us to compute the raw power ${P}_{\mathrm{raw}}(k)$ of each spectrum for a linearly spaced set of modes from the fundamental mode given by the length of the spectral chunk to the Nyquist limit. We subtract the noise power P_noise(k) from this raw power and divide the difference by the window function W_R resulting from finite spectroscopic resolution R, following the fast Fourier transform (FFT) method described in PD+13:

$$\begin{eqnarray}&&{P}_{\mathrm{data}}(k)=\left\langle \displaystyle \frac{{P}_{\mathrm{raw}}(k)-{P}_{\mathrm{noise}}(k)}{{W}_{R}^{2}(k,R,{\rm{\Delta }}v)}\right\rangle ,\end{eqnarray} \tag{ 4 }$$

with

$$\begin{eqnarray}&&{W}_{R}(k,R,{\rm{\Delta }}v)=\exp \left(-\displaystyle \frac{1}{2}{({kR})}^{2}\right)\displaystyle \frac{\sin (k{\rm{\Delta }}v/2)}{(k{\rm{\Delta }}v/2)},\end{eqnarray} \tag{ 5 }$$

and Δv refers to the pixel scale of the spectra in velocity units.

The noise power P_noise is measured by creating 100 realizations of Gaussian random noise generated from the 1σ error vector of each quasar spectrum. The resolution assumed for the window function correction is the nominal slit resolution of the spectrograph, which is different from the actual spectral resolution of the data, which also depends on the seeing of the observations. For the typical resolution in our data set (a resolving power of 50,000 or equivalently $R\,=2.55\,\mathrm{km}\,{{\rm{s}}}^{-1}$), we get ${W}_{R}^{2}(0.1\,{\rm{s}}\,{\mathrm{km}}^{-1})\sim 0.94$. Given this small correction, even ∼20% error in our knowledge of the resolution (due to the unknown seeing) would only lead to a small ≲4% correction of the power at $k\lt 0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$ (for further discussion of this point, see Section 4.2 and Appendix A). Therefore, the resolution uncertainty does not significantly affect our measurement of the k modes we consider.

We chose⁹ to use the same logarithmically spaced k bins as in McDonald et al. (2000) for our analysis. Note that this is different from the linear spacing adopted by PD+13 and Iršič et al. (2017b).

We adopt the same Fourier normalization convention as the BOSS measurement (PD+13), such that the variance in flux contrast is ${\sigma }_{{\delta }_{F}}^{2}={\sigma }_{F}^{2}/{\bar{F}}^{2}={\int }_{-\infty }^{\infty }P(k){dk}/(2\pi )$. Note that this differs from the conventions used by some older high-resolution measurements (see Appendix C), leading to additional normalization factors being needed when comparing to those results.

As described in Section 2.2, we masked out parts of the data in real space due to metal contamination. This is a different approach than the one used by PD+13, who estimated the metal power from transitions with $\lambda \gt {\lambda }_{\mathrm{Ly}\alpha ,\mathrm{rest}}$ by measuring the power redward of the forest in lower-redshift quasar spectra. The power measured in this way was then subtracted from the measurement, but this method can never account for transitions with $\lambda \lt {\lambda }_{\mathrm{Ly}\alpha ,\mathrm{rest}}$, which always fall into the forest (e.g., the Si iii λ1206 line that leads to the Si iii correlation feature in those measurements; see McDonald et al. 2006).

Masking spectra in configuration space with a window function W_m leads to the measured P_masked for each spectrum being effectively a convolution of the true power P_true with the square of the Fourier transform ${\tilde{W}}_{m}$ of this window function:

$$\begin{eqnarray}&&{P}_{\mathrm{masked}}={P}_{\mathrm{true}}\,\ast \,{{\tilde{W}}_{m}}^{2}.\end{eqnarray} \tag{ 6 }$$

Determining the true power thus requires deconvolving the window function or adopting a different technique to measure the power taking into account the windowing (like the minimum variance estimator used in McDonald et al. 2006). In this work, we opt to use a simpler approach by generating forward models of the data with and without the windowing applied to be able to determine the effect this window function has on our data. In the end, we correct our data for those effects by dividing the measurement by a correction function based on these models. In the remainder of this section, we describe how our forward models are generated, how we estimate the data covariance matrix based on those models, and how we fit the data using those models.

We use simulations of the Lyα forest for two reasons. First, we need to simulate the effect of noise, resolution, and the window function due to masking on the power spectrum. Second, we want to connect the information encoded in the power spectrum to a thermal history of the IGM (parameterized by the temperature at mean density T₀, the slope of the temperature–density relation γ, and the pressure smoothing scale λ_P). For this purpose, one generally needs to run hydrodynamical simulations with different thermal histories. However, these are computationally expensive, and at least for the first point we only need a model that is flexible enough to provide a good fit to the observed power spectra, but it need not necessarily provide the correct thermal parameters. Because DM-only simulations of the Lyα forest are more flexible and computationally inexpensive to generate, for all the forward modeling in this paper, we use approximate DM-only simulations (a single box with different thermal parameters generated in postprocessing) with a seminumerical approach to paint on the thermal state of the IGM (Hui & Gnedin 1997; Croft et al. 1998; Gnedin & Hui 1998; Meiksin & White 2001; Gnedin et al. 2003; Rorai et al. 2013) This fast simulation scheme allows us to generate a grid of ∼500 combinations of thermal parameters in a reasonable amount of time.

These DM-based simulations are, however, not sufficiently accurate to infer the thermal state of the IGM, as they produce significant biases in thermal parameters when fitted to mock data based on hydrodynamical simulations (see, e.g., Sorini et al. 2016 for a detailed comparison between hydrodynamical and DM-only simulations). They do, however, provide a good fit to the data (compare lines and same-color error bars in Figure 6, which we will discuss in more detail later), and we therefore use them to correct for the window function in our measurement as well as test our analysis procedure and our fast power spectrum emulator (interpolation scheme; see Section 3.6). For inference of IGM thermal parameters, a grid of hydrodynamical simulations will be used in a companion paper (M. Walther et al. 2017b, in preparation).

For the DM-only simulations, we use an updated version of the TreePM code described in White (2002) to simulate the evolution of DM particles from initial conditions at z = 150 up to z = 1.8. We use a simulation with L_box = 30 h⁻¹ Mpc and 2048³ particles with a Plummer equivalent smoothing of 1.2 h⁻¹ Kpc based on a Planck Collaboration et al. (2014) cosmology with Ω_m = 0.30851, Ω_bh² = 0.022161, h = 0.6777, n_s = 0.9611, and σ₈ = 0.8288. We do not include uncertainties in cosmological parameters in our models because cosmic microwave background (CMB) measurement errors on these parameters (Hinshaw et al. 2013; Planck Collaboration et al. 2016b) are much smaller than current constraints on thermal parameters.

The model is generated for snapshots with z between 1.8 and 3.4 with a separation of 0.2, the same as our power spectrum measurements, and provides a DM density and velocity field. However, the relevant quantities for absorption in the IGM are baryonic density and temperature fields. The results of previous hydrodynamic simulations suggested a computation of the relevant baryonic fields using scaling relations on the DM quantities (Hui & Gnedin 1997; Gnedin & Hui 1998; Gnedin et al. 2003). This basically consists of smoothing the DM density field with a Gaussian kernel to mimic pressure support as well as rescaling the densities to get the right Ω_b. Temperatures are then introduced by applying the power-law temperature–density relation from Equation (1) to the density field (see Rorai et al. 2013, Section 2.2 for the exact procedure).

Given that the UVB is not known perfectly, we created a sequence of models with different mean transmissions $\bar{F}$ spanning an ∼5σ range around the current observational constraints by Becker et al. (2013), Faucher-Giguère et al. (2008b), and Kirkman et al. (2005). This was done by rescaling the optical depths of the full set of skewers to match the desired $\bar{F}$. In total, we have a parameter grid of ∼500 different thermal parameter combinations (T₀, λ_P, γ), with each of those evaluated for five different values of $\bar{F}$.

As the observed spectra are much longer than the simulation box, we first divided each spectral chunk of Δz = 0.2 into regions smaller than our box and assigned a random simulated model skewer to each of the pieces. The skewers were than assumed to fall on the respective position of the spectrum and truncated to have the same length as the spectral chunk. A model of a single real spectrum therefore consists of four to eight simulation skewers. We generated ∼15,000 skewers of the same length as the simulation box for each parameter combination, so we have the equivalent of ≳1900 mock spectra (or equivalently ≳100 times our whole data set) available at each set of parameters. This step is required to enable us to add the wavelength-dependent noise and masks to the models.

While the overall mean flux of the box is a free parameter of our models, we renormalize the flux in each pixel of the skewer again to account for the slight redshift evolution of the mean flux along the skewer with respect to the mean flux of the simulation snapshot. To do this, the fluxes are converted back to optical depth by $\tau =-\mathrm{log}(F)$. We then use the best-fit relation to the τ evolution by Becker et al. (2013, see Equation (3)) to compute the fractional change in τ_eff between box redshift and the individual pixel redshifts and rescaled τ at each pixel with the corresponding value. After this we convert back to fluxes. The same mean flux evolution function is then later taken out in the power spectrum computation when we, analogous to our procedure on the data, divide by the mean flux in order to compute the flux contrast. We then convolved these spectra with the respective instrument resolution and interpolated them onto the same wavelength grid as the observed data. The results of these steps are mock Lyα forest spectra with the same noise properties, spectral coverage, and masking as our data. We henceforth call these the "forward models."

We also compute the power spectrum of the same number of model skewers without any noise, masking, or degradation of resolution, which we compare to the power spectrum of our mocks to validate our power spectrum pipeline and to determine the window function correction. We will henceforth refer to these as "perfect models."

To fit the BOSS data, there is no need to create full forward models, because all the imperfections our forward model approach treats are already accounted for by the BOSS power spectrum pipeline, and therefore we simply compare the BOSS data to the perfect models.

In addition to measurements of power spectra, our statistical analysis requires a covariance matrix. The full covariance matrix consists of ${n}_{\mathrm{bins}}\cdot ({n}_{\mathrm{bins}}+1)/2$ independent entries where n_bins ≃ 15 is the number of band powers of the wavenumber k. For a data set like ours that consists of only n_QSO ∼ 10 quasars, this estimate will be extremely noisy. McDonald et al. (2000) did tests on bootstrapped covariance matrices based on spectra that were subdivided into five spectral chunks (which then were treated as independent data in the same redshift bin) and concluded that there was still significant statistical uncertainties on the estimates of individual covariance matrix entries. Also, Iršič et al. (2017a) tested the bootstrap covariance estimation technique on models with different amounts of skewers. When using only 100 model skewers, which is comparable to the size of their data set, noise in their correlation matrix is still clearly visible. They also provide a correlation matrix of their measurement that looks similar to this model estimate. As our typical redshift bin has less data available than theirs and as our data are additionally masked for metal absorption, we do not believe that we can estimate reliable covariance matrices from our data (see Figure 5 and discussion below). To circumvent this problem, we use a hybrid approach (in a similar way as Lidz et al. 2006) and only measure the diagonal values ${\sigma }_{\mathrm{data}}^{2}$ of the covariance from the data itself while computing the off-diagonal correlations R_ij from simulations (see Section 3.3) for which we can obtain sufficient statistics:

$$\begin{eqnarray}&&{R}_{\mathrm{ij}}=\displaystyle \frac{\langle ({P}_{{\rm{m}},\mathrm{single}}({k}_{{\rm{i}}})-{P}_{{\rm{m}}}({k}_{{\rm{i}}}))({P}_{{\rm{m}},\mathrm{single}}({k}_{{\rm{j}}})-{P}_{{\rm{m}}}({k}_{{\rm{j}}}))\rangle }{{\sigma }_{{\rm{m}}}({k}_{{\rm{i}}}){\sigma }_{{\rm{m}}}({k}_{{\rm{j}}})},\end{eqnarray} \tag{ 7 }$$

where ${P}_{{\rm{m}},\mathrm{single}}$ is the power estimated from a single model skewer, P_m is its mean over all model skewers, σ_m its standard deviation, and the average $\langle ...\rangle$ is performed over all model skewers. The covariance matrix is then computed as

$$\begin{eqnarray}&&{C}_{\mathrm{ij}}={\sigma }_{\mathrm{data}}({k}_{{\rm{i}}}){\sigma }_{\mathrm{data}}({k}_{{\rm{j}}}){R}_{\mathrm{ij}}.\end{eqnarray} \tag{ 8 }$$

Zoom In Zoom Out Reset image size

To estimate the data uncertainties ${\sigma }_{\mathrm{data}}$, we use a bootstrapping technique to resample the data. For this, we draw 1000 random sets of quasars with replacement from those contributing to any given redshift bin, and for each we compute the power spectrum using Equation (4). For the correlation matrix R, we use the same mock spectra generated with our forward-modeling procedure (the "forward models" in Section 3.4) to determine the correlation matrix according to Equation (7). No bootstrapping was performed, since for the mocks we have access to different forward-modeling realizations of the same data set. Note that we have a correlation matrix for each model, because the shape of the power spectrum impacts the correlations, leading to stronger correlations for modes on scales smaller than the power spectrum cutoff.

In Figure 5 we show the correlation matrix for one of our redshift bins determined from the model with parameters closest to our best fit at this redshift. The typical correlation between neighboring bins is ∼15% and decreases strongly for bins farther apart. For comparison, we also compute the correlation matrix from our data set at this redshift as well. This correlation matrix is noisy and is not used anywhere else in our analysis but is shown here for the sake of illustration. However, one sees that it qualitatively agrees well with our simulated correlation matrix, validating our approach. We also show the correlations of a perfect model, and comparing to the forward model we can see that the additional masking added significant correlations on the 20% level between nonneighboring k bins.

We also tested our bootstrap estimation of the diagonal elements of the covariance using random simulated data sets with a size comparable to our measurement data set that were drawn from our full set of models without replacement. We found that the variance of the power spectrum determined from the ensemble of simulated data sets was in good agreement with the bootstrap estimate obtained from a single mock data set. Therefore we are confident that our bootstrap estimates of the diagonal elements of the covariance are converged and reflect the actual uncertainties of the power spectrum.

To be able to fit the power spectrum measurement, we need to create a model for the power in the full region of interest for thermal parameters (T₀, γ, λ_P) and the mean transmission of the forest $\bar{F}$. As simulations are relatively expensive, there is no way to run a full simulation for every combination of thermal parameters at which the power spectrum needs to be evaluated during the fitting process. Therefore we generate a grid of simulations with different thermal parameters (using our seminumerical DM-only simulations) and adopt a fast emulation procedure to compute the power spectrum for parameters between the grid points (in a way similar to what is used for the cosmic calibration framework by Heitmann et al. 2006, 2009, 2013; Habib et al. 2007). Following the procedure described in Section 3.3, we generated Lyα forest skewers for a grid of thermal parameters (T₀, γ, λ_P) and mean flux values $\bar{F}$. We computed the power spectrum from perfect skewers generated from these models, as well as for mock data run through the forward-modeling pipeline described in Section 3.4.

For both sets of power spectra (perfect and forward modeled), we perform a principal component analysis (PCA) of the full set of model power spectra. As a result, the power at any point $\theta =\{{T}_{0},\gamma ,{\lambda }_{P},\bar{F}\}$ in our model grid can be written as

$$\begin{eqnarray}&&P(k,\theta )=\displaystyle \sum _{i}{\omega }_{i}(\theta ){{\rm{\Phi }}}_{i}(k),\end{eqnarray} \tag{ 9 }$$

where both the eigenvectors ${{\rm{\Phi }}}_{i}(k)$ and the PCA weights ${\omega }_{i}(\theta )$ are a result of the PCA decomposition. Note that the PCA weights depend on the grid points in our model parameter space, and given a suitable interpolation scheme, they can be used to generate model predictions at any point in this space. Along these lines, and following Habib et al. (2007), the PCA weights ω_i are then used as an input to train a Gaussian-process interpolation scheme. From this Gaussian process, we can later evaluate new weights for parameter combinations that lie between our original grid points, allowing us to evaluate the power at any location. For more details on this approach, see a companion paper to this analysis (M. Walther et al. 2017a, in preparation) or, for example, Rorai et al. (2013) and Habib et al. (2007). The Gaussian-process interpolation uses a squared exponential kernel with smoothing lengths chosen to be larger than the separation between model grid points of the parameter space. To speed up computations, only the first nine principal components are actually used for the analysis. We found that the additional errors due to discarding the higher-order components is less than 1%.

In general, the same approach can be used to emulate models from hydrodynamical simulations, but this means that a separate hydrodynamical simulation must be run for each parameter combination in the grid, which is far more costly than the approach we chose for this work. Nevertheless, we will use an emulator based on full hydrodynamic simulations to infer thermal parameters from our power spectrum measurements in a companion paper (M. Walther et al. 2017b, in preparation).

To explore the parameter space and fit the measured data power spectra from both BOSS and high-resolution data sets, we use a Bayesian Markov chain Monte Carlo (MCMC) approach with a Gaussian multivariate likelihood:

$$\begin{eqnarray}{ \mathcal L } & \equiv & P(\mathrm{data}| \mathrm{model})\\ & \propto & \displaystyle \prod _{\mathrm{data}\,\mathrm{sets}}\displaystyle \frac{1}{\sqrt{\det (C)}}\exp \left(-\displaystyle \frac{{{\boldsymbol{\Delta }}}^{{\rm{T}}}{C}^{-1}{\boldsymbol{\Delta }}}{2}\right)\\ {\boldsymbol{\Delta }} & = & {{\boldsymbol{P}}}_{\mathrm{data}}-{{\boldsymbol{P}}}_{\mathrm{model}},\end{eqnarray} \tag{ 10 }$$

where ${{\boldsymbol{P}}}_{\mathrm{data}}$ and C are the power spectrum and covariance matrix obtained for any data set, and the product is over all data sets considered (i.e., BOSS and our high-resolution measurement, but including more than two data sets would just add factors to this product).

For the high-resolution data set, ${C}^{-1}$ is estimated using our hybrid approach (see Section 3.5) for each model on the parameter grid and therefore is a model-dependent quantity. Inside the likelihood computation, we use nearest-neighbor interpolation in model parameter space to obtain an inverse covariance estimate for each combination of model parameters (which is then just the matrix at the closest model grid point). For the covariance matrix of the BOSS data set, we use the tabulated values from PD+13, which are directly measured from the data and are therefore model independent. The model power ${\boldsymbol{P}}$_model is determined for any parameter combination using our emulation procedure. Note that we actually have two distinct emulators, one for the high-resolution data (using the full forward-modeling procedure) and one for the BOSS data (using our perfect simulation skewers).

It is well known that correlated absorption of hydrogen Lyα and Si iii at $1206.5\,\mathring{\rm A}$ leads to a bump in the Lyα forest flux correlation function $\xi ({\rm{\Delta }}v)$ at ${\rm{\Delta }}v=2271\,\mathrm{km}\,{{\rm{s}}}^{-1}$ (McDonald et al. 2006; PD+13).¹⁰ This imprints wiggles on the power spectrum with separation ${\rm{\Delta }}k=2\pi /{\rm{\Delta }}v=0.0028\,{\rm{s}}\,{\mathrm{km}}^{-1}$. Following McDonald et al. (2006), we model this contamination with a multiplicative correction to the power

$$\begin{eqnarray}&&P(k)=(1+{a}_{\mathrm{Si}\,{\rm{III}}}^{2}+2a\cos ({\rm{\Delta }}v\,k)){P}_{{\rm{H}}{\rm{I}}}(k),\end{eqnarray} \tag{ 11 }$$

with ${a}_{\mathrm{Si}{\rm{III}}}$ being a free nuisance parameter for the strength of the correlation. In previous works, this was typically expressed as ${a}_{\mathrm{Si}{\rm{III}}}={f}_{\mathrm{Si}{\rm{III}}}/(1-\bar{F})$, with ${f}_{\mathrm{Si}{\rm{III}}}$ being a redshift-independent quantity that was fit using the entire data set. We adopt this same parameterization but opt to fit for a unique value of ${f}_{\mathrm{Si}{\rm{III}}}$ at each redshift. Therefore we have five free parameters ${T}_{0},\gamma ,{\lambda }_{J},\bar{F},{f}_{\mathrm{Si}{\rm{III}}}$ to fit to our power spectrum measurements in each redshift bin.

We used the publicly available MCMC package "emcee" (Foreman-Mackey et al. 2013) based on an affine-invariant ensemble sampler (Goodman & Weare 2010) to sample the posterior distribution:

$$\begin{eqnarray}&&P(\mathrm{model}| \mathrm{data})=\displaystyle \frac{P(\mathrm{data}| \mathrm{model})P(\mathrm{model})}{P(\mathrm{data})}\end{eqnarray} \tag{ 12 }$$

with $P(\mathrm{data}| \mathrm{model})$ being the combined likelihood of both data sets for a given set of model parameters, and $P(\mathrm{model})$ being the prior on that combination of parameters. We assume that the priors on parameters are independent, so $P(\mathrm{model})$ is just the product of the individual priors for each parameter. For the thermal parameters, we adopt flat priors on γ in the range 0.5 < γ < 2.1, $\mathrm{log}({\lambda }_{J})$ in the range $22\,\mathrm{ckpc}\lt {\lambda }_{J}\lt 150\,\mathrm{ckpc}$, and $\mathrm{log}({T}_{0})$ in the range $3000\,{\rm{K}}\lt {T}_{0}\lt 20000\,{\rm{K}}$. For the mean flux, we adopt Gaussian priors based on the most recent measurements of Becker et al. (2013), Faucher-Giguère et al. (2008b), and Kirkman et al. (2005) (depending on redshift) for $\bar{F}$. For the Si iii correlations, we used a Gaussian prior defined by the best-fit value in ${f}_{\mathrm{Si}{\rm{III}}}$ and the 1σ region of Palanque-Delabrouille et al. (2013). All priors for all redshifts are listed in Table 4. Note that as long as a good fit to the data can be obtained (i.e., the model is sufficiently flexible), the exact values of thermal parameters (and therefore also the exact priors chosen) do not matter in the context of this work as the fits are only used to correct out the window function, as we discuss in the next subsection. In a companion paper doing thermal parameter analysis based on hydrodynamical simulations, we will loosen the priors on both $\bar{F}$ and ${f}_{\mathrm{Si}{\rm{III}}}$ to not bias our results. Fits were always performed for individual redshifts, and no correlation between thermal properties at different redshifts was assumed.

Table 4.  Priors Used for Each Fitting Parameter

Parameter	Type of prior	Lower limit	Upper limit	μ	σ
$\mathrm{log}({T}_{0}/{\rm{K}})$	flat	3.48	4.30	⋯	⋯
γ	flat	0.5	2.1	⋯	⋯
$\mathrm{log}({\lambda }_{J}/\mathrm{ckpc})$	flat	1.34	2.18	⋯	⋯
${f}_{\mathrm{Si}{\rm{III}}}$	Gaussian	−0.002	0.018	0.008	0.001
$\bar{F}(z=1.8)$	Gaussian	0.871	0.931	0.901	0.006
$\bar{F}(z=2.0)$	Gaussian	0.785	0.976	0.881	0.019
$\bar{F}(z=2.2)$	Gaussian	0.818	0.921	0.869	0.010
$\bar{F}(z=2.4)$	Gaussian	0.766	0.859	0.812	0.009
$\bar{F}(z=2.6)$	Gaussian	0.723	0.813	0.768	0.009
$\bar{F}(z=2.8)$	Gaussian	0.683	0.771	0.727	0.009
$\bar{F}(z=3.0)$	Gaussian	0.640	0.724	0.682	0.008
$\bar{F}(z=3.2)$	Gaussian	0.581	0.661	0.621	0.008
$\bar{F}(z=3.4)$	Gaussian	0.528	0.602	0.565	0.007

Note. Priors in $\bar{F}$ are based on Kirkman et al. (2005) for z = 1.8, Faucher-Giguère et al. (2008b) for z = 2, and Becker et al. (2013) for the higher redshifts.

Download table as:  ASCIITypeset image

The impact of masking on the power depends on the underlying shape of the power spectrum. We use the fits to the underlying power based on our DM models to determine the impact of the window function on the shape of the power spectrum. We emphasize here that we use this window-function-corrected power mainly for visualization purposes only as the covariance in the measurement is subtly changed by the correction. While we do provide an approximate covariance matrix in Appendix D, ideally a full forward model of the power spectrum should be used when doing parameter estimations.

Our measurement, emulation, and window function procedures are illustrated in Figure 6. On the left side we show a comparison between our raw-data power spectrum (i.e., before window function correction) and the BOSS measurement as points. After applying our fitting procedure, we obtain an MCMC chain from which we can draw parameter combinations ${\rm{\Theta }}=\{{T}_{0},{\lambda }_{P},\gamma ,\bar{F},{f}_{\mathrm{Si}{\rm{III}}}\}$ that are compatible with the data. Feeding random draws from this chain into our emulator routines for both the forward model and the perfect model produces the black and purple bands. We can see that these provide good fits to our data set and the BOSS data set, respectively. For a single draw Θ from the posterior, we can then measure a window function correction f_window using both emulators as

$$\begin{eqnarray}&&{f}_{\mathrm{window}}(k)=\displaystyle \frac{{P}_{\mathrm{forward}\mathrm{model}}(k,{\rm{\Theta }})}{{P}_{\mathrm{perfect}\mathrm{model}}(k,{\rm{\Theta }})}.\end{eqnarray} \tag{ 13 }$$

The gray bands on the right side of Figure 6 show the resulting 16% and 84% quantiles of this quantity. The dominant effect of windowing our data is a strong increase in power on the small-scale ($k\gtrsim 0.07\,{\rm{s}}\,{\mathrm{km}}^{-1}$) end of the measurement.

Zoom In Zoom Out Reset image size

To extract the underlying (i.e., window function corrected) power spectrum from our raw measurement, we proceed as follows. We generate 1000 random draws from a multivariate normal distribution with a mean given by the raw P(k) measurement (which includes the effect of the window function) and a covariance matrix defined in Equation (7) (based on both the raw measurement and the best-fit model). We then obtain the window function correction ${f}_{\mathrm{window}}(k)$ for 1000 draws from the posterior. By multiplying each P(k) draw with each ${f}_{\mathrm{window}}(k)$ draw, we obtain samples of the distribution describing the window-corrected measurement. We take the mean of these samples to represent our window-corrected power, and the covariance of these samples gives the respective window-corrected covariance. We show the resulting correlation matrix applying this procedure at z = 2.8 in Figure 5 (lower right panel). Compared to the uncorrected forward model, we can see additional correlation at the smallest scales ($k\gtrsim 0.07\,{\rm{s}}\,{\mathrm{km}}^{-1}$) where the influence of the window correction is strongest. The window-function-corrected measurements and correlation matrices are tabulated in Appendix D and can be used for comparison with other data sets as well as model fitting. We also provide MCMC chains of f_window based on our forward models to facilitate reproduction of our measurements.

In previous work, the effect of metals on the small-scale power spectrum has been largely ignored. While McDonald et al. (2000) computed the power spectrum on a data set with and without metals masked and also considered the effect of masking on the power spectrum, they did not combine both results to see the net effect of metal contamination. However, they did note that, especially for $k\gt 0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$, the effect of metals as well as noise becomes too strong for their measurement to be usable. Motivated by this conclusion, later measurements by Croft et al. (2002), Kim et al. (2004), and Viel et al. (2008, 2013) also ignored these small-scale (larger k) modes. For lower resolution measurements using SDSS (McDonald et al. 2006) or BOSS (PD+13) spectra, these modes are well above the resolution limit of the data, and this has not been an issue. Instead, for those measurements, the power in metals is estimated from lower redshift data (using a spectral range redward of the Lyα forest). As this procedure cannot treat metal lines that are always inside the forest, large-scale correlations between Si iii and Lyα are the dominant effect of metal contamination in this case.

In Figure 7 we show a comparison between our power spectrum measurement applying the metal masking procedure described in Section 2.2 and performing the analysis without masking metals. Both measurements have been noise subtracted following the discussion in Section 3 and corrected for their respective window functions according to Section 3.8. We also show the BOSS measurement as a comparison for z > 2.2 where those measurements exist.

Zoom In Zoom Out Reset image size

It is clear that, particularly at small scales ($k\gtrsim 0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$), metal lines significantly contribute to the measured power spectrum. This increase in small-scale power in general leads to an underestimation of the small-scale thermal cutoff, and naively fitting this metal-contaminated power would result in lower overall IGM temperatures (i.e., more small-scale structure and hence less thermal broadening or pressure smoothing). For all further analysis, we therefore use the metal-masked power (black dots in the figure), and our model fitting is conservatively restricted to modes with $k\lt 0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$, where the impact of metal line contamination is relatively weak, and hence our metal-masked power is relatively insensitive to the fact that we may not have masked all of the metals (see discussion in Section 2.2).

Figure 7 also indicates that the impact of metal line contamination is not restricted to small scales. Multiple ionic metal-line transitions are typically associated with a given absorber redshift, and because the velocity separations are thousands of km s⁻¹, these large-scale velocity correlations can impact the power spectrum on large scales (low k) as well. Masking metals in the data at, for example, z = 2.6 decreases this effect, lowering the power spectrum, and therefore increases the agreement with the PD+13 measurement.

To summarize, metal line contamination of the power spectrum does not significantly impact our results, given that we mask the metals and conservatively restrict our power spectrum fits to low $k\lt 0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$, where the impact of residual (i.e., metals we missed in our masking) metal line contamination should be negligible. In principle, a more careful treatment of metal lines (such as by improved masking, or forward modeling the metals, or subtracting the red-side metal-line power) could allow one to access even smaller scales (higher k), although this would also require a very careful treatment of the noise and instrumental resolution, whose effect also increases for smaller scales (larger k).

In Figure 8 we show our new metal- and window-corrected measurement of the high-resolution power spectrum compared to the BOSS measurement from PD+13. Note that different power spectrum bins are correlated, and the error bars only reflect the diagonal elements of the covariance matrix. Where both measurements overlap ($k\lesssim 0.02\,{\rm{s}}\,{\mathrm{km}}^{-1}$, 2.2 < z < 3.4), we find good agreement with the BOSS power spectrum (which has a typical accuracy of ∼2%) within our errors at modes $k\gt 0.01\,{\rm{s}}\,{\mathrm{km}}^{-1}$, and the agreement is generally good for larger scale (lower k) modes as well. The discrepancies at low k are most likely due to the fact that our measurements are rather noisy for small wavenumbers.¹¹

Zoom In Zoom Out Reset image size

We also compare to the recent results of Iršič et al. (2017a) and Yèche et al. (2017) based on the XQ-100 data set (López et al. 2016). Specifically, Iršič et al. (2017a) measured the power at 3.0 ≤ z ≤ 4.2, and their redshift bins at z = 3.0, 3.2, and 3.4 match those used in our analysis (Yèche et al. 2017 chose to use broader bins, and their z = 3.2 bin can be compared to our results). Those are compared in Figure 9.

Zoom In Zoom Out Reset image size

While the agreement between our measurements at z = 3 is good, at z = 3.2 and z = 3.4 Iršič et al. (2017a) measures ∼10% to 50% higher power than we do at small scales $k\gtrsim 0.02\,{\rm{s}}\,{\mathrm{km}}^{-1}$, which is clearly statistically significant given the error bars for the respective data sets. Note, however, that this disagreement is restricted to high k, and we agree well with Iršič et al. (2017a) at intermediate k in all redshift bins. While it is difficult to be certain about the source of this discrepancy, given the different methodologies used for measuring the power, the k dependence of this disagreement provides a very important clue. Note that X-SHOOTER data are significantly lower in resolution than the data set presented here, and the resolution corrections become significant at higher k.

Possible uncertainties of spectroscopic resolution can come from several sources: (1) As X-SHOOTER is a slit spectrograph, its resolution is seeing dependent, so the seeing itself needs to be accurately determined to get an accurate resolution. (2) Seeing changes between different observations and correcting assuming a constant resolution across the data set might bias the measurement (see Appendix A for more details on those points). (3) The UVB slit resolution quoted for X-SHOOTER might be significantly underestimated (see Appendix B). Because of those problems, a correction based on the nominal slit resolution can overproduce small-scale power in the data analysis. However, Iršič et al. (2017a) already used a higher resolution than provided in the XQ-100 data release paper when performing his corrections.¹² Using the Iršič et al. (2017a) value for the resolution leads to agreement between both XQ-100 power spectrum analyses by Iršič et al. (2017a) and the independent determination of the power spectrum from the same data set by Yèche et al. (2017). Yèche et al. (2017) determined the spectral resolution by assuming the XQ-100 resolution and estimating the seeing on an object-by-object basis, analyzing the width of the Lyα forest in the spatial direction of the 2D spectra. While this approach should remove the sensitivity toward seeing, it cannot tackle a potentially misestimated slit resolution or changes of this quantity along the spectral arm. Yèche et al. (2017) also chose to not combine the different spectral arms. Thereby, there should not be a strong change of resolution when reaching the overlap region between both spectral arms. Thus, inside each redshift bin, their data should be more homogeneous regarding resolution. Nevertheless, in the end, both measurements seem to give essentially the same result.

The influence of resolution errors on the resolution correction factor W_R² (see Equation (5)) can be found by simple error propagation,

$$\begin{eqnarray}{\rm{\Delta }}\mathrm{ln}{W}_{R}^{2} & = & -{k}^{2}{\rm{\Delta }}({R}^{2})\\ & \simeq & -2{k}^{2}{R}^{2}{\rm{\Delta }}\mathrm{ln}(R),\end{eqnarray} \tag{ 14 }$$

and propagates to an error on the estimated power spectrum P as

$$\begin{eqnarray}&&{\rm{\Delta }}\mathrm{ln}P=-{\rm{\Delta }}\mathrm{ln}{W}_{R}^{2}.\end{eqnarray} \tag{ 15 }$$

Assuming the nominal resolving power of 4350 (according to López et al. 2016, listing a slightly higher value than the ESO instrument description) using the X-SHOOTER 1'' slit on the UVB arm as a worst-case scenario, a spectrum with 10% higher resolution than assumed when performing the correction will lead to a ∼45% (28% assuming a resolving power of 5350) overestimate of the power in the resolution-corrected measurement. As our own power spectrum measurement is based on ∼10 times higher resolution data, a comparable error in the knowledge of the resolution will have a ∼100 times smaller effect (so ∼0.4% at $k=0.05\,{\rm{s}}\,{\mathrm{km}}^{-1}$ and ∼1.5% at $k=0.1\,{\rm{s}}\,{\mathrm{km}}^{-1}$). Therefore, a lack of knowledge of the precise resolution of the spectrograph, a significant concern for the X-SHOOTER measurement, can be safely ignored for our study.

To determine the influence of possible errors in the resolution estimates due to points (1) and (2), we divided out the resolution correction from the Iršič et al. (2017a) power spectrum measurement and corrected using different assumptions. First, we used the nominal slit resolution R = 4350 from the X-SHOOTER spectrograph and generated Gaussian-distributed seeing values (with a mean of 075 and FWHM of 02, similar to the distribution shown in Yèche et al. 2017) for each, and we computed the resolution correction W_R² according to Equation (5) and obtained the mean correction, which we then used as the new resolution correction. This gives results that are basically identical(that we do not show in our comparison figure for clarity) to the original Iršič et al. (2017a) measurement, showing that the seeing estimate used for their measurement is in agreement with the distribution determined by Yèche et al. (2017) and cannot be the reason for the disagreement with our measurement. In addition, we also generate a measurement assuming a higher slit resolution of R = 5350 (due to point (3) and in agreement with measurements based on calibration spectra; see Appendix B) and otherwise performing the same analysis. This comparison is shown as red diamonds in Figure 9, and we can see that the agreement between our high-resolution measurement and XQ-100 at z = 3.2 and z = 3.4 is good in this case. However, the agreement in the z = 3.0 bin without assuming a different X-SHOOTER resolution correction is unclear, but might be due to possible variations in the resolution between echelle orders. Note that for the other spectral arms (that cover the Lyα forest for z > 3.6), this is a less severe problem as data from those are intrinsically higher resolution. Also note that, for those arms, one can in principle obtain the resolution of the science observation from the width of the telluric absorption lines (see Bosman et al. 2017), but those are rare in the UVB arm. This can also explain the agreement between XQ-100 measurements and older high-resolution data by McDonald et al. (2000) at z ∼ 3.8 (but as the redshift bins of both measurements are significantly different, this comparison is tricky) and Viel et al. (2013) at z = 4.2.

Because of the severe impact the resolution correction can have on the XQ-100 power spectrum measurement, we caution against using the smallest scales ($k\gt 0.02\,{\rm{s}}\,{\mathrm{km}}^{-1}$) of this measurement (at least for the lowest redshift bins) for parameter studies without additional treatments for spectroscopic resolution. This is especially true for measurements that rely on the accurate determination of the power spectrum cutoff, like determining the thermal state of the IGM or the nature of DM (e.g., Baur et al. 2017; Iršič et al. 2017b, 2017c). Although to be fair, for the latter most of the sensitivity comes from the higher redshift (z ≳ 4) bins where the resolution of X-SHOOTER is higher and additionally high-resolution (R ≃ 50,000) data from Viel et al. (2013) are used.

In Figure 4 we can see that our metal removal and masking correction techniques do not change the data hugely (the difference is covered by our error bars) at the redshifts and scales where we disagree with Iršič et al. (2017a). Additionally, the changes due to metal masking do not exhibit the same shape as the discrepancy. Finally, we also checked the raw power spectra not corrected for any masking and could not get a small-scale power as high as in the Iršič et al. (2017a) result. We are therefore confident that metal masking and window correction are not the reason for the discrepancies.

Our measurements probe the small-scale cutoff in the power (0.02 ≲ k ≲ 0.1) in all redshift bins with a typical precision of 5%–15%. The position of this cutoff is still at far larger scales then the expected cutoff, due to the spectroscopic resolution of our data and the observed cutoff in the power results from thermal or pressure broadening of the absorption lines (or, e.g., warm/fuzzy DM). In the next subsection, we also compare to previous high-resolution measurements to make sure our measurement agrees with existing data on small scales as well.

Previous power spectrum measurements based on high-resolution data were obtained by different groups using redshift bins of different size and location. They also differed in the Fourier normalization convention used (leading to factors of 2 between some measurements), the field of which the power is computed (flux contrast or transmission leading to additional factors of ${\bar{F}}^{2}$ in the power), as well as whether metals were masked and noise was subtracted. We therefore show comparisons at the quoted redshifts for the old high-resolution measurements and linearly interpolate our results to those redshifts and renormalize the different measurements to the power spectrum convention that we use (see Appendix C).

For authors who chose to study the F field instead of δ_F, there is ambiguity in the mean flux. While the mean flux of the IGM is clearly the same, the mean flux of the data is sensitive to where the continuum is placed. It is well known that hand continuum fits to high-resolution spectra are biased low (Faucher-Giguère et al. 2008a), and this systematic effect is a bigger issue at higher redshift. If one measures the power spectrum of F and the continuum fits are biased low, then the power will be biased high. McDonald et al. (2000) provide measurements of their mean flux, and we can therefore easily correct this effect, whereas Croft et al. (2002) does not. Thus a direct comparison to Croft's measurements is not straightforward, but we do our best by simply assuming their continua are unbiased and multiplying in the mean flux of the IGM measured by Becker et al. (2013) at the respective redshifts of their measurements. The differences between power spectrum conventions clearly limit the precision at which our measurements will agree with previous work.

The comparison between high-resolution measurements is shown in Figure 10. While overall agreement is good considering the different approaches, some comparisons show disagreement. The strongest mismatch is with Croft et al. (2002) at z = 2.13 on all scales. At similar redshifts we agree with both Kim et al. (2004) and PD+13, which hints at an incorrect mean flux or a problem with the Croft et al. (2002) measurement at this redshift (which is not part of their fiducial sample). We can also see that the shape of the different measurements on scales smaller (larger k) than the cutoff matches between the four high-resolution data sets for scales 0.01 ≲ k ≲ 0.08. On smaller scales, the difference in treatment of metals and S/N of the data sets as well as removal of noise power can probably explain the deviations between these measurements. We do note that our z = 3.0 bin exhibits a cutoff at slightly smaller scales (larger k) than the one of McDonald et al. (2000), but agrees with Croft et al. (2002) at essentially the same redshift. This shows that there are clear limitations of comparisons to the previous measurements from high-resolution data sets due to the different conventions.

Zoom In Zoom Out Reset image size

To summarize, we find reasonable agreement between our measurements and previous analyses. We will use our new result for parameter estimations in future works.