Detecting Neutrino Mass by Combining Matter Clustering, Halos, and Voids

High-precision measurements of large-scale structure from upcoming cosmological surveys, such as the Dark Energy Spectroscopic Instrument (DESI), ¹⁷ Euclid, ¹⁸ Prime Focus Spectrograph, ¹⁹ Roman Space Telescope, ²⁰ Vera Rubin Observatory, ²¹ Square Kilometre Array, ²² and SPHEREx, ²³ are expected to revolutionize our understanding of fundamental physics, for example, by measuring neutrino mass. To fully realize the potential of these surveys, an urgent task is to determine the key observables that can maximize the scientific return. For Gaussian density fields, the answer is well known: the power spectrum, or equivalently, the correlation function, is the statistic that completely characterizes the field. Therefore, on large scales and at high redshift, where the density fluctuation in the universe resembles a Gaussian field, the power spectrum encapsulates all the information.

However, at low redshift and on small scales, nonlinear gravitational evolution moves information from the power spectrum into higher-order moments. It is currently ill-understood which observable(s) will allow retrieval of the maximum information in the nonlinear regime. For instance, it has been shown that for non-Gaussian fields, all clustering information may not be embedded in the infinite N-point statistics (Carron 2011, 2012). Since the number of modes increases rapidly by going to small scales, it is expected that the amount of information will also increase by considering observables in the mildly to fully nonlinear regime. While the amount of information, at least for some parameters, may saturate in the power spectrum (Rimes & Hamilton 2005; Villaescusa-Navarro et al. 2020) (see however Blot et al. 2016), many authors have shown that other statistics contain complementary information (see, e.g., Takada & Jain 2004; Sefusatti et al. 2006; Bergé et al. 2010; Kayo et al. 2013; Schaan et al. 2014; Liu et al. 2015a, 2015b; Kacprzak et al. 2016; Martinet et al. 2018; Shan et al. 2018; Allys et al. 2020; Dai et al. 2020; Hahn et al. 2020; Uhlemann et al. 2020; Banerjee & Abel 2021; Gualdi et al. 2021; Hahn & Villaescusa-Navarro 2021; Harnois-Déraps et al. 2021; Massara et al. 2021).

In this paper we quantify the information embedded in the nonlinear matter power spectrum, the halo mass function (HMF), and the void size function (VSF). We apply the Fisher formalism using a subset of the Quijote simulations (Villaescusa-Navarro et al. 2020), comprising of 23,000 N-body simulations for 16 different cosmologies spanning six cosmological parameters: Ω_m , Ω_b , h, n_s , σ₈, and M_ν . We study the information that these probes contain individually and when combined together, showing how the combination of these three statistics breaks degeneracies among the cosmological parameters, in turn setting very tight constraints. We consider the effects of both the autocorrelation for each probe and the cross correlation between different probes when computing the total information content. A simpler, theoretical, treatment combining cluster and void abundances has been studied by Sahlén (2019).

Of particular interest in this work are constraints on the sum of the neutrino masses M_ν ≡ ∑_ν m_ν . The first evidence for neutrino mass came from oscillation experiments (Fukuda et al. 1998; Ahmad et al. 2002; Araki et al. 2005; Ahn et al. 2006; An et al. 2012), which measured the difference in the squares of the masses of the three neutrino mass eigenstates. The best-fit results obtained from a joint analysis of oscillation experiments are ${\rm{\Delta }}{m}_{21}^{2}\equiv {m}_{2}^{2}-{m}_{1}^{2}\simeq 7.55\times {10}^{-5}\,{\mathrm{eV}}^{2}$ from solar neutrinos, and $| {\rm{\Delta }}{m}_{31}^{2}| \equiv | {m}_{3}^{2}-{m}_{1}^{2}| \simeq 2.50\times {10}^{-3}{\mathrm{eV}}^{2}$ from atmospheric neutrinos (de Salas et al. 2018). Since atmospheric neutrino experiments only probe the magnitude of the mass difference, there are two possibilities for the neutrino mass hierarchy: ${\rm{\Delta }}{m}_{31}^{2}\gt 0$, known as the normal hierarchy, or ${\rm{\Delta }}{m}_{31}^{2}\lt 0$, known as the inverted hierarchy. This gives a lower bound on the sum of the neutrino masses of M_ν ≳ 0.06eV for the normal hierarchy, or M_ν ≳ 0.1 eV for the inverted hierarchy. The current tightest upper bound on the effective electron neutrino mass from particle experiments is obtained by the KATRIN β-decay experiment, ${m}_{{\nu }_{e}}^{\mathrm{eff}}\lesssim 1.1\,\mathrm{eV}$ (Aker et al. 2019). ²⁴

Neutrinos also play an important role in the universe's history, as the presence of massive neutrinos both shifts the time of matter-radiation equality and suppresses the growth of structure on small scales. Measuring these effects enables determination of neutrino mass via cosmology, providing a complementary probe to particle physics (Doroshkevich et al. 1981; Hu et al. 1998; Eisenstein & Hu 1999; Lesgourgues et al. 2013). While the effects of neutrinos on linear (i.e., relatively large) scales are well understood theoretically, understanding the effects on nonlinear (i.e., relatively small) scales is an active field of research. There are numerous approaches to obtain theoretical predictions of the nonlinear effects of neutrinos, with varying computational efficiency (see, e.g., Saito et al. 2008; Brandbyge & Hannestad 2009, 2010; Shoji & Komatsu 2010; Viel et al. 2010; Bird et al. 2012; Ali-Haïmoud & Bird 2013; Costanzi et al. 2013; Castorina et al. 2014; Villaescusa-Navarro et al. 2014; Castorina et al. 2015; Archidiacono & Hannestad 2016; Banerjee & Dalal 2016; Carbone et al. 2016; Upadhye et al. 2016; Adamek et al. 2017; Emberson et al. 2017; Inman & Pen 2017; Senatore & Zaldarriaga 2017; Yu et al. 2017; Banerjee et al. 2018; Bird et al. 2018; Liu et al. 2018; Villaescusa-Navarro et al. 2018; Dakin et al. 2019; Chen et al. 2020, 2021; Bayer et al. 2021).

The current best constraints on M_ν arise by considering the cosmic microwave background (CMB) and combining it with other cosmological probes. Assuming a ΛCDM cosmological model, the upper bound on the neutrino mass from the Planck 2018 CMB temperature and polarization data is M_ν < 0.26 eV (95% CL) (Planck Collaboration et al. 2020). When combined with baryonic acoustic oscillations (BAOs) a more stringent bound of M_ν < 0.13 eV (95% CL) is obtained. Further combining with CMB lensing gives M_ν < 0.12 eV (95% CL).

A major limiting factor of current cosmological constraints is that CMB experiments measure the combined quantity A_s e^−2τ , where A_s is the amplitude of scalar perturbations and τ is the optical depth of reionization. Hence, accurate determination of τ is imperative to obtaining tight constraints when combining CMB with clustering/lensing (Allison et al. 2015; Liu et al. 2016; Archidiacono et al. 2017; Yu et al. 2018; Brinckmann et al. 2019). Most upcoming ground-based CMB experiments, such as Simons Observatory and CMB-S4, will not observe scales larger than ℓ ∼ 30, and will therefore be unable to directly constrain τ (Abazajian et al. 2016). Planck currently provides the best constraint of τ = 0.054 ± 0.007, with large improvements expected from the ongoing Cosmology Large Angular Scale Surveyor experiment (Watts et al. 2018) and the upcoming LiteBIRD (Hazumi et al. 2012) space mission. Furthermore, future radio 21 cm and, e.g., near-infrared/optical galaxy observations will provide new information on the optical depth which would also help improve the constraints form the CMB Liu et al. (2016); Brinckmann et al. (2019).

Before significant progress will be made in measuring τ, improved measurements of M_ν are expected from galaxy surveys such as DESI, the Legacy Survey of Space and Time, and Euclid. These surveys will measure fluctuations on nonlinear scales with unprecedented precision. There is thus much motivation to explore other probes of neutrino mass, beyond the traditional two-point clustering. By adding probes such as the halo and void abundances, we demonstrate that it is possible to break the strong degeneracy between M_ν and σ₈ usually seen in two-point clustering constraints (see, e.g., Villaescusa-Navarro et al. 2018). In turn, this gives tight constraints on neutrino mass, and in fact all cosmological parameters, potentially without the need for including CMB priors. In addition to improved constraints, having multiple independent probes of neutrino masses will allow for more robust controls of systematics.

The paper is organized as follows. We first review the Quijote simulations in Section 2. The Fisher formalism used to quantify the information content on the different observables is described in Section 3. We explain how the matter power spectrum, HMF, and VSF are obtained in Section 4. We show the results of our analysis in Section 5. Finally, we conclude in Section 6.

We quantify the information content of different cosmological observables using the Fisher matrix formalism. We model the observables using the Quijote simulations (Villaescusa-Navarro et al. 2020), a set of 23,000 N-body simulations that at a given redshift contain about 8 trillion (8 × 10¹²) particles over a total combined volume of 44,100 ${\left({h}^{-1}{\rm{G}}{\rm{p}}{\rm{c}}\right)}^{3}$. Each simulation considers a box of size 1 ${\left({h}^{-1}\,\mathrm{Gpc}\right)}^{3}$. The simulation subset used in this work spans a total of 16 different cosmological models that have been designed to evaluate the two ingredients required to compute the Fisher matrix: (1) the covariance matrix of the observables and (2) the derivatives of the observables with respect to the cosmological parameters. Despite their larger computational cost than analytic approaches (e.g., perturbation theory or the halo model), numerical simulations are more accurate into the fully nonlinear regime and rely on fewer assumptions and approximations.

We consider six cosmological parameters: Ω_m , Ω_b , h, n_s , σ₈, and M_ν . The set of cosmological parameters is shown in Table 1. To evaluate the covariance matrix, we use the 15,000 simulations of the fiducial cosmology. We compute the derivatives by considering simulations where only one cosmological parameter is varied, with all others fixed. We use 1000 simulations (500 pairs) for each derivative, with the exception of neutrino mass, where we use 1500 (see below).

Table 1. Characteristics of the Subset of the Quijote Simulations Used in this Work

Quijote Simulations
Name	Ωm	Ωb	h	ns	σ8	Mν (eV)	ICs	Realizations
Fiducial	0.3175	0.049	0.6711	0.9624	0.834	0.0	2LPT	15,000
Fiducial ZA	0.3175	0.049	0.6711	0.9624	0.834	0.0	Zel'dovich	500
${{\rm{\Omega }}}_{m}^{+}$	0.3275	0.049	0.6711	0.9624	0.834	0.0	2LPT	500
${{\rm{\Omega }}}_{m}^{-}$	0.3075	0.049	0.6711	0.9624	0.834	0.0	2LPT	500
${{\rm{\Omega }}}_{b}^{++}$	0.3175	0.051	0.6711	0.9624	0.834	0.0	2LPT	500
${{\rm{\Omega }}}_{b}^{--}$	0.3175	0.047	0.6711	0.9624	0.834	0.0	2LPT	500
h+	0.3175	0.049	0.6911	0.9624	0.834	0.0	2LPT	500
h−	0.3175	0.049	0.6511	0.9624	0.834	0.0	2LPT	500
${n}_{s}^{+}$	0.3175	0.049	0.6711	0.9824	0.834	0.0	2LPT	500
${n}_{s}^{-}$	0.3175	0.049	0.6711	0.9424	0.834	0.0	2LPT	500
${\sigma }_{8}^{+}$	0.3175	0.049	0.6711	0.9624	0.849	0.0	2LPT	500
${\sigma }_{8}^{-}$	0.3175	0.049	0.6711	0.9624	0.819	0.0	2LPT	500
${M}_{\nu }^{+}$	0.3175	0.049	0.6711	0.9624	0.834	0.1	Zel'dovich	500
${M}_{\nu }^{++}$	0.3175	0.049	0.6711	0.9624	0.834	0.2	Zel'dovich	500
${M}_{\nu }^{+++}$	0.3175	0.049	0.6711	0.9624	0.834	0.4	Zel'dovich	500

Note. The fiducial cosmology contains 15,000 simulations, that are used to compute the covariance matrix. In the other cosmological models, one parameter is varied at a time, and these simulations are used to compute the numerical derivatives. The ICs of all simulations were generated at z = 127 using 2LPT, except for the simulations with massive neutrinos and a copy of the fiducial cosmology, where the Zel'dovich approximation is used (see main text for further details). All realizations follow the evolution of 512³ cold dark matter (CDM; + 512³ neutrino) particles in a box of size 1 h⁻¹Gpc down to z = 0, with a gravitational softening length 50 h⁻¹kpc. For massive neutrino simulations, we assume three degenerate neutrino masses.

Download table as:  ASCIITypeset image

The initial conditions (ICs) were generated in all cases at z = 127 using second-order Lagrangian perturbation theory (2LPT) for simulations with massless neutrinos, by rescaling the z = 0 matter power spectrum using the scale-independent growth factor from linear theory. Because the 2LPT formalism has not yet been developed to account for massive neutrinos, the ICs for massive neutrino cosmologies adopt the Zel'dovich approximation with scale-dependent growth factors and rates, following Zennaro et al. (2017). For this reason there is also a "Fiducial (ZA)" class of simulations, which is identical to the fiducial simulations but with Zel'dovich ICs to match the M_ν simulations (see Villaescusa-Navarro et al. 2020, for further details); this enables accurate computation of derivatives with respect to M_ν . Note that in the full Quijote simulations there are two sets of Ω_b cosmologies; we use the ${{\rm{\Omega }}}_{b}^{++}$ and ${{\rm{\Omega }}}_{b}^{--}$ set too obtain smoother derivatives.

All simulations follow the evolution of 512³ dark matter particles down to z = 0. The simulations with massive neutrinos also contain 512³ neutrino particles. The gravitational force tree for neutrinos is turned on at z = 9. The gravitational softening for both dark matter and neutrinos is 50 h⁻¹kpc (1/40 of the mean interparticle distance). In this work, we consider redshift z = 0 only.

We use the Fisher matrix formalism (Tegmark et al. 1997; Heavens et al. 2007; Heavens 2009; Verde 2010) to calculate the information embedded in the nonlinear matter power spectrum, the HMF and the VSF, individually and when combined. The Fisher matrix is defined as

$$\begin{eqnarray}&&{F}_{{ij}}=-\left\langle \displaystyle \frac{{\partial }^{2}\mathrm{log}{ \mathcal L }}{\partial {\theta }_{i}\partial {\theta }_{j}}\right\rangle \,,\end{eqnarray} \tag{ 1 }$$

where ${ \mathcal L }$ is the likelihood and θ is the vector representing the parameters of the model (Fisher 1925). Under the assumption that the region around the maximum of the likelihood can be approximated as a multivariate normal distribution, one can write the Fisher matrix as

$$\begin{eqnarray}\begin{array}{rcl}{F}_{{ij}} & = & \displaystyle \frac{1}{2}\left[\displaystyle \frac{\partial {\boldsymbol{O}}}{\partial {\theta }_{i}}{C}^{-1}\displaystyle \frac{\partial {{\boldsymbol{O}}}^{T}}{\partial {\theta }_{j}}+\displaystyle \frac{\partial {\boldsymbol{O}}}{\partial {\theta }_{j}}{C}^{-1}\displaystyle \frac{\partial {{\boldsymbol{O}}}^{T}}{\partial {\theta }_{i}}\right]\\ & & +\displaystyle \frac{1}{2}\mathrm{Tr}\left[{C}^{-1}\displaystyle \frac{\partial C}{\partial {\theta }_{i}}{C}^{-1}\displaystyle \frac{\partial C}{\partial {\theta }_{j}}\right]\,,\end{array}\end{eqnarray} \tag{ 2 }$$

where O is the vector with the values of the observables and C is the covariance matrix. In order to avoid underestimating the errors, we follow Carron (2013) and neglect the dependence of the covariance on the cosmological parameters, by setting the last term of Equation (2) to zero. This is necessary when assuming a Gaussian likelihood. Note that we use Greek (Latin) characters to index observables (parameters).

In this work, the observables and parameters are given by

$$\begin{eqnarray*}&&\begin{array}{l}{\boldsymbol{O}}=\{{P}_{m}({k}_{1}),\ldots ,{P}_{m}({k}_{A}),{ \mathcal H }({M}_{1}),\ldots ,{ \mathcal H }({M}_{B}),\\ \quad { \mathcal V }({R}_{1}),\ldots ,{ \mathcal V }({R}_{D})\},\\ {\boldsymbol{\theta }}=\{{{\rm{\Omega }}}_{m},{{\rm{\Omega }}}_{b},h,{n}_{s},{\sigma }_{8},{M}_{\nu }\},\end{array}\end{eqnarray*}$$

respectively, where P_m (k) is the matter power spectrum at wavenumber k, ${ \mathcal H }(M)$ is the HMF at mass M, and ${ \mathcal V }(R)$ is the VSF at radius R. Note there are a total of A, B, and D bins for the matter power spectrum, the HMF, and the VSF respectively, giving a total dimensionality of A + B + D.

We quantify the information content by considering the marginalized error on the cosmological parameters,

$$\begin{eqnarray}&&\sigma ({\theta }_{i})\equiv \sqrt{{\left({F}^{-1}\right)}_{{ii}}}\,,\end{eqnarray} \tag{ 3 }$$

which is a lower bound.

3.1. Covariance Matrix

We estimate the covariance matrix using the N_cov = 15,000 simulations of the fiducial cosmology as

$$\begin{eqnarray}&&{C}_{\alpha \beta }=\left\langle \left({O}_{\alpha }-\langle {O}_{\alpha }\rangle \right)\left({O}_{\beta }-\langle {O}_{\beta }\rangle \right)\right\rangle ,\end{eqnarray} \tag{ 4 }$$

where 〈〉 denotes the mean over simulations. This is the largest number of simulations used for covariance estimation to date. We have verified that our combined results are converged even with half of the simulations. We show the results of our convergence tests in Appendix A.

3.2. Derivatives

For the cosmological parameters Ω_m , Ω_b , h, n_s , and σ₈, we approximate the derivatives using a central difference scheme centered on the fiducial cosmology,

$$\begin{eqnarray}&&\displaystyle \frac{\partial {\boldsymbol{O}}}{\partial {\theta }_{i}}\simeq \displaystyle \frac{{\boldsymbol{O}}({\theta }_{i}+\delta {\theta }_{i})-{\boldsymbol{O}}({\theta }_{i}-\delta {\theta }_{i})}{2{\theta }_{i}}.\end{eqnarray} \tag{ 5 }$$

Note that only the value of the ith cosmological parameter is perturbed about its fiducial value, θ_i , while the values of all other parameters are held fixed. The error of this approximation is ${ \mathcal O }(\delta {\theta }_{i}^{2})$.

For neutrinos we cannot use Equation (5) because the fiducial model has massless neutrinos, so O (θ_i − δ θ_i ) would correspond to a cosmology with negative neutrino mass. We thus compute the derivatives for neutrinos using a second-order forward difference scheme,

$$\begin{eqnarray}&&\begin{array}{l}\displaystyle \frac{\partial {\boldsymbol{O}}}{\partial {M}_{\nu }}\\ \simeq \,\displaystyle \frac{-3{\boldsymbol{O}}({M}_{\nu }+2\delta {M}_{\nu })+4{\boldsymbol{O}}({M}_{\nu }+\delta {M}_{\nu })-3{\boldsymbol{O}}({M}_{\nu })}{2\delta {M}_{\nu }},\end{array}\end{eqnarray} \tag{ 6 }$$

which has error ${ \mathcal O }(\delta {M}_{\nu }^{2})$. We exclusively use the ${M}_{\nu }^{++}$ and ${M}_{\nu }^{+++}$ cosmologies in Equation (6) throughout this work.

We use a total of N_der = 1000 (500+500) simulations to compute derivatives when using Equation (5), and 1500 when using Equation (6). In Appendix A we show that our results are robust and converged with this number of simulations. We also give evidence of robustness with respect to the choice of finite difference scheme for M_ν .

In this section we outline the cosmological observables considered in this work: the matter power spectrum, the HMF, and the VSF.

4.1. Matter Power Spectrum

The first observable we study is the matter power spectrum. For each realization, the density field is computed by depositing particle masses to a regular grid using the cloud-in-cell mass assignment scheme. In simulations with massive neutrinos we consider both CDM and neutrino particles when constructing the density field. The density contrast field, $\delta ({\boldsymbol{x}})=\rho ({\boldsymbol{x}})/\bar{\rho }-1$, is then Fourier transformed and the power spectrum is computing by averaging ∣δ( k )∣² over spherical bins in ∣k∣. The size of each bin is equal to the fundamental frequency, 2π/L, where L = 1 h⁻¹Gpc is the simulation box size.

A grid with 1024³ cells is used, which is large enough to avoid aliasing effects on the scales of interest for this work. In our analysis we consider wavenumbers up to ${k}_{\max }=0.5\,h\,{\mathrm{Mpc}}^{-1}$, using 79 bins. This choice of ${k}_{\max }$ is based on the fact that the clustering of the simulations is converged at this scale for this mass resolution (see Villaescusa-Navarro et al. 2020). We will however show that using a larger ${k}_{\max }$ would likely lead to even tighter constraints than the ones we report. We show the power spectrum for the fiducial cosmology in Figure 1.

Zoom In Zoom Out Reset image size

4.2. HMF

The second observable we consider is the HMF. Dark matter halos are identified using the friends-of-friends (FoF) algorithm (Davis et al. 1985), with a linking length b = 0.2. The halo finder considers only the dark matter distribution, as the contribution of neutrinos to the total mass of a halo is expected to be negligible (Villaescusa-Navarro et al. 2011; Ichiki & Takada 2012; Villaescusa-Navarro et al. 2013; LoVerde & Zaldarriaga 2014).

The HMF is defined as the comoving number density of halos per unit of (log) halo mass, ${dn}/d\mathrm{ln}M$. The mass of a halo is estimated as

$$\begin{eqnarray}&&M={{Nm}}_{p},\end{eqnarray} \tag{ 7 }$$

where N is the number of dark matter particles in the halo and m_p is the mass of a single dark matter particle. Note that in the Quijote simulations, there are only dark matter and neutrino particles, i.e., dark matter particles represent the CDM+baryon fluid. The mass of a dark matter particle is thus normalized according to Ω_cb , such that

$$\begin{eqnarray}&&{m}_{p}=\displaystyle \frac{V{\rho }_{c}}{{N}_{p}}{{\rm{\Omega }}}_{{cb}}=\displaystyle \frac{V{\rho }_{c}}{{N}_{p}}\left({{\rm{\Omega }}}_{m}-\displaystyle \frac{{M}_{\nu }}{93.14{h}^{2}}\right),\end{eqnarray} \tag{ 8 }$$

where V = L³ is the simulation volume, N_p is the total number of dark matter particles in the simulation, and ρ_c is the universe's critical energy density at z = 0. Thus m_p = m_p (Ω_m , M_ν ) is a cosmology dependent quantity, which induces noise when computing the derivatives of the HMF with respect to Ω_m or M_ν in a fixed mass bin. This is because it is the number of dark matter particles that is the fundamental constituent of the halo mass: a halo with a given number of particles will lie in the same number bin for all cosmologies, whereas it may lie in a different mass bin depending on the value of m_p . This noise can thus be avoided by instead working with bins of fixed particle number by considering the derivative of the comoving number density of halos per unit (log) number of particles, ${dn}/d\mathrm{ln}N$. One can then transform these derivatives in bins of fixed N to derivatives in bins of fixed M to obtain the derivatives of the HMF.

Using the shorthand ${ \mathcal H }$ to denote the HMF, we now derive this transformation. In practice, one measures the HMF for a fixed cosmology, thus working in logarithmic bins gives

$$\begin{eqnarray}&&{ \mathcal H }:= \displaystyle \frac{{dn}}{d\mathrm{ln}M}=\displaystyle \frac{{dn}}{d\mathrm{ln}N},\end{eqnarray} \tag{ 9 }$$

where it is understood that the derivative is taken with fixed cosmological parameters, θ . Explicitly, one can think of the HMF as a function of the cosmological parameters and halo mass, ${ \mathcal H }({\boldsymbol{\theta }},M)$, or the cosmological parameters and number of particles, ${ \mathcal H }({\boldsymbol{\theta }},N)$. Thus, the derivative of the HMF with respect to one of the cosmological parameters, θ, while holding all other cosmological parameters, ${/}\!\!\!{\theta }$, fixed can be written as

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{{/}\!\!\!{\theta }}={\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{M,{/}\!\!\!{\theta }}+{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \mathrm{ln}M}\right)}_{{\boldsymbol{\theta }}}{\left(\displaystyle \frac{\partial \mathrm{ln}M}{\partial \theta }\right)}_{{/}\!\!\!{\theta }},\end{eqnarray} \tag{ 10 }$$

or

$$\begin{eqnarray}&&{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{{/}\!\!\!{\theta }}={\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{N,{/}\!\!\!{\theta }}+{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \mathrm{ln}N}\right)}_{{\boldsymbol{\theta }}}{\left(\displaystyle \frac{\partial \mathrm{ln}N}{\partial \theta }\right)}_{{/}\!\!\!{\theta }}.\end{eqnarray} \tag{ 11 }$$

Equating these two equations and rearranging gives

$$\begin{eqnarray}&&\begin{array}{l}{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{M,{/}\!\!\!{\theta }}={\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{N,{/}\!\!\!{\theta }}+{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \mathrm{ln}N}\right)}_{{\boldsymbol{\theta }}}{\left[\displaystyle \frac{\partial \mathrm{ln}N}{\partial \theta }-\displaystyle \frac{\partial \mathrm{ln}M}{\partial \theta }\right]}_{{/}\!\!\!{\theta }}\\ \quad =\,{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \theta }\right)}_{N,{/}\!\!\!{\theta }}-{\left(\displaystyle \frac{\partial { \mathcal H }}{\partial \mathrm{ln}N}\right)}_{{\boldsymbol{\theta }}}{\left(\displaystyle \frac{\partial \mathrm{ln}{m}_{p}}{\partial \theta }\right)}_{{/}\!\!\!{\theta }},\end{array}\end{eqnarray} \tag{ 12 }$$

where Equation (7) was used in the final step.

The cosmology dependence of m_p takes effect in the final term of Equation (12). There is only a difference between the fixed N and fixed M derivative of the HMF when m_p depends on θ, i.e., when θ ∈ {Ω_m , M_ν }. Using Equation (8), one finds that

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{ln}{m}_{p}}{\partial {{\rm{\Omega }}}_{m}}=\displaystyle \frac{1}{{{\rm{\Omega }}}_{m}},\end{eqnarray} \tag{ 13 }$$

$$\begin{eqnarray}&&\displaystyle \frac{\partial \mathrm{ln}{m}_{p}}{\partial {M}_{\nu }}=\displaystyle \frac{1}{{{\rm{\Omega }}}_{{cb}}93.14{h}^{2}},\end{eqnarray} \tag{ 14 }$$

where it is understood that all cosmological parameters apart from the one in the derivative are held fixed at their fiducial values.

Thus, our procedure to compute derivatives of the HMF using Equation (12) is as follows. We first bin the number of halos according to the number of dark matter particles they contain. We then compute the derivatives for each fixed-N bin using the equations from Section 3.2, yielding the first term on the right-hand side of Equation (12). This will be sufficient for all cosmological parameters except for Ω_m and M_ν , as these require a correction term to transform to fixed-M bins due to the variation of m_p . The $\partial { \mathcal H }/\partial \mathrm{ln}N$ term can be computed via spline interpolation or by using finite difference methods between the bins of the HMF of the fiducial cosmology. We have confirmed the stability of both approaches. Finally, the derivative of $\mathrm{ln}{m}_{p}$ with respect to θ is computed using Equations (13) and 14 evaluated at the fiducial values.

We consider halos with a number of dark matter particles between 30 and 7000, using 15 logarithmically spaced bins. The corresponding halo mass range is approximately 2.0 × 10¹³ to 4.6 × 10¹⁵ h⁻¹ M_⊙. As with the matter power spectrum, this choice of binning and cuts is made to ensure convergence of the derivatives based on the resolution and number of the simulations available. Hence, using more bins and/or a larger mass range would likely lead to stronger constraints than we report. We show the HMF for the fiducial cosmology in Figure 2.

Zoom In Zoom Out Reset image size

4.3. VSF

We identify voids in the underlying matter field using a spherical void finding algorithm developed by Banerjee & Dalal (2016), which we now outline. We use a grid of resolution 768³ to look for voids—this is slightly finer than the CDM grid resolution of 512³ to enable detection of small voids. The density contrast field is then smoothed with a top-hat filter over a large-scale, R = 53.4 h⁻¹Mpc, which is a multiple of the grid spacing and is chosen to be bigger than the size of the largest void. Next, minima that are smaller than the threshold δ_th = − 0.7 in the smoothed field are considered as voids with radius R, unless they overlap with existing voids. This procedure is then performed iteratively while decrementing R by the grid spacing. In this work we use a threshold of δ_th = − 0.7, but have checked that results are similar for δ_th = − 0.5.

The VSF is then computed as the comoving number density of voids per unit of radius, denoted $d\tilde{n}/{dR}$. Unlike the HMF, the VSF is not prone to the changes in particle mass, since the void finder operates directly in the same unit as the VSF. The range of void sizes is limited by our resolution and the size of our simulated volume. Having found the voids, we apply radius cuts of ${R}_{\min }=10.4$ and ${R}_{\max }=29.9$ h⁻¹ Mpc, corresponding to 15 bins linear in R. As with the matter power spectrum and the HMF, this choice of binning and cuts is made to ensure convergence of the derivatives based on the resolution and number of the simulations available. Hence, using more bins and/or a larger range of void sizes may lead to stronger constraints than we report. We show the VSF for the fiducial cosmology in Figure 3.

Zoom In Zoom Out Reset image size

Investigation of the VSF, and void abundances, is a rich field that has shown promising theoretical work to match mocks (see, e.g., Platen et al. 2008; Bos et al. 2012; Sutter et al. 2012; Jennings et al. 2013; Pisani et al. 2015; Paillas et al. 2017; Contarini et al. 2019; Sahlén 2019; Verza et al. 2019).

In this section we present the main results of this work.

5.1. Full Covariance of the Probes

In Figure 4 we show the correlation matrix, defined as $\mathrm{Corr}({O}_{\alpha },{O}_{\beta }):= {C}_{\alpha \beta }/\sqrt{{C}_{\alpha \alpha }{C}_{\beta \beta }}$, where C_{α β} is the covariance matrix (Equation (4)). First we discuss the correlations for each individual probe (autocorrelations). For the matter power spectrum (bottom-left region of Figure 4), we observe some well-known structures: the covariance is almost diagonal on large scales, while mode-coupling induces significant off-diagonal correlations on small scales. For the HMF (central region of Figure 4), the covariance matrix is almost diagonal, with some small correlations between the different mass bins; the correlations are negative for heavy halos, but are positive for the lightest halos considered in this work. The covariance of the VSF (top-right region of Figure 4) is also almost diagonal, with the abundance of different void sizes slightly anticorrelated with nearby bins due to conservation of volume.

Zoom In Zoom Out Reset image size

Next, we consider the correlations between different probes (cross correlations). The HMF shows an interesting correlation pattern with the matter power spectrum: the abundance of the more (less) massive halos shows a ∼20% correlation (anticorrelation) with small scales of the matter power spectrum. Similar trends are seen between halos and large scales of the matter power spectrum, albeit at a weaker level. On the other hand, voids can be seen to be somewhat independent of both the matter power spectrum and halos, as their cross correlation is ≲5% for all scales and masses.

As discussed in Section 3, we combine the covariance matrix with the numerically computed derivatives to calculate the Fisher matrix. The numerical derivatives and related numerical convergence tests are shown in Appendix A.

5.2. Cosmological Constraints

We show the two-dimensional (2D) 68% and 95% confidence intervals obtained from our Fisher analysis for each individual probe, and the combination of all probes, in Figure 5. The constraints on the parameters are not generally tight when considering any of three probes alone, because we adopt a conservative survey volume of $1{\left({h}^{-1}\,\mathrm{Gpc}\right)}^{3}$, which is significantly smaller than what is achievable by DESI, $\sim {10}^{2}{\left({h}^{-1}\,\mathrm{Gpc}\right)}^{3}$.

Zoom In Zoom Out Reset image size

The three probes show different degeneracies and are sensitive to each parameter at different levels. For example, the HMF provides a relatively tight constraint on Ω_m when compared to the other two probes, as the HMF depends nonlinearly on and is highly sensitive to Ω_m (see, e.g., Haiman et al. 2001). The VSF provides weaker constraints than the other two probes on almost all parameters, except for n_s compared to P_m . Naively, this is not surprising, considering the relatively smaller range of scales being probed by the VSF compared to the matter power spectrum. More information could probably be retrieved by using other void-related observables, such as the void-matter correlation function.

Because the degeneracies between parameters are often very different for each probe, it is expected that combining the probes will break the degeneracies and in turn yield significantly tighter constraints on the cosmological parameters than the individual probes do. Indeed, the black ellipses in Figure 5 show the tight constraints obtained by combining the three probes. We emphasize that these constraints account for all the correlations between the different observables, i.e. , by using the full covariance matrix of Figure 4.

The benefit of combining the three probes is particularly well demonstrated in the M_ν –σ₈ plane. Because the combined constraints are too small to be visible in Figure 5, we zoom in on this plane in Figure 6. We find that, despite not being as powerful tools as P_m in constraining M_ν , the HMF and VSF both show degeneracies in different directions from that of P_m , which guarantees that constraints on the neutrino masses will be largely reduced by combining the three probes. In turn this helps break the well-known M_ν –σ₈ degeneracy for the matter power spectrum. We note that the area of these confidence contours, particularly for the HMF, can potentially be reduced by increasing the bin boundaries and/or by fine-tuning the binning schemes. Our choice of binning is restricted by our simulation resolution. We leave these investigations to future works.

Zoom In Zoom Out Reset image size

For a direct comparison to the usual constraints expected from the matter power spectrum, we show the 1D marginalized errors (Equation (3)) from different combinations of the probes with P_m in Figure 7. We study how the errors vary with the cutoff scale ${k}_{\max }$. Combining P_m with either the HMF, VSF, or both, can achieve a significant level of improvement on all 6 parameters. The combination with the HMF is typically more beneficial than the combination with the VSF. The only exception is for M_ν , where the VSF is the better probe to combine with P_m .

Zoom In Zoom Out Reset image size

While the constraints from P_m alone saturate at around ${k}_{\max }=0.2\,h\,{\mathrm{Mpc}}^{-1}$ for all parameters, the combined constraints for M_ν (and Ω_m ) continue to improve beyond ${k}_{\max }=0.5\,h\,{\mathrm{Mpc}}^{-1}$. This can be explained by the breaking of degeneracies when combining probes. It was shown in Figure 5 of Villaescusa-Navarro et al. (2020) that increasing ${k}_{\max }$ beyond 0.2 h Mpc⁻¹ leads to a squeezing along the semiminor axes (i.e., the most constraining direction) for the P_m ellipses. While this squeezing has little effect on the marginalized error on M_ν from P_m alone, its effects are manifest when combined with other probes with misaligned contours, resulting in significant tightening of constraints. Even though the numerical resolution of the Quijote simulations prevent us from confidently investigating beyond ${k}_{\max }=0.5\,h\,{\mathrm{Mpc}}^{-1}$, our results hint that even tighter constraints could be achieved by including smaller scales.

In Table 2 we list the errors for ${k}_{\max }=0.5\,h\,{\mathrm{Mpc}}^{-1}$ using different probe combinations. We list the constraints obtained by combining all three probes while (1) only using the diagonals of the covariance matrix (diag), (2) only considering auto-covariance (auto), and (3) considering the full covariance (full). We find that using only the diagonal components of the covariance matrix, effectively ignoring both the correlation between the probes and between different bins of the same probe, leads to a factor of 1.7 increase on the error on the neutrino mass. Using only block cross correlations, i.e., ignoring the correlation between the probes, leads to a factor of 1.2 increase on the error on the neutrino mass. Therefore, to obtain the tightest constraints, it is crucial to model the full covariance matrix. It is interesting to note that when considering the matter power spectrum alone, correlations cause an increase in errors due to the positive correlation between different scales (see Figure 4). However, it is the complex correlation structure, notably the anticorrelations, introduced by considering the HMF and VSF that leads to a reduction in error, both for the HMF and VSF individually, and in turn when combining all probes. The association of anticorrelation with the tightening of constraints was also pointed out by Chartier et al. (2021).

Table 2. Marginalized Errors of Cosmological Parameters for ${k}_{\max }=0.5\,h\,{\mathrm{Mpc}}^{-1}$ Using Different Probe Combinations

Marginalized Fisher Constraints
Probe(s)	Ωm	Ωb	h	ns	σ8	Mν (eV)
Pm	0.098	0.039	0.51	0.50	0.014	0.77
HMF	0.034	0.042	0.28	0.12	0.082	1.6
VSF	0.31	0.12	1.3	0.42	0.083	1.1
Pm + HMF	0.00077	0.0089	0.076	0.034	0.0016	0.061
Pm + VSF	0.016	0.011	0.12	0.074	0.0018	0.025
HMF + VSF	0.0063	0.037	0.23	0.10	0.0069	0.096
Pm + HMF + VSF (diag)	0.0015	0.0088	0.066	0.028	0.00061	0.031
Pm + HMF + VSF (auto)	0.0015	0.0086	0.071	0.033	0.0016	0.025
Pm + HMF + VSF (full)	0.00071	0.0084	0.064	0.025	0.0015	0.018
Multiplicative improvement	137	5	8	20	10	43

Note. We list the constraints obtained by combining all three probes while (1) only using the diagonals of the covariance matrix (diag), (2) only considering auto-covariance (auto), and (3) considering the full covariance (full). We highlight in bold the full constraints on the sum of the neutrino masses. We also list the multiplicative improvement in the constraints from the full covariance compared to those from P_m alone.

Download table as:  ASCIITypeset image

In Table 2 we quantify the improvement of the combined constraints compared to those achieved from P_m alone. We find the improvements to be a factor of 137, 5, 8, 20, 10, and 43, for Ω_m , Ω_b , h, n_s , σ₈, and M_ν , respectively. Thus, we achieve 43 times tighter constraints on neutrino mass by combining all three probes. Specifically, the marginalized errors on M_ν are 0.77 eV (P_m alone) and 0.018 eV (P_m +HMF+VSF). We provide an additional plot in Appendix B to show the confidence ellipses when combining only two of the probes at a time.

Upcoming galaxy surveys will map large volumes of the universe at low redshifts, with the potential to drastically improve our understanding of the underlying cosmological model. With the unprecedentedly precision achievable by these surveys, it is expected that a very large amount of cosmological (and astrophysical) information will lie in the mildly to fully nonlinear regime, where analytic methods are often intractable. It remains an open question which observable(s) will lead to the tightest bounds on the cosmological parameters.

In this paper, we use the Quijote simulations, based on the Fisher formalism, to quantify the information content embedded in the nonlinear matter power spectrum, the HMF, and the VSF, both individually and when combined, at z = 0. We find that the HMF and VSF have different degeneracies to each other and to the matter power spectrum, particularly in the M_ν –σ₈ plane (Figures 5 & 6). In terms of measuring neutrino mass, we find the VSF to be the more complementary probe to combine with the matter power spectrum. This is consistent with findings that void properties are particularly sensitive to matter components that are less clustered, such as neutrinos (Massara et al. 2015; Kreisch et al. 2019).

By combining the nonlinear matter power spectrum (${k}_{\max }=0.5\,h\,{\mathrm{Mpc}}^{-1}$), with the HMF (M ≳ 2 × 10¹³ h⁻¹ M_⊙), and the VSF (R ≥ 10.4 h⁻¹Mpc), we achieve significantly tighter constraints on the cosmological parameters compared to P_m alone (Figure 7). In particular, we find that with a volume of just $1\,{\left({h}^{-1}\,\mathrm{Gpc}\right)}^{3}$, the error on the sum of neutrino masses from the combined probes is at the 0.018 eV level, compared to 0.77 eV from the matter power spectrum alone—a factor of 43 improvement. We emphasize that this value mainly demonstrates the information content in the late-time statistics, and they are not forecasts for any particular survey.

Also of particular interest is the factor 137 improvement in the error on Ω_m . This is driven by the information in the HMF, and gives a marginalized error of σ(Ω_m ) = 7.1 × 10⁻⁴, which is almost 100 times smaller than the error obtained from a joint large-scale structure analysis by DES Y1 (σ(Ω_m ) ≈ 0.04, To et al. 2021), and 8 times smaller than Planck 2018, (σ(Ω_m ) ≈ 5.6 × 10⁻³ (TT,TE,EE+lowE+lensing+BAO), Planck Collaboration et al. 2020). In addition, we found σ(h) = 0.064 by combining the three probes, which is 8 times tighter than the constraints from the matter power spectrum alone. This could provide a new angle to investigate the Hubble tension.

There are several caveats in this work. First, we assumed perfect knowledge of the three-dimensional spatial distribution of the underlying matter field in real space. However, in reality, one observes either tracers of the matter field in redshift space, or the projected matter field through lensing. Therefore, additional links must be made to bridge the galaxy–matter connection and the 2D lensing–3D matter distribution gaps. This effect is also relevant for voids: in this work we considered voids in the 3D matter field, which is not something current surveys are able to observe directly. Detecting voids in the matter field from photometric (2D lensing) data has been considered in works such as Pollina et al. (2019) and Davies et al. (2021). Alternatively, one can measure voids in the 3D halo field (see, e.g., Nadathur 2016; Contarini et al. 2019). If we were to instead have considered voids in the 3D CDM field, the combined error on M_ν slightly degrades to 0.025 eV. However, considering voids in the CDM field versus halo field can lead to nontrivial differences in void properties, which might increase or decrease constraints (Kreisch et al. 2019). We will consider voids in the halo field in a future work.

A further note regarding voids is that there are various conventions when it comes to defining voids (see, e.g., Platen et al. 2007; Sutter et al. 2015). It would thus be interesting further work to consider how the choice of void finder impacts constraints. A different void finder may be able to extract additional information compared to the spherical void finder applied here.

Another limitation of this work is that our simulations consider only gravitational interactions and hence ignore baryonic effects, which can impact the small-scale matter distribution. This is particularly relevant for both clustering and halos (see, e.g., Villaescusa-Navarro et al. 2021; Cromer et al. 2021; Debackere et al. 2021, and references therein) while it is expected that baryons have a lower impact on voids (Paillas et al. 2017). Furthermore, halo clustering is influenced by various properties, such as spin, concentration, and velocity anisotropy, which have not been considered in this work (see, e.g., Wechsler et al. 2006; Gao & White 2007; Faltenbacher & White 2010; Lacerna & Padilla 2012; Lacerna et al. 2014; Paranjape et al. 2018; Shi & Sheth 2018).

Additionally, we have neglected supersample covariance (Takada & Hu 2013; Li et al. 2014), which could modify the errors reported in this work.

We also note that the constraints obtained here may be overly conservative due to the limited number and resolution of simulations available. First, this means that the number of bins used are likely suboptimal. Second, applying more aggressive bounds on the observables, e.g., a higher ${k}_{\max }$, a larger halo mass range, or a larger void size range, would likely also reduce the combined constraints. Third, we only considered a single redshift, z = 0: in practice, surveys measure z > 0 where the universe is more linear and the constraints will thus be weaker, however, combining multiple redshifts could tighten the constraints as found in works such as Liu & Madhavacheril (2019). Fourth, we considered a volume of only $1\,{\left({h}^{-1}\,\mathrm{Gpc}\right)}^{3}$, whereas surveys such as Euclid and DESI will cover volumes of around ${10}^{2}\,{\left({h}^{-1}\,\mathrm{Gpc}\right)}^{3}$, so, conservatively, the error on the parameters will shrink by a factor of $1/\sqrt{{10}^{2}}=0.1$. Fifth, we have only considered three probes; using the same observations, one can derive other statistics such as the bispectrum, void profile, and BAO, which could be combined with the statistics considered here to further break degeneracies. Finally, considering redshift space distortions would also tighten constraints as neutrinos are distinguishable from CDM via their higher thermal velocity.

We have demonstrated that combining multiple probes of cosmological structure using their full covariance matrix provides remarkably tight constraints on the cosmological parameters, and helps extract much additional information from small scales. In particular, we have shown that there is, in principle, sufficient information to measure the sum of the neutrino masses at the minimum mass of 0.06 eV. Our results are in good agreement with Sahlén 2019 who found that combining halo and void abundances can yield ${ \mathcal O }(0.01\,\mathrm{eV})$ constraints on the neutrino mass. This approach opens a promising pathway to measure neutrino mass, potentially without relying on CMB-based measurements which require accurate knowledge of the optical depth, τ. In addition, comparing constraints from different combinations of observables, e.g., CMB+P_m and P_m +HMF+VSF, will help identify systematic issues and provide robust evidence for any discovery. We thus hope our work will motivate galaxy survey collaborations to build the simulations and analytic tools necessary to implement this approach on upcoming observational data.

We thank Ravi Sheth for useful conversations on the early stages of this project, and Alice Pisani for fruitful discussion regarding voids. The Quijote simulations can be found at https://github.com/franciscovillaescusa/Quijote-simulations. The analysis of the simulations has made use of the Pylians libraries, publicly available at https://github.com/franciscovillaescusa/Pylians3. The work of D.N.S., B.D.W., and S.H. is supported by the Simons Foundation. L.V. acknowledges support from the European Union Horizon 2020 research and innovation program ERC (BePreSySe, grant agreement 725327) and MDM-2014-0369 of ICCUB (Unidad de Excelencia Maria de Maeztu). M.V. is supported by INFN PD51 Indark grant and by the agreement ASI-INAF n.2017-14-H.0.

In this section we verify the stability of our results to reduction in the number of simulations used to compute the covariance matrix and derivatives. In Figure 8 we show the derivatives of the matter power spectrum (top), HMF (middle), and VSF (bottom) with respect to the cosmological parameters when using a different number of realizations. For the matter power spectrum, the derivatives are already converged when the mean values for each model are computed with 300 realizations. Results are slightly noisier for the HMF and the VSF, but still sufficiently converged by 500 realizations.

Zoom In Zoom Out Reset image size

Next, we comment on the convergence of our simulated results with theory. The convergence of the matter power spectrum in Quijote has been thoroughly tested (Aviles & Banerjee 2020; Hahn et al. 2020; Villaescusa-Navarro et al. 2020). For the VSF there is no theoretical formula accurate enough to compute derivatives, but we have checked results are robust to the parameters used in the void finder. Therefore, we only compare our measured HMF to theoretical predictions. For the HMF, we plot the theoretical predictions of Sheth–Tormen (ST; Sheth & Tormen 1999, 2002) and Tinker (Tinker et al. 2008). We use the prescription of Costanzi et al. (2013) in the case of massive neutrino cosmologies by replacing Ω_m → Ω_cb and P_m → P_cb as neutrinos have negligible contribution to halo mass. There is good agreements between these predictions and Quijote. We have also checked that there is good agreement for different choice of step size (not shown). Note that these theoretical formulae provide a guideline rather than exact predictions, as they were fitted to simpler simulations or calibrated on spherical overdensity halos, as opposed to FoF here.

In Figure 9 we show the convergence of the Fisher matrix elements with respect to the number of realizations used to compute the covariance, N_cov, and derivatives N_der. We consider the Fisher matrix components for P_m (red), the HMF (blue), the VSF (green), and the combined probes (black). The gray bands corresponds to the ±5% interval. While there is some noise in the σ₈ component of the Fisher matrix for P_m as function of N_cov, good convergence is achieved by 15,000. Likewise the Fisher matrix is well converged as a function of N_der. Crucially, the Fisher matrix elements for the combined probes (black) all show good convergence. Note that when combining probes we scale the power spectrum by a factor of 10⁻¹⁰ to ensure that the condition number of the covariance matrix is sufficiently low for accurate inversion.

Zoom In Zoom Out Reset image size

Finally, we comment on the choice of finite difference scheme used to compute the derivative of probes with respect to M_ν . Throughout the paper we used Equation (6) with δ M_ν = 0.2 eV, thus making use of simulations with M_ν = 0, 0.2, and 0.4 eV. Using this scheme we found the full combined constraint on M_ν is 0.018 eV, as shown in Table 2. To illustrate robustness to this choice of finite difference scheme, we also performed the analysis using Equation (6) with δ M_ν = 0.1 eV and found it to give an identical constraint of 0.018 eV. Additionally, we tried a forward difference scheme between M_ν = 0 and 0.1 eV, which also gave identical constraints. Hence, the results are consistent with the choice of finite difference scheme. We do also note that since the joint constraints on the parameters given in Table 2 are smaller than the step sizes used to compute derivatives, it would be interesting to investigate the effect of smaller step sizes on the joint constraints. This would reduce the error in the numerical derivatives, and thus may slightly modify the joint constraints.

Given these results, we believe that our conclusions are robust against potential numerical systematics. We note again that our bin configuration has been chosen with these results in mind, to ensure sufficiently converged derivatives and Fisher matrix components, but in principle one could consider more bins over a wider range to potentially obtain tighter constraints.

Figure 10 shows the 2D Fisher contours to illustrate the effects of only combining two out of the three probes at a time. In most cases, the constraints obtained by combining the HMF with the VSF are the weakest, indicating that it is important to use the information from the nonlinear matter power spectrum to break degeneracies.

Zoom In Zoom Out Reset image size