Simulating the effect of high column density absorbers on the one-dimensional Lyman-alpha forest flux power spectrum

TL;DR

This paper investigates how high column density hydrogen absorbers affect the 1D Lyman-alpha forest flux power spectrum using hydrodynamical simulations, providing models to improve cosmological analyses.

Contribution

It introduces the first simulation-based analysis of high column density absorbers' impact on the Lyman-alpha flux power spectrum and offers templates for residual contamination modeling.

Findings

High column density absorbers bias the flux power spectrum through damping wings.

Provided templates as a function of column density and redshift for modeling.

Fitted functions enable straightforward incorporation into data analysis pipelines.

Abstract

We measure the effect of high column density absorbing systems of neutral hydrogen (HI) on the one-dimensional (1D) Lyman-alpha forest flux power spectrum using cosmological hydrodynamical simulations from the Illustris project. High column density absorbers (which we define to be those with HI column densities $N(\mathrm{HI}) > 1.6 \times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$) cause broadened absorption lines with characteristic damping wings. These damping wings bias the 1D Lyman-alpha forest flux power spectrum by causing absorption in quasar spectra away from the location of the absorber itself. We investigate the effect of high column density absorbers on the Lyman-alpha forest using hydrodynamical simulations for the first time. We provide templates as a function of column density and redshift, allowing the flexibility to accurately model residual contamination, i.e., if an analysis selectively clips out the largest damping wings. This flexibility will improve cosmological parameter estimation, e.g., allowing more accurate measurement of the shape of the power spectrum, with implications for cosmological models containing massive neutrinos or a running of the spectral index. We provide fitting functions to reproduce these results so that they can be incorporated straightforwardly into a data analysis pipeline.

Figures

Tables

Table 1. Table 1: The neutral hydrogen ( Hi ) column density limits [ N ​ ( Hi ) min , N ​ ( Hi ) max ] 𝑁 subscript Hi min 𝑁 subscript Hi max [N(\textsc{Hi})_{\mathrm{min}},N(\textsc{Hi})_{\mathrm{max}}] that define the categories of absorbing systems used in this work. The columns on the right show the percentage of spectra (at each redshift z 𝑧 z that is considered) in our ( 106.5 ​ Mpc ) 3 superscript 106.5 Mpc 3 (106.5\,\mathrm{Mpc})^{3} simulation box (Vogelsberger et al., 2014b ; Nelson et al., 2015 , Illustris-1) where the highest-density system belongs to a given category.

Absorber category	$N{\left(\text{Hi}\right)}_{\min }$	$N{\left(\text{Hi}\right)}_{\max }$	% of spectra in ${\left(106.5\mathrm{Mpc}\right)}^3$ simulation at
	[$\mathrm{atoms}{\mathrm{cm}}^{-2}$]		$z=2.00$	$z=2.44$	$z=3.01$	$z=3.49$	$z=4.43$
Lyman-alpha forest	0	$1.6\times {10}^{17}$	77.7	69.6	57.4	45.7	22.0
LLS	$1.6\times {10}^{17}$	$1\times {10}^{19}$	10.6	14.9	21.8	27.0	36.6
Sub-DLA	$1\times {10}^{19}$	$2\times {10}^{20}$	5.9	8.1	11.4	14.3	20.1
Small DLA	$2\times {10}^{20}$	$1\times {10}^{21}$	3.1	4.1	5.5	7.8	12.8
Large DLA	$1\times {10}^{21}$	$\infty$	2.7	3.3	3.9	5.2	8.5

Table 2. Table 2: Best-fit values of the parameters in our templates for the bias of spectra contaminated by high column density absorbers on the one-dimensional Lyman-alpha forest flux power spectrum. The template parameters are defined in Eqs. ( 6 ) and ( 7 ). Values are shown for each high column density absorber category. The definitions of the different categories of high column density absorber are given in Table 1 .

Absorber category	Template parameter values
	$a_0$	$a_1$	$b_0$	$b_1$	$c_0$	$c_1$
LLS	2.2001	0.0134	36.449	-0.0674	0.9849	-0.0631
Sub-DLA	1.5083	0.0994	81.388	-0.2287	0.8667	0.0196
Small DLA	1.1415	0.0937	162.95	0.0126	0.6572	0.1169
Large DLA	0.8633	0.2943	429.58	-0.4964	0.3339	0.4653

Equations

$$P^{1D}(k_{||},z)=\int P^{3D}(k_{||},\text{mn@boldsymbol}{k}_{\perp },z)\frac{\mathrm{d}\text{mn@boldsymbol}{k}_{\perp }}{(2\pi {)}^2},$$

$$P_i^{1D}(k_{||},z)={⟨\frac{P_{\mathrm{Raw}}^{1D}(k_{||},z)}{W^2(k_{||},\Delta v,R)}⟩}_i,$$

$$W(k_{||},\Delta v,R)=\exp \left(-\frac{1}{2}(k_{||}R{)}^2\right)\times \frac{\sin (k_{||}\Delta v/2)}{k_{||}\Delta v/2},$$

$$P_{\mathrm{Total}}^{1D}(k_{||},z)=\sum _i{\alpha }_i(z)P_i^{1D}(k_{||},z),$$

$$P_{\mathrm{Total}}^{1D}(k_{||},z)=P_{\mathrm{Forest}}^{1D}(k_{||},z)\left[{\alpha }_{\mathrm{Forest}}(z)+\sum _{i\ne \mathrm{Forest}}{\alpha }_i(z)\frac{P_i^{1D}(k_{||},z)}{P_{\mathrm{Forest}}^{1D}(k_{||},z)}\right].$$

$$\frac{P_i^{1D}(k_{||},z)}{P_{\mathrm{Forest}}^{1D}(k_{||},z)}={\left(\frac{1+z}{1+z_0}\right)}^{-3.55}\times \frac{1}{(a(z)e^{b(z)k_{||}}-1{)}^2}+c(z),$$

$$a(z)=a_0{\left(\frac{1+z}{1+z_0}\right)}^{a_1};b(z)=b_0{\left(\frac{1+z}{1+z_0}\right)}^{b_1};c(z)=c_0{\left(\frac{1+z}{1+z_0}\right)}^{c_1};$$

$$\begin{array}{rl}{P}_{\mathrm{Total}}^{1D}(k_{||},z)=P_{\mathrm{Forest}}^{1D}& (k_{||},z)[{\alpha }_0(z)\\ & +\sum _{i\ne \mathrm{Forest}}{\alpha }_i(z){\left(\frac{1+z}{1+z_0}\right)}^{-3.55}\times \frac{1}{(a(z)e^{b(z)k_{||}}-1{)}^2}].\end{array}$$

$$\tau (\lambda )=N\alpha (\lambda )=N\frac{\sqrt{\pi }{e}^2}{4\pi {ϵ}_0{m}_{\mathrm{e}}{c}^2}\frac{f{\lambda }_{\mathrm{t}}^2}{\Delta {\lambda }_{\mathrm{D}}}u(x,y),$$

$$\Delta {\lambda }_{\mathrm{D}}=\frac{{\lambda }_{\mathrm{t}}}{c}{\left(\frac{2k_{\mathrm{B}}T}{m_{\mathrm{ion}}}\right)}^{\frac{1}{2}}.$$

$$w(z)=e^{-z^2}\mathrm{erfc}(-iz)=u(x,y)+iv(x,y),$$

$$x(\lambda )=\frac{\lambda -{\lambda }_{\mathrm{c}}}{\Delta {\lambda }_{\mathrm{D}}};y=\frac{\Gamma {\lambda }_{\mathrm{t}}^2}{4\pi c}\frac{1}{\Delta {\lambda }_{\mathrm{D}}},$$

Full text

1 Introduction

The Lyman-alpha forest (a series of neutral hydrogen absorption lines in the spectra of quasars) is a uniquely powerful probe of the clustering of matter at redshifts from about $z=2$ to $z=6$ (Croft et al., 1998; McDonald et al., 2000, 2005b; Viel et al., 2013; Palanque-Delabrouille et al., 2013; Iršič et al., 2016) and from sub-Mpc to hundreds of Mpc scales. The one-dimensional (1D) Lyman-alpha forest flux power spectrum (along the line of sight) is particularly sensitive to small-scale clustering in the quasi-linear regime and provides important constraints on extended cosmological models that suppress small-scale power (Seljak et al., 2005; Palanque-Delabrouille et al., 2015; Iršič et al., 2017b; Yeche et al., 2017; Iršič et al., 2017a; Armengaud et al., 2017), notably those containing massive neutrinos and warm dark matter. This small-scale information complements the larger scales probed by the angular power spectrum of the cosmic microwave background (CMB). For example, the best upper limit on the sum of neutrino masses (Palanque-Delabrouille et al., 2015) comes from combining CMB data from the Planck Collaboration (Planck Collaboration et al., 2016) with the 1D Lyman-alpha forest power spectrum as measured from Sloan Digital Sky Survey (SDSS)-III/Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 9 (DR9) quasar spectra (Eisenstein et al., 2011; Dawson et al., 2013; Palanque-Delabrouille et al., 2013).

Future surveys like the Dark Energy Spectroscopic Instrument (DESI; DESI Collaboration et al., 2016a, b) will further improve constraints on extended cosmological models. Font-Ribera et al. (2014) forecast one-sigma errors on a DESI measurement of the sum of neutrino masses to be 0.017 eV111This is the full forecasted constraint considering a combination of Planck CMB data, DESI broadband galaxy power spectrum, DESI broadband Lyman-alpha forest flux power spectrum and $\sim 100$ high-resolution Lyman-alpha forest quasar spectra.. Considering that the lower limit on the sum of neutrino masses from neutrino oscillation experiments is 0.06 eV (Gonzalez-Garcia et al., 2014; Forero et al., 2014; Esteban et al., 2017), this would constitute at least a three-sigma detection. Furthermore, the 1D Lyman-alpha forest flux power spectrum probes the primordial power spectrum on the smallest currently accessible scales, $k\sim 4$ Mpc*-1*. Including Lyman-alpha forest data will improve constraints on the running of the spectral index (which quantifies deviations from a pure power-law spectrum) by a factor of two, reaching one-sigma errors of $\pm 0.002$ (Font-Ribera et al., 2014). This would provide new insights into early universe physics, potentially ruling out classes of models of inflation. Importantly, it will also provide a unique independent cross-check at small scales of the primordial power spectrum shape inferred from CMB measurements at large scales.

Achieving these limits requires marginalisation over the uncertain impact of a number of astrophysical effects on the 1D Lyman-alpha forest power spectrum. In particular, this includes broadened absorption features from high column density absorbers. High column density absorbers are usually classified as either damped Lyman-alpha absorbers (DLAs), with column densities $N(\textsc{Hi})$ exceeding $2\times 10^{20}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$ (Wolfe et al., 1986), or Lyman-limit systems (LLS), which correspond to $2\times 10^{20}\,\mathrm{atoms}\,\mathrm{cm}^{-2}>N(\textsc{Hi})>1.6\times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$. Both types of system produce broad damping wings which extend large distances in redshift space. If not accounted for, they will bias cosmological parameter estimation from the Lyman-alpha forest. The systems are formed at peaks of the underlying density distribution; consequently, they cluster more strongly than the forest itself (Font-Ribera et al., 2012).

To remove the bias induced by damped absorbers, one can fit a model for their effect on power spectra. The most widely used approach (McDonald et al., 2005a) is now more than a decade old. Although this model was adequate for the data available at the time, future surveys will be substantially more constraining and therefore demand tighter control over systematics. Furthermore, there have been significant improvements in theoretical modelling of these systems (e. g., Pontzen et al., 2008; Bird et al., 2015). An updated model for the effects of high column density absorbers is therefore both timely and essential in order to achieve the forecasted cosmological constraints from future surveys.

Different column densities correspond to gas at different physical densities, so that simulations suitable for modelling the forest are often not suited to reproducing high column density systems. The Lyman-alpha forest is largely insensitive to the physics of galaxy formation since it is sourced by gas at below mean density; the primary uncertainties arise from cosmological parameters and the thermal history of the intergalactic medium. Conversely, high column density absorbers arise largely from regions within or around galaxies and are thus very sensitive to the physics of galaxy formation and less sensitive to large-scale cosmology. It is consequently essential to model the effect of high column density absorbers using simulations which include detailed galaxy formation physics and can thus reproduce their characteristics and statistics.

In Lyman-alpha forest studies, damping wings are sometimes “clipped” (i. e., removed or masked) from quasar spectra (e. g., see Lee et al., 2013, for details of the process for BOSS DR9 spectra). However, not all damping wings are identified and many will remain in the spectra, especially in noisier spectra where they are harder to spot and for lower-density absorbers (i. e., LLS) which have narrower wings. Therefore, in the final cosmological parameter estimation from the 1D Lyman-alpha forest power spectrum, the effect of residual high column density absorbers is modelled as a multiplicative scale-dependent bias of the power spectrum with an amplitude (reflecting the level of residual contamination) that is fitted and marginalised (Palanque-Delabrouille et al., 2015). The functional form of this model (i. e., its scale and redshift dependence) is based on the measurements made in McDonald et al. (2005a).

McDonald et al. (2005a) investigated the effect with lognormal model mock quasar spectra (i. e., generated without hydrodynamical simulations; details of their generation are given in McDonald et al., 2006), since the numerical simulations available at the time were not large enough to generate spectra encompassing the full width of damping wings. They then probe the effect of high column density absorbers on the Lyman-alpha forest by inserting damping wings in mock spectra at the peaks of the lognormal field, based on the observationally-determined column density distribution function (CDDF). They find a systematic effect on the observed 1D Lyman-alpha forest power spectrum that is maximised on scales corresponding to the width of a damped system and which has negligible redshift evolution (considering three redshift slices at $z=[2.2,3.2,4.2]$). They provide a single template to fit their bias measurement, including the effect of all LLS and DLAs together. However, as discussed above, in current data analysis pipelines, damping wings are removed from quasar spectra in a way that preferentially removes higher density systems. Therefore, when the template is used in parameter inference, it may not correctly model the bias of the residual contamination, which will have a different CDDF to the total — the clipping of the survey spectra changes the survey CDDF. The bias will have a different scale-dependence (not just amplitude), since this is driven by the distribution of the widths of damping wings remaining in quasar spectra.

In this work, we investigate the effect of high column density absorbers on the 1D Lyman-alpha forest power spectrum as a function of their column density and redshift using hydrodynamical simulations of galaxy formation from the Illustris project (Vogelsberger et al., 2014b; Nelson et al., 2015). Comparison to relevant observations has shown that Illustris reproduces the observed CDDF and spatial clustering of high-density systems (Vogelsberger et al., 2014b; Bird et al., 2014, see § 3.1 for more details) at the 95% confidence level. Spectra are generated from this simulation, then separated into categories according to the maximum column density within each spectrum (see § 2 for more details). We measure the 1D flux power spectrum of each of these types of spectrum and measure the (multiplicative) bias of each type compared to the power spectrum of the Lyman-alpha forest alone. We make this measurement at multiple redshifts and so probe the redshift evolution of this effect.

We discuss high column density absorbers in more detail in § 2. In § 3, our methodology in going from hydrodynamical simulations to measurements of the 1D flux power spectrum is explained. We present our main results in § 4. These results are discussed in § 5 and in § 6, we present the templates that we have fitted to our measurements. Finally, conclusions are drawn in § 7.

2 Damped Lyman-alpha absorbers and Lyman-limit systems

High column density absorbers are regions of neutral hydrogen (Hi) gas that are above a column density threshold of $N(\textsc{Hi})>1.6\times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$. By contrast, lower column density absorbers form the Lyman-alpha forest. The absorption lines formed by high column density absorbers are broadened, forming damping wings and hence absorption in the spectrum away from the location of the absorbing gas. The damping wings have a characteristic Voigt profile, which is a convolution of a Gaussian profile (caused by Doppler broadening) and a Lorentzian profile (caused by natural or collision broadening). The width of these wings in velocity space increases with the column density of the absorbing system. High column density absorbers are then usually classified as either damped Lyman-alpha absorbers (DLAs), whose damping wings are considered significantly broadened and which correspond to $N(\textsc{Hi})>2\times 10^{20}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$ (Wolfe et al., 1986); or Lyman-limit systems (LLS), which correspond to column densities in the range $2\times 10^{20}\,\mathrm{atoms}\,\mathrm{cm}^{-2}>N(\textsc{Hi})>1.6\times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$.

In this work, we aim to investigate the effect of high column density absorbers (and especially their damping wings) on the one-dimensional Lyman-alpha forest flux power spectrum, as a function of their column density (and redshift). We therefore use a more refined classification of high column density absorbers based on their column densities, in particular accounting for the fact that higher density LLS do have wide damping wings. Table 1 shows the column density limits that define our categories, as well as the percentage of simulated spectra (see § 3.1) where the highest-density system is a given type and hence is the main contaminant. The overall percentage of spectra contaminated by high column density absorbers (LLS, sub-DLAs, small and large DLAs) increases with redshift because the Hi CDDF increases at higher densities at higher redshifts, but always there are more LLS than DLAs.

3 Method

We first outline the method we have used and then explain the steps in more detail in the following subsections (§ 3.1 to 3.3).

1. (1).

We use a cosmological hydrodynamical simulation from the Illustris project (Vogelsberger et al., 2014b; Nelson et al., 2015) and generate mock spectra on a grid (562 500 in total, each at a velocity resolution of $25\,\mathrm{km}\,\mathrm{s}^{-1}$ and with a typical length of $\simeq 8\,000\,\mathrm{km}\,\mathrm{s}^{-1}$). We repeat this for a number of redshift slices at which the Lyman-alpha forest is observed ($z=[2.00,2.44,3.01,3.49,4.43]$). (See § 3.1.) 2. (2).

For each redshift slice, we separate the spectra according to the highest column density system within that spectrum using the absorber categories defined in Table 1. For each absorber category (and the total set of spectra), we measure the one-dimensional (1D) flux power spectrum (i. e., along the line of sight, integrating over transverse directions) using a fast Fourier transform (FFT). (See § 3.2.) 3. (3).

We then measure the (multiplicative) bias of the flux power spectra from each category relative to the 1D flux power spectrum of the Lyman-alpha forest, as a function of absorber type (i. e., maximum column density) and redshift (see § 3.3). We fit parametric models to these bias measurements and provide these templates in § 6.

3.1 Hydrodynamical simulations and mock spectra

Our main results make use of snapshots from the highest-resolution (in terms of both dark matter particles and hydrodynamical cells) cosmological hydrodynamical simulation from the Illustris project (Vogelsberger et al., 2014b; Nelson et al., 2015, Illustris-1222The simulation we use is publically available at http://www.illustris-project.org/data.). The simulation adopts the following cosmological parameters: $\Omega_{\mathrm{m}}=0.2726$, $\Omega_{\Lambda}=0.7274$, $\Omega_{\mathrm{b}}=0.0456$, $\sigma_{8}=0.809$, $n_{\mathrm{s}}=0.963$ and $H_{0}=100\,h\,\mathrm{km}\,\mathrm{s}^{-1}\,\mathrm{Mpc}^{-1}$, where $h=0.704$ (Vogelsberger et al., 2014a). The box has a volume in comoving units of $(106.5\,\mathrm{Mpc})^{3}$ and we consider snapshots at redshifts $z=[2.00,2.44,3.01,3.49,4.43]$.

The Illustris simulations use the moving mesh code AREPO (Springel, 2010). The galaxy formation physics implemented is of relevance to dense regions of neutral hydrogen gas, and therefore we describe it briefly here. The subgrid models include prescriptions for supernova (Springel & Hernquist, 2003; Vogelsberger et al., 2013) and active galactic nuclei (AGN) (Springel et al., 2005; Sijacki et al., 2007) feedback (Bird et al. 2014 showed that the properties of DLAs are quite insensitive to the details of AGN feedback); radiative cooling; star formation and metal enrichment of gas. Self-shielding is implemented as a correction to the photoionization rate, which is a function of hydrogen density and gas temperature. The potential ionising effect of local stellar radiation within the most dense absorbers (i. e., large DLAs) (e. g., Fumagalli et al., 2011) is neglected. Pontzen et al. (2010) found this effect to be negligible and accurate calculations in any case require physics on parsec scales, well below the resolution of the simulation (it can then be viewed as part of the unresolved physics included in the above feedback prescriptions). More details of these models are given in Vogelsberger et al. (2013); Bird et al. (2014). Gravitational interactions are computed using the TreePM approach (Springel, 2005).

We require that these simulations accurately reproduce the necessary statistics of high column density absorbers that are observed in surveys. As a means of quantifying this, we can first consider the CDDF of neutral hydrogen over relevant column densities ($N(\textsc{Hi})>1.6\times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$). Vogelsberger et al. (2014b) make a comparison of the CDDF as produced by Illustris centered at $z=3$ to the distribution observed in a number of surveys [Prochaska et al. (2010) for LLS; Zafar et al. (2013) for sub-DLAs; Noterdaeme et al. (2009) for DLAs]. In particular, the distributions are consistent with the feature in the CDDF around the DLA threshold, where the distribution rises, being reproduced well (the results of Bird et al. 2017 from SDSS-III DR12 spectra are also consistent for DLAs). Bird et al. (2014) showed that the AREPO code with the above hydrodynamical models can produce values of the DLA halo bias (at $z=2.3$) which are in agreement with measured values from real surveys (Font-Ribera et al., 2012), indicating that the clustering of high column density absorbers is well reproduced. Bird et al. (2015) compared the distribution function of velocity widths of low ionization metal absorbers associated with DLAs as produced by the simulations at $z=3$ to the distribution observed in Neeleman et al. (2013). The data points are within the 68% confidence interval of the simulated distribution. This suggests that the simulations are reproducing the kinematics, and thus the host halo distribution, of high column density absorbers. One potential caveat is that these simulations are found to produce too high a total incidence rate of DLAs when compared to observations (Noterdaeme et al., 2012) at $z=2$ (Bird et al., 2014). However, the overall incidence rate is absorbed into a normalisation that must in any case be allowed to float during analysis of clipped spectra (as discussed in § 5).

For each snapshot, we generate mock spectra containing only the Lyman-alpha absorption line (i. e., with a rest wavelength of $1215.67\textup{\AA}$) from neutral hydrogen. We do this on a square grid of 562 500 spectra, in the plane perpendicular to a direction that we define as the line of sight. Each spectrum extends the full length of the simulation box with periodic boundary conditions, giving a size in velocity (or “redshift”) space of $\{7111,7501,8000,8420,9199\}\,\mathrm{km}\,\mathrm{s}^{-1}$ respectively at $z=[2.00,2.44,3.01,3.49,4.43]$333We convert the comoving length of the box to a proper velocity by the Hubble law.. We first measure the optical depth $\tau$ in velocity bins of size $25\,\mathrm{km}\,\mathrm{s}^{-1}$ along the spectrum444For comparison, BOSS DR9 spectra are binned at a velocity resolution of $69.02\,\mathrm{km}\,\mathrm{s}^{-1}$ (Lee et al., 2013).. We further convolve our spectra with a Gaussian kernel of $\mathrm{FWHM}=8\,\mathrm{km}\,\mathrm{s}^{-1}$, setting the simulated spectrographic resolution. We then calculate the transmitted flux $\mathcal{F}=e^{-\tau}$. In this way, the spectra we have constructed are insensitive to contamination from other absorption (or emission) lines, estimation of the emitted quasar continuum (which here is effectively set to unity) or instrumental noise. In each spectrum pixel, we are also able to measure the column density (integrated along the line of sight in each bin) of neutral hydrogen, which we use in measuring the maximum density systems in each spectrum (§ 3.2).

3.2 One-dimensional flux power spectrum

We separate our spectra into the absorber categories (Lyman-alpha forest, LLS, sub-DLAs, small and large DLAs) defined in Table 1 according to the maximum column density system within each spectrum. We search for the highest column density integrated over any four neighbouring velocity bins; this amounts to a comoving length along the line of sight of $\{1.50,1.42,1.33,1.27,1.16\}\,\mathrm{Mpc}$ respectively at $z=[2.00,2.44,3.01,3.49,4.43]$. The categorisation is insensitive to the number of neighbouring velocity bins that we use, as the boundaries between categories differ by orders of magnitude in column density. Moreover, the method is efficient in identifying high column density absorbers since they are vastly more dense than the surrounding gas forming the Lyman-alpha forest555We have explicitly tested the impact of doubling or halving the number of neighbouring velocity bins we use on the 1D flux power spectra we measure in each absorber category. We find that the maximum absolute difference in any power spectrum bin is a negligible 0.2%.. We have chosen a length that is much larger than the most extensive DLAs as found by recent studies (Krogager et al., 2012) and so we are sure to integrate over the full length of any high column density absorbers. Our definition of high column density absorbers includes blends, where a number of smaller, lower column density systems have been added together. In this way, we have associated with each spectrum the most dominant absorbing system and in the case where high column density absorbers are identified, these are the main contamination to the spectrum through their associated damping wings. The percentage of spectra in each absorber category at each redshift slice is given in Table 1.

We measure the 1D flux power spectrum of all the spectra and each absorber category at each redshift slice. The 1D power spectrum $P^{\mathrm{1D}}(k_{||},z)$ is defined as the integral of the three-dimensional (3D) power spectrum $P^{\mathrm{3D}}(k_{||},\mn@boldsymbol{k_{\perp}},z)$ over directions perpendicular to the line of sight:

[TABLE]

where the wavevector $\mn@boldsymbol{k}=[k_{||},\mn@boldsymbol{k_{\perp}}]$ is conjugate to velocities in real space and so is measured in units of inverse velocity (e. g., $\mathrm{s}\,\mathrm{km}^{-1}$). We also use the convention of absorbing the $2\pi$ into the conjugate variable666I. e., we define the Fourier transform as $\delta(k)=\int\delta(x)e^{-ikx}\mathrm{d}x$..

To measure $P^{\mathrm{1D}}$ for an individual line of sight, we first calculate the fluctuation in each velocity $v_{||}$ bin $\delta_{\mathcal{F}}(v_{||})=\frac{\mathcal{F}(v_{||})}{\langle\mathcal{F}\rangle}-1$, where $\langle\mathcal{F}\rangle$ is the average flux over all spectra at each redshift (Croft et al., 1998). We calculate the 1D Fourier transform along the line of sight $\hat{\delta}_{\mathcal{F}}(k_{||})$ using a fast Fourier transform (FFT)-based method since we have evenly-spaced velocity bins. We then estimate the 1D flux power spectrum for each sightline $P_{\mathrm{Raw}}^{\mathrm{1D}}(k_{||})=|\hat{\delta}_{\mathcal{F}}(k_{||})|^{2}$. Finally, we estimate the 1D flux power spectrum in Eq. (1) for each absorber category $i$ by (e. g., Palanque-Delabrouille et al., 2013)

[TABLE]

where we explicitly indicate that the raw 1D power spectra depend on redshift $z$. The average is taken over spectra of a given category (or all spectra for the total power spectrum) at each redshift slice. The window function $W(k_{||},\Delta v,R)$ that is divided out arises from the binning in velocity space ($\Delta v$) and the simulated spectrographic resolution $R$:

[TABLE]

where $\Delta v=25\,\mathrm{km}\,\mathrm{s}^{-1}$ and $R=3.40\,\mathrm{km}\,\mathrm{s}^{-1}$ (not to be confused with the spectrographic resolving power; see § 3.1). We then have an estimate of the 1D flux power spectrum for each absorber category of spectra at each redshift slice.

3.3 Modelling the effect of high column density absorbers

The total 1D flux power spectrum of a set of spectra $P_{\mathrm{Total}}^{\mathrm{1D}}(k_{||},z)$ can be expressed as a weighted sum of the 1D flux power spectra calculated in Eq. (2) for each absorber category $i$:

[TABLE]

where $\alpha_{i}(z)$ are the fraction of spectra in each absorber category at each redshift (as given in Table 1 for our simulated ensemble of spectra). In a real survey, $\alpha_{i}(z)$ may change from their raw values due to the attempt to clip (i. e., remove) high column density absorbers discussed in § 1. We can rearrange Eq. (4) to isolate the 1D flux power spectrum of the Lyman-alpha forest alone:

[TABLE]

In this way, we have isolated the effect of spectra containing high column density absorbers on the 1D flux power spectrum of the Lyman-alpha forest as a multiplicative bias (i. e., the terms in square brackets)777We could simplify this form further by asserting the fact that $\sum_{i}\alpha_{i}(z)=1$ to remove the parameter $\alpha_{\mathrm{Forest}}(z)$, but it is useful to keep this form as we explain in § 6.. This matches the general form of modelling this effect in previous studies, as explained in Palanque-Delabrouille et al. (2015) (based on the results in McDonald et al. 2005a), but now additionally probing the bias as a function of column density (i. e., by using the different absorber categories). We discuss in more detail in § 5 our motivations for using this particular form of the bias (as opposed to e. g., an additive bias). Using the 1D flux power spectra we have calculated in § 3.2, we are able to measure the fractions in Eq. (5) ($P_{i}^{\mathrm{1D}}(k_{||},z)/P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$) and we present the results in § 4.

4 Results

Figure 1 shows the 1D flux power spectra of different subsets of sightlines that we have measured from our simulations [see § 3.2 and in particular Eq. (2)] at redshift $z=2.00$. The different subsets shown are: the total as would be measured if no distinction between different types of spectra was made; spectra containing only Lyman-alpha forest (i. e., the ensemble that is uncontaminated by high column density absorbers); and spectra contaminated by different categories of high column density absorber, as defined in Table 1. We first note that the total 1D flux power spectrum at any redshift slice can be reconstructed as a weighted sum of the other 1D flux power spectra for each absorber category at that redshift (see § 3.3 and in particular Eq. (4)). The weights are the fraction of spectra in each category (the values we measure for our simulated ensemble are given in Table 1). We can estimate the fractional ($1\sigma$) statistical error on each power spectrum data-point as $1/\sqrt{N_{i}}$, where $N_{i}$ is the number of input modes (i. e., simulated spectra) per data-point $i$. This assumes that data-points and input modes are independent. This is largest for the large DLA power spectrum at $z=2.00$, which has 15,188 input simulated spectra giving an error of 0.81%, and smallest for the forest power spectrum at $z=2.00$, which has 437,063 input simulated spectra, giving an error of 0.15%. All the other uncertainties for each measured power spectrum range in-between these values and can be computed from Table 1.

The total power spectrum deviates from the Lyman-alpha forest power spectrum at all redshifts, showing there is a bias from contamination of spectra by high column density absorbers. This bias can be deconstructed as a function of column density by looking at the power spectra of different absorber categories. The power spectra of high column density absorbers have a distinctive increase on large scales (small $k_{||}$). This is caused by self-correlations across the width of damping wings, which (as discussed in § 2) can be modelled by a Voigt profile (a convolution of a Gaussian and a Lorentzian). Therefore, the power spectrum of high column density absorbers (on large scales) is connected to the Fourier transform of a Voigt profile. This increases for higher column density systems since there is more line broadening, and starts on larger scales for higher column density systems since the damping wings are wider. (See Appendix A for more analysis and discussion of the Voigt profile model.) On small scales, all the power spectra converge to a scaled version of the Lyman-alpha forest flux power spectrum. This reflects the fact that contaminated spectra do contain some uncontaminated spectral pixels. The amplitude of the small-scale power spectrum reflects the fraction of spectra that is uncontaminated, increasing for lower-column density systems since their damping wings are narrower. There is some sensitivity to the length of our simulated spectra, which primarily manifests in our results as the amplitude of the small-scale residual Lyman-alpha forest power in the contaminated power spectra. This is because longer simulated contaminated spectra would have a larger fraction with residual Lyman-alpha forest. This is discussed further and explicitly modelled such that this effect is removed in § 6.

Figure 2 shows 1D flux power spectra as in Fig. 1, but for more of the redshift slices that we consider ($z=[2.00,2.44,3.49,4.43]$), for spectra containing only Lyman-alpha forest and spectra contaminated by large DLAs. The Lyman-alpha forest flux power spectrum has the expected shape, amplitude and redshift evolution, matching observations (e. g., Palanque-Delabrouille et al., 2013) and reflecting the fact that it is an integral of a (biased) matter power spectrum. A peculiarity of the Lyman-alpha forest flux power spectrum is that its amplitude increases with redshift (unlike the linear matter power spectrum); this is because neutral hydrogen is more abundant at higher redshift and so there is more absorption in quasar spectra (i. e., the Lyman-alpha forest becomes a more negatively biased tracer of the matter distribution). By contrast, it can be seen that the large-scale correlations associated with the large DLAs are converging to a single point as redshift changes. This reflects the fact that these correlations arise from individual damping wings, which do not evolve with redshift.

Figure 3 shows the same 1D flux power spectra as in Figs. 1 and 2, but now as ratios between the flux power of spectra contaminated by high column density absorbers and the flux power of spectra containing only Lyman-alpha forest, for $z=2.00$ and $z=4.43$. These ratios are the quantities to which we fit our templates (see § 6) as part of our bias model (see § 3.3). Plotted in this form, it is clear that the large-scale corrections associated with damping wings increase with column density of the damped system. The corrections also decrease with increasing redshift because the Lyman-alpha forest flux power spectrum (on the denominator of the ratio) increases with redshift. On small scales, the ratios converge to a constant value, which reflects the fraction of a line of sight that is uncontaminated (see above). The redshift evolution of this constant value is driven by the transformation from comoving to velocity space: spectra are longer in velocity space at higher redshift (despite being drawn from the same comoving length of the simulation). Conversely, the width of damping wings does not change with redshift (for a given column density) because this width just arises from the physical processes within the hydrogen gas (rather than cosmological evolution). Therefore, the fraction of spectra uncontaminated by the damping wings increases with redshift.

5 Discussion

We first discuss and summarise the results we have presented in § 4. Using our measurements from cosmological hydrodynamical simulations, we have been able to confirm and characterise the effect of high column density absorbers on the 1D Lyman-alpha forest flux power spectrum as a function of column density, scale and redshift. There are distinctive large-scale correlations across the widths of individual damping wings (a “one-halo” term) arising from high column density absorbers that are seen to bias the 1D flux power spectrum of a set of spectra, relative to the power spectrum of the Lyman-alpha forest alone (Fig. 1). These correlations persist for all the high column density absorbing systems that we identify (i. e., for all column densities $N(\textsc{Hi})>1.6\times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$). Our results can be further understood by relating the shape and amplitude of the large-scale power spectrum of spectra contaminated by high column density absorbers to the Fourier transform of the Voigt profile that is normally used to model damping wings (due to the combination of physical effects that broaden absorption lines; see Appendix A). We find evidence in our simulation results that the 1D flux power spectrum of high column density absorbers does not evolve with redshift (Fig. 2). This reflects the fact that the Voigt profiles of damping wings depend only on column density (i. e., the physical processes within high column density absorbing regions) and not redshift (i. e., cosmological evolution) (see Eq. (9)).

The most recent previous investigation into the effect of high column density absorbers on the Lyman-alpha forest was performed by McDonald et al. (2005a) (see also Croft et al., 1999; Viel et al., 2004). These authors measured a single bias function for the 1D Lyman-alpha forest flux power spectrum (at each redshift they consider) that includes the combined effect of all high column density absorbers (i. e., all LLS and DLAs). Our results are qualitatively similar to those of the previous study; however, by investigating different absorber categories based on column density ranges (Table 1), we have shown that the form of the bias as a function of wavenumber depends strongly on column density.

This will have implications for any parameter inference from the 1D flux power spectrum. For instance, the analysis by Palanque-Delabrouille et al. (2015) uses a single multiplicative bias model for the Lyman-alpha forest flux power spectrum based on the results in McDonald et al. (2005a)888The model is reported in Palanque-Delabrouille et al. (2015) as $1-0.2\,\alpha_{\mathrm{DLA}}\,[1/(15000.0\,k_{||}-8.9)+0.018]$, where $\alpha_{\mathrm{DLA}}$ is the free amplitude.. The model has a free amplitude that is allowed to vary (reflecting the level of contamination in a given survey) and is then marginalised. The shape of this model is therefore based on the observed CDDF of high column density absorbers. However, as discussed in § 1, in the measurement of the 1D flux power spectrum, high column density absorbers in the quasar spectra are clipped out in the hope of removing noise (Lee et al., 2013; Palanque-Delabrouille et al., 2013). This process changes the CDDF of high column density absorbers by preferentially removing higher column density systems which are easier to spot in the noisy spectra. Hence, the shape of the bias from residual high column density absorbers is different (as we have shown in § 4) and the model used by Palanque-Delabrouille et al. (2013) may not be flexible enough to account for this, especially at the level of accuracy required by future surveys. Our measurements provide a set of templates for the effect of different absorber categories as a function of column density. By using our templates as part of the model in Eq. (5), it is now possible to more accurately characterise the bias of the residual contamination. We also find evidence for redshift dependence of the fractional effect of high column density absorbers on the forest power spectrum (driven by the changing amplitude of the forest power spectrum), which is also not included in the current model, but is incorporated into our templates. Fits allowing incorporation of our new results into future pipelines are described in § 6.

We now discuss our motivations for some of the choices made in our analysis. We have chosen to present our main results as the 1D flux power spectra of different sets of simulated spectra, where we have categorised spectra according to the maximum column density system within each spectrum. This means that we are measuring the power spectra of ensembles of spectra that are contaminated to similar extents, rather than the flux power spectra of high column density absorbers alone. Furthermore, a consequence of this categorisation is that within the spectra of a given category, there may be high column density absorbers of a lower column density (e. g., there may be LLS in the large-DLA category of spectra). In the first instance, this does not affect our results because the power spectrum measurements we have made (§ 4) and the templates that we construct (§ 6) include the effect of this possible additional lower column density contamination. A subtlety arises because the amount of additional lower column density contamination will be partly sensitive to the length of simulated spectra, since longer spectra have a greater chance of being contaminated. However, the damping wings of the highest column density systems already produce zero transmitted flux over a significant fraction of the length of our simulation box, so that the presence of possible additional high column density absorbers will make very little difference in any case. We tested this by inserting an LLS into a spectrum contaminated by a large DLA, which reduced the total transmitted flux by 0.07%. By carrying out this insertion test with a “control” scenario without the additional contamination, we are able to show that this subtlety will have negligible impact on our conclusions and the validity of our templates.

Finally, we comment on the particular form of our preferred model for the bias of high column density absorbers to the 1D Lyman-alpha forest flux power spectrum (as shown in Eq. 5). We model the bias as a multiplicative correction, rather than e. g., an additive form. First, this matches the form of the currently-used model (as shown in Palanque-Delabrouille et al., 2015). Moreover, an additive form would require either the separation of high column density absorbers and the Lyman-alpha forest in the simulated spectra or a complete physical understanding of how the two components interact at the ensemble level. The former is not trivial for our analysis since we are not inserting high column density absorbers (as previous studies have done), but are simultaneously simulating the low and high column density regions of gas. We avoid the latter due to any remaining physical uncertainties and instead form a parametric multiplicative model based on our simulated results (see § 6).

6 Templates for the effect of high column density absorbers

To aid incorporation in future pipelines, we have produced fits to the biases induced by contaminants in our different column density bins. The parametric form of our templates is

[TABLE]

where

[TABLE]

and the pivot redshift $z_{0}=2.00$. $[a_{0},a_{1},b_{0},b_{1},c_{0},c_{1}]$ are free parameters that we fit simultaneously in $k_{||}$ and $z$ space for each absorber category $i\in$ {LLS, sub-DLA, small DLA, large DLA}. We fit using the Levenberg-Marquardt algorithm (Levenberg, 1944; Marquardt, 1963)999We were able to further validate our modelling by initially fitting using a subset of the available redshift slices and using this preliminary fit to predict the results at $z=3.01$. We found the model to accurately predict the results at this intermediate redshift, acting as a form of successful blind test for our model. Our final best-fit parameters use all available data..

Figure 4 shows the result of these fits (dashed lines) with the raw ratios measured from the simulation (solid lines); the corresponding parameter values are given in Table 2. These can be used to reconstruct a final model for the bias of spectra containing high column density absorbers by using Eq. (5). The model described by Eq. (6) characterises the results we have measured in our simulations and through Eq. (7) allows interpolation of our results to intermediate redshifts that we have not explicitly probed. (Use of the model outside the limits of redshift and scale we have considered would constitute an extrapolation, but this should not be necessary since our measurements bracket the main redshifts and scales of interest to Lyman-alpha forest studies.) No strong physical meaning should be attached to its terms, although we can motivate the first term on the right-hand side of Eq. (6) as being the (reciprocal of the) main term of the redshift evolution of $P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$ as found by Palanque-Delabrouille et al. (2013) (using a maximum likelihood estimator). In this way, the parametric form isolates the redshift evolution from $P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$ and then fits the residual redshift evolution using the terms in Eq. (7). The best-fit values of the exponents in Eq. (7) (as given in Table 2) are small, indicating that most of the redshift evolution can indeed be ascribed to the expected cosmological evolution of $P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$.

Our results are dependent on the length of our simulated spectra. This manifests in the value of the constant that the ratios $P_{i}^{\mathrm{1D}}(k_{||},z)/P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$ have at high $k_{||}$, which is set by the fraction of the length of contaminated spectra which are unaffected by damping wings and contain only Lyman-alpha forest. Since the incidence rates of high column density absorbers are such that one per contaminated spectrum is most likely, a longer spectrum will have a larger fraction that is uncontaminated, causing the constant value at high $k_{||}$ to rise with spectrum length. However, in an analysis of observational data this will be absorbed into a free parameter. We have used a parametric form for our templates such that all this dependency is measured by the term $c(z)$101010It can then be understood why we do not factor out the redshift evolution of $P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$, as we do for the first term on the right-hand side of Eq. (6).. By inserting Eq. (6) into Eq. (5), it can be seen that the term $c(z)$ is degenerate with $\alpha_{\mathrm{Forest}}(z)$ and hence these terms can be combined and allowed to vary. It follows that the full parametric form of our model for the effect of high column density absorbers on the 1D Lyman-alpha forest flux power spectrum is

[TABLE]

When using this model in inference from the 1D Lyman-alpha forest power spectrum $P_{\mathrm{Forest}}^{\mathrm{1D}}(k_{||},z)$, it will be necessary to vary five free parameters $\alpha_{0}$ and $\alpha_{i}$, where $i$ indexes each high column density absorber category. In this way, the column density, scale and redshift dependence of the effect of high column density absorbers is fully determined by our templates, while the relative impact of each absorber category is fitted since this is specific to the survey at hand, as well as the details of any clipping of damping wings that changes the survey CDDF. (See § 5 for more discussion of these details.) Note that the parameter $\alpha_{0}$ is degenerate with factors that rescale the mean flux and could be omitted in an end-to-end analysis.

Figure 5 compares the model we have constructed to the existing model presented in Palanque-Delabrouille et al. (2015) and based on the results in McDonald et al. (2005a). There is broad agreement between the existing model and our model for the total contamination of high column density absorbers, although our model is less steep in its scale dependence. We also show our model applied to a possible “residual” contamination, i. e., under the assumption that all DLAs are identified and clipped out in an analysis, leaving only contamination from LLS and sub-DLAs (e. g., as assumed by Bautista et al. 2017). The model for this lower column density residual contamination has a shallow scale dependence that the model of McDonald et al. (2005a) is unable to characterise. The use of our more flexible model will avoid potential biases due to mischaracterisation of the scale dependence of the residual contamination, thus improving estimation of cosmological effects such as massive neutrinos or the tilt of the primordial power spectrum.

We now discuss the prior probability distributions that can be adopted for $\alpha_{i}(z)$ in any inference using the model we have presented. The $\alpha_{i}(z)$ are technically not independent parameters, but are each related to integrals of the Hi CDDF for a particular survey over the appropriate column density ranges (and absorption distance per sightline). The effect of spectrum clipping which changes the survey CDDF can be modelled by applying a weighting function to the CDDF, which down-weights higher column densities, which are easier to spot and remove. If one wanted to reduce the dimensionality of these nuisance parameters, in particular in redshift space, they could be replaced by a parameterisation which quantifies deviations from the expected redshift evolution of the CDDF with only one or two parameters (rather than a parameter for each redshift bin considered). We leave the details of the construction of prior distributions to individual analyses, since the precise considerations will be survey-specific.

To conclude this section, we present a summary of the steps required to incorporate our final model for the effect of high column denisty absorbers into future 1D Lyman-alpha forest analyses:

• •

Our model describes the effect of quasar spectra contaminated by high column density absorbers as a multiplicative bias to the 1D Lyman-alpha forest flux power spectrum, as given by Eq. (8). It can therefore be incorporated into a pipeline at the stage of flux power spectrum interpretation to marginalise over effects of these absorbers.

• •

The free parameters are $\alpha_{i}(z)$, where $i$ indexes different categories of high column density absorber (as given in Table 1). Our model is of use to any Lyman-alpha forest survey that contains spectra which may be contaminated by high column density absorbers (both LLS and DLAs). The relative impacts of different categories of high column density absorbers will be determined in the estimation of posterior distributions of these nuisance parameters. While normalisation is necessarily floating, the model fully specifies the scale, column density and redshift dependence of the effect of high column density absorbers, using the results we have measured from hydrodynamical simulations.

• •

In a survey which does not clip its quasar spectra, strong priors can be given for the free parameters of our model, based on the expected or measured Hi CDDF.

• •

In a survey which does clip its quasar spectra in an attempt to remove high column density absorbers (and therefore changes the survey CDDF), strong priors can still be given for our model parameters, assuming a model can be constructed for the effect of the clipping process on the CDDF. This will constitute some re-weighting of the CDDF.

• •

In order to reduce the dimensionality of our nuisance parameters, rather than having a separate parameter for each redshift bin in a given analysis, one could parameterise the redshift evolution by a simple deviation from the CDDF with only one or two numbers.

7 Conclusions

We have used a cosmological hydrodynamical simulation (Illustris; Vogelsberger et al., 2014b; Nelson et al., 2015) to investigate the effect of high column density absorbing systems of neutral hydrogen and their associated damping wings on the 1D Lyman-alpha forest flux power spectrum. We find that the effect of high column density absorbers on the Lyman-alpha forest flux power spectrum is a strong function of column density. Accounting for this change in scale-dependence with column density will remove a source of bias in cosmological inference from the Lyman-alpha forest. Previous models (Palanque-Delabrouille et al., 2015) combine the effect of all high column density absorbers together (i. e., all neutral hydrogen column densities $N(\textsc{Hi})>1.6\times 10^{17}\,\mathrm{atoms}\,\mathrm{cm}^{-2}$) based on the column density distribution function (CDDF) in the raw spectra (McDonald et al., 2006). However, the damping wings of some high column density absorbers are clipped out in the final analysis (Lee et al., 2013), which preferentially removes higher density systems (because they are easier to spot) and changes the column density distribution in the residual contamination. Our results apply for both clipped and unclipped survey spectra, since we separately model the effect for different column densities of the dominant absorber, allowing us to accurately account for the contamination in the 1D flux power spectrum. We discuss in § 6 the practicalities of employing our model in future analyses.

The shape and amplitude of the distortions in the power spectrum due to a damped absorber depend on its column density because they are driven by the width of the damping wings; i. e., the dominant effect is a “one-halo” term. We defer investigation of potential “two-halo” terms to future work where we measure the effect of high column density absorbers on the 3D Lyman-alpha forest flux power spectrum.

We anticipate that our model will help realise forecasted cosmological constraints from upcoming surveys like DESI. E. g., Font-Ribera et al. (2014) forecast that DESI will have the constraining power to make a $\sim$ three-sigma detection of the sum of neutrino masses (in combination with Planck CMB data); and they show the power of the 1D Lyman-alpha forest power spectrum in probing the primordial power spectrum, e. g., halving the one-sigma error on the running of the spectral index, with implications for inflationary models. It will be necessary to use the models we have presented here, alongside carefully constructed priors on the residual CDDF, to remove degeneracies between the effect of high column density absorbers and cosmological effects.

Acknowledgements

KKR, SB, HVP and BL thank the organisers of the COSMO21 symposium in 2016, where this project was conceived. KKR was supported by the Science and Technology Facilities Council (STFC). SB was supported by NASA through Einstein Postdoctoral Fellowship Award Number PF5-160133. HVP was partially supported by the European Research Council (ERC) under the European Community’s Seventh Framework Programme (FP7/2007-2013)/ERC grant agreement number 306478-CosmicDawn. AP was supported by the Royal Society. AFR was supported by an STFC Ernest Rutherford Fellowship, grant reference ST/N003853/1. BL was supported by NASA through Einstein Postdoctoral Fellowship Award Number PF6-170154.

Appendix A One-dimensional flux power spectrum of a Voigt profile

As discussed in § 2, the broadened absorption lines of high column density absorbers are usually modelled by a Voigt profile. A Voigt profile is a convolution of a Lorentzian profile and a Gaussian profile. It therefore appropriately models the combination of the main physical processes that broaden atomic transition lines: the Lorentzian profile from e. g., natural or collisional broadening and the Gaussian profile from e. g., Doppler broadening. The optical depth as a function of wavelength $\tau(\lambda)$ is the product of the line-of-sight column density $N$ and the atomic absorption coefficient $\alpha(\lambda)$ (e. g., Humlicek, 1979)111111Eq. (9) is valid in SI units.:

[TABLE]

where the fundamental physical constants have their usual meaning, $f$ is the oscillator strength of the atomic transition, $\lambda_{\mathrm{t}}$ is the transition wavelength and the Doppler wavelength “shift” associated with a gas of temperature $T$ for an ion of mass $m_{\mathrm{ion}}$,

[TABLE]

$u(x,y)$ is an unnormalised form of the Voigt function (the normalisation is already expressed in the pre-factors of Eq. (9)), specifically the real part of the Faddeeva function:

[TABLE]

where $\mathrm{erfc}(x)$ is the complementary error function and $z=x+iy$. $x$ and $y$ are respectively the wavelength difference from the line centre $\lambda_{\mathrm{c}}$ and the natural width of the transition, in units of the Doppler shift:

[TABLE]

where $\Gamma$ is the damping constant of the transition, i. e., the inverse of the time scale for the electron to remain in the upper level of the transition in the vacuum. For the Lyman-alpha transition, $f=0.4164$, $\lambda_{\mathrm{t}}=1215.67\textup{\AA}$, $m_{\mathrm{ion}}=m_{\mathrm{proton}}$ and $\Gamma=6.265\times 10^{8}\,\mathrm{Hz}$ (Morton, 2004). For the column densities that we consider, we assume a gas temperature $T\approx 10^{4}K$. In order to calculate the 1D flux power spectrum arising from these Voigt profiles, the same procedure is followed as in § 3.2, i. e., we form flux spectra and carry out a Fourier transform. We transform from wavelengths to velocities by $\Delta v/c=\Delta\lambda/\lambda$.

Figure 6 shows the 1D flux power spectra of Voigt profiles as given by Eq. (9) for the Lyman-alpha absorption line for three different column densities of neutral hydrogen $N(\textsc{Hi})=[10^{19},10^{20},10^{21}]\,\mathrm{atoms}\,\mathrm{cm}^{-2}$, spanning the column densities for LLS and DLAs. This figure should be compared with Fig. 1 in § 4, which shows the 1D flux power spectra we have measured in the hydrodynamical simulations. The trends in Fig. 6 broadly support the arguments made in § 5, relating the large-scale power spectrum of simulated spectra contaminated by high column density absorbers to the power spectrum of relevant Voigt profiles. The shape of the large-scale power spectrum of the Voigt profiles is similar in amplitude and scale-dependence as the excesses on large scales for the 1D flux power spectra of simulated spectra in high column density absorber categories. Moreover, these excesses get steeper, increase in amplitude and become prominent on larger scales for higher column densities, both in the simulated and analytic spectra. This reflects the fact that a higher column density means wider damping wings and so correlations on larger scales. In the analytic power spectra in Fig. 6, we observe oscillations in the power spectrum on smaller scales that rapidly decrease in amplitude. These are not observed in the fully-simulated power spectra since the oscillations are orders of magnitude lower in amplitude than the flux power spectrum of residual Lyman-alpha forest (see Fig. 1). Furthermore, in our results, we are effectively averaging over a number of column densities in each column density bin (or absorber category) that we consider; this will have the additional effect of averaging out these smaller-scale oscillations in the power spectrum to form a smoother scaling.