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
DDraft version February 11, 2021
Preprint typeset using L A TEX style emulateapj v. 01/23/15
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 , , Berkeley Center for Cosmological Physics, University of California, 341 Campbell Hall, Berkeley, CA 94720, USA Department of Physics, University of California, 366 LeConte Hall, Berkeley, CA 94720, USA Department of Astrophysical Sciences, Princeton University, Peyton Hall, Princeton NJ 08544-0010, USA Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, 10010, New York, NY, USA Waterloo Centre for Astrophysics, University of Waterloo, 200 University Ave W, Waterloo, ON N2L 3G1, Canada Lawrence Berkeley National Laboratory, 1 Cyclotron Road, Berkeley, CA 94720, USA Institut de Ci`encies del Cosmos, University of Barcelona, ICCUB, Barcelona 08028, Spain Institucio Catalana de Recerca i Estudis Avancats, Passeig Lluis Companys 23, Barcelona 08010, Spain Institut d’Astrophysique de Paris, UMR 7095, CNRS, 98 bis boulevard Arago, F-75014 Paris, France Institut Lagrange de Paris, Sorbonne Universites, 98 bis Boulevard Arago, 75014 Paris, France SISSA, Via Bonomea 265, 34136 Trieste, Italy INFN - Istituto Nazionale di Fisica Nucleare, Sezione di Trieste, Via Bonomea 265, 34136 Trieste, Italy INAF - Osservatorio Astronomico di Trieste, via Tiepolo 11, I-34143 Trieste, Italy IFPU - Institute for Fundamental Physics of the Universe, Via Beirut 2, 34014 Trieste, Italy and Department of Physics, Carnegie Mellon University, Pittsburgh, PA 15213, USA
Draft version February 11, 2021
ABSTRACTWe quantify the information content of the non-linear matter power spectrum, the halo massfunction, and the void size function, using the
Quijote
N-body simulations. We find that thesethree statistics exhibit very different degeneracies amongst the cosmological parameters, and thusthe combination of all three probes enables the breaking of degeneracies, in turn yielding remarkablytight 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 comparedto a correlation-free analysis. The multiplicative improvement of the constraints on the cosmologicalparameters obtained by combining all three probes compared to using the power spectrum aloneare: 137, 5, 8, 20, 10, and 43, for Ω m , Ω b , h , n s , σ , and M ν , respectively. The marginalized erroron the sum of the neutrino masses is σ ( M ν ) = 0 . h − Gpc) ,using k max = 0 . h Mpc − , and without CMB priors. We note that this error is an underestimateinsomuch as we do not consider super-sample covariance, baryonic effects, and realistic surveynoises and systematics. On the other hand, it is an overestimate insomuch as our cuts and binningare suboptimal due to restrictions imposed by the simulation resolution. Given upcoming galaxysurveys will observe volumes spanning ∼
100 ( h − Gpc) , this presents a promising new avenue tomeasure neutrino mass without being restricted by the need for accurate knowledge of the opticaldepth, which is required for CMB-based measurements. Furthermore, the improved constraints onother cosmological parameters, notably Ω m , may also be competitive with CMB-based measurements. Keywords: neutrinos, cosmological parameters, large-scale structure of Universe, methods: numerical INTRODUCTION
High-precision measurements of large-scale structurefrom upcoming cosmological surveys, such as DESI , Eu-clid , PFS , Roman Space Telescope , Vera Rubin Ob-servatory , and SKA , are expected to revolutionize ourunderstanding of fundamental physics, for example bymeasuring neutrino mass. To fully realize the potentialof these surveys, an urgent task is to determine the keyobservables that can maximize the scientific return. ForGaussian density fields, the answer is well known: thepower spectrum, or equivalently, the correlation func- * [email protected] † fvillaescusa@flatironinstitute.org https://pfs.ipmu.jp/index.html https://wfirst.gsfc.nasa.gov/index.html tion, is the statistic that completely characterizes thefield. Therefore, on large scales and at high redshift,where the density fluctuation in the Universe resemblesa Gaussian field, the power spectrum encapsulates all theinformation.However, at low redshift and on small scales, non-linear gravitational evolution moves information from thepower spectrum into higher-order moments. It is cur-rently ill-understood which observable(s) will allow re-trieval of the maximum information in the non-linearregime. For instance, it has been shown that for non-Gaussian fields, all clustering information may not beembedded in the infinite N-point statistics (Carron 2011,2012). Since the number of modes increases rapidly bygoing to small scales, it is expected that the amountof information will also increase by considering observ-ables in the mildly to fully non-linear regime. Whilethe amount of information, at least for some parameters,may saturate in the power spectrum (Rimes & Hamil- a r X i v : . [ a s t r o - ph . C O ] F e b Bayer, Villaescusa-Navarro et al. ton 2005; Villaescusa-Navarro et al. 2020b) (see howeverBlot et al. 2016), many authors have shown that otherstatistics contain complementary information (see e.g.Takada & Jain 2004; Sefusatti et al. 2006; Berg´e et al.2010; Kayo et al. 2013; Schaan et al. 2014; Liu et al.2015a,b; Kacprzak et al. 2016; Shan et al. 2018; Mar-tinet et al. 2018; Hahn et al. 2020; Hahn & Villaescusa-Navarro 2020; Dai et al. 2020; Uhlemann et al. 2020;Allys et al. 2020; Gualdi et al. 2020; Harnois-D´eraps et al.2020; Banerjee & Abel 2021; Massara et al. 2021).In this paper we quantify the information embeddedin the non-linear matter power spectrum, the halo massfunction, and the void size function. We apply theFisher formalism using a subset of the
Quijote simula-tions (Villaescusa-Navarro et al. 2020b), comprising of23,000 N-body simulations for 16 different cosmologiesspanning six cosmological parameters: Ω m , Ω b , h , n s , σ , and M ν . We study the information that these probescontain individually and when combined together, show-ing how the combination of these three statistics breaksdegeneracies amongst the cosmological parameters, inturn setting very tight constraints. We consider the ef-fects of both the auto-correlation for each probe and thecross-correlation between different probes when comput-ing the total information content. A simpler, theoret-ical, treatment combining cluster and void abundanceshas been studied by Sahl´en (2019).Of particular interest in this work are constraints onthe sum of the neutrino masses M ν ≡ (cid:80) ν m ν . The firstevidence for neutrino mass came from oscillation exper-iments (Fukuda et al. 1998; Ahmad et al. 2002; Arakiet al. 2005; Ahn et al. 2006; An et al. 2012), which mea-sured the difference in the squares of the masses of thethree neutrino mass eigenstates. The best fit results ob-tained from a joint analysis of oscillation experiments are∆ m ≡ m − m (cid:39) . × − eV from solar neutri-nos, and | ∆ m | ≡ | m − m | (cid:39) . × − eV fromatmospheric neutrinos (de Salas et al. 2018). Since at-mospheric neutrino experiments only probe the magni-tude of the mass difference, there are two possibilities forthe neutrino mass hierarchy: ∆ m >
0, known as thenormal hierarchy, or ∆ m <
0, known as the invertedhierarchy. This gives a lower bound on the sum of theneutrino masses of M ν (cid:38) . M ν (cid:38) . β -decay experiment, m eff ν e (cid:46) . .Neutrinos also play an important role in the Universe’shistory, as the presence of massive neutrinos both shiftsthe time of matter-radiation equality and suppresses thegrowth of structure on small scales. Measuring these ef-fects enables determination of neutrino mass via cosmol-ogy, 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 ofneutrinos on linear (i.e. relatively large) scales are well Single β -decay experiments do not directly measure the neu-trino mass sum, but rather the effective mass of electron neutrinos.In the quasi-degenerate regime where the eigenmasses m i > . i = 1 , , m eff ν e ≈ / M ν . understood theoretically, understanding the effects onnon-linear (i.e. relatively small) scales is an active field ofresearch. There are numerous approaches to obtain the-oretical predictions of the non-linear effects of neutrinos,with varying computational efficiency (see e.g. Saito et al.2008; Brandbyge & Hannestad 2009, 2010; Shoji & Ko-matsu 2010; Viel et al. 2010; Ali-Ha¨ımoud & Bird 2012;Bird et al. 2012, 2018; Costanzi et al. 2013; Villaescusa-Navarro et al. 2014, 2018; Castorina et al. 2014, 2015;Archidiacono & Hannestad 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; Liu et al. 2018; Dakin et al. 2019; Bayeret al. 2021).The current best constraints on M ν arise by consider-ing the cosmic microwave background (CMB) and com-bining it with other cosmological probes. Assuming a Λcold dark matter (ΛCDM) cosmological model, the upperbound on the neutrino mass from the Planck 2018 CMBtemperature and polarization data is M ν < . M ν < . M ν < . A s e − τ , where A s is the amplitude of scalarperturbations and τ is the optical depth of reioniza-tion. Hence, accurate determination of τ is imperative toobtaining tight constraints when combining CMB withclustering/lensing (Allison et al. 2015; Liu et al. 2016;Archidiacono et al. 2017; Brinckmann et al. 2019). Mostupcoming ground-based CMB experiments, such as Si-mons Observatory and CMB-S4, will not observe scaleslarger than (cid:96) ∼
30, and will therefore be unable to di-rectly constrain τ (Abazajian et al. 2016). Planck cur-rently provides the best constraint of τ = 0 . ± . τ ,improved measurements of M ν are expected from galaxysurveys such as DESI, LSST, and Euclid. These surveyswill measure fluctuations on non-linear scales with un-precedented precision. There is thus much motivation toexplore other probes of neutrino mass, beyond the tra-ditional 2-point clustering. By adding probes such asthe halo and void abundances, we demonstrate that it ispossible to break the strong degeneracy between M ν and σ usually seen in 2-point clustering constraints (see e.g.Villaescusa-Navarro et al. 2018). In turn, this gives tightconstraints on neutrino mass, and in fact all cosmologicalparameters, potentially without the need for includingCMB priors. In addition to improved constraints, hav-ing multiple independent probes of neutrino masses willallow for more robust controls of systematics.The paper is organized as follows. We first review the Quijote simulations in Section 2. The Fisher formalismused to quantify the information content on the differentobservables is described in Section 3. We explain how thematter power spectrum, halo mass function, and voidsize function are obtained in Section 4. We show theresults of our analysis in Section 5. Finally, we conclude etecting neutrino mass by combining matter clustering, halos, and voids SIMULATIONS
We quantify the information content of different cos-mological observables using the Fisher matrix formalism.We model the observables using the
Quijote simulations(Villaescusa-Navarro et al. 2020b), a set of 23,000 N-bodysimulations that at a given redshift contain about 8 tril-lion (8 × ) particles over a total combined volumeof 44,100 ( h − Gpc) . Each simulation considers a boxof size 1 ( h − Gpc) . The simulation subset used in thiswork spans a total of 16 different cosmological modelsthat have been designed to evaluate the two ingredientsrequired to compute the Fisher matrix: 1) the covariancematrix of the observables and 2) the derivatives of theobservables with respect to the cosmological parameters.Despite their larger computational cost than analytic ap-proaches (e.g. perturbation theory or the halo model),numerical simulations are more accurate into the fullynon-linear regime and rely on fewer assumptions and ap-proximations.We consider 6 cosmological parameters: Ω m , Ω b , h , n s , σ , and M ν . The set of cosmological parameters isshown in Table 1. To evaluate the covariance matrix,we use the 15,000 simulations of the fiducial cosmology.We compute the derivatives by considering simulationswhere only one cosmological parameter is varied, withall others fixed. We use 1,000 simulations (500 pairs)for each derivative, with the exception of neutrino mass,where we use 1,500 (see below).The initial conditions (ICs) were generated in all casesat z = 127 using second-order Lagrangian PerturbationTheory (2LPT) for simulations with massless neutrinos,by rescaling the z = 0 matter power spectrum using thescale-independent growth factor from linear theory. Be-cause the 2LPT formalism has not yet been developedto account for massive neutrinos, the ICs for massiveneutrino cosmologies adopt the Zel’dovich approximationwith scale-dependent growth factors and rates, followingZennaro et al. (2017). For this reason there is also a‘Fiducial (ZA)’ class of simulations, which is identical tothe fiducial simulations but with Zel’dovich initial con-ditions to match the M ν simulations (see Villaescusa-Navarro et al. 2020b, for further details); this enablesaccurate computation of derivatives with respect to M ν .Note that in the full Quijote simulations there are twosets of Ω b cosmologies; we use the Ω ++ b and Ω −− b set tooobtain smoother derivatives.All simulations follow the evolution of 512 dark mat-ter particles down to z = 0. The simulations withmassive neutrinos also contain 512 neutrino particles.The gravitational force tree for neutrinos is turned onat z = 9. The gravitational softening for both darkmatter and neutrinos is 50 h − kpc (1 /
40 of the meaninter-particle distance). In this work, we consider red-shift z = 0 only. 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 non-linearmatter power spectrum, the halo mass function and thevoid size function, individually and when combined. The Fisher matrix is defined as F ij = − (cid:28) ∂ log L ∂θ i ∂θ j (cid:29) , (1)where L is the likelihood and (cid:126)θ is the vector representingthe parameters of the model (Fisher 1925). Under theassumption that the region around the maximum of thelikelihood can be approximated as a multivariate normaldistribution, one can write the Fisher matrix as F ij = 12 (cid:34) ∂ (cid:126)O∂θ i C − ∂ (cid:126)O T ∂θ j + ∂ (cid:126)O∂θ j C − ∂ (cid:126)O T ∂θ i (cid:35) + 12 Tr (cid:20) C − ∂C∂θ i C − ∂C∂θ j (cid:21) , (2)where (cid:126)O is the vector with the values of the observablesand C is the covariance matrix. In order to avoid un-derestimating the errors, we follow Carron (2013) andneglect the dependence of the covariance on the cosmo-logical parameters, by setting the last term of Eq. 2 tozero. This is necessary when assuming a Gaussian like-lihood. Note that we use Greek (Latin) characters toindex observables (parameters).In this work, the observables and parameters are givenby (cid:126)O = { P m ( k ) , ..., P m ( k A ) , H ( M ) , ..., H ( M B ) ,....................................... V ( R ) , ..., V ( R D ) } ,(cid:126)θ = { Ω m , Ω b , h, n s , σ , M ν } respectively, where P m ( k ) is the matter power spectrumat wavenumber k , H ( M ) is the halo mass function atmass M , and V ( R ) is the void size function at radius R .Note there are a total of A , B , and D bins for the matterpower spectrum, the halo mass function, and the voidsize function respectively, giving a total dimensionalityof A + B + D .We quantify the information content by considering themarginalized error on the cosmological parameters, σ ( θ i ) ≡ (cid:112) ( F − ) ii , (3)which is a lower bound. Covariance matrix
We estimate the covariance matrix using the N cov =15,000 simulations of the fiducial cosmology as C αβ = (cid:104) ( O α − (cid:104) O α (cid:105) ) ( O β − (cid:104) O β (cid:105) ) (cid:105) , (4)where (cid:104)(cid:105) denotes the mean over simulations. This is thelargest number of simulations used for covariance estima-tion to date. We have verified that our combined resultsare converged even with half of the simulations. We showthe results of our convergence tests in Appendix A. Derivatives
For the cosmological parameters Ω m , Ω b , h , n s , and σ ,we approximate the derivatives using a central differencescheme centered on the fiducial cosmology, ∂ (cid:126)O∂θ i (cid:39) (cid:126)O ( θ i + δθ i ) − (cid:126)O ( θ i − δθ i )2 θ i . (5) Bayer, Villaescusa-Navarro et al.
Qujiote SimulationsName Ω m Ω b h n s σ M ν (eV) ICs RealizationsFiducial 0.3175 0.049 0.6711 0.9624 0.834 0.0 2LPT 15,000Fiducial ZA 0.3175 0.049 0.6711 0.9624 0.834 0.0 Zel’dovich 500Ω + m − m ++ b −− b h + h − n + s n − s σ +8 σ − M + ν M ++ ν M +++ ν Table 1
Characteristics of the subset of the
Quijote simulations used in this work. The fiducial cosmology contains 15,000 simulations, that areused to compute the covariance matrix. In the other cosmological models, one parameter is varied at a time, and these simulations areused to compute the numerical derivatives. The initial conditions of all simulations were generated at z = 127 using 2LPT, except for thesimulations with massive neutrinos and a copy of the fiducial cosmology, where the Zel’dovich approximation is used (see main text forfurther details). All realizations follow the evolution of 512 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. Note that only the value of the i th cosmological param-eter is perturbed about its fiducial value, θ i , while thevalues of all other parameters are held fixed. The errorof this approximation is O ( δθ i ).For neutrinos we cannot use Eq. 5 because the fidu-cial model has massless neutrinos, so (cid:126)O ( θ i − δθ i ) wouldcorrespond to a cosmology with negative neutrino mass.We thus compute the derivatives for neutrinos using asecond-order forward difference scheme, ∂ (cid:126)O∂M ν (cid:39) − (cid:126)O ( M ν + 2 δM ν ) + 4 (cid:126)O ( M ν + δM ν ) − (cid:126)O ( M ν )2 δM ν , (6)which has error O ( δM ν ). We exclusively use the M ++ ν and M +++ ν cosmologies in Eq. 6 throughout this work,but we have checked that our results are robust to thechoice of finite difference scheme and massive-neutrinocosmologies.We use a total of N der = 1,000 (500+500) simulationsto compute derivatives when using Eq. 5, and 1,500 whenusing Eq. 6. In Appendix A we show that our results arerobust and converged with this number of simulations. COSMOLOGICAL PROBES
In this section we outline the cosmological observablesconsidered in this work: the matter power spectrum, thehalo mass function, and the void size function.
Matter power spectrum
The first observable we study is the matter power spec-trum. For each realization, the density field is computedby depositing particle masses to a regular grid using thecloud-in-cell mass assignment scheme. In simulations with massive neutrinos we consider both CDM and neu-trino particles when constructing the density field. Thedensity contrast field, δ ( (cid:126)x ) = ρ ( (cid:126)x ) / ¯ ρ −
1, is then Fouriertransformed and the power spectrum is computing byaveraging | δ ( (cid:126)k ) | over spherical bins in | k | . The size ofeach 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 enoughto avoid aliasing effects on the scales of interest for thiswork. In our analysis we consider wavenumbers up to k max = 0 . h Mpc − , using 79 bins. This choice of k max is based on the fact that the clustering of the simulationsis converged at this scale for this mass resolution (seeVillaescusa-Navarro et al. 2020b). We will however showthat using a larger k max would likely lead to even tighterconstraints than the ones we report. Halo mass function
The second observable we consider is the halo massfunction (HMF). Dark matter halos are identified usingthe Friends-of-Friends algorithm (Davis et al. 1985), witha linking length b = 0 .
2. The halo finder considers onlythe dark matter distribution, as the contribution of neu-trinos to the total mass of a halo is expected to be neg-ligible (Villaescusa-Navarro et al. 2011, 2013; Ichiki &Takada 2012; LoVerde & Zaldarriaga 2014).The halo mass function is defined as the comovingnumber 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 , (7)where N is the number of dark matter particles in thehalo and m p is the mass of a single dark matter particle. etecting neutrino mass by combining matter clustering, halos, and voids Quijote simulations, there are only darkmatter and neutrino particles, i.e. dark matter particlesrepresent the CDM+baryon fluid. The mass of a darkmatter particle is thus normalized according to Ω cb , suchthat m p = V ρ c N p Ω cb = V ρ c N p (cid:18) Ω m − M ν . h (cid:19) , (8)where V = L is the simulation volume, N p is the totalnumber 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 ofthe HMF with respect to Ω m or M ν in a fixed mass bin.This is because it is the number of dark matter particlesthat is the fundamental constituent of the halo mass: ahalo with a given number of particles will lie in the same number bin for all cosmologies, whereas it may lie in adifferent mass bin depending on the value of m p . Thisnoise can thus be avoided by instead working with binsof fixed particle number by considering the derivativeof the comoving number density of halos per unit (log)number of particles, dn/d ln N . One can then transformthese derivatives in bins of fixed N to derivatives in binsof fixed M to obtain the derivatives of the halo massfunction.Using the shorthand H to denote the halo mass func-tion, we now derive this transformation. In practice, onemeasures the halo mass function for a fixed cosmology,thus working in logarithmic bins gives H := dnd ln M = dnd ln N , (9)where it is understood that the derivative is taken withfixed cosmological parameters, (cid:126)θ . Explicitly, one canthink of the halo mass function as a function of the cos-mological parameters and halo mass, H ( (cid:126)θ, M ), or the cos-mological parameters and number of particles, H ( (cid:126)θ, N ).Thus the derivative of the HMF with respect to one ofthe cosmological parameters, θ , while holding all othercosmological parameters, /θ , fixed can be written as (cid:18) ∂ H ∂θ (cid:19) /θ = (cid:18) ∂ H ∂θ (cid:19) M,/θ + (cid:18) ∂ H ∂ ln M (cid:19) (cid:126)θ (cid:18) ∂ ln M∂θ (cid:19) /θ , (10)or (cid:18) ∂ H ∂θ (cid:19) /θ = (cid:18) ∂ H ∂θ (cid:19) N,/θ + (cid:18) ∂ H ∂ ln N (cid:19) (cid:126)θ (cid:18) ∂ ln N∂θ (cid:19) /θ . (11)Equating these two equations and rearranging gives (cid:18) ∂ H ∂θ (cid:19) M,/θ = (cid:18) ∂ H ∂θ (cid:19) N,/θ + (cid:18) ∂ H ∂ ln N (cid:19) (cid:126)θ (cid:20) ∂ ln N∂θ − ∂ ln M∂θ (cid:21) /θ = (cid:18) ∂ H ∂θ (cid:19) N,/θ − (cid:18) ∂ H ∂ ln N (cid:19) (cid:126)θ (cid:18) ∂ ln m p ∂θ (cid:19) /θ , (12)where Eq. 7 was used in the final step.The cosmology dependence of m p takes effect in thefinal term of Eq. 12. There is only a difference betweenthe fixed N and fixed M derivative of the HMF when m p depends on θ , i.e. when θ ∈ { Ω m , M ν } . Using Eq. 8, one finds that ∂ ln m p ∂ Ω m = 1Ω m , (13) ∂ ln m p ∂M ν = 1Ω cb . h , (14)where it is understood that all cosmological parametersapart from the one in the derivative are held fixed attheir fiducial values.Thus our procedure to compute derivatives of the HMFusing Eq. 12 is as follows. We first bin the number ofhalos according to the number of dark matter particlesthey contain. We then compute the derivatives for eachfixed- N bin using the equations from Section 3.2, yieldingthe first term on the RHS of Eq. 12. This will be sufficientfor all cosmological parameters except for Ω m and M ν , asthese require a correction term to transform to fixed- M bins due to the variation of m p . The ∂ H /∂ ln N term canbe computed via spline interpolation or by using finitedifference methods between the bins of the halo massfunction of the fiducial cosmology. We have confirmedthe stability of both approaches. Finally, the derivativeof ln m p with respect to θ is computed using Eqs. 13 and14 evaluated at the fiducial values.We consider halos with a number of dark matter par-ticles between 30 and 7,000, using 15 logarithmicallyspaced bins. The corresponding halo mass range is ap-proximately 2 . × to 4 . × h − M (cid:12) . As with thematter power spectrum, this choice of binning and cutsis made to ensure convergence of the derivatives basedon the resolution and number of the simulations avail-able. Hence, using more bins and/or a larger mass rangewould likely lead to stronger constraints than we report. Void size function
We identify voids in the underlying matter field using aspherical void finding algorithm developed by Banerjee &Dalal (2016). First, the density contrast field is smoothedwith a top-hat filter over a large-scale R =53.4 h − Mpc.Next, minima that are smaller than the threshold δ th = − . R , unless they overlap with existing voids. Thisprocedure is repeated for increasingly smaller values of R . We consider voids of radius between R min = 10 . R max = 29 . h − Mpc using 15 linear bins.The void size function (VSF) is then computed as thecomoving number density of voids within a given radiusinterval per unit of radius. Unlike the halo mass function,the VSF is not prone to the changes in particle mass,since the void finder operates directly in the same unitas the VSF. The range of void sizes is limited by ourresolution and the size of our simulated volume. Hence,using more bins and/or a larger range of void sizes wouldlikely lead to stronger constraints than we report. In thiswork we use a threshold of δ th = − .
7, but have checkedthat results are similar for δ th = − . RESULTS
In this section we present the main results of this work.
Full covariance of the probes
In Fig. 1 we show the correlation matrix, defined asCorr( O α , O β ) := C αβ / (cid:112) C αα C ββ , where C αβ is the co- Bayer, Villaescusa-Navarro et al. P m HMF
VSF P m H M F V S F Correlation matrix : P m + HMF + VSF C o rr ( O α , O β ) Figure 1.
Correlation matrix for the matter power spectrum ( P m ,with 72 linear bins and k max = 0 . h Mpc − ), the halo mass func-tion (HMF, 15 log bins between 2 . × and 4 . × h − M (cid:12) ),and the void size function (VSF, 15 linear bins between 10.4 and29.9 h − Mpc), from bottom-left to top-right. Bin values increasefrom left to right for each probe. While the HMF shows clear off-block correlation with P m , the VSF is somewhat independent fromboth P m and the HMF. variance matrix (Eq. 4). First we discuss the correla-tions for each individual probe (auto-correlations). Forthe matter power spectrum (bottom-left region of Fig. 1),we observe some well-known structures: the covariance isalmost diagonal on large scales, while mode-coupling in-duces significant off-diagonal correlations on small scales.For the halo mass function (central region of Fig. 1), thecovariance matrix is almost diagonal, with some smallcorrelations between the different mass bins; the corre-lations are negative for heavy halos, but are positive forthe lightest halos considered in this work. The covari-ance of the void size function (top-right region of Fig. 1)is also almost diagonal, with the abundance of differentvoid sizes slightly anti-correlated with nearby bins.Next, we consider the correlations between differentprobes (cross-correlations). The halo mass functionshows an interesting correlation pattern with the matterpower spectrum: the abundance of the more (less) mas-sive halos shows a ∼
20% correlation (anti-correlation)with small scales of the matter power spectrum. Simi-lar trends are seen between halos and large scales of thematter power spectrum, albeit at a weaker level. On theother hand, voids can be seen to be somewhat indepen-dent of both the matter power spectrum and halos, astheir cross-correlation is (cid:46)
5% for all scales and masses.As discussed in Section 3, we combine the covariancematrix with the numerically computed derivatives to cal-culate the Fisher matrix. The numerical derivatives andrelated numerical convergence tests are shown in Ap-pendix A.
Cosmological constraints
We show the two-dimensional (2D) 68% and 95% confi-dence intervals obtained from our Fisher analysis for eachindividual probe, and the combination of all probes, inFig. 2. The constraints on the parameters are not gener- ally tight when considering any of three probes alone,because we adopt a conservative survey volume of 1( h − Gpc) , which is significantly smaller than what isachievable by DESI, ∼ ( h − Gpc) .The three probes show different degeneracies and aresensitive to each parameter at different levels. Forexample, the halo mass function provides a relativelytight constraint on Ω m when compared to the other twoprobes, as the halo mass function depends non-linearlyon and is highly sensitive to Ω m (see e.g. Haiman et al.2001). The void size function provides weaker constraintsthan the other two probes on almost all parameters, ex-cept for n s compared to P m . Naively, this is not sur-prising, considering the relatively smaller range of scalesbeing probed by the void size function compared to thematter power spectrum. More information could proba-bly be retrieved by using other void-related observables,such as the void-matter correlation function.Because the degeneracies between parameters are oftenvery different for each probe, it is expected that combin-ing the probes will break the degeneracies and in turnyield significantly tighter constraints on the cosmologi-cal parameters than the individual probes do. Indeed,the black ellipses in Fig. 2 show the tight constraints ob-tained by combining the three probes. We emphasizethat these constraints account for all the correlations be-tween the different observables, i.e. by using the full co-variance matrix of Fig. 1.The benefit of combining the three probes is particu-larly well demonstrated in the M ν – σ plane. Because thecombined constraints are too small to be visible in Fig. 2,we zoom in on this plane in Fig. 3. We find that, despitenot being as powerful tools as P m in constraining M ν ,the HMF and VSF both show degeneracies in differentdirections from that of P m , which guaranties that con-straints on the neutrino masses will be largely reducedby combining the three probes. In turn this helps breakthe well-known M ν – σ degeneracy for the matter powerspectrum. We note that the area of these confidencecontours, particularly for the HMF, can potentially bereduced by increasing the bin boundaries and/or by finetuning the binning schemes. Our choice of binning isrestricted by our simulation resolution. We leave theseinvestigations to future works.For a direct comparison to the usual constraints ex-pected from the matter power spectrum, we show the1D marginalized errors (Eq. 3) from different combina-tions of the probes with P m in Fig. 4. We study howthe errors vary with the cutoff scale k max . Combining P m with either the HMF, VSF, or both, can achieve asignificant level of improvement on all 6 parameters. Thecombination with the HMF is typically more beneficialthan the combination with the VSF. The only exceptionis for M ν , where the VSF is the better probe to combinewith P m .While the constraints from P m alone saturate ataround k max = 0 . h Mpc − for all parameters, the com-bined constraints for M ν (and Ω m ) continue to improvebeyond k max = 0 . h Mpc − . This can be explained bythe breaking of degeneracies when combining probes.It was shown in Fig. 5 of Villaescusa-Navarro et al.(2020b) that increasing k max beyond 0 . h Mpc − leadsto a squeezing along the semi-minor axes (i.e, the most etecting neutrino mass by combining matter clustering, halos, and voids Ω b h n s σ Ω m M ν ( e V ) Ω b h n s σ P m : k max = 0 . h Mpc − HMF : M ∈ (2 . × , . × ) h − M fl VSF : R ∈ (10 . , . h − Mpc P m + HMF + VSF Figure 2.
68% (darker shades) and 95% (lighter shades) confidence contours for the cosmological parameters for the non-linear matterpower spectrum ( P m , red), the halo mass function (HMF, blue), and the void size function (VSF, green). Due to the often differentdegeneracies of each probe, we obtain significantly tighter constraints when combining the three probes (black). constraining direction) for the P m ellipses. While thissqueezing has little effect on the marginalized error on M ν from P m alone, its effects are manifest when com-bined with other probes with misaligned contours, re-sulting in significant tightening of constraints. Eventhough the numerical resolution of the Quijote simula-tions prevent us from confidently investigating beyond k max = 0 . h Mpc − , our results hint that even tighterconstraints could be achieved by including smaller scales.In Table 2 we list the errors for k max = 0 . h Mpc − us-ing different probe combinations. We list the constraintsobtained by combining all 3 probes while: 1) only using the diagonals of the covariance matrix (diag), 2) onlyconsidering auto-covariance (auto), and 3) consideringthe full covariance (full). We find that using only the di-agonal components of the covariance matrix, effectivelyignoring both the correlation between the probes and be-tween different bins of the same probe, leads to a factorof 1.7 increase on the error on the neutrino mass. Usingonly block cross-correlations, i.e. ignoring the correlationbetween the probes, leads to a factor of 1.2 increase onthe error on the neutrino mass. Therefore, to obtain thetightest constraints, it is crucial to model the full covari-ance matrix. It is interesting to note that when consider- Bayer, Villaescusa-Navarro et al. σ M ν ( e V ) P m HMFVSFALL
Figure 3.
The M ν – σ plane from Fig. 2. We inset a zoom-in ofthe contour obtained by combining all three probes. The marginal-ized error on M ν from P m alone is 0 . . ∼
43 improvement. ing the matter power spectrum alone, correlations causean increase in errors due to the positive correlation be-tween different scales (see Fig. 1). However, it is the com-plex correlation structure, notably the anti-correlations,introduced by considering the HMF and VSF that leadsto a reduction in error, both for the HMF and VSF in-dividually, and in turn when combining all probes. Theassociation of anti-correlation with the tightening of con-straints was also pointed out by Chartier et al. (2020).In Table 2 we quantify the improvement of the com-bined 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 ν , re-spectively. Thus we achieve 43 times tighter constraintson neutrino mass by combining all three probes. Specif-ically, the marginalized errors on M ν are 0 . P m alone) and 0 . P m +HMF+VSF). We provide anadditional plot in Appendix B to show the confidence el-lipses when combining only two of the probes at a time. DISCUSSION AND CONCLUSIONS
Upcoming galaxy surveys will map large volumes ofthe Universe at low redshifts, with the potential to dras-tically improve our understanding of the underlying cos-mological model. With the unprecedentedly precisionachievable by these surveys, it is expected that a verylarge amount of cosmological (and astrophysical) infor-mation will lie in the mildly to fully non-linear regime,where analytic methods are often intractable. It remainsan open question which observable(s) will lead to thetightest bounds on the cosmological parameters.In this paper, we use the
Quijote simulations, based onthe Fisher formalism, to quantify the information contentembedded in the non-linear matter power spectrum, thehalo mass function, and the void size function, both in-dividually and when combined, at z = 0. We find that the HMF and VSF have different degeneracies to eachother and to the matter power spectrum, particularly inthe M ν – σ plane (Figs. 2 & 3). In terms of measuringneutrino mass, we find the void size function to be themore complementary probe to combine with the matterpower spectrum. This is consistent with findings thatvoid properties are particularly sensitive to matter com-ponents that are less clustered, such as neutrinos (Mas-sara et al. 2015; Kreisch et al. 2019).By combining the non-linear matter power spectrum( k max = 0 . h Mpc − ), with the halo mass function( M (cid:38) × h − M (cid:12) ), and the void size function( R (cid:62) . h − Mpc), we achieve significantly tighter con-straints on the cosmological parameters compared to P m alone (Fig. 4). In particular, we find that with a vol-ume of just 1 ( h − Gpc) , the error on the sum of neu-trino masses from the combined probes is at the 0 . . m . This is driven bythe information in the HMF, and gives a marginal-ized error of σ (Ω m ) = 7 . × − , which is 8times smaller than the error obtained by Planck 2018(TT,TE,EE+lowE+lensing+BAO), 5 . × − (PlanckCollaboration et al. 2018). In addition, we found σ ( h ) =0 .
064 by combining the three probes, which is 8 timestighter than the constraints from the matter power spec-trum alone. This could provide a new angle to investigatethe Hubble tension.There are several caveats in this work. Firstly, we as-sumed perfect knowledge of the three-dimensional spa-tial distribution of the underlying matter field in real-space. However, in reality, one observes either tracersof the matter field in redshift-space, or projected mat-ter fields through lensing. Therefore, additional linksmust be made to bridge the galaxy–matter connectionand the 2D lensing–3D matter distribution gaps. In ad-dition, our simulations consider only gravitational inter-actions and hence ignore baryonic effects that can impactthe small-scale matter distribution (Villaescusa-Navarroet al. 2020a). We also neglected super-sample covari-ance. These effects will increase our quoted errors, andpotentially change the relative improvements we discov-ered from combining the probes. We plan to addressthese issues in future works.On the other hand, the constraints obtained here maybe overly conservative due to the limited number andresolution of simulations available. Firstly, this meansthat the number of bins used are likely suboptimal. Sec-ond, applying more aggressive bounds on the observ-ables, e.g. a higher k max , a larger halo mass range, or alarger void size range, would likely also reduce the com-bined constraints. Third, we only considered a single red-shift z = 0: combining multiple redshifts should furthertighten the constraints. Fourth, we considered a volumeof only 1 ( h − Gpc) , whereas surveys such as Euclid andDESI will cover volumes of around 10 ( h − Gpc) , so,conservatively, the error on the parameters will shrinkby a factor of 1 / √ = 0 .
1. Finally, we have only con- etecting neutrino mass by combining matter clustering, halos, and voids -3 -2 -1 σ ( Ω m ) -2 -1 σ ( Ω b ) -1 σ ( h ) P m P m + H P m + V P m + H + V k max [ h Mpc − ] -1 σ ( n s ) k max [ h Mpc − ] -2 σ ( σ ) k max [ h Mpc − ] -1 σ ( M ν ( e V )) Figure 4.
The 1D marginalized error for each of the cosmological parameters as a function of k max . We consider 4 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 (cid:39) . h Mpc − , the combined constraints for M ν (and Ω m ) continue to improve until k max = 0 . h Mpc − , and likely beyond. Marginalized Fisher ConstraintsProbe(s) Ω m Ω b h n s σ M ν (eV) P m P m + HMF 0.00077 0.0089 0.076 0.034 0.0016 0.061 P m + VSF 0.016 0.011 0.12 0.074 0.0018 0.025HMF + VSF 0.0063 0.037 0.23 0.10 0.0069 0.096 P m + HMF + VSF (diag) 0.0015 0.0088 0.066 0.028 0.00061 0.031 P m + HMF + VSF (auto) 0.0015 0.0086 0.071 0.033 0.0016 0.025 P m + HMF + VSF (full) 0.00071 0.0084 0.064 0.025 0.0015 Multiplicative improvement 137 5 8 20 10 43
Table 2
Marginalized errors of cosmological parameters for k max = 0 . h Mpc − using different probe combinations. Note, we list the constraintsobtained by combining all 3 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 also list the multiplicative improvement in the constraints from the full covariancecompared to those from P m alone. sidered three probes; using the same observations, onecan derive other statistics such as the void profile, BAO,and redshift space distortions, which could be combinedwith the statistics considered here to further break de-generacies.We have demonstrated that combining multiple probesof cosmological structure using their full covariance ma-trix provides remarkably tight constraints on the cos-mological parameters. Our approach opens a promisingpathway to measure neutrino mass, potentially withoutrelying on CMB-based measurements which require ac-curate knowledge of the optical depth, τ . In addition,comparing constraints from different combinations of ob- servables, e.g. CMB+ P m and P m +HMF+VSF, will helpidentify systematic issues and provide robust evidencefor any discovery. We have shown that there is sufficientinformation in upcoming galaxy surveys to measure theneutrino mass sum at the minimum mass 0 . ACKNOWLEDGEMENTS
We thank Ravi Sheth for useful conversations onthe early stages of this project. The
Quijote Bayer, Villaescusa-Navarro et al. simulations can be found at https://github.com/franciscovillaescusa/Quijote-simulations . Theanalysis of the simulations has made use of the
Pylians libraries, publicly available at https://github.com/franciscovillaescusa/Pylians3 . The work of DS,BW, and SH is supported by the Simons Foundation. LV acknowledges support from the European Union Horizon2020 research and innovation program ERC (BePreSySe,grant agreement 725327) and MDM-2014-0369 of ICCUB(Unidad de Excelencia Maria de Maeztu). MV is sup-ported by INFN PD51 Indark grant and by the agree-ment 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 thecovariance matrix and derivatives. In Fig. 5 we show the derivatives of the matter power spectrum (top), halo massfunction (middle), and void size function (bottom) with respect to the cosmological parameters when using a differentnumber of realizations. For the matter power spectrum, the derivatives are already converged when the mean valuesfor each model are computed with 300 realizations. Results are slightly noisier for the halo mass function and the voidsize function, but still sufficiently converged by 500 realizations.Next, we comment on the convergence of our simulated results with theory. The convergence of the matter powerspectrum in
Quijote has been thoroughly tested (Villaescusa-Navarro et al. 2020b; Aviles & Banerjee 2020; Hahn et al.2020). For the void size function there is no theoretical formula accurate enough to compute derivatives, but we havechecked results are robust to the parameters used in the void finder. Therefore, we only compare our measured HMFto theoretical predictions. For the HMF, we plot the theoretical predictions of Sheth-Tormen (ST) (Sheth & Tormen1999; Sheth & Tormen 2002) and Tinker (Tinker et al. 2008). We use the prescription of Costanzi et al. (2013) in thecase of massive neutrino cosmologies by replacing Ω m → Ω cb and P m → P cb as neutrinos have negligible contributionto halo mass. There is good agreements between these predictions and Quijote . We have also checked that there isgood agreement for different choice of step size (not shown). Note that these theoretical formulae provide a guidelinerather than exact predictions, as they were fitted to simpler simulations or calibrated on spherical overdensity halos,as opposed to FoF here.In Fig. 6 we show the convergence of the Fisher matrix elements with respect to the number of realizations usedto compute the covariance, N cov , and derivatives N der . We consider the Fisher matrix components for P m (red), theHMF (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 isachieved by 15000. Likewise the Fisher matrix is well converged as a function of N der . Crucially, the Fisher matrixelements for the combined probes (black) all show good convergence. Note that when combining probes we scale thepower spectrum by a factor of 10 − to ensure that the condition number of the covariance matrix is sufficiently lowfor accurate inversion.Given these results, we believe that our conclusions are robust against potential numerical systematics. We noteagain that our bin configuration has been chosen with these results in mind, to ensure sufficiently converged derivativesand Fisher matrix components, but in principle one could consider more bins over a wider range to potentially obtaintighter constraints. B. COMBINING TWO PROBES AT A TIME
Fig. 7 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 halo mass function with the void size function are the weakest,indicating that it is important to use the information from the non-linear matter power spectrum to break degeneracies. etecting neutrino mass by combining matter clustering, halos, and voids P m / Ω m
300 realizations400 realizations500 realizations P m / Ω b P m / h -2 -1 k [ h Mpc − ] P m / n s -2 -1 k [ h Mpc − ] P m / σ -2 -1 k [ h Mpc − ] × P m / M ν -7 -6 -5 -4 -3 H M F / Ω m STTinker
300 realizations400 realizations500 realizations -7 -6 -5 -4 -3 H M F / Ω b -7 -6 -5 -4 -3 H M F / h M [ h − M fl ] -7 -6 -5 -4 -3 H M F / n s M [ h − M fl ] -7 -6 -5 -4 -3 H M F / σ M [ h − M fl ] -7 -6 -5 -4 -3 × H M F / M ν -10 -9 -8 -7 -6 -5 V S F / Ω m
300 realizations400 realizations500 realizations -10 -9 -8 -7 -6 -5 V S F / Ω b -10 -9 -8 -7 -6 -5 V S F / h R [ h − Mpc] -10 -9 -8 -7 -6 -5 V S F / n s R [ h − Mpc] -10 -9 -8 -7 -6 -5 V S F / σ R [ h − Mpc] -10 -9 -8 -7 -6 -5 × V S F / M ν Figure 5.
Derivatives of the matter power spectrum (top), halo mass function (middle) and void size function (bottom) with respect tothe different cosmological parameters at z = 0. We show results when the mean values are estimated using 300 (red), 400 (blue) and 500realizations (black). Solid/dashed lines indicate that the value of the derivative is positive/negative. While the derivatives for the matterpower spectrum are well converged already with 300 realizations, more simulations are required for halos and voids. Bayer, Villaescusa-Navarro et al. N cov Ω m , Ω m Ω m , Ω b Ω m , h Ω m , n s Ω m , σ Ω m , M ν Ω b , Ω b Ω b , h Ω b , n s Ω b , σ Ω b , M ν h , hh , n s h , σ h , M ν n s , n s n s , σ n s , M ν σ , σ σ , M ν M ν , M ν F ij ( N cov ) /F ij ( N cov = 15000) P m HMF VSF ALL
250 300 350 400 450 500 N der Ω m , Ω m Ω m , Ω b Ω m , h Ω m , n s Ω m , σ Ω m , M ν Ω b , Ω b Ω b , h Ω b , n s Ω b , σ Ω b , M ν h , hh , n s h , σ h , M ν n s , n s n s , σ n s , M ν σ , σ σ , M ν M ν , M ν F ij ( N der ) /F ij ( N der = 500) P m HMF VSF ALL
Figure 6.
Left:
Convergence of all Fisher matrix components as a function of number of simulations used to compute the covariancematrix, 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 computederivatives, N der . Each line shows the ratio between the Fisher matrix elements computed using N der simulations and 500 simulationsfor 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 σ componentof 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 asa function of N der . Crucially, the Fisher matrix elements for the combined probes (black) all show good convergence. etecting neutrino mass by combining matter clustering, halos, and voids Ω b h n s σ .
30 0 .
31 0 .
32 0 .
33 0 . Ω m M ν ( e V ) .
000 0 .
025 0 .
050 0 .
075 0 . Ω b . . . . h . . . n s .
825 0 .
830 0 .
835 0 . σ P m + HMF P m + VSFHMF + VSF P m + HMF + VSF Figure 7.
Same as Fig. 2 but for pair combinations of the probes: power spectrum + halo mass function (magenta), power spectrum +voids size function (yellow), halo mass function + void size function (cyan) and power spectrum + halo mass function + void size function(black). Bayer, Villaescusa-Navarro et al.
REFERENCES
Abazajian, K. N. et al. 2016, CMB-S4 Science Book, FirstEdition, [arXiv:1610.02743]Adamek, J., Durrer, R., & Kunz, M. 2017, Journal of Cosmologyand Astroparticle Physics, 2017, 004–004Ahmad, Q. R., et al. 2002, Phys. Rev. Lett., 89, 011301,[arXiv:nucl-ex/0204008]Ahn, M. H. et al. 2006, Phys. Rev. D, 74, 072003Aker, M. et al. 2019, Physical Review Letters, 123Ali-Ha¨ımoud, Y., & Bird, S. 2012, Monthly Notices of the RoyalAstronomical Society, 428, 3375–3389Allison, R., Caucal, P., Calabrese, E., Dunkley, J., & Louis, T.2015, Phys. Rev. D, 92, 123535Allys, E., Marchand, T., Cardoso, J. F., Villaescusa-Navarro, F.,Ho, S., & Mallat, S. 2020, Phys. Rev. D, 102, 103506,[arXiv:2006.06298]An, F. P. et al. 2012, Phys. Rev. Lett., 108, 171803Araki, T., et al. 2005, Phys. Rev. Lett., 94, 081801,[arXiv:hep-ex/0406035]Archidiacono, M., Brinckmann, T., Lesgourgues, J., & Poulin, V.2017, Journal of Cosmology and Astroparticle Physics, 2017,052Archidiacono, M., & Hannestad, S. 2016, Journal of Cosmologyand Astroparticle Physics, 2016, 018–018Aviles, A., & Banerjee, A. 2020, Journal of Cosmology andAstroparticle Physics, 2020, 034–034Banerjee, A., & Abel, T. 2021, MNRAS, 500, 5479,[arXiv:2007.13342]Banerjee, A., & Dalal, N. 2016, J. Cosmology Astropart. Phys.,11, 015, [arXiv:1606.06167]Bayer, A. E., Banerjee, A., & Feng, Y. 2021, Journal ofCosmology and Astroparticle Physics, 2021, 016Berg´e, J., Amara, A., & R´efr´egier, A. 2010, ApJ, 712, 992,[arXiv:0909.0529]Bird, S., Ali-Ha¨ımoud, Y., Feng, Y., & Liu, J. 2018, MonthlyNotices of the Royal Astronomical Society, 481, 1486,[arXiv:https://academic.oup.com/mnras/article-pdf/481/2/1486/25716087/sty2376.pdf]Bird, S., Viel, M., & Haehnelt, M. G. 2012, Monthly Notices ofthe Royal Astronomical Society, 420, 2551–2561Blot, L., Corasaniti, P. S., Amendola, L., & Kitching, T. D. 2016,MNRAS, 458, 4462, [arXiv:1512.05383]Brandbyge, J., & Hannestad, S. 2009, Journal of Cosmology andAstroparticle Physics, 2009, 002–002——. 2010, Journal of Cosmology and Astroparticle Physics,2010, 021–021Brinckmann, T., Hooper, D. C., Archidiacono, M., Lesgourgues,J., & Sprenger, T. 2019, Journal of Cosmology andAstroparticle Physics, 2019, 059–059Carbone, C., Petkova, M., & Dolag, K. 2016, Journal ofCosmology and Astroparticle Physics, 2016, 034–034Carron, J. 2011, ApJ, 738, 86, [arXiv:1105.4467]Carron, J. 2012, Physical Review Letters, 108Carron, J. 2013, A&A, 551, A88, [arXiv:1204.4724]Castorina, E., Carbone, C., Bel, J., Sefusatti, E., & Dolag, K.2015, Journal of Cosmology and Astroparticle Physics, 2015,043–043Castorina, E., Sefusatti, E., Sheth, R. K., Villaescusa-Navarro, F.,& Viel, M. 2014, Journal of Cosmology and AstroparticlePhysics, 2014, 049–049Chartier, N., Wandelt, B., Akrami, Y., & Villaescusa-Navarro, F.2020, CARPool: fast, accurate computation of large-scalestructure statistics by pairing costly and cheap cosmologicalsimulations, [arXiv:2009.08970]Collaboration, P. et al. 2018, Planck 2018 results. VI.Cosmological parameters, [arXiv:1807.06209]Costanzi, M., Villaescusa-Navarro, F., Viel, M., Xia, J.-Q.,Borgani, S., Castorina, E., & Sefusatti, E. 2013, Journal ofCosmology and Astroparticle Physics, 2013, 012–012Dai, J.-P., Verde, L., & Xia, J.-Q. 2020, J. Cosmology Astropart.Phys., 2020, 007, [arXiv:2002.09904]Dakin, J., Brandbyge, J., Hannestad, S., HaugbØlle, T., & Tram,T. 2019, Journal of Cosmology and Astroparticle Physics, 2019,052Davis, 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, Phys. Lett., B782, 633, [arXiv:1708.01186]Doroshkevich, A. G., Khlopov, M. Y., Sunyaev, R. A., Szalay,A. S., & Zeldovich, Y. B. 1981, Annals of the New YorkAcademy of Sciences, 375, 32Eisenstein, D. J., & Hu, W. 1999, The Astrophysical Journal,511, 5–15Emberson, J. D. et al. 2017, Research in Astronomy andAstrophysics, 17, 085Fisher, R. A. 1925, Mathematical Proceedings of the CambridgePhilosophical Society, 22, 700–725Fukuda, Y. et al. 1998, Physical Review Letters, 81, 1562–1567Gualdi, D., Novell, S., Gil-Mar´ın, H., & Verde, L. 2020, arXive-prints, arXiv:2009.02290, [arXiv:2009.02290]Hahn, C., & Villaescusa-Navarro, F. 2020, Constraining M ν withthe Bispectrum II: The Total Information Content of theGalaxy Bispectrum, [arXiv:2012.02200]Hahn, C., Villaescusa-Navarro, F., Castorina, E., & Scoccimarro,R. 2020, Journal of Cosmology and Astroparticle Physics, 2020,040–040Haiman, Z., Mohr, J. J., & Holder, G. P. 2001, ApJ, 553, 545,[arXiv:astro-ph/0002336]Harnois-D´eraps, J., Martinet, N., Castro, T., Dolag, K., Giblin,B., Heymans, C., Hildebrandt, H., & Xia, Q. 2020, arXive-prints, arXiv:2012.02777, [arXiv:2012.02777]Hazumi, M. et al. 2012, in Space Telescopes and Instrumentation2012: Optical, Infrared, and Millimeter Wave, ed. M. C.Clampin, G. G. Fazio, H. A. MacEwen, & J. M. O. Jr., Vol.8442, International Society for Optics and Photonics (SPIE),451 – 459Heavens, A. 2009, arXiv e-prints, arXiv:0906.0664,[arXiv:0906.0664]Heavens, A. F., Kitching, T. D., & Verde, L. 2007, MNRAS, 380,1029, [arXiv:astro-ph/0703191]Hu, W., Eisenstein, D. J., & Tegmark, M. 1998, Physical ReviewLetters, 80, 5255–5258Ichiki, K., & Takada, M. 2012, Phys. Rev. D, 85, 063521,[arXiv:1108.4688]Inman, D., & Pen, U.-L. 2017, Physical Review D, 95Kacprzak, T. et al. 2016, ArXiv e-prints, [arXiv:1603.05040]Kayo, I., Takada, M., & Jain, B. 2013, MNRAS, 429, 344,[arXiv:1207.6322]Kreisch, C. D., Pisani, A., Carbone, C., Liu, J., Hawken, A. J.,Massara, E., Spergel, D. N., & Wandelt, B. D. 2019, MNRAS,488, 4413, [arXiv:1808.07464]Lesgourgues, J., Mangano, G., Miele, G., & Pastor, S. 2013,Recent times: neutrinos and structure formation (CambridgeUniversity Press), 273–347Liu, A., Pritchard, J. R., Allison, R., Parsons, A. R., Seljak, U.c. v., & Sherwin, B. D. 2016, Phys. Rev. D, 93, 043013Liu, J., Bird, S., Matilla, J. M. Z., Hill, J. C., Haiman, Z.,Madhavacheril, M. S., Spergel, D. N., & Petri, A. 2018, Journalof Cosmology and Astroparticle Physics, 2018Liu, J., Petri, A., Haiman, Z., Hui, L., Kratochvil, J. M., & May,M. 2015a, Phys. Rev. D, 91, 063507, [arXiv:1412.0757]Liu, X. et al. 2015b, MNRAS, 450, 2888, [arXiv:1412.3683]LoVerde, M., & Zaldarriaga, M. 2014, Phys. Rev. D, 89, 063502,[arXiv:1310.6459]Martinet, N. et al. 2018, MNRAS, 474, 712, [arXiv:1709.07678]Massara, E., Villaescusa-Navarro, F., Ho, S., Dalal, N., & Spergel,D. N. 2021, Phys. Rev. Lett., 126, 011301, [arXiv:2001.11024]Massara, E., Villaescusa-Navarro, F., Viel, M., & Sutter, P. M.2015, J. Cosmology Astropart. Phys., 2015, 018,[arXiv:1506.03088]Planck Collaboration et al. 2018, ArXiv e-prints,[arXiv:1807.06209]Rimes, C. D., & Hamilton, A. J. S. 2005, MNRAS, 360, L82,[arXiv:astro-ph/0502081]Sahl´en, M. 2019, Physical Review D, 99Saito, S., Takada, M., & Taruya, A. 2008, Physical ReviewLetters, 100Schaan, E., Takada, M., & Spergel, D. N. 2014, Phys. Rev. D, 90,123523, [arXiv:1406.3330]Sefusatti, E., Crocce, M., Pueblas, S., & Scoccimarro, R. 2006,Phys. Rev. D, 74, 023522, [arXiv:astro-ph/0604505] etecting neutrino mass by combining matter clustering, halos, and voids15