Neutrino signatures on the high-transmission regions of the Lyman |$\boldsymbol {\alpha }$| forest

ABSTRACT

We quantify the impact of massive neutrinos on the statistics of low-density regions in the intergalactic medium as probed by the Lyman α forest at redshifts z = 2.2–4. Based on mock but realistic quasar (QSO) spectra extracted from hydrodynamic simulations with cold dark matter, baryons and neutrinos, we find that the probability distribution of weak Lyman α absorption features, as sampled by Lyman α flux regions at high transmissivity, is strongly affected by the presence of massive neutrinos. We show that systematic errors affecting the Lyman α forest reduce but do not erase the neutrino signal. Using the Fisher matrix formalism, we conclude that the sum of the neutrino masses can be measured, using the method proposed in this paper, with a precision smaller than 0.4 eV using a catalogue of 200 high-resolution (signal-to-noise ratio ∼100) QSO spectra. This number reduces to 0.27 eV by making use of reasonable priors in the other parameters that also affect the statistics of the high-transitivity regions of the Lyman α forest. The constraints obtained with this method can be combined with independent bounds from the cosmic microwave background, large-scale structures and measurements of the matter power spectrum from the Lyman α forest to produce tighter upper limits on the sum of the masses of the neutrinos.

neutrinos, intergalactic medium, quasars: absorption lines, cosmology: theory, large-scale structure of Universe

1 INTRODUCTION

Neutrino oscillation experiments revealed that neutrinos are not massless particles. Since then a major effort has been dedicated to measure or constrain neutrino masses. Current laboratory bounds constrain the electron neutrino mass to |$m_{\nu _{\rm e}}<2.05$| eV (Lobashev 2003; Kraus et al. 2005). Cosmological bounds for the sum of all neutrino masses are still significantly stronger: constraints from 7-year Wilkinson Microwave Anisotropy Probe (WMAP7) alone yield |$\Sigma _im_{\nu _i}<1.3$| eV (Komatsu et al. 2009), while combined with large-scale structure (LSS) measurements they constraint the mass to |$\Sigma _im_{\nu _i}<0.3$| eV (Wang et al. 2005; Gonzalez-Garcia, Maltoni & Salvado 2010; Reid et al. 2010; Thomas, Abdalla & Lahav 2010; de Putter et al. 2012; Xia et al. 2012). The tightest 2σ upper limit of |$\Sigma _im_{\nu _i}<0.17$| eV is obtained by combining cosmic microwave background (CMB) results, LSS and Lyman α forest (Seljak, Slosar & McDonald 2006) data sets (see Abazajian et al. 2011, for a summary of current and future neutrino mass constrains). Among all the different observables the Lyman α forest is particularly constraining since it probes structures over a wide range of redshift, in a mildly non-linear regime and at small scales where the neutrino signature is present (Lesgourgues & Pastor 2006; Rauch 1998; Meiksin 2009).

The dynamics of cosmological neutrinos is very different from that of the dominant cold dark matter (CDM) component. The large velocity dispersion of neutrinos suppresses their power spectrum of density fluctuations at small scales, making the shape of the total power spectrum a potential probe of neutrino masses.

Previous studies have addressed the role of neutrinos in dark matter haloes (Singh & Ma 2003; Ringwald & Wong 2004; Villaescusa-Navarro et al. 2011), LSS (Ma & Bertschinger 1994; Brandbyge et al. 2008, 2010; Brandbyge & Hannestad 2009, 2010; Marulli et al. 2011; Bird, Viel & Haehnelt 2012) and the intergalactic medium (IGM; Viel, Haehnelt & Springel 2010), using both linear theory and N-body/hydrodynamic techniques for the non-linear regime. It has been shown that on scales of 1–10 h⁻¹ Mpc the non-linear suppression is redshift, scale and mass dependent in a way that is different from a naive extrapolation of linear theory.

In this paper we study the effect of massive neutrinos on the properties of low-density regions or voids in the IGM. Neutrinos have only a mild effect on dark matter haloes (Singh & Ma 2003; Ringwald & Wong 2004; Villaescusa-Navarro et al. 2011), since their large velocity dispersion prevents their clustering on small scales. In contrast, we find that the impact of neutrinos on void properties is much stronger. Voids are relatively empty regions with |$\delta =\rho _{\rm m}/\bar{\rho }_{\rm m}-1$| ranging from almost −1 in their cores to ∼−0.7 at radii 10–20 Mpc at z = 0 (Colberg et al. 2008). By solving the dynamical equations for an isolated spherical top-hat underdense perturbation, we find that neutrinos modify the evolution of underdense regions by making them smaller and denser. Neutrinos contribute to the interior mass of the underdense region delaying the rate at which CDM is being evacuated from its interior and slowing down the velocity of the shell surrounding it. We find that the linearly extrapolated density contrast when the underdense region enters into its non-linear phase decreases by ∼10 per cent for neutrinos with |$\Sigma _im_{\nu _i}\sim 1$| eV. Using the analytic model presented in Sheth & van de Weygaert (2004) we find that the statistics of voids depend on both σ₈ and |$\Sigma _i m_{\nu _i}$|⁠. Lyman α voids and their dependence on other cosmological parameters have been investigated in Viel, Colberg & Kim (2008). Here we focus on the dependence of void properties on the sum of the neutrino masses. We consider the Lyman α signature of low-density regions, and introduce a new and simple statistical tool (note that a similar observable was already studied in Miralda-Escudé et al. 1996; Fan et al. 2002) that samples most of the IGM volume and appears to be highly sensitive to neutrino masses.

2 NUMERICAL METHOD

Our mock quasar spectra are based on cosmological simulations run with the TreePM-SPH code gadget-3 (Springel 2005). The code has been extended to include neutrinos either by solving their potential on the mesh or by representing them as discrete particles (Viel et al. 2010). Here, we use primarily the first implementation and refer the reader to Brandbyge et al. (2008), Brandbyge & Hannestad (2009) and Marulli et al. (2011) for a critical comparison of the two methods and also for comparison with another new method (Ali-Haïmoud & Bird 2012). Our simulations consist of 2 × 512³ CDM plus gas particles sampling a periodic box of 512 h⁻¹ Mpc. We adopt a flat ΛCDM background with cosmological parameters Ω_CDM + Ω_ν = 0.25, |$\Omega _\Lambda =0.7$|⁠, Ω_b = 0.05, h = 0.7 and n_s = 1. We consider three degenerates neutrino species with a total mass of |$\Sigma _i m_{\nu _i}=0.0,\ 0.3$| and 0.6 eV. The initial power spectra of most of our simulations, produced with camb,¹ are normalized for all neutrino masses at a wavenumber 2 × 10⁻³h Mpc⁻¹, corresponding to the scale constrained by CMB data. This produces different values of σ₈ = 0.877, 0.806 and 0.732 at z = 0 for the models with |$\Sigma _i m_{\nu _i}=0.0,\ 0.3$| and 0.6 eV, respectively (see Fig. 1). Our initial conditions are generated at z = 49. A summary of the different simulations we carried out is shown in Table 1.

Figure 1.

Linear matter power spectra for different neutrino masses at z = 0 (upper lines) and z = 3 (bottom lines). The inner panel shows the linear (solid lines) and non-linear (dotted lines) matter power spectrum at z = 3 normalized to the case without neutrinos. On scales k < 0.03 h Mpc⁻¹ changes in the non-linear power spectra are driven by the differences in the linear power spectra.

Table 1.

Summary of the N-body/hydrodynamic parameters of the simulations. The cosmological parameters are the same for all simulations and are given in the text. Ω_M = Ω_cdm + Ω_b + Ω_ν is kept constant. N_CDM, N_b and N_ν correspond to the number of CDM, baryon and neutrino particles, respectively. All the simulations except P0 and P6 are based on the Fourier space implementation of neutrinos (see text).

Name	|$\Sigma _i m_{\nu _i}$| (eV)	σ8 (z = 0)	|$N_\mathrm{CDM}^{1/3}$|	|$N_\mathrm{b}^{1/3}$|	|$N_\nu ^{1/3}$|
S0	0.0	0.877	512	512	0
S0+	0.0	0.928	512	512	0
S0−	0.0	0.828	512	512	0
S3	0.3	0.948	512	512	0
S3+	0.3	0.877	512	512	0
S3−	0.3	0.807	512	512	0
S5	0.5	0.755	512	512	0
S6	0.6	0.732	512	512	0
S6+	0.6	0.877	512	512	0
S7	0.7	0.709	512	512	0
LR0	0.0	0.877	448	448	0
LR6	0.6	0.732	448	448	0
P0	0.0	0.877	512	512	512
P6	0.6	0.732	512	512	512

For each simulation we consider snapshots at redshifts z = 2.2 and 4 that bracket the range of interest for the observed Lyman α forest in quasar spectra from ground-based telescopes. For each snapshot we sample 4500 random line-of-sights (RLOSs) uniformly distributed along each x, y or z direction. For each RLOS we extract the baryon density contrast |$\rho _{\rm b}(r)/\overline{\rho }_{\rm b}$| and the peculiar velocity V_p(r) along the line-of-sight and then compute the transmitted flux e^−τ(u) in redshift space (with u in km s⁻¹), where τ is the Lyman α optical depth, by using the fluctuating Gunn–Peterson approximation:

\begin{equation} \tau (u)=A\int _{-\infty }^{+\infty }\mathrm{d}x\, \delta [u-x-V_{\rm p}(x)]\ \bigg (\frac{\rho _{\rm b}(x)}{\overline{\rho }_{\rm b}}\bigg )^{1.6}, \end{equation}

(1)

where x = H(z)r/(1 + z) is the redshift-space coordinate and A is a factor that depends on the global thermal history of the IGM (Croft et al. 2002):

\begin{eqnarray} A&=&0.433\bigg (\frac{1+z}{3.5}\bigg )^6\bigg (\frac{\Omega _{\rm b}h^2}{0.02}\bigg )^2\bigg (\frac{0.65}{h}\bigg )\bigg (\frac{3.68H_0}{H(z)}\bigg )\nonumber \\ &&\times \bigg (\frac{1.5\times 10^{-12}\,\mathrm{s}^{-1}}{\Gamma _{\rm H\,{\small I}}}\bigg )\bigg (\frac{6000\,{\rm K}}{T_0}\bigg )^{0.7}, \end{eqnarray}

(2)

with |$\Gamma _{\rm H\,{\small I}}$| being the hydrogen photoionization rate. The power-law index in the scaling with |$\rho _{\rm b}/\overline{\rho }_{\rm b}$| arises from the equation of state for the IGM temperature, |$T=T_0(\rho _{\rm b}/\overline{\rho }_{\rm b})^\alpha$| (Hui & Gnedin 1997), with α ≈ 0.6. In all our calculations we adopt T₀ = 10⁴ K and choose |$\Gamma _{\rm H\,{\small I}}$| such that the mean flux over the whole set of RLOS reproduce the observed mean flux at redshift z (Miralda-Escudé et al. 1996) |$\langle F \rangle ={\rm e}^{-\tau _{\rm eff}(z)}$| with τ_eff(z) = 0.0023(1 + z)^3.65 (Kim et al. 2007). We neglect the effects of thermal broadening. Finally, we smooth the flux over a scale of 1 h⁻¹ Mpc, which is larger than the Jeans length, to avoid sensitivity to substructure below the Jeans scale, which is affected by numerical resolution and astrophysical processes (e.g. feedback from galactic winds).

Fig. 2 shows the baryon density contrast, |$\rho _{\rm b}/\overline{\rho }_{\rm b}$|⁠, and peculiar velocity, V_p, extracted along a RLOS as a function of the comoving coordinate r together with the corresponding transmitted flux F = e^−τ in redshift space, plotted in terms of the mean flux at redshift z.

Figure 2.

Real-space distribution of the baryon density contrast, |$\rho _{\rm b}/\overline{\rho }_{\rm b}$| (top panel), and the peculiar velocity, V_p (middle), along a RLOS. In the bottom panel we plot the transmitted flux F = e^−τ, in units of the mean flux 〈F〉, in redshift space.

3 ANALYSIS OF THE SIMULATIONS

We focus our analysis on the statistical properties of low-density regions that are expected to produce weak absorption features. We define as void ‘region’ a continuous domain in the transmitted flux profile which remains always above a given threshold. The higher the threshold, the lower the absorption in that region. For each RLOS we extract the transmitted flux from |$\rho _{\rm b}/\overline{\rho }_{\rm b}$| and V_p and count the number of regions above the selected threshold. This results in a statistical estimate of the low absorption contribution to the Lyman α signal, and allows us to quantify the impact of neutrinos on those regions.

In Fig. 3 we plot the probability distribution function (PDF) for the number of regions per path length of 100 h⁻¹ Mpc² above a threshold of F/〈F〉 = 1.14 at redshift z = 2.2 (top) and at redshift z = 4.0 for a threshold F/〈F〉 = 1.70 (bottom) for three different neutrino masses, |$\Sigma _i m_{\nu _i}=0.0,\ 0.3,\ 0.6$| eV. The two numbers above are chosen taken into account two competing effects: on one side, the larger the value of F/〈F〉, the larger the differences between the models. On the other side, the number of regions above the threshold drops rapidly as F/〈F〉 increases, requiring a larger QSO spectra catalogue to obtain converged results. By choosing the numbers above we make sure that differences are large enough having converged results. We have verified that these PDFs do not change if we increase the number of RLOS, i.e. our statistical sample of RLOS is large enough to reliably measure the PDF. Fig. 3 shows that the neutrino mass has a significant impact on the mean of the distributions. In Fig. 4 we plot the mean of the distributions, i.e. the average number of regions per path length of 100 h⁻¹ Mpc above a given threshold, as a function of the threshold for the three different neutrino masses (⁠|$\Sigma _i m_{\nu _i}=0.0,\ 0.3,\ 0.6$| eV) at redshift z = 2.2 (top) and z = 4.0 (bottom). This shows clearly that the higher the threshold, the larger the differences between the various neutrino cosmologies. This is the expected neutrino signature as we discuss below.

Figure 3.

PDF for the number of regions per path length of 100 h⁻¹ Mpc above a threshold of F/〈F〉 = 1.14 (top), 1.70 (bottom) as a function of |$\Sigma _i m_{\nu _i}$| and σ₈ at z = 2.2 (top) and z = 4 (bottom). The PDFs have long tails with a very low probability that extend up to 10–12. The σ₈ − Ω_ν degeneracy is not perfect and can be broken by studying the spectra at different redshifts.

We explicitly checked that relative differences between our neutrino models are numerically converged against mass and spatial resolution. Furthermore, we used the neutrino particle implementation (simulations P0 and P6 in Table 1) and found the same trends in the neutrino signature as with the grid method, although relative differences between the different models are even slightly larger when we use the particle implementation. This is due to the fact that non-linear neutrino effects, such as phase mixing, are only properly captured by using the particle implementation. However, we note that the grid implementation in the mildly non-linear Lyman α regime is fully justified since non-linear neutrinos effects should not be particularly important at those redshifts and at k < 1 h Mpc⁻¹.

4 SYSTEMATIC ERRORS

The high-transmissivity regions of the Lyman α forest are prone to systematic errors such as those induced by the continuum fitting procedure and the signal-to-noise ratio (S/N) of the QSO spectrum. We have investigated whether these effects are able to spoil the neutrino signal. In particular, we have studied how the presence of systematic errors affects the differences between models (in terms of the average number of regions above the threshold at z = 2.2), and their impact on the sensitivity to Ω_ν of a catalogue containing 200 QSO spectra.

We first focus our attention to the case without systematic errors. In the subplot of Fig. 4 we show the average number of regions per path length of 100 h⁻¹ Mpc as a function of the threshold at z = 2.2 normalized to the neutrinoless model. The black error bars show the 90 per cent (interior tick marks) and 99 per cent (exterior tick marks) confident intervals for a mock catalogue consisting of 200 RLOS taken from the simulation with |$(\Sigma _i m_{\nu _i}=0.0$| eV, σ₈ = 0.877). We find that with a catalogue consisting of 200 QSO spectra (not affected by systematic errors) we can rule out models |$(\Sigma _i m_{\nu _i}=0.3$| eV, σ₈ = 0.806) and |$(\Sigma _i m_{\nu _i}=0.6$| eV, σ₈ = 0.732) with a high significance. Distinguishing models with the same σ₈ would require a larger QSO spectra catalogue and combining results at different redshifts in order to disentangle the different redshift evolution of the models.

Figure 4.

Average number of regions per path length of 100 h⁻¹ Mpc as a function of flux threshold at redshift z = 2.2 (top) and z = 4 (bottom) for different neutrino masses and σ₈. The subplot in the upper panel shows the ratio between models with |$\Sigma _i m_{\nu _i}\ne 0.0$| and the model with |$\Sigma _i m_{\nu _i}=0.0$|⁠. The black error bars indicate the 90 per cent (interior tick marks) and 99 per cent (exterior tick marks) confident intervals for a mock catalogue consisting of 200 RLOS taken from the simulation with |$(\Sigma _i m_{\nu _i}=0.0$| eV, σ₈ = 0.877). Models with |$\Sigma _i m_{\nu _i}=0.3,\ 0.6$| and σ₈ = 0.806, 0.732, respectively, can be ruled out with a high significance by using a catalogue of 200 QSO spectra.

In order to study how systematic errors impact on our results, we have created a realistic mock catalogue of high-resolution QSO spectra. We have mimicked the continuum fitting errors by rescaling each flux pixel by the quantity F_i/F_max (F_max being the maximum value of the transmitted flux along the spectrum) as done by McDonald et al. (2000); we considered a S/N = 100 which is reasonable for Ultraviolet and Visual Echelle Spectrograph (UVES)/Very Large Telescope (VLT) QSO spectra. Note that this treatment of the systematic errors induced by the continuum fitting is rather simplistic and conservative, and ideally one would like to put more refined models for the QSO continuum and generate a set of mock spectra that would be as close as possible to the observed one (e.g. by including the redshift distribution of the sources), however, this is beyond the scope of the present paper.

Using this catalogue, we have repeated the analysis described above and we show the results in Fig. 5. Two things can be pointed out from this figure. On one side, we find that the mean number of regions above a flux threshold is strongly affected by the Gaussian noise and by the bin size in the QSO spectra. This dependence can be easily understood considering that the Gaussian noise on a pixel can divide a single region above a threshold into two, and also by the fact that spurious regions will appear because the Gaussian noise can increase the value of the transmitted flux (ending up with a final value above the threshold) on a pixel which is below the threshold. On the other side, it turns out that, as expected, the differences between models become smaller as the QSO catalogue S/N drops. However, the remaining differences between models points out that neutrino effects are not erased by the presence of the systematics discussed in this section. The subplot of Fig. 5 shows the ratio between the different models and the model with |$\Sigma _i m_{\nu _i}=0.0$|⁠. The black error bars indicate the 90 per cent (interior tick marks) and 99 per cent (exterior tick marks) confident intervals for a mock catalogue consisting of 200 high-resolution QSO spectra extracted from the simulation with |$(\Sigma _i m_{\nu _i}=0.0$| eV, σ₈ = 0.877). We find that with a high-resolution catalogue (S/N = 100) consisting of 200 QSO spectra, models with |$\Sigma _i m_{\nu _i}=0.0$| and σ₈ = 0.877 and |$\Sigma _i m_{\nu _i}=0.6$| and σ₈ = 0.732 can be distinguished at a very high significance. The models with |$\Sigma _i m_{\nu _i}=0.0$| and σ₈ = 0.877 and |$\Sigma _i m_{\nu _i}=0.3$| and σ₈ = 0.806 can be distinguished with a lower significance, while in order to distinguish models with the same σ₈, a catalogue with either more QSO spectra or higher S/N value is needed.

Figure 5.

Effects of the Gaussian noise in the QSO spectra and the continuum fitting procedure on the properties of the high-transmission regions of the Lyman α forest. We create high-resolution (S/N = 100) mock QSO spectra mimicking the continuum fitting errors by using the prescription of McDonald et al. (2000). In the figure we plot the average number of regions per path length of 100 h⁻¹ Mpc as a function of threshold at redshifts z = 2.2 (top) and z = 4.0 (bottom) for different neutrino masses and σ₈. The subplot shows the ratio between models with |$\Sigma _i m_{\nu _i}\ne 0.0$| and the model with |$\Sigma _i m_{\nu _i}=0.0$| at z = 2.2. The black error bars indicate the 90 per cent (interior tick marks) and 99 per cent (exterior tick marks) confident intervals for a mock catalogue consisting of 200 high-resolution QSO spectra taken from the simulation with |$(\Sigma _i m_{\nu _i}=0.0$| eV, σ₈ = 0.877). Systematic errors do not erase the neutrino signal, however, the differences between models become smaller. At z = 2.2, only the model with |$\Sigma _i m_{\nu _i}=0.6$| eV and σ₈ = 0.732 can be ruled out with a high significance using the catalogue considered in this example.

In addition to the systematic errors arising from the noise in the spectra and the continuum fitting procedure, the presence of metal lines should not affect strongly our findings since in high-resolution spectra these lines can be identified and metal-free regions can be conservatively used in the analysis and, by smoothing the transmitted flux over a region which is typically ∼1 com. h⁻¹ Mpc (roughly 20 times larger than the typical width of a metal line), we should be less sensitive to these contaminants.

5 SENSITIVITY TO NEUTRINO MASSES USING THE FISHER MATRIX FORMALISM

In this section we quantify the sensitivity to Ω_ν of the observable described in this paper.

From Figs 3–5 it is clear that both Ω_ν and σ₈ impact on the statistics of the high-transmission regions of the Lyman α forest. Variations in the mean flux in the thermal and ionization history of the IGM are also expected to impact on the Lyman α properties of large size voids (Viel et al. 2008). For real data, the S/N is not known with infinite precision (a 10 per cent error could be a reasonable conservative assumption), and as can be seen when comparing Figs 4 and 5, variations in it are expected to strongly impact on the statistics of low-density regions. The scale over which we smooth the transmitted flux spectra will also impact on our results. However, we have full control on this scale that we can set to an any particular value. Therefore we do not need to consider this parameter in this analysis.

We use the Fisher matrix formalism to study how the above parameters affect the statistics of the high-transmissivity regions of the Lyman α forest and with which error can Ω_ν be constrained by using the method proposed in this paper. The parameters, |$\boldsymbol {p}=(p_1,p_2,p_3,p_4)$|⁠, we use in the analysis are the sum of the neutrino masses, |$\Sigma _i m_{\nu _i}$|⁠, the catalogue mean transmitted flux, 〈F〉, the parameter α present in the IGM temperature–density relation, |$T=T_0(\rho _{\rm b}/\overline{\rho }_{\rm b})^{\alpha }$|⁠, and σ₈. Note that the parameters T₀ and |$\Gamma _{\rm H\,{\small I}}$| affect the properties of the Lyman α forest only through 〈F〉, and therefore, we do not need to include them in the analysis. For the realistic catalogue, the one that incorporates the systematic errors, we also need to consider a further parameter, p₅, which corresponds to the catalogue S/N.

We carry out the Fisher matrix analysis for two different catalogues. In the first one we consider a catalogue consisting in 200 QSO spectra with no continuum fitting errors and with an S/N equal to infinite. This catalogue, although unrealistic, help us to determine the tightest constrains on Ω_ν we can achieve by using the method here proposed. For the second one, we use a 200 QSO spectra catalogue with S/N = 100. The continuum fitting errors are mimicked using the prescription of McDonald et al. (2000) (see Section 4). In both catalogues, we smooth the transmitted flux spectra over a scale of 1 h⁻¹ com. Mpc, which is comparable to the Jeans scale. This smoothing is intended in order to be less sensitive to the substructure below this scale and to the noise level properties and can be applied on both simulations and real spectra in exactly the same way.

Our observables, f_b, with b = 0, 1, 2, …, 29, correspond to the average number of regions per path length of 100 h⁻¹ Mpc above a threshold equal to F_b = 0.90 + 0.09b/29. We assume flat priors on the parameters, and therefore, the posterior distribution is equal to the likelihood. The likelihood associated with each model will be given by

\begin{equation} \mathcal {L}\propto {\rm {exp}}\left [-\frac{1}{2} \displaystyle \sum _{a=0}^{29} \displaystyle \sum _{b=0}^{29} \left(f_a-\widetilde{f}_a(\boldsymbol {p})\right)\bf{C}^{-1}_{ab} \left(f_b-\widetilde{f}_b(\boldsymbol {p})\right)\right ], \end{equation}

(3)

where |$\widetilde{f}_b(\boldsymbol {p})$| is the theoretical prediction for f_b for a model with parameters p_i; |$\bf{C}$| is the covariance matrix:

\begin{equation} \bf{C}_{ab}=\langle (f_a-\widetilde{f}_a) (f_b-\widetilde{f}_b)\rangle . \end{equation}

(4)

Let |$\boldsymbol {p}^0$| be the values that maximize the likelihood. Around |$\boldsymbol {p}^0$| we can expand |$\rm {ln}\ \mathcal {L}$| in a Taylor series as

\begin{equation} {\rm {ln}}\, \mathcal {L}={\rm {ln}}\, \mathcal {L}(\boldsymbol {p}^0)+\frac{1}{2}\displaystyle \sum _{ij}(p_i-p_i^0)\frac{\mathrm{\partial} ^2{\rm {ln}}\, \mathcal {L}}{\mathrm{\partial} p_i \mathrm{\partial} p_j}\left(p_j-p_j^0\right)+\cdots . \end{equation}

(5)

Note that because the likelihood has a maximum in |$\boldsymbol {p}^0$|⁠, |$\mathrm{\partial} {\rm {ln}}\, \mathcal {L}/\mathrm{\partial} p_i(\boldsymbol {p}^0)=0$|⁠. The Fisher matrix is defined as

\begin{equation} F_{ij}=\left\langle \frac{\mathrm{\partial} ^2{\rm {ln}}\, \mathcal {L}}{\mathrm{\partial} p_i \mathrm{\partial} p_j}\right\rangle , \end{equation}

(6)

and the marginalized error in the parameter p_i satisfies the relation

\begin{equation} \sigma _{p_i}\geqslant \sqrt{\bf{F}_{ii}^{-1}}. \end{equation}

(7)

If the data are distributed according to a Gaussian distribution, the Fisher matrix can be written in the following way (Tegmark, Taylor & Heavens 1997; Heavens 2009):

\begin{eqnarray} \bf{F}_{ij}=\frac{1}{2}{\rm {Tr}}[\bf{C}^{-1}\mathrm{\partial} _i\bf{CC}^{-1}\mathrm{\partial} _j\bf{C}&+&\bf{C}^{-1}\mathrm{\partial} _i \widetilde{f}_b(\boldsymbol {p}) \mathrm{\partial} _j (\widetilde{f}_b(\boldsymbol {p}))^{\rm T} \end{eqnarray}

(8)

\begin{eqnarray} &&\qquad\qquad\qquad\qquad\qquad\,\,\,\,\,+\,\bf{C}^{-1}\mathrm{\partial} _j \widetilde{f}_b(\boldsymbol {p})\mathrm{\partial} _i (\widetilde{f}_b(\boldsymbol {p}))^{\rm T}], \end{eqnarray}

(9)

\begin{equation} \bf{F}_{ij}=\frac{1}{2}{\rm {Tr}}[\bf{C}^{-1}\mathrm{\partial} _i\bf{CC}^{-1}\mathrm{\partial} _j\bf{C}]+\displaystyle \sum _{a=0}^{29} \displaystyle \sum _{b=0}^{29} \frac{\mathrm{\partial} \widetilde{f}_a}{\mathrm{\partial} _{p_i}}\bf{C}^{-1}_{ab}\frac{\mathrm{\partial} \widetilde{f}_b}{\mathrm{\partial} _{p_j}}, \end{equation}

(10)

where ∂_i = ∂/∂p_i and Tr is the matrix trace. We find that the first term in the previous equation contributes very little to the Fisher matrix. This happens because the errors (the covariance matrix) depend very weakly of the parameters |$\boldsymbol {p}$|⁠. For that reason, we have neglected that term in our analysis.³

The values of the parameters for our fiducial model are |$p_1^0=\Sigma _i m_{\nu _i}=0.3$| eV, |$p_2^0=\langle F \rangle =0.8517$|⁠, |$p_3^0=\alpha =0.5714$| and |$p_4^0=\sigma _8=0.877$|⁠. For the realistic catalogue, the fifth parameter has a fiducial value equal to |$p_5^0={\rm S/N}=100$|⁠. We first computed the Fisher matrix for a catalogue with a S/N = ∞, i.e. for the ideal case, unaffected by systematic errors, that we consider in the body of this paper. Without using priors on the parameters, we find that the marginalized error in |$\Sigma _i m_{\nu _i}$| is roughly given |$0.30\sqrt{200/N}$| eV, where N is the number of QSO spectra in the catalogue. In the left-hand column of Fig. 6 we show the contour plots at 1σ (blue) and 2σ (red) for a catalogue consisting of 200 QSO spectra with S/N = ∞.

Figure 6.

Contour plots at 1σ (blue) and 2σ (red) showing the correlation between |$\Sigma _i m_{\nu _i}$| and α, 〈F〉, σ₈ and S/N for a catalogue consisting in 200 QSO spectra for three different situations: a catalogue with S/N = ∞ and without continuum fitting errors, assuming no priors on the parameters (left-hand column), a catalogue with S/N = 100 and without assuming priors on the parameters (middle column) and a catalogue with S/N = 100 assuming priors of 10 per cent in the value of α and 3 per cent in the value of σ₈ (right-hand column). The black points show the position of the fiducial model.

We have also studied the sensitivity of our method making one further assumption: we suppose that the amplitude of the matter power spectra is fixed on large scales. By making that assumption, Ω_ν and σ₈ are no longer independent parameters. We have carried out the Fisher matrix analysis taking as parameters |$\boldsymbol {p}=(\Sigma _i m_{\nu _i},\langle F \rangle , \alpha )$| and we found that the marginalized error in |$\Sigma _i m_{\nu _i}$| is |${\sim }0.25\sqrt{200/N}$| eV. Although we have one parameter less than in the previous case, the marginalized error in Ω_ν is not significantly reduced. This happens because we find a strong correlation between the parameters |$\Sigma _i m_{\nu _i}$| and α. This correlation is the main source of error when computing the marginalized error in |$\Sigma _i m_{\nu _i}$|⁠. We note that this correlation is expected and has a physical meaning: a larger value of α makes the low-density IGM colder and thereby would result in a larger amount of neutral hydrogen. In such a case the voids will contain more matter and this is the same effect that can be achieved in a universe with a larger value of |$\Sigma _i m_{\nu _i}$|⁠. We note that this huge degeneracy between |$\Sigma _i m_{\nu _i}$| and α is broken once we introduce σ₈ in the Fisher matrix. This can be understood by looking at Fig. 4: the mean number of regions above the threshold varies differently depending on whether σ₈ is fixed (compare red and green lines) or the amplitude of the power spectra is fixed on large scales (compare red, blue and magenta lines).

By repeating the Fisher matrix analysis with the 200 high-resolution QSO spectra catalogue (S/N ∼ 100), without assuming priors on the parameters, we find that the marginalized error in |$\Sigma _i m_{\nu _i}$| is given by |$0.4\sqrt{200/N}$| eV. In the middle column of Fig. 6 we show the contour plots for a catalogue consisting in 200 high-resolution QSO spectra.

In order to investigate the tightest constraints on |$\Sigma _i m_{\nu _i}$| that we can achieve with this catalogue, we impose realistic priors on the parameters. We assume that 〈F〉, α and the S/N are known within a 10 per cent error and that σ₈ is known within 3 per cent. The previous error intervals are at 1σ. We repeat the Fisher matrix analysis and we find a marginalized 1σ error on |$\Sigma _i m_{\nu _i}$| equal to 0.27 eV for a catalogue consisting in 200 high-resolution QSO spectra. The right-hand column of Fig. 6 shows the contour plots between |$\Sigma _i m_{\nu _i}$| and the other parameters.

Finally, we repeat the whole analysis for a different fiducial model with |$p_1^0=\Sigma _i m_{\nu _i}=0.0$| eV (the values of the other parameters are the same as those of the above fiducial model) and we find that the marginalized error in Ω_ν and the correlation with the other parameters are in very well agreement with those obtained for the 0.3 eV above fiducial model. This reinforces our assessment that the contribution of the first term in equation (10) to the whole Fisher matrix is negligible.

6 DISCUSSION AND CONCLUSIONS

We have studied in detail the sensitivity of the threshold crossing statistics (see e.g. Miralda-Escudé et al. 1996; Fan et al. 2002) to the masses of the neutrinos. We have focused our study on the low-density regions of the Lyman α forest in quasar spectra. Those regions correspond to the innermost parts of non-linear voids. We find that the number of regions above a given threshold in the flux is strongly affected by the masses of the neutrinos. The changes between different models are due to two factors: the change in amplitude and slope in the linear power spectrum driven by neutrinos and non-linear effects associated with CDM and neutrinos (note that neutrinos modify the non-linear evolution of the CDM distribution). The inner panel of Fig. 1 shows the linear (solid lines) and non-linear (dotted lines) versions of the power spectrum at z = 3 normalized to the case without neutrinos. Whereas the modification on large scales (k < 0.03 h Mpc⁻¹) is due to the change in the linear power spectrum, we find that on smaller spatial scales the non-linear effects dominate.

By creating realistic mock catalogues of high-resolution QSO spectra, we have found that systematic errors, such as the noise in the spectrum of the continuum fitting procedure, reduce the differences between different models (in terms of the average number of regions above the threshold as a function of the threshold) but do not erase the neutrino signal in the high-transitivity regions of the Lyman α forest. We have used the Fisher matrix formalism to forecast the errors associated with a measurement of Ω_ν using the method proposed in this paper. We conclude that |$\Sigma _i m_{\nu _i}$| can be measured with an error of 0.4 eV (1σ) using a catalogue of 200 high-resolution QSO spectra. By assuming that the values of the parameters 〈F〉, α and the S/N are known with a 10 per cent error and that the value of σ₈ is known with a 3 per cent error, we find that the |$\Sigma _i m_{\nu _i}$| can be measured with a precision of 0.27 eV (1σ) by using 200 high-resolution QSO spectra. This is not an unreasonable number of high-resolution spectra and is already available to the community. This accuracy is achievable given the current efforts to measure the IGM thermal state by using the flux PDF and power spectrum, wavelets and the line-width distribution from Voigt profile fitting.

We note that the statistics we have studied here to measure the neutrino masses implicitly contain more information than the one that can be extracted from measurements of the power spectrum. The threshold crossing statistics is analogous to the genus curve used to characterize the topology of the three-dimensional galaxy distribution (Fan et al. 2002). The amplitude of the genus per unit of volume depends only on the second moment of the power spectrum (Hamilton, Gott & Weinberg 1986) if the field is Gaussian. For non-Gaussian fields (as the ones we are studying here) the amplitude of the genus curve contains information of higher order correlation functions (see e.g. Zhang, Springel & Yang 2010).

The aim of this method is to put a new and independent constrain on Ω_ν using the Lyman α forest. The results found with this method can be combined with other cosmological measurements such as the CMB or LSS to improve current bounds on the masses of the neutrinos. Previous studies (Viel et al. 2009) have demonstrated the advantages of using the flux PDF or void statistics rather than the power spectrum or bispectrum to distinguish cosmological models with small differences (in this case authors investigated non-Gaussianities).

We have also studied another set of statistics which carry more information than the one compressed in the matter power spectrum. One of these statistics,⁴C₂(R, F_th), is widely used in material science (Torquato, Beasley & Chiew 1988; Jiao, Stillinger & Torquato 2009), and has been recently applied to the Lyman α forest (Lee & Spergel 2011). For the Lyman α forest, C₂(R, F_th), known as the cluster function, is defined as the probability function of finding a pair of pixels in the same phase, belonging to the same region, separated by a distance R. The pixels in the transmitted flux spectrum are assigned to two different phases: phase 1 if the value of the transmitted flux is larger than F_th and phase 2 otherwise. We have studied a statistics directly related to C₂(R, F_th): the PDF of the sizes of the connected regions above a given threshold in the transmitted flux spectrum. We have also focused on the high-transitivity regions of the Lyman α forest, restricting our study to the phase 1, to z = 2.2 and to values of F_th larger than 0.9. We have found that neutrino masses leave an imprint in this statistics, although it is smaller than the one we have presented in this paper. By using the Fisher matrix formalism, we conclude that this method could distinguish neutrino masses with a precision about 0.4 eV using 200 QSO spectra unaffected by systematic errors. This value should be compared with the error equal to 0.3 eV that we obtained with our method for the same number of QSO spectra. For this reason, we believe that the statistics we have presented in this paper is one of the most suitable observables, containing higher order information, to place new independent bounds on the masses of the neutrinos, even if the 1σ error bar is less competitive than that obtained by other cosmological probes.

Given that the error in |$\Sigma _i m_{\nu _i}$| can be significantly reduced by adding priors to the parameters, measuring independently the IGM thermal history and/or using QSO spectra with S/N larger than 100, we conclude that in the near future, a large (but not unreasonable) number of high-resolution QSO spectra could provide a relatively tight, new and independent constraint on neutrino masses which will be complementary to that provided by other large-scale structure probes.

ACKNOWLEDGMENTS

We thank the referee for the very helpful report. We thank Enzo Branchini, Olga Mena, Jordi Miralda-Escudé, Carlos Peña-Garay, Khee-Gan Lee, Jonathan Pritchard and Jesús Zavala for discussions and helpful comments on the manuscript and Volker Springel for providing us with gadget-3. The simulations were run on the Darwin Supercomputer Centre HPCF in Cambridge (UK) and at the research computing centre of Harvard University. FV-N was supported by the CSIC under the JAE-pre-doc program when this work started. This work was supported in part by grants: ASI/AAE, INFN/PD51, ASI/EUCLID, PRIN-INAF 2009, PRIN-MIUR (for MVi), as well as NSF grant AST-0907890 and NASA grants NNX08AL43G and NNA09DB30A (for AL). MVi and FV-N are supported by the ERC Starting Grant ‘CosmoIGM’.

1

http://camb.info/

2

Non-integer numbers are due to the path-length normalization.

3

We have explicitly checked that the presence of this term do not change our results.

4

They are called threshold probability functions.