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

Adrian E. Bayer, Francisco Villaescusa-Navarro, Elena Massara, Jia Liu, David N. Spergel, Licia Verde, Benjamin Wandelt, Matteo Viel, Shirley Ho

► To cite this version:

Adrian E. Bayer, Francisco Villaescusa-Navarro, Elena Massara, Jia Liu, David N. Spergel, et al. Detecting Neutrino Mass by Combining Matter Clustering, Halos, and Voids. The Astrophysical Journal, 2021, 919 (1), pp.24. 10.3847/1538-4357/ac0e91. hal-03157122

HAL Id: hal-03157122

https://hal.science/hal-03157122v1

Submitted on 6 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

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

Adrian E. Bayer$^{1,2,16}$, Francisco Villaescusa-Navarro$^{3,4}$, Elena Massara$^{4,5}$, Jia Liu$^{1,2,6}$, David N. Spergel$^{3,4}$, Licia Verde$^{7,8}$, Benjamin D. Wandelt$^{4,9,10}$, Matteo Viel$^{1,12,13,14}$, and Shirley Ho$^{3,4,15}$

$^{1}$ Berkeley Center for Cosmological Physics, University of California, 341 Campbell Hall, Berkeley, CA 94720, USA; abayer@berkeley.edu

$^{2}$ Department of Physics, University of California, 366 LeConte Hall, Berkeley, CA 94720, USA

$^{3}$ Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton, NJ 08544-0010, USA; fvillaescusa@flatironinstitute.org

$^{4}$ Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

$^{5}$ Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada

$^{6}$ Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA

$^{7}$ Institut de Ciències del Cosmos, University of Barcelona, ICCUB, Barcelona E-08028, Spain

$^{8}$ Institutució Catalana de Recerca i Estudis Avancats, Passeig Lluís Companys 23, Barcelona E-08010, Spain

$^{9}$ Institut d'Astrophysique de Paris, UMR 7095, CNRS, 98 bis Boulevard Arago, F-75014 Paris, France

$^{10}$ Institut Lagrange de Paris, Sorbonne Universités, 98 bis Boulevard Arago, F-75014 Paris, France

$^{11}$ SISSA, Via Bonomea 265, I-34136 Trieste, Italy

$^{12}$ INFN—Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, Via Bonomea 265, I-34136 Trieste, Italy

$^{13}$ INAF—Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34143 Trieste, Italy

$^{14}$ IFPU—Institute for Fundamental Physics of the Universe, Via Beirut 2, I-34014 Trieste, Italy

$^{15}$ Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA

Received 2021 March 22; revised 2021 June 21; accepted 2021 June 23; published 2021 September 20

## Abstract

We quantify the information content of the nonlinear matter power spectrum, the halo mass function, and the void size function, using the *Quijote N*-body simulations. We find that these three statistics exhibit very different degeneracies among the cosmological parameters, and thus the combination of all three probes enables the breaking of degeneracies, in turn yielding remarkably tight constraints. We perform a Fisher analysis using the full covariance matrix, including all auto- and cross correlations, finding that this increases the information content for neutrino mass compared to a correlation-free analysis. The multiplicative improvement of the constraints on the cosmological parameters obtained by combining all three probes compared to using the power spectrum alone are: 137, 5, 8, 20, 10, and 43, for $\Omega_m$, $\Omega_b$, $h$, $n_s$, $\sigma_8$, and $M_{\odot}$, respectively. The marginalized error on the sum of the neutrino masses is $\sigma(M_{\odot}) = 0.018$ eV for a cosmological volume of $l(h^{-1} \text{ Gpc})^3$, using $k_{\text{max}} = 0.5~h~\text{Mpc}^{-1}$, and without cosmic microwave background (CMB) priors. We note that this error is an underestimate insomuch as we do not consider super-sample covariance, baryonic effects, and realistic survey noises and systematics. On the other hand, it is an overestimate insomuch as our cuts and binning are suboptimal due to restrictions imposed by the simulation resolution. Given upcoming galaxy surveys will observe volumes spanning $\sim 100(h^{-1} \text{ Gpc})^3$, this presents a promising new avenue to measure neutrino mass without being restricted by the need for accurate knowledge of the optical depth, which is required for CMB-based measurements. Furthermore, the improved constraints on other cosmological parameters, notably $\Omega_m$, may also be competitive with CMB-based measurements.

*Unified Astronomy Thesaurus concepts: Cosmological neutrinos (338); Neutrino masses (1102); Large-scale structure of the universe (902); N-body simulations (1083); Voids (1779); Galaxy dark matter halos (1880); Cosmological parameters (339); Cosmological parameters from large-scale structure (340)*

## 1. Introduction

High-precision measurements of large-scale structure from upcoming cosmological surveys, such as the Dark Energy Spectroscopic Instrument (DESI),$^{17}$ Euclid,$^{18}$ Prime Focus Spectrograph,$^{19}$ Roman Space Telescope,$^{20}$ Vera Rubin Observatory,$^{21}$ Square Kilometre Array,$^{22}$ and SPHERE,$^{23}$ 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

<footnote>$^{16}$ Corresponding author.</footnote>

<footnote>$^{17}$ https://www.desi.lbl.gov</footnote>

<footnote>$^{18}$ https://www.euclid-ec.org</footnote>

<footnote>$^{19}$ https://pfs.ipmu.jp/index.html</footnote>

<footnote>$^{20}$ https://wfirst.gsfc.nasa.gov/index.html</footnote>

<footnote>$^{21}$ https://www.lsst.org</footnote>

<footnote>$^{22}$ https://www.skatelescope.org</footnote>

<footnote>$^{23}$ https://www.jpl.nasa.gov/missions/spherex</footnote>
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: $\Omega_m$, $\Omega_b$, $h$, $n_s$, $\sigma_8$, and $M_\nu$. 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_\nu \equiv \sum_i m_\nu$. 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 $\Delta m_{21}^2 \equiv m_2^2 - m_1^2 \simeq 7.55 \times 10^{-5} \text{ eV}^2$ from solar neutrinos, and $|\Delta m_{31}^2| \equiv |m_3^2 - m_1^2| \simeq 2.50 \times 10^{-3} \text{ 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: $\Delta m_{31}^2 > 0$, known as the normal hierarchy, or $\Delta m_{31}^2 < 0$, known as the inverted hierarchy. This gives a lower bound on the sum of the neutrino masses of $M_\nu \gtrsim 0.06 \text{ eV}$ for the normal hierarchy, or $M_\nu \gtrsim 0.1 \text{ eV}$ for the inverted hierarchy. The current tightest upper bound on the effective electron neutrino mass from particle experiments is obtained by the KATRIN $\beta$-decay experiment, $m_{\nu_e}^{\text{eff}} \lesssim 1.1 \text{ eV}$ (Aker et al. 2019).<fnref n="24" />

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_\nu$ arise by considering the cosmic microwave background (CMB) and combining it with other cosmological probes. Assuming a $\Lambda$CDM cosmological model, the upper bound on the neutrino mass from the Planck 2018 CMB temperature and polarization data is $M_\nu < 0.26 \text{ eV}$ (95% CL) (Planck Collaboration et al. 2020). When combined with baryonic acoustic oscillations (BAOs) a more stringent bound of $M_\nu < 0.13 \text{ eV}$ (95% CL) is obtained. Further combining with CMB lensing gives $M_\nu < 0.12 \text{ eV}$ (95% CL).

A major limiting factor of current cosmological constraints is that CMB experiments measure the combined quantity $A_s e^{-2\ell}$, where $A_s$ is the amplitude of scalar perturbations and $\ell$ is the optical depth of reionization. Hence, accurate determination of $\tau$ 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 $\ell \sim 30$, and will therefore be unable to directly constrain $\tau$ (Abazajian et al. 2016). Planck currently provides the best constraint of $\tau = 0.054 \pm 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 from the CMB Liu et al. (2016); Brinckmann et al. (2019).

Before significant progress will be made in measuring $\tau$, improved measurements of $M_\nu$ 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_\nu$ and $\sigma_8$ 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

<footnote n="24">²⁴ Single $\beta$-decay experiments do not directly measure the neutrino mass sum, but rather the effective mass of electron neutrinos. In the quasi-degenerate regime where the eigenmasses $m_i > 0.2 \text{ eV}$ ($i = 1, 2, 3$), the three eigenmasses are the same to better than 3%, and hence $m_{\nu_e}^{\text{eff}} \approx 1/3 M_\nu$.</footnote>
<table><caption>**Table 1**

Characteristics of the Subset of the *Quijote* Simulations Used in this Work</caption><thead><tr><th rowspan="2">Name</th><th colspan="7">Quijote Simulations</th></tr><tr><th>Ω<sub>m</sub></th><th>Ω<sub>b</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th><th>M<sub>ν</sub> (eV)</th><th>ICs</th><th>Realizations</th></tr></thead><tbody><tr><td>Fiducial</td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>15,000</td></tr><tr><td>Fiducial ZA</td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td>0.0</td><td>Zel'dovich</td><td>500</td></tr><tr><td>Ω<sup>+</sup><sub>m</sub></td><td><u>0.3275</u></td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>Ω<sub>m</sub></td><td><u>0.3075</u></td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>Ω<sup>++</sup><sub>b</sub></td><td>0.3175</td><td><u>0.051</u></td><td>0.6711</td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>Ω<sub>b</sub><sup>-</sup></td><td>0.3175</td><td><u>0.047</u></td><td>0.6711</td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>h<sup>+</sup></td><td>0.3175</td><td>0.049</td><td><u>0.6911</u></td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>h<sup>-</sup></td><td>0.3175</td><td>0.049</td><td><u>0.6511</u></td><td>0.9624</td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>n<sub>s</sub><sup>+</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td><u>0.9824</u></td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>n<sub>s</sub><sup>-</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td><u>0.9424</u></td><td>0.834</td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>σ<sub>8</sub><sup>+</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td><u>0.849</u></td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>σ<sub>8</sub><sup>-</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td><u>0.819</u></td><td>0.0</td><td>2LPT</td><td>500</td></tr><tr><td>M<sub>ν</sub><sup>+</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td><u>0.1</u></td><td>Zel'dovich</td><td>500</td></tr><tr><td>M<sub>ν</sub><sup>++</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td><u>0.2</u></td><td>Zel'dovich</td><td>500</td></tr><tr><td>M<sub>ν</sub><sup>+++</sup></td><td>0.3175</td><td>0.049</td><td>0.6711</td><td>0.9624</td><td>0.834</td><td><u>0.4</u></td><td>Zel'dovich</td><td>500</td></tr></tbody></table>

**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^3$ cold dark matter (CDM; + $512^3$ neutrino) particles in a box of size $1 h^{-1}$Gpc down to $z = 0$, with a gravitational softening length $50 h^{-1}$kpc. For massive neutrino simulations, we assume three degenerate neutrino masses.

the results of our analysis in Section 5. Finally, we conclude in Section 6.

## 2. Simulations

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 \times 10^{12}$) particles over a total combined volume of 44,100 ($h^{-1}$Gpc)$^3$. Each simulation considers a box of size $1 (h^{-1}$ Gpc)$^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: $\Omega_m$, $\Omega_b$, $h$, $n_s$, $\sigma_8$, and $M_\nu$. 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).

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_\nu$ simulations (see Villaescusa-Navarro et al. 2020, for further details); this enables accurate computation of derivatives with respect to $M_\nu$. Note that in the full *Quijote* simulations there are two sets of $\Omega_b$ cosmologies; we use the $\Omega_b^{++}$ and $\Omega_b^{--}$ set too obtain smoother derivatives.

All simulations follow the evolution of $512^3$ dark matter particles down to $z = 0$. The simulations with massive neutrinos also contain $512^3$ 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^{-1}$kpc (1/40 of the mean interparticle distance). In this work, we consider redshift $z = 0$ only.

## 3. Fisher Information

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

$$F_{ij} = - \left\langle \frac{\partial^2 \log \mathcal{L}}{\partial \theta_i \partial \theta_j} \right\rangle, \quad (1)$$

where $\mathcal{L}$ is the likelihood and $\theta$ 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

$$F_{ij} = \frac{1}{2} \left[ \frac{\partial \mathcal{O}}{\partial \theta_i} C^{-1} \frac{\partial \mathcal{O}^T}{\partial \theta_j} + \frac{\partial \mathcal{O}}{\partial \theta_j} C^{-1} \frac{\partial \mathcal{O}^T}{\partial \theta_i} \right] + \frac{1}{2} \mathrm{Tr} \left[ C^{-1} \frac{\partial C}{\partial \theta_i} C^{-1} \frac{\partial C}{\partial \theta_j} \right], \quad (2)$$

where $\mathcal{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{aligned} \mathcal{O} &= \{P_m(k_1), \dots, P_m(k_A), \mathcal{H}(M_1), \dots, \mathcal{H}(M_B), \\ & \qquad \mathcal{V}(R_1), \dots, \mathcal{V}(R_D)\}, \\ \theta &= \{\Omega_m, \Omega_b, h, n_s, \sigma_8, M_\nu\}, \end{aligned}$$

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,

$$\sigma(\theta_i) \equiv \sqrt{(F^{-1})_{ii}}, \quad (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

$$C_{\alpha\beta} = \langle (O_{\alpha} - \langle O_{\alpha} \rangle)(O_{\beta} - \langle O_{\beta} \rangle) \rangle, \quad (4)$$

where $\langle \rangle$ 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 $\Omega_m$, $\Omega_b$, $h$, $n_s$, and $\sigma_8$, we approximate the derivatives using a central difference scheme centered on the fiducial cosmology,

$$\frac{\partial \mathcal{O}}{\partial \theta_i} \approx \frac{\mathcal{O}(\theta_i + \delta\theta_i) - \mathcal{O}(\theta_i - \delta\theta_i)}{2\theta_i}. \quad (5)$$

Note that only the value of the ith cosmological parameter is perturbed about its fiducial value, $\theta_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 $\mathcal{O}(\theta_i - \delta\theta_i)$ would correspond to a cosmology with negative neutrino mass. We thus compute the derivatives for neutrinos using a second-order

forward difference scheme,

$$\frac{\partial \mathcal{O}}{\partial M_{\nu}} \approx \frac{-3\mathcal{O}(M_{\nu} + 2\delta M_{\nu}) + 4\mathcal{O}(M_{\nu} + \delta M_{\nu}) - 3\mathcal{O}(M_{\nu})}{2\delta M_{\nu}}, \quad (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_{\nu}$.

## 4. Cosmological Probes

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(\mathbf{x}) = \rho(\mathbf{x})/\bar{\rho} - 1$, is then Fourier transformed and the power spectrum is computed by averaging $|\delta(\mathbf{k})|^2$ over spherical bins in $|k|$. The size of each bin is equal to the fundamental frequency, $2\pi/L$, where $L = 1 h^{-1} \text{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 \text{ 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.

### 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 \ln M$. The mass of a halo is estimated as

$$M = N m_p, \quad (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
<figure><img src="image_3.png" /><figcaption>Figure 1. The matter power spectrum for the fiducial cosmology.</figcaption></figure>

fluid. The mass of a dark matter particle is thus normalized according to $\Omega_{cb}$, such that

$$m_p = \frac{V\rho_c}{N_p}\Omega_{cb} = \frac{V\rho_c}{N_p}\left(\Omega_m - \frac{M_\nu}{93.14h^2}\right), \quad (8)$$

where $V=L^3$ is the simulation volume, $N_p$ is the total number of dark matter particles in the simulation, and $\rho_c$ is the universe's critical energy density at $z=0$. Thus $m_p=m_p(\Omega_m, M_\nu)$ is a cosmology dependent quantity, which induces noise when computing the derivatives of the HMF with respect to $\Omega_m$ or $M_\nu$ 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\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

$$\mathcal{H} := \frac{dn}{d \ln M} = \frac{dn}{d \ln N}, \quad (9)$$

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

all other cosmological parameters, $\theta$, fixed can be written as

$$\left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{\theta} = \left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{M, \theta} + \left(\frac{\partial \mathcal{H}}{\partial \ln M}\right)_{\theta}\left(\frac{\partial \ln M}{\partial \theta}\right)_{\theta}, \quad (10)$$

or

$$\left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{\theta} = \left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{N, \theta} + \left(\frac{\partial \mathcal{H}}{\partial \ln N}\right)_{\theta}\left(\frac{\partial \ln N}{\partial \theta}\right)_{\theta}. \quad (11)$$

Equating these two equations and rearranging gives

$$\begin{aligned} \left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{M, \theta} &= \left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{N, \theta} + \left(\frac{\partial \mathcal{H}}{\partial \ln N}\right)_{\theta}\left[\frac{\partial \ln N}{\partial \theta} - \frac{\partial \ln M}{\partial \theta}\right]_{\theta} \\ &= \left(\frac{\partial \mathcal{H}}{\partial \theta}\right)_{N, \theta} - \left(\frac{\partial \mathcal{H}}{\partial \ln N}\right)_{\theta}\left(\frac{\partial \ln m_p}{\partial \theta}\right)_{\theta}, \end{aligned} \quad (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 $\theta$, i.e., when $\theta \in \{\Omega_m, M_\nu\}$. Using Equation (8), one finds that

$$\frac{\partial \ln m_p}{\partial \Omega_m} = \frac{1}{\Omega_m}, \quad (13)$$

$$\frac{\partial \ln m_p}{\partial M_\nu} = \frac{1}{\Omega_{cb} 93.14 h^2}, \quad (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 $\Omega_m$ and $M_\nu$, as these require a correction term to transform to fixed-$M$ bins due to the variation of $m_p$. The $\partial \mathcal{H}/\partial \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 $\ln m_p$ with respect to $\theta$ 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 \times 10^{13}$ to $4.6 \times 10^{15} h^{-1} M_\odot$. 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.

### 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^3$ to
<figure><img src="image_4.png" /><figcaption>**Figure 2.** The HMF for the fiducial cosmology.</figcaption></figure>

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^{-1}\mathrm{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 $\delta_{\text{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 $\delta_{\text{th}} = -0.7$, but have checked that results are similar for $\delta_{\text{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_{\text{min}} = 10.4$ and $R_{\text{max}} = 29.9 \ h^{-1}\mathrm{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.

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).

## 5. Results

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 $\text{Corr}(O_{\alpha}, O_{\beta}) := C_{\alpha\beta} / \sqrt{C_{\alpha\alpha} C_{\beta\beta}}$, where $C_{\alpha\beta}$ is the covariance matrix (Equation (4)). First we discuss the correlations for each

<figure><img src="image_5.png" /><figcaption>**Figure 3.** The VSF for the fiducial cosmology.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Correlation matrix for the matter power spectrum ($P_m$, with 72 linear bins and $k_{\text{max}} = 0.5 \ h \ \mathrm{Mpc}^{-1}$), the HMF (15 log bins between $2.0 \times 10^{13}$ and $4.6 \times 10^{16} \ h^{-1}\mathrm{Mpc}$), and the VSF (15 linear bins between 10.4 and 29.9 $h^{-1}\mathrm{Mpc}$), from bottom left to top right. Bin values increase from left to right for each probe. While the HMF shows clear off-block correlation with $P_m$, the VSF is somewhat independent from both $P_m$ and the HMF.</figcaption></figure>

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.
<figure><img src="image_7.png" /><figcaption>**Figure 5.** 68% (darker shades) and 95% (lighter shades) confidence contours for the cosmological parameters for the nonlinear matter power spectrum ($P_m$, red), the HMF (blue), and the VSF (green). Due to the often different degeneracies of each probe, we obtain significantly tighter constraints when combining the three probes (black). We note that some contours extend into unphysical regions ($\Omega_b < 0$, $h < 0$, $M_\nu < 0$): this is just a result of the Gaussian approximation associated with a Fisher analysis.</figcaption></figure>

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
<figure><img src="image_8.png" /><figcaption>**Figure 6.** The $M_{\nu-\sigma_8}$ plane from Figure 5. We inset a zoom-in of the contour obtained by combining all three probes. The marginalized error on $M_{\nu}$ from $P_m$ alone is 0.77 eV, while the error after combining all three probes is 0.018 eV, corresponding to a factor ~43 improvement.</figcaption></figure>

constraints on the parameters are not generally tight when considering any of three probes alone, because we adopt a conservative survey volume of $1(h^{-1} \text{ Gpc})^3$, which is significantly smaller than what is achievable by DESI, $\sim 10^2(h^{-1} \text{ Gpc})^3$.

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 $\Omega_m$ when compared to the other two probes, as the HMF depends nonlinearly on and is highly sensitive to $\Omega_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_{\nu-\sigma_8}$ 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_{\nu}$, 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_{\nu-\sigma_8}$ 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.

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_{\text{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_{\nu}$, where the VSF is the better probe to combine with $P_m$.

While the constraints from $P_m$ alone saturate at around $k_{\text{max}} = 0.2 \text{ h Mpc}^{-1}$ for all parameters, the combined constraints for $M_{\nu}$ (and $\Omega_m$) continue to improve beyond $k_{\text{max}} = 0.5 \text{ h 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_{\text{max}}$ beyond $0.2 \text{ h Mpc}^{-1}$ 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_{\nu}$ 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_{\text{max}} = 0.5 \text{ h 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_{\text{max}} = 0.5 \text{ h 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).

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 $\Omega_m$, $\Omega_b$, $h$, $n_s$, $\sigma_8$, and $M_{\nu}$, respectively. Thus, we achieve 43 times tighter constraints on neutrino mass by combining all three probes. Specifically, the marginalized errors on $M_{\nu}$ are 0.77 eV ($P_m$ alone) and 0.018 eV ($P_m+\text{HMF}+\text{VSF}$). We provide an additional plot in Appendix B to show the confidence ellipses when combining only two of the probes at a time.
<figure><img src="image_9.png" /><figcaption>**Figure 7.** The 1D marginalized error for each of the cosmological parameters as a function of $k_{max}$. We consider four scenarios: $P_m$ alone (red), $P_m$ + HMF (magenta), $P_m$ + VSF (yellow), and $P_m$ + HMF + VSF (black). While the constraints from $P_m$ alone saturate at $k_{max} \approx 0.2\ h\ Mpc^{-1}$, the combined constraints for $M_{\nu}$ (and $\Omega_m$) continue to improve until $k_{max} = 0.5\ h\ Mpc^{-1}$, and likely beyond.</figcaption></figure>

<table><caption>**Table 2**
Marginalized Errors of Cosmological Parameters for $k_{max} = 0.5\ h\ Mpc^{-1}$ Using Different Probe Combinations</caption><thead><tr><th rowspan="2">Probe(s)</th><th colspan="6">Marginalized Fisher Constraints</th></tr><tr><th>Ω<sub>m</sub></th><th>Ω<sub>b</sub></th><th>h</th><th>n<sub>s</sub></th><th>σ<sub>8</sub></th><th>M<sub>ν</sub> (eV)</th></tr></thead><tbody><tr><td>P<sub>m</sub></td><td>0.098</td><td>0.039</td><td>0.51</td><td>0.50</td><td>0.014</td><td>0.77</td></tr><tr><td>HMF</td><td>0.034</td><td>0.042</td><td>0.28</td><td>0.12</td><td>0.082</td><td>1.6</td></tr><tr><td>VSF</td><td>0.31</td><td>0.12</td><td>1.3</td><td>0.42</td><td>0.083</td><td>1.1</td></tr><tr><td>P<sub>m</sub> + HMF</td><td>0.00077</td><td>0.0089</td><td>0.076</td><td>0.034</td><td>0.0016</td><td>0.061</td></tr><tr><td>P<sub>m</sub> + VSF</td><td>0.016</td><td>0.011</td><td>0.12</td><td>0.074</td><td>0.0018</td><td>0.025</td></tr><tr><td>HMF + VSF</td><td>0.0063</td><td>0.037</td><td>0.23</td><td>0.10</td><td>0.0069</td><td>0.096</td></tr><tr><td>P<sub>m</sub> + HMF + VSF (diag)</td><td>0.0015</td><td>0.0088</td><td>0.066</td><td>0.028</td><td>0.00061</td><td>0.031</td></tr><tr><td>P<sub>m</sub> + HMF + VSF (auto)</td><td>0.0015</td><td>0.0086</td><td>0.071</td><td>0.033</td><td>0.0016</td><td>0.025</td></tr><tr><td>P<sub>m</sub> + HMF + VSF (full)</td><td>0.00071</td><td>0.0084</td><td>0.064</td><td>0.025</td><td>0.0015</td><td><strong>0.018</strong></td></tr><tr><td>Multiplicative improvement</td><td>137</td><td>5</td><td>8</td><td>20</td><td>10</td><td>43</td></tr></tbody></table>

**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.

## 6. Discussion and Conclusions

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_{\nu}-\sigma_8$ 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\ Mpc^{-1}$), with the HMF ($M \gtrsim 2 \times 10^{13}\ h^{-1}M_{\odot}$), and the VSF ($R \gtrsim 10.4\ h^{-1}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 ($h^{-1}$ Gpc)$^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 $\Omega_m$. This is driven by the information in the HMF, and gives a marginalized error of $\sigma(\Omega_m) = 7.1 \times 10^{-4}$, which is almost 100 times smaller than the error obtained from a joint large-scale structure analysis by DES Y1 ($\sigma(\Omega_m) \approx 0.04$, To et al. 2021), and 8 times smaller than Planck 2018, ($\sigma(\Omega_m) \approx 5.6 \times 10^{-3}$ (TT,TE, EE+lowE+lensing+BAO), Planck Collaboration et al. 2020). In addition, we found $\sigma(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_v$ 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 (h^{-1} \text{ Gpc})^3$, whereas surveys such as Euclid and DESI will cover volumes of around $10^2 (h^{-1} \text{ Gpc})^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 \text{ 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, $\tau$. 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.

# Appendix A
## Robustness of Results to Numerical Systematics

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),
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Derivatives of the matter power spectrum (top), HMF (middle), and VSF (bottom) with respect to the different cosmological parameters at $z=0$. We show results when the mean values are estimated using 300 (red), 400 (blue), and 500 realizations (black). Solid/dashed lines indicate that the value of the derivative is positive/negative. While the derivatives for the matter power spectrum are well converged already with 300 realizations, more simulations are required for halos and voids.</figcaption></figure>

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.
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
$\Omega_m \rightarrow \Omega_{cb}$ and $P_m \rightarrow P_{cb}$ as neutrinos have negligible contrib-
ution 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
σ_8 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.

Finally, we comment on the choice of finite difference
scheme used to compute the derivative of probes with respect
to $M_\nu$. Throughout the paper we used Equation (6) with
$\delta M_\nu = 0.2$ eV, thus making use of simulations with $M_\nu = 0,
0.2$, and $0.4$ eV. Using this scheme we found the full combined
constraint on $M_\nu$ 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 $\delta M_\nu = 0.1$ eV
and found it to give an identical constraint of 0.018 eV.
Additionally, we tried a forward difference scheme between
$M_\nu = 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><img src="image_11.png" /><figcaption>**Figure 9.** Left: Convergence of all Fisher matrix components as a function of number of simulations used to compute the covariance matrix, $N_{cov}$. Each line shows the ratio between the Fisher matrix elements computed using $N_{cov}$ simulations and 15,000 simulations (as used in the paper). Right: Convergence of all Fisher matrix components as a function of number of simulations used to compute derivatives, $N_{der}$. Each line shows the ratio between the Fisher matrix elements computed using $N_{der}$ simulations and 500 simulations for each cosmology (as used in the paper). In both cases, we plot the Fisher matrix components for $P_m$ (red), the HMF (blue), the VSF (green), and the combined probes (black). The gray bands correspond to the ±5% interval. While there is some noise in the $\sigma_8$ component of the Fisher matrix for $P_m$ as a 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.</figcaption></figure>
**Appendix B**
**Combining Two Probes at a Time**

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.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Same as Figure 5 but for pair combinations of the probes: power spectrum + HMF (magenta), power spectrum + VSF (yellow), HMF + VSF (cyan), and power spectrum + HMF + VSF (black).</figcaption></figure>
ORCID iDs

Adrian E. Bayer https://orcid.org/0000-0002-3568-3900

Francisco Villaescusa-Navarro https://orcid.org/0000-0002-4816-0455

Elena Massara https://orcid.org/0000-0002-0637-8042

Jia Liu https://orcid.org/0000-0001-8219-1995

David N. Spergel https://orcid.org/0000-0002-5151-0006

Licia Verde https://orcid.org/0000-0003-2601-8770

Benjamin D. Wandelt https://orcid.org/0000-0002-5854-8269

Matteo Viel https://orcid.org/0000-0002-2642-5707

References

Abazajian, K. N., Adshead, P., Ahmed, Z., et al. 2016, arXiv:1610.02743

Adamek, J., Durrer, R., & Kunz, M. 2017, JCAP, 2017, 004

Ahmad, Q. R., Allen, R. C., Andersen, T. C., et al. 2002, PhRvL, 89, 011301

Ahn, M. H., Aliu, E., Andringa, S., et al. 2006, PhRvD, 74, 072003

Aker, M., Altenmüller, K., Arenz, M., et al. 2019, PhRvL, 123, 221802

Ali-Haïmoud, Y., & Bird, S. 2013, MNRAS, 428, 3375

Allison, R., Caucal, P., Calabrese, E., Dunkley, J., & Louis, T. 2015, PhRvD, 92, 123535

Allys, E., Marchand, T., Cardoso, J. F., et al. 2020, PhRvD, 102, 103506

An, F. P., Bai, J. Z., Balantekin, B., et al. 2012, PhRvL, 108, 171803

Araki, T., Eguchi, K., Enomoto, S., et al. 2005, PhRvL, 94, 081801

Archidiacono, M., Brinckmann, T., Lesgourgues, J., & Poulin, V. 2017, JCAP, 2017, 052

Archidiacono, M., & Hannestad, S. 2016, JCAP, 2016, 018

Aviles, A., & Banerjee, A. 2020, JCAP, 2020, 034

Banerjee, A., & Abel, T. 2021, MNRAS, 500, 5479

Banerjee, A., & Dalal, N. 2016, JCAP, 11, 015

Banerjee, A., Powell, D., Abel, T., & Villaescusa-Navarro, F. 2018, JCAP, 2018, 028

Bayer, A. E., Banerjee, A., & Feng, Y. 2021, JCAP, 2021, 016

Bergé, J., Amara, A., & Réfrégier, A. 2010, ApJ, 712, 992

Bird, S., Ali-Haïmoud, Y., Feng, Y., & Liu, J. 2018, MNRAS, 481, 1486

Bird, S., Viel, M., & Haehnelt, M. G. 2012, MNRAS, 420, 2551

Blot, L., Corasaniti, P. S., Amendola, L., & Kitching, T. D. 2016, MNRAS, 458, 4462

Bos, E. G. P., van de Weygaert, R., Dolag, K., & Pettorino, V. 2012, MNRAS, 426, 440

Brandbyge, J., & Hannestad, S. 2009, JCAP, 2009, 002

Brandbyge, J., & Hannestad, S. 2010, JCAP, 2010, 021

Brinckmann, T., Hooper, D. S., Archidiacono, M., Lesgourgues, J., & Sprenger, T. 2019, JCAP, 2019, 059

Carbone, C., Petkova, M., & Dolag, K. 2016, JCAP, 2016, 034

Carron, J. 2011, ApJ, 738, 86

Carron, J. 2012, PhRvL, 108, 071301

Carron, J. 2013, A&A, 551, A88

Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K. 2015, JCAP, 2015, 043

Castorina, E., Sefusatti, E., Sheth, R. K., Villaescusa-Navarro, F., & Viel, M. 2014, JCAP, 2014, 049

Chartier, N., Wandelt, B., Akrami, Y., & Villaescusa-Navarro, F. 2021, MNRAS, 503, 1897

Chen, J. Z., Upadhye, A., & Wong, Y. Y. Y. 2020, arXiv:2011.12504

Chen, J. Z., Upadhye, A., & Wong, Y. Y. Y. 2021, JCAP, 2021, 065

Contarini, S., Ronconi, T., Marulli, F., et al. 2019, MNRAS, 488, 3526

Costanzi, M., Villaescusa-Navarro, F., Viel, M., et al. 2013, JCAP, 2013, 012

Cromer, D., Battaglia, N., Miyatake, H., & Simet, M. 2021, arXiv:2104.06925

Dai, J.-P., Verde, L., & Xia, J.-Q. 2020, JCAP, 2020, 007

Dakin, J., Brandbyge, J., Hannestad, S., Haugbølle, T., & Tram, T. 2019, JCAP, 2019, 052

Davies, C. T., Cautun, M., Giblin, B., et al. 2021, MNRAS, 507, 2267

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371
de Salas, P. F., Forero, D. V., Ternes, C. A., Tortola, M., & Valle, J. W. F.
2018, PhLB, 782, 633

Debackere, S. N. B., Schaye, J., & Hoekstra, H. 2021, MNRAS, 505, 593

Doroshkevich, A. G., Khlopov, M. Y., Sunyaev, R. A., Szalay, A. S., & Zeldovich, Y. B. 1981, NYASA, 375, 32

Eisenstein, D. J., & Hu, W. 1999, ApJ, 511, 5

Emberson, J. D., Yu, H.-R., Inman, D., et al. 2017, RAA, 17, 085

Faltenbacher, A., & White, S. D. M. 2010, ApJ, 708, 469

Fisher, R. A. 1925, PCPS, 22, 700

Fukuda, Y., Hayakawa, T., Ichihara, E., et al. 1998, PhRvL, 81, 1562

Gao, L., & White, S. D. M. 2007, MNRAS, 377, L5

Gualdi, D., Novell, S., Gil-Marín, H., & Verde, L. 2021, JCAP, 2021, 015

Hahn, C., & Villaescusa-Navarro, F. 2021, JCAP, 2021, 029

Hahn, C., Villaescusa-Navarro, F., Castorina, E., & Scoccimarro, R. 2020, JCAP, 2020, 040

Haiman, Z., Mohr, J. J., & Holder, G. P. 2001, ApJ, 553, 545

Harnois-Déraps, J., Martinet, N., Castro, T., et al. 2021, MNRAS, 506, 1623

Hazumi, M., Borrill, J., Chinone, Y., et al. 2012, Proc. SPIE, 8442, 844219

Heavens, A. 2009, arXiv:0906.0664

Heavens, A. F., Kitching, T. D., & Verde, L. 2007, MNRAS, 380, 1029

Hu, W., Eisenstein, D. J., & Tegmark, M. 1998, PhRvL, 80, 5255

Ichiki, K., & Takada, M. 2012, PhRvD, 85, 063521

Inman, D., & Pen, U.-L. 2017, PhRvD, 95, 063535

Jennings, E., Li, Y., & Hu, W. 2013, MNRAS, 434, 2167

Kacprzak, T., Kirk, D., Friedrich, O., et al. 2016, MNRAS, 463, 3653

Kayo, I., Takada, M., & Jain, B. 2013, MNRAS, 429, 344

Kreisch, C. D., Pisani, A., Carbone, C., et al. 2019, MNRAS, 488, 4413

Lacerna, I., & Padilla, N. 2012, MNRAS, 426, L26

Lacerna, I., Rodríguez-Puebla, A., Avila-Reese, V., & Hernández-Toledo, H. M. 2014, ApJ, 788, 29

Lesgourgues, J., Mangano, G., Miele, G., & Pastor, S. 2013, Recent Times: Neutrinos and Structure Formation (Cambridge: Cambridge Univ. Press), 273

Li, Y., Hu, W., & Takada, M. 2014, PhRvD, 89, 083519

Liu, A., Pritchard, J. R., Allison, R., et al. 2016, PhRvD, 93, 043013

Liu, J., Bird, S., Matilla, J. M. Z., et al. 2018, JCAP, 2018

Liu, J., & Madhavacheril, M. S. 2019, PhRvD, 99, 083508

Liu, J., Petri, A., Haiman, Z., et al. 2015a, PhRvD, 91, 063507

Liu, X., Pan, C., Li, R., et al. 2015b, MNRAS, 450, 2888

LoVerde, M., & Zaldarriaga, M. 2014, PhRvD, 89, 063502

Martinet, N., & Schneider, P., Hildebrandt, H., et al. 2018, MNRAS, 474, 712

Massara, E., Villaescusa-Navarro, F., Ho, S., Dalal, N., & Spergel, D. N. 2021, PhRvL, 126, 011301

Massara, E., Villaescusa-Navarro, F., Viel, M., & Sutter, P. M. 2015, JCAP, 2015, 018

Nadathur, S. 2016, MNRAS, 461, 358

Paillas, E., Lagos, C. D. P., Padilla, N., et al. 2017, MNRAS, 470, 4434

Paranjape, A., Hahn, O., & Sheth, R. K. 2018, MNRAS, 476, 5442

Pisani, A., Sutter, P. M., Hamaus, N., et al. 2015, PhRvD, 92, 083531

Planck Collaboration, Aghani, N., Akrami, Y., et al. 2020, A&A, 641, A6

Platen, E., Van De Weygaert, R., & Jones, B. J. T. 2007, MNRAS, 380, 551

Platen, E., van de Weygaert, R., & Jones, B. J. T. 2008, MNRAS, 387, 128

Pollina, G., Hamaus, N., Paech, K., et al. 2019, MNRAS, 487, 2836

Rimes, C. D., & Hamilton, A. J. S. 2005, MNRAS, 360, L82

Sahlén, M. 2019, PhRvD, 99, 063525

Saito, S., Takada, M., & Taruya, A. 2008, PhRvL, 100, 191301

Schaan, E., Takada, M., & Spergel, D. N. 2014, PhRvD, 90, 123523

Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2006, PhRvD, 74, 023522

Senatore, L., & Zaldarriaga, M. 2017, arXiv:1707.04698

Shan, H., Liu, X., Hildebrandt, H., et al. 2018, MNRAS, 474, 1116

Sheth, R. K., & Tormen, G. 1999, MNRAS, 308, 119

Sheth, R. K., & Tormen, G. 2002, MNRAS, 329, 61

Shi, J., & Sheth, R. K. 2018, MNRAS, 473, 2486

Shoji, M., & Komatsu, E. 2010, arXiv:1003.0942

Sutter, P. M., Lavaux, G., Hamaus, N., et al. 2015, A&C, 9, 1

Sutter, P. M., Lavaux, G., Wandelt, B. D., & Weinberg, D. H. 2012, ApJ, 761, 187

Takada, M., & Hu, W. 2013, PhRvD, 87, 123504

Takada, M., & Jain, B. 2004, MNRAS, 348, 897

Tegmark, M., Taylor, A. N., & Heavens, A. F. 1997, ApJ, 480, 22

Tinker, J., Kravtsov, A. V., Klypin, A., et al. 2008, ApJ, 688, 709

To, C., Krause, E., Rozo, E., et al. 2021, PhRvL, 126, 141301

Uhlemann, C., Friedrich, O., Villaescusa-Navarro, F., Banerjee, A., & Codis, S. 2020, MNRAS, 495, 4006

Upadhye, A., Kwan, J., Pope, A., et al. 2016, PhRvD, 93, 063515

Verde, L. 2010, Lectures on Cosmology: Accelerated Expansion of the Universe (Berlin: Springer), 147

Verza, G., Pisani, A., Carbone, C., Hamaus, N., & Guzzo, L. 2019, JCAP, 2019, 040

Viel, M., Haehnelt, M. G., & Springel, V. 2010, JCAP, 6, 015

Villaescusa-Navarro, F., Anglés-Alcázar, D., Genel, S., et al. 2021, ApJ, 915, 71

Villaescusa-Navarro, F., Banerjee, A., Dalal, N., et al. 2018, ApJ, 861, 53
Villaescusa-Navarro, F., Bird, S., Peña-Garay, C., & Viel, M. 2013, JCAP, 3, 019

Villaescusa-Navarro, F., Hahn, C., Massara, E., et al. 2020, ApJS, 250, 2

Villaescusa-Navarro, F., Marulli, F., Viel, M., et al. 2014, JCAP, 2014, 011

Villaescusa-Navarro, F., Miralda-Escudé, J., Peña-Garay, C., & Quilis, V. 2011, JCAP, 6, 027

Watts, D. J., Wang, B., Ali, A., et al. 2018, ApJ, 863, 121

Wechsler, R. H., Zentner, A. R., Bullock, J. S., Kravtsov, A. V., & Allgood, B. 2006, ApJ, 652, 71

Yu, B., Knight, R. Z., Sherwin, B. D., et al. 2018, arXiv:1809.02120

Yu, H.-R., Emberson, J. D., Derek, I., et al. 2017, NatAs, 1, 143

Zennaro, M., Bel, J., Villaescusa-Navarro, F., et al. 2017, MNRAS, 466, 3244