# Neutrino Mass Constraint from an Implicit Likelihood Analysis of BOSS Voids

Leander Thiele<sup>1,2,3</sup> <sup>1</sup>ID, Elena Massara<sup>4,5</sup>, Alice Pisani<sup>6,7,8</sup> <sup>1</sup>ID, ChangHoon Hahn<sup>8</sup> <sup>1</sup>ID, David N. Spergel<sup>7</sup> <sup>1</sup>ID, Shirley Ho<sup>7</sup>, and Benjamin Wandelt<sup>7,9</sup> <sup>1</sup>ID

<sup>1</sup>Department of Physics, Princeton University, Princeton, NJ 08544, USA; leander.thiele@ipmu.jp

<sup>2</sup>Kavli Institute for the Physics and Mathematics of the Universe (WPI), The University of Tokyo Institutes for Advanced Study (UTIAS), The University of Tokyo, 5-1-5 Kashiwanoha, Kashiwa-shi, Chiba 277-8583, Japan

<sup>3</sup>Center for Data-Driven Discovery (CD3), Kavli IPMU (WPI), UTIAS, The University of Tokyo, Kashiwa-shi, Chiba 277-8583, Japan

<sup>4</sup>Waterloo Centre for Astrophysics, University of Waterloo, 200 University Avenue W, Waterloo, ON N2L 3G1, Canada

<sup>5</sup>Department of Physics and Astronomy, University of Waterloo, 200 University Avenue W, Waterloo, ON N2L 3G1, Canada

<sup>6</sup>The Cooper Union for the Advancement of Science and Art, 41 Cooper Square, New York, NY 10003, USA

<sup>7</sup>Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

<sup>8</sup>Department of Astrophysical Sciences, Princeton University, Princeton, NJ 08544, USA

<sup>9</sup>Institut d'Astrophysique de Paris (IAP), UMR 7095, CNRS, Sorbonne Université, Paris, France

Received 2023 August 14; revised 2024 March 29; accepted 2024 April 23; published 2024 July 2

## Abstract

Cosmic voids identified in the spatial distribution of galaxies provide complementary information to two-point statistics. In particular, constraints on the neutrino mass sum, $\sum m_{\nu}$, promise to benefit from the inclusion of void statistics. We perform inference on the CMASS Northern Galactic Cap sample of the third-generation Sloan Digital Sky Survey/Baryon Oscillation Spectroscopic Survey (BOSS) with the aim of constraining $\sum m_{\nu}$. We utilize the void size function, the void-galaxy cross power spectrum, and the galaxy auto power spectrum. To extract constraints from these summary statistics we use a simulation-based approach, specifically implicit likelihood inference. We populate approximate gravity-only, particle-neutrino cosmological simulations with an expressive halo occupation distribution model. With a conservative scale cut of $k_{max} = 0.15 h \text{ Mpc}^{-1}$ and a Planck-inspired Lambda cold dark matter prior, we find upper bounds on $\sum m_{\nu}$ of 0.43 and 0.35 eV from the galaxy auto power spectrum and the full data vector, respectively (95% credible interval). Relaxing the scale cut to $k_{max} = 0.2 h \text{ Mpc}^{-1}$, we find 0.28 and 0.32 eV. Generally, void statistics appear to add modest but nonzero information beyond the auto power spectrum. We also substantiate the usual assumption that the void size function is Poisson distributed.

Unified Astronomy Thesaurus concepts: Redshift surveys (1378); Voids (1779); Cosmological neutrinos (338); Neutrino masses (1102); N-body simulations (1083); Neural networks (1933)

## 1. Introduction

The Universe's ability to provide glimpses into experimentally inaccessible conditions has a long history, including the deduction of the laws of gravity and the discovery of helium. In the present day, cosmology offers a unique view on the properties of neutrinos, which are among the last unknowns in the standard model of particle physics. First evidence for a nonzero neutrino mass sum, $\sum m_{\nu}$, came from the solar neutrino problem (Bahcall & Davis Jr., 1976; Wolfenstein 1978; Mikheyev & Smirnov 1985). Subsequently, oscillation experiments provided proof that neutrinos must have mass (Super-Kamiokande Collaboration et al. 1998; SNO Collaboration et al. 2002; KamLAND Collaboration et al. 2005; K2K Collaboration et al. 2006; Daya Bay Collaboration et al. 2012) and established lower bounds of 0.06 and 0.1 eV in the normal and inverted hierarchy, respectively. The terrestrial experiment KATRIN currently sets an upper bound of $\ge 0.8$ eV (KATRIN Collaboration et al. 2022).<sup>10</sup> However, the strongest upper bounds are already provided by cosmological data, the

primary cosmic microwave background (CMB) alone giving 0.38 eV (Planck Collaboration et al. 2020), for example. It will be one of cosmology's primary goals in the coming decade to tighten this bound and eventually detect neutrino mass.

One of the natural regimes to look at to constrain $\sum m_{\nu}$ is extremely underdense regions, cosmic voids (Icke 1984; Pisani et al. 2019; Moresco et al. 2022; Schuster et al. 2023). As the cold dark matter (CDM) flows out of the voids and into filaments and clusters, neutrinos are more smoothly distributed. Thus, the neutrino-to-CDM ratio is higher in the voids and lower in the clusters. These qualitative considerations have spawned considerable theoretical interest in the use of void properties to constrain $\sum m_{\nu}$. This includes simulated data vector-level investigations (Massara et al. 2015; Banerjee & Dalal 2016; Kreisch et al. 2019; Schuster et al. 2019; Contarini et al. 2021; Verza et al. 2023) as well as forecasts (Sahlén 2019; Bayer et al. 2021b; Kreisch et al. 2022). The forecasts find promising error bars on $\sum m_{\nu}$, albeit under simplifying assumptions. Voids may be the first regime in which nonlinear signatures of massive neutrinos will be observed (also see Hotinli et al. 2023).

Being large objects, voids had to wait for the era of relatively deep, large-volume surveys with approximately uniform selection function to be statistically usable. While the original detections focused on individual objects (Gregory & Thompson 1978; Jöeveer et al. 1978; Tully & Fisher 1978;

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

<footnote n="10"><sup>10</sup> The KATRIN bound is on $m_{\nu}^2 \equiv \sum_{\nu} |U_{e\nu}^{\text{PMNS}}|^2 m_{\nu}^2$, so it only equals a bound on $\sum m_{\nu}$ for a special, experimentally excluded, choice of the PMNS matrix and the mass hierarchy. In general, the bound is weaker.</footnote>
Kirshner et al. 1981; de Lapparent et al. 1986), which already contain cosmological information (Sahlén et al. 2016), we are now able to utilize catalogs of hundreds and thousands of voids (Hoyle & Vogeley 2004; Sutter et al. 2012a, 2014a; Pan et al. 2012; Nadathur 2016; Mao et al. 2017b) to perform precision cosmology with void shapes (Sutter et al. 2012b, 2014b; Hamaus et al. 2016, 2017, 2020; Mao et al. 2017a; Nadathur et al. 2020; Aubert et al. 2022; Woodfinden et al. 2023) and sizes (Contarini et al. 2023, 2024).

In this work, we use voids identified in the CMASS sample of the third-generation Sloan Digital Sky Survey (SDSS-III) Baryon Oscillation Spectroscopic Survey (BOSS; SDSS Collaboration et al. 2000; SDSS Collaboration et al. 2011, 2013) to place constraints on $\sum m_\nu$. The void statistics we consider are the void size function (VSF) and the void-galaxy cross power spectrum. We combine these with the usual galaxy auto power spectrum multipoles which by themselves already place a tight upper bound on $\sum m_\nu$ (through the suppression of matter power below the neutrino free-streaming scale; e.g., Ivanov et al. 2020a; Semenaite et al. 2023).

Since voids can be considered antihalos (Stopyra et al. 2021), a popular model for the VSF descends from Press-Schechter theory and the excursion set formalism (Press & Schechter 1974; Bond et al. 1991; Sheth & Tormen 1999), with slight modifications (Sheth & van de Weygaert 2004; Paranjape et al. 2012, 2013; Jennings et al. 2013; Pisani et al. 2015a; Ronconi & Marulli 2017; Contarini et al. 2019, 2022; Ronconi et al. 2019; Verza et al. 2019; Pelliciari et al. 2023). While the void-galaxy correlation function can be used for cosmological purposes without explicit knowledge of the void profile through the Alcock-Paczynski test and redshift space distortions (Alcock & Paczynski 1979; Lavaux & Wandelt 2012; Pisani et al. 2014), modeling of the profile itself has also been considered (Padilla et al. 2005; Paz et al. 2013; Hamaus et al. 2014a, 2014b; Massara & Sheth 2018; Kreisch et al. 2022).

However, all analytic approaches to modeling void statistics are problematic for our purposes. First, it is difficult to construct a consistent galaxy bias model across the different statistics comprising the data vector. Second, the calibration of analytic models typically did not utilize large simulations with varied neutrino mass. Third, existing models typically apply only to an aggressively cleaned subset of the entire void catalog, potentially leading to appreciable losses in constraining power.

Therefore, we choose to work in a simulation-based framework. Our simulations are based on particle neutrino, approximate gravity-only FastPM (Feng et al. 2016; Bayer et al. 2021a) realizations, in which we place galaxies through an expressive halo occupation distribution (HOD) model (Berlind & Weinberg 2002; Cooray & Sheth 2002; Wechsler & Tinker 2018). We then postprocess the galaxy catalogs to generate light cones incorporating survey realism.

The likelihood analysis with these simulations is a nontrivial problem. A popular approach is to build emulators of the mean data vector and perform the analysis under the assumption of a usually Gaussian likelihood where the covariance matrix is estimated from simulations. However, this approach turns out to be challenging for our problem. First, constructing an emulator in a 17-dimensional space (six cosmology, 11 HOD) is quite difficult, especially given a feasible number of simulations. Second, the assumption of a Gaussian likelihood is wrong. We demonstrate in Section 4.5 that the VSF is very

close to Poisson distributed (as long as bins are chosen wide enough, as would be naively expected), but modeling its covariance with the void-galaxy cross power spectrum and the galaxy auto power spectrum is difficult.

For these reasons, we opt for an implicit-likelihood<sup>11</sup> approach (Cranmer et al. 2020). This formalism uses neural networks to approximate functions that can be converted into posteriors. In general structure, this work is therefore similar to the SIMBIG papers (Hahn et al. 2022, 2023), but it differs in almost all details (statistics, simulations, objective, HOD, and code). The resulting complementarity will therefore be useful to assess the state of implicit likelihood inference (ILI) in galaxy clustering cosmology.

The rest of this paper is structured as follows. Section 2 describes our simulation pipeline. Section 3 contains details on the data vector and the inference procedure. Section 4 collects our results and their interpretation. We conclude in Section 5. The appendices contain additional material as well as information about the data and code availability.

## 2. Simulations

### 2.1. Cosmological Prior

Since our objective is $\sum m_\nu$, we place a tight prior on Lambda cold dark matter (ΛCDM). For this, we use the posterior from the Planck (Planck Collaboration et al. 2020) primary CMB analysis.<sup>12</sup> Specifically, we use the chains run with fixed $\sum m_\nu$ and measure the mean and covariance matrix in the five “CMB parameters,” $\omega_c$, $\omega_b$, $\log A_s$, $n_s$, and $\theta_{\text{MC}}$. In these parameters the posterior is close to Gaussian and we approximate it as such. To ensure the robustness of our conclusions, we inflate the Planck error bars on the cosmological parameters by a factor of 2. For $\sum m_\nu$, we choose a flat prior between 0 and 0.6 eV, the upper boundary being motivated by preliminary tests (see Section 2.3) in which we established a sensitivity of the order $\sigma \sim 0.2$ eV. We assume three neutrino species with degenerate masses.

Of course, the primary CMB's information leads to some correlation between $\sum m_\nu$ and the CMB parameters. This correlation is not included in our prior. However, given the sensitivity of the data used (compared to Planck), these residual correlations have a relatively small effect. For example, projecting the $\sum m_\nu-\omega_c$ correlation in the Planck posterior to our upper prior boundary of 0.6 eV, we obtain a shift $\Delta\omega_c/\sigma_{\text{EFT}} \sim 0.25$ where $\sigma_{\text{EFT}}$ is the error bar obtained from the EFTofLSS analysis of BOSS (Philcox & Ivanov 2022; of which we only use a subset). From these considerations, it also follows that our results do not depend strongly on the precise choice of ΛCDM prior.

We draw from the cosmological prior using an open quasi-random sequence. In contrast to popular sampling methods such as Latin hypercube or Sobol, an open sequence does not require initial knowledge of the total number of samples. Our sequence is constructed by taking integer multiples of a vector whose elements are powers of a generalized golden ratio.<sup>13</sup> In hindsight it turns out that our results are insensitive to taking out random subsets from the sequence (see Appendix F), indicating that a simple pseudorandom sampling would have

<footnote n="11"><sup>11</sup> Likelihood free and simulation based.</footnote>

<footnote n="12"><sup>12</sup> plikHM-TTIEEE-low1-lowE</footnote>

<footnote n="13"><sup>13</sup> https://tinyurl.com/4726n999</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** The cosmological parameter values sampled. We adopt a prior on ΛCDM that has the same shape and twice the size (in σ) as the Planck posterior. We sample this prior with an open quasi-random sequence designed to yield low-variance estimates of integrals. The color scale correlates with $\sum m_\nu$. The red marker indicates the position of our fiducial simulations.</figcaption></figure>

been sufficient. This is likely due to the aforementioned compactness of the prior compared to the data's sensitivity.

## 2.2. Cosmological Simulations

We run 127 simulations with varied cosmologies and 69 at a fiducial cosmology, illustrated in Figure 1.<fnref n="14" /> We choose the fiducial cosmology close to the mean of the ΛCDM prior, with $\sum m_\nu = 0.1$ eV. The cosmo-varied simulations share the random seed. This is because the decision to adopt an implicit-likelihood approach was only made after we encountered severe challenges with the conventional approach, as described in the introduction. However, as we shall see below, the simulations are large enough that sufficient quasi-independent data vectors can be generated.

After generating a cosmological parameter vector in the CMB parameters we replace $\theta_{\text{MC}}$ by the Hubble constant using CAMB (Lewis et al. 2000; Howlett et al. 2012). We then produce power spectra at $z=99$ using CLASS (Blas et al. 2011) and REPS (Zennaro et al. 2016, 2017).

We run particle neutrino gravity-only simulations using the approximate FastPM solver. We choose a box size of $2.5 h^{-1}$ Gpc with $2800^3$ CDM particles, leading to a minimum resolved halo mass $\sim 1.3 \times 10^{12} h^{-1} M_{\odot}$ that is sufficient for the CMASS luminous red galaxy sample. With regard to the neutrino options in FastPM we follow Bayer et al. (2021a).<fnref n="15" /> In particular, we follow their recommendation to increase the number of early time steps for larger neutrino masses. Specifically, at $\sum m_\nu = 0$ we take seven logarithmic steps (in scale factor) between $z=99$ and $z=19$, followed by 12 linear steps until $z = 0.68$. Afterwards, we take 20 steps until $z = 0.44$.

during each of which we write a snapshot of CDM particles to disk. At higher neutrino masses we insert up to 12 additional logarithmic steps before $z=79$. The seemingly large number of 20 snapshots was established in our preliminary tests during which we saw slight differences between 10 and 20 snapshots (in the summary statistics considered and for a single FastPM run) and thus decided to err on the side of caution. Note that a recent study supports the accuracy of FastPM for a task similar to ours, albeit with some important differences (Variu et al. 2024).

## 2.3. Galaxies

We fine tuned the method to populate CDM snapshots with galaxies through preliminary tests. Specifically, we performed global optimization over a large and partly discrete space of HODs considering two objectives: (1) the power spectrum multipoles of galaxies placed in Quijote simulations (Villaescusa-Navarro et al. 2020) and (2) the VSF of the CMASS data. The first step was intended to identify degrees of freedom that are necessary to correct for potential approximation errors in FastPM, while the second test was primarily meant to test our simulations' fidelity. Discrete degrees of freedom included in the optimization are the choice of halo catalog, the assembly bias variable, whether and where to include redshift dependence, and whether to include velocity biases. The continuous degrees of freedom were used to establish intuition for reasonable prior ranges. During optimization, we try to fit the summary statistics described above by varying the HOD parameters and use the best-fit cost (as well as the cost profile) to inform subsequent analysis choices. The optimization problems were solved with optuna (Akiba et al. 2019). In the following, we briefly describe the HOD model resulting from these preliminary tests. More detail can be found in Appendix A.

We identify halos in the CDM snapshots using Rockstar (Behroozi et al. 2012, 2013), which in our preliminary tests performed better than the friends-of-friends finder shipped with FastPM.

Then, galaxies are assigned stochastically to halos using an HOD. Besides the usual five-dimensional model parameterized by $M_{\text{min}}$, $\sigma_{\text{log}M}$, $M_0$, $M_1$, and $\alpha$ (Zheng et al. 2005), we introduce six additional degrees of freedom.

Although galaxy assembly bias (Zhu et al. 2006; Pujol & Gaztañaga 2014; Zentner et al. 2014) has been argued to be not necessary to describe the clustering of CMASS galaxies (Reid et al. 2014; Lin et al. 2016; Kobayashi et al. 2022), we decide to be conservative by adding galaxy assembly bias parameterized by $P_1$ and $a_{\text{bias}}$. Furthermore, we add velocity bias (Berlind et al. 2003; Yoshikawa et al. 2003; van den Bosch et al. 2005; Guo et al. 2015) parameterized by $\eta_{\text{cen}}$ and $\eta_{\text{sat}}$. Finally, we introduce redshift dependence to $M_{\text{min}}$ and $M_1$, parameterized by $\mu(M_{\text{min}})$ and $\mu(M_1)$, respectively. One advantage of having these slopes as free parameters is that we know them to be relatively close to zero, enabling useful sanity checks on any posteriors.

The resulting 11-dimensional HOD parameterization is quite similar to, e.g., Hahn et al. (2022, 2023) and Zhai et al. (2023).

We populate the cosmo-varied simulations with galaxies according to the HOD parameters drawn from the priors given in Appendix A. For the simulations at the fiducial cosmology we only populate with a single HOD. We choose the HOD parameters used for the fiducial mocks based on preliminary

<footnote n="14">¹⁴ We attempted 130 and 71, respectively, but some runs failed. As discussed before, the fact that due to the failures we do not sample the quasi-random sequence strictly sequentially does not affect our results.</footnote>

<footnote n="15">¹⁵ every_ncdm = 4, n_side = 4, n_shell = 10.</footnote>
inference runs using the VSF only. It turns out that these
parameters are not very close to best fit when considering the
entire data vector; generating new mocks closer to the best-fit
point may increase the efficiency of the compression step
described below. However, the difference in HOD parameters
cannot cause biases since the fiducial mocks are not used in
constructing the likelihood.

## 2.4. Light Cones

We use the cuboid remapping code BoxRemap<fnref n="16" /> (Carlson & White 2010) to deform our simulated cubes to the CMASS Northern Galactic Cap (NGC) geometry. It turns out that there are two possible choices of remapping and we use both (as part of the augmentation scheme discussed below).<fnref n="17" />

When projecting galaxies onto the light cone, we extrapolate
their positions from the snapshots using the host halo velocities
(using the stochastic galaxy velocities would weaken the
correlation between centrals and satellites). The resulting
corrections are small thanks to the large number of available
snapshots.

After mapping galaxies to the light cone, we apply all
angular masks and approximately mimic fiber collisions using
the procedure described in Hahn et al. (2022).

In contrast to some other works, we downsample the galaxy field predicted by the HOD to the data's density $n(z)$ (only, of course, if the simulation contains more galaxies at the given redshift). This downsampling is performed iteratively in conjunction with the implementation of fiber collisions so that both are self-consistent. Our implementation performs any necessary downsampling regardless of host halo properties; future work could take a prediction for stellar mass into account.

# 3. Inference

## 3.1. Data Vectors

We use the NGC part of the CMASS sample. The Southern Galactic Cap (SGC) is smaller and we do not expect dramatic improvements from its inclusion. Since our focus is on better understanding the impact of void statistics, rather than the tightest possible bounds on neutrino mass, we ignore the SGC for simplicity. Similarly, we do not include the LOWZ sample; its smaller volume makes it less suitable for void science.

We cut galaxies into the redshift interval $0.42 < z < 0.7$ and map them to comoving space using a fixed $\Omega_m = 0.3439$. Voids are identified using the VIDE code (Sutter et al. 2015) which is based on ZOBOV (Neyrinck 2008) and works by Voronoi tessellating the galaxies and then applying a watershed algorithm to find contiguous density minima. We use the “untrimmed” catalog computed by VIDE as it does not require arbitrary assumptions.

While many different void finders exist (e.g., Colberg et al. 2008; van de Weygaert & Schaap 2009; Cautun et al. 2018), prior work suggests that shape-agnostic void finders such as VIDE yield voids with better constraining power on $\sum m_\nu$ than spherical finders (Kreisch et al. 2022). Future work could

investigate the influence of void definition on signal-to-noise
ratio.

Galaxy auto power spectra $P_l^{gg}(k)$ and void-galaxy cross power spectra $P_l^{vg}(k)$ are computed using `nbodykit` (Hand et al. 2018, 2019) and `pypower`,<fnref n="18" /> reducing variance with Feldman-Kaiser-Peacock weights (Feldman et al. 1994) and correcting for observational systematics using the provided weights (except for fiber collisions, of course; Reid et al. 2016). We only utilize the systematics weights when computing $P^{gg}$. In the case of $P^{vg}$, we find no significant change in posteriors when using the galaxy weights, consistent with Hamaus et al. (2017). In the case of void identification, there is no guarantee that the obvious method to incorporate the systematics weights would yield cleaner voids. Given the relatively large void sizes considered in this work, we do not expect significant contamination by unmodeled survey systematics, but suggest that this point may warrant future work. Galaxy randoms are taken from the public catalogs. Void randoms are constructed by taking a large catalog of voids from many different mocks to choose angular positions and constructing a kernel density estimator in redshift matched to the specific void catalog. This procedure ensures that the randoms are consistent with a given cut in void radius since voids of different sizes have somewhat different angular distributions due to the survey mask.

We consider the VSF in 32 linearly spaced effective radius bins between 30 and $80 h^{-1}$ Mpc. The minimum radius cut is well above the mean tracer separation and thus we expect contamination by Poisson voids to be small (Pisani et al. 2015b; Cousinou et al. 2019). We split voids for the VSF into two redshift bins, separated at $z = 0.53$. This splits the CMASS sample approximately equally.

We perform analyses with $k_{max} = 0.15$ and $0.2$.<fnref n="19" /> We consider $k_{max} = 0.15$ the conservative baseline choice but $k_{max} = 0.2$ is still expected to be reliably modeled by FastPM as well as the halo model (e.g., Kobayashi et al. 2022). We do not use power spectra on scales larger than $k_{min} = 0.01$ (Ivanov et al. 2020b). Since our theoretical model is simulation based, we do not deconvolve the survey window function. This means that there is a small level of contamination by Fourier modes outside the $k$-range considered, but we assess this effect to be negligible.

In our baseline analysis we only use the monopole $P_0^{gg}$ of the galaxy auto power spectrum multipoles. This choice was made based on the limited information content of the quadrupole (from EFTofLSS posteriors) and with the aim of simplicity. We discuss the effect of including the quadrupole below. For $P^{vg}$ we use both the monopole and the quadrupole.

It is worth noting that we opt to use reciprocal space void-galaxy cross power spectra $P^{vg}$ instead of the more popular configuration space correlation function. The correlation function has the primary advantage that one can rescale its argument on a void-by-void basis by the respective radius and thus sharpen the resulting void profile (this is also possible in reciprocal space but computationally expensive; future work could explore this point). We believe, however, that the mixing of Fourier modes in the correlation function could lead to problems with approximate solvers like FastPM whose domain of validity is better localized in reciprocal space. In order to optimize the signal-to-noise ratio, we consider $P^{vg}$

<footnote n="16">¹⁶ http://mwhite.berkeley.edu/BoxRemap</footnote>

<footnote n="17">¹⁷ The two cuboids have side lengths {1.412, 1, 0.7071} and {1.7321, 0.8165, 0.7071} (in units such that $2.5 h^{-1}$ Gpc ≡ 1). They are generated by the matrices {{1, 1, 0}, {0, 0, 1}, {1, 0, 0}} and {{1, 1, 1}, {1, 0, 0}, {0, 1, 0}}, respectively.</footnote>

<footnote n="18">¹⁸ https://pypower.readthedocs.io</footnote>

<footnote n="19">¹⁹ For brevity, we implicitly take all wavenumbers in units of $h$ Mpc$^{-1}$.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Illustration of the data vector considered. In the top panel, the vertical axis has been scaled such that the different data vector components are well visible and the power spectra are plotted as kP(k). The solid lines display best-fit mocks, considering separate parts of the data vector (the notation $P_{0,2}^{gg}$ denotes the combination of monopole and quadrupole of the galaxy auto power spectrum, etc.). The $\chi^2$ was computed under the approximation of a Gaussian likelihood and the mocks were averaged over eight augmentations (no interpolation/emulation was performed). The fact that, e.g., the model with the best-fit VSF still reproduces the other parts of the data vector reasonably well, as well as the reduced $\chi^2$ values close to unity, indicate our mocks' high fidelity. Note that we do not use $k < 0.01$ and $k > 0.15$ in our baseline analysis and these scales were not included in identifying the best-fit models shown here. The bottom panel shows the residuals in units of estimated standard deviations. For better visibility, the y-axis has linear and logarithmic scales with the transition indicated by the dotted lines.</figcaption></figure>

computed with three different choices of minimum void radius, 30, 40, and 50 $h^{-1}$ Mpc.

An illustration of the data vector is given in Figure 2. We emphasize that the “best-fit” curves shown there are simply picked from the available simulations, implying that the $\chi^2$ values and residuals cannot be interpreted in the usual way.

## 3.2. Implicit Likelihood Inference

As discussed in the Introduction, the standard emulator-based approach is difficult in 17 dimensions. The main reason is that the training objective for an emulator does not directly map to the ultimate goal of accurate posteriors, implying that the optimum needs to be very sharp (which requires many simulations).

Combined with the unknown likelihood function, we believe ILI to be the appropriate tool for our task.

We opt for neural ratio estimation (NRE; Cranmer et al. 2015; Hermans et al. 2019; Durkan et al. 2020; Delaunoy et al. 2022; Miller et al. 2022), which recasts inference as a classification problem. The choice of an amortized instead of a sequential method<fnref n="20" /> was made based on the hierarchical structure of our simulations; we then opt for NRE because of its simplicity. In its original and simplest form, NRE works with pairs of parameter vectors $\theta$, $\theta'$ drawn from the prior $p(\theta)$. We then consider a data vector $x$ drawn from the likelihood $p(x|\theta)$, where the simulation process described above approximates this draw. A neural network $f$ maps the pairs $(x, \theta)$, $(x, \theta')$ to scalars $y, y'$. If we now choose a classification loss function $L(y, y')$, e.g., binary cross entropy

$$L(y, y') = -\log(y) - \log(1 - y'), \quad (1)$$

it is easy to show that the functional optimization problem

$$f^* = \underset{f}{\operatorname{argmin}} \int d\theta d\theta' dx p(\theta)p(\theta')p(x|\theta)L, \quad (2)$$

has the solution

$$\frac{p(x|\theta)}{p(x)} = \frac{f^*(x, \theta)}{1 - f^*(x, \theta)}. \quad (3)$$

In other words, a neural network trained to distinguish between samples from $p(x, \theta)$ and samples from $p(x)p(\theta)$ approximates the likelihood-to-evidence ratio at optimum. Posteriors can then be obtained through usual Monte Carlo Markov Chain sampling, which we perform with emcee (Foreman-Mackey et al. 2013a, 2013b). In practice, this general idea of approximating $p(x|\theta)/p(x)$ through a classifier works better in the multiclass version “NRE-B” (Durkan et al. 2020). We use the implementation provided in the sbi package (Tejero-Cantero et al. 2020).

In the above, it is actually not necessary for the parameter vectors $\theta$ to be drawn independently from the prior $p(\theta)$. In fact, all that is required is that a sum over the simulated parameter vectors approximates the integral in Equation (2). For this reason, it is correct for us to populate each of the 127 cosm varied simulations with multiple draws from the HOD prior (~230). For each HOD draw we compute eight augmentations as described below, yielding ~1.7 × 10⁵ training samples.

The ILI framework allows implicit marginalization over nuisance parameters. This is one of its primary benefits in high-dimensional parameter spaces.<fnref n="21" /> In principle, we could take $\theta = \{\sum m_i\}$ as one dimensional. In practice, it is likely better to include a subset of the nuisance parameters in $\theta$. This is because we have intuition for the posteriors expected for some nuisance parameters and thus making them explicit allows

<footnote n="20">²⁰ Sequential ILI considers the observed data to sequentially narrow down the space in which to run simulations, thus being more efficient usually. On the other hand, amortized ILI does not specialize on the data at hand; this enables important validations like coverage tests (see Figure 6).</footnote>

<footnote n="21">²¹ Consider computation of a high-dimensional integral over $f(x)$ given samples $f(x_i)$. Interpolating $f(x)$ using these samples and then performing quadrature is a more difficult problem than using the Monte Carlo estimator.</footnote>
useful checks. We opt to include log $M_{\min}$ and $\mu(M_{\min})$ in the
parameter vector. For the former we know that the data should
provide a constraint considerably tighter than the prior, while
for the latter we expect a result close to zero. The posterior on
$\sum m_\nu$ is unaffected by this choice of $\theta$ and the extra
computational cost in training and sampling the neural network
is marginal when making two nuisance parameters explicit
(however, making all nuisance parameters explicit would
complicate the training unnecessarily).

We parameterize the classifier $f$ as a residual neural network. Hyperparameter optimization was performed considering the loss on a validation set of 13 cosmologies (i.e., ~2.4 × 10⁴ mocks).<fnref n="22" /> We converged at a relatively large network with 1.7 × 10⁷ trainable parameters but high dropout rates.

High-dimensional data vectors $x$ are often problematic for ILI, our problem being no exception. This necessitates a compression step before the data vector is passed to the neural network. Since we expect our likelihood to be close to Gaussian/Poissonian, we use the linear score compression MOPED (Heavens et al. 2000) to 17 compressed statistics. Indeed, MOPED is locally optimal both for a pure Gaussian and a pure Poissonian likelihood. We also experimented with the nuisance-hardened generalization (Alsing & Wandelt 2019) to 5 and 10 dimensions, obtaining consistent but slightly wider posteriors. In order to construct the MOPED compression matrix, estimates for the covariance and derivatives of the data vector are required. We construct the covariance matrix from our fiducial mocks using the usual estimator. For the derivatives, we generate ~10⁵ additional mocks (10³ parameter vectors, each with 96 augmentations) in a small ball around the fiducial model. We then perform linear regression and read off the derivatives. Simple tests indicate that the dependence on parameters is close to linear in the region considered.

### 3.3. Augmentations

As discussed before, our cosmo-varied simulations share the random seed. This fact ostensibly makes them unsuitable for the ILI approach discussed before since the integral in Equation (2) requires a sampling of initial conditions.

However, as we shall discuss in this section, it is possible to generate many quasi-independent realizations from a single simulation. As mentioned before, we do not require independent identically distributed realizations, so this is in fact enough to approximate Equation (2) with sufficient accuracy.

For a single 2.5 $h^{-1}$ Gpc simulation box populated with galaxies, we can take the product of the following transformations: two cuboid remappings, eight reflections, and six axis transpositions. This results in 96 augmentations. In principle, many more augmentations can be generated through translations, but we expect these to be more correlated.

The crucial question now is whether these 96 augmentations approximate the distribution over initial conditions. We can answer this question by considering our fiducial simulations, which have 69 different random seeds. Given the fiducial parameter vector, we generate a matrix $D_{\mu a}$ whose elements are data vector valued and where $\mu = 1...69$, $a = 1...96$ index the initial conditions and augmentations, respectively. We can perform statistical tests by computing marginals over $\mu$ and $a$

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Checks for statistical independence of augmentations, in the compressed space. The data labeled “augments” are obtained by marginalizing over augmentations, while those labeled “ICs” are obtained by marginalizing over the initial conditions. Top: ratio of average standard deviations, with the difference with unity not exceeding two percent. This check can also be performed with the uncompressed data vector, leading to similar conclusions. Bottom: distributions of log-likelihoods of covariance matrices under the Wishart distribution with our fiducial covariance matrix. The distributions have large overlap, again indicating that the augmentations are very close to statistical independence.</figcaption></figure>

individually or jointly. In order to simplify the statistical interpretation, we restrict $a$ to 69 randomly chosen indices.

In the upper panel of Figure 3, we compare the diagonals of covariance matrices in the MOPED compressed space. We see that the standard deviations are almost identical for marginalization over $\mu$ and $a$. This test can also be performed for the uncompressed data vector, yielding consistent results and no systematic differences between different summary statistics or scales.

In the lower panel, we perform a test considering the entire content of the covariance matrices. We construct the Wishart distribution given the covariance matrix $C_{\text{joint}}$ obtained by marginalizing over $\mu$ and $a$ jointly and then compute the log-likelihood of the individually marginalized covariance matrices. If these covariance matrices were drawn from the Wishart distribution sourced by $C_{\text{joint}}$, their log-likelihoods would be distributed as indicated by the green line. We see that the distributions are somewhat different but still have large overlap. In conclusion, the 96 augmentations reproduce the distribution over initial conditions reasonably well. Since 96 ≫ 1, we expect the augmentations to provide a good approximation to the integral in Equation (2).

Why does this augmentation procedure work? First, our simulation boxes are about 5.7× larger than the survey volume. Second, the augmentations alter the redshift direction. Third, galaxies and the survey mask interact. Fourth, galaxies are captured at different times so their peculiar motions alter their real space positions. All these points need to be seen relative to the specific survey and simulation configuration; the described augmentation procedure is certainly not expected to work universally.

## 4. Results

In this section, we first present our main posteriors on $\sum m_\nu$ from the CMASS NGC data, taking various combinations of

<footnote n="22">²² It is actually important to separate training and validation sets by cosmologies. Initial trials that mixed simulations exhibited hidden overfitting.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cumulative posteriors on $\sum m_\nu$ from different data vector combinations, with our baseline choice of $k_{\text{max}} = 0.15$. The addition of void–galaxy cross power spectra improves the constraint compared to the galaxy auto power spectrum, whereas the VSF has negligible impact.</figcaption></figure>

the summary statistics VSF $N_\nu$, void–galaxy cross power spectrum $P^{\text{vg}}$, and galaxy auto power spectrum $P^{\text{gg}}$.

We present most of our posteriors in their cumulative form. This is because at the current level of precision, no neutrino mass detection is expected and upper bounds are the main objective. The cumulative posterior is the most direct visualization of upper bounds. In all plots we include a diagonal dashed line indicating the prior.

In the following, we will occasionally compare with results obtained with the EFTofLSS (Ivanov et al. 2020b; d’Amico et al. 2020; Chen et al. 2022). The EFTofLSS allows for the analysis of the full-shape galaxy auto power spectrum (as well as other statistics we will not consider here). We use the window-less full-shape likelihood (Philcox et al. 2020; Philcox 2021; Philcox & Ivanov 2022)<fnref n="23" /> and the CLASS-PT code (Chudaykin et al. 2020).<fnref n="24" /> We restrict the data included in the likelihood to the NGC high-$z$ sample, approximately equal to the data we use for our analysis. Furthermore, we impose the same $\Lambda$CDM prior while keeping the nuisance parameter priors equal to those implemented in the public likelihood code. In any comparison with our results we use identical $k_{\text{max}}$. Likewise, we usually only use the monopole $P_0^{\text{gg}}$, consistent with our simulation-based analysis. The likelihood part termed “Alcock–Paczynski” in the EFTofLSS likelihood is included, since our method also effectively includes this term. On the other hand, we do not include the baryon acoustic oscillations (BAO) reconstruction or real space likelihoods.

We emphasize that a comparison between the EFTofLSS and HOD methods is beyond the scope of this work. Therefore, we will use the EFTofLSS posteriors to provide intuition, show that our posteriors are at least qualitatively reasonable, and for an interesting observation about the quadrupole $P_2^{\text{gg}}$ later on.

### 4.1. Neutrino Mass Posterior

In Figure 4, we show the baseline posterior on $\sum m_\nu$, with $k_{\text{max}} = 0.15$. The galaxy auto power spectrum gives a 95% credible interval constraint of $\sum m_\nu < 0.43$ eV. Upon inclusion

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Effect of increasing $k_{\text{max}}$ to 0.2. The posterior shrinks, as expected, but remains broader than the EFTofLSS result. The impact of adding voids to the data vector is now reversed.</figcaption></figure>

of the VSF, the posterior broadens somewhat. Including the void–galaxy cross power spectrum tightens the posterior to $\sum m_\nu < 0.35$ eV, an ~20% improvement. Further adding the VSF does not lead to any appreciable change. Posteriors are generally wider than the EFTofLSS result. We suspect that the dominant effect driving the differences between the EFTofLSS posteriors and our simulation-based posteriors stems from the choice of priors. The HOD prior used in this work likely contains some redundancies, especially since we downsample to the observed $n(z)$. Prior volume effects also exist in the EFTofLSS (e.g., Simon et al. 2023), but are likely subdominant in this case.

In Figure 5, we show a similar set of posteriors obtained with $k_{\text{max}} = 0.2$. We believe that our simulated model should still have a high level of fidelity at these somewhat smaller scales. We observe that including smaller scales tightens the posterior, as expected. However, adding void statistics to $P^{\text{gg}}$ now slightly broadens the posterior. Most of the remainder of this section will be devoted to better understanding the observations from Figures 4 and 5.

### 4.2. Validation

Any simulation-based, and especially implicit-likelihood, inference necessitates rigorous validation of the simulated model, the likelihood approximation, and the resulting posteriors. In this section, we present three tests verifying different aspects of our pipeline.

First, in Figure 6, we present the usual coverage (or $q-q$) plot (Hermans et al. 2021). For this diagnostic, we perform inference on mocks drawn from the prior; in particular, we use the ~2.4 × 10⁴ validation mocks discussed before. We use the $N_\nu + P^{\text{vg}} + P^{\text{gg}}$ likelihood with $k_{\text{max}} = 0.15$. For each resulting chain, we compute the marginal distributions of the explicit parameters and then the confidence level at which the true input parameter is located. In Figure 6, we show the cumulative histograms of these confidence levels. If the posterior is well-calibrated, these cumulative distribution functions (CDFs) should coincide with the diagonal. As can be seen, for all parameters considered this is the case. This

<footnote n="23">²³ https://github.com/oliverphilcox/full_shape_likelihoods</footnote>

<footnote n="24">²⁴ https://github.com/michalychforever/CLASS-PT</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Coverage (q-q) plot, demonstrating that when performing inference on mocks the posteriors are well calibrated.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Effect of adding the quadrupole $P_2^{gg}$ to the data vector. Consistent with the EFTofLSS result, the $\sum m_\nu$ posterior broadens slightly. This effect is likely a statistical fluctuation. However, it provides a useful cross-check for how well our simulations model redshift space distortions.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 8.** Comparison of the data posterior with posteriors obtained on fiducial mocks, using the $N_\nu + P_{0,2}^{vg} + P_0^{gg}$ data vector. Note that the HOD used in generating these mocks is somewhat different from the best fit. While the data posterior at $k_{\text{max}} = 0.15$ is slightly unusual, at $k_{\text{max}} = 0.20$ it becomes more typical (compared to the distribution of mock posteriors).</figcaption></figure>

diagnostic is a powerful internal consistency check and verifies that the neural network is well trained.

Second, in Figure 7, we show an interesting observation concerning the galaxy auto power spectrum quadrupole $P_2^{gg}$. As discussed before, this summary statistic has limited constraining power and we do not use it for our main posteriors. As can be seen in Figure 7, adding the quadrupole to the data vector slightly broadens the posterior. This happens consistently in our analysis and in the EFTofLSS. We believe that this observation increases confidence in the validity of our simulation model, in particular the modeling of redshift space distortions.

Third, in Figure 8 we compare the data posteriors with posteriors obtained by running inference on randomly chosen

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Effect of considering voids above different radius cuts. Voids with radii between 40 and $50 \ h^{-1} \ \text{Mpc}$ contribute the majority of the observed tightening of the posterior relative to the $P^{gg}$-only result.</figcaption></figure>

mocks generated at the fiducial point. We remind the reader that the fiducial HOD is rather far from best fit, which somewhat complicates the interpretation. We observe that at $k_{\text{max}} = 0.15$ the data posterior is tighter than most of the mock ones. If the cosmological simulations were to blame for this, the naive expectation would be for the discrepancy to become more severe as smaller scales are included in the analysis. However, this appears not to be the case: at $k_{\text{max}} = 0.20$ the data posterior becomes more typical. We conclude that even though we observe hints of differences between the data and simulations, the evidence is not conclusive and the data posterior could well be consistent with the observed distribution. It should also be noted that the real $\sum m_\nu$ may be less than the choice 0.1 eV with which the fiducials were run, potentially leading to a tighter data posterior.

### 4.3. Broadening of Posteriors

One peculiar observation is that the inclusion of void statistics can broaden the posterior on $\sum m_\nu$. We do not fully understand this phenomenon and can only provide some suggestive results. These are more comprehensively described in Appendix B; here we only provide a summary.

We observe similar broadening on fiducial mocks and thus propose that we are in fact observing a generic phenomenon. Therefore, we suggest the possibility that void statistics are most effective at constraining the neutrino mass sum from below. A further test using artificially enlarged volumes supports this theory.

For a potential physical explanation, we consider the free-streaming length. At $z = 0.5$, $\lambda_{\text{fs}} = 90 \ h^{-1} \ \text{Mpc}$ ($0.3 \ \text{eV}/\sum m_\nu$) for degenerate masses. This length scale is comparable to the diameter of the voids that seem to contribute most (see Section 4.4). Thus, it may be that $\sum m_\nu$ at the upper end of the posterior is “invisible” to voids. However, we identify voids using tracers of small-scale fluctuations, so the full picture is much more complicated and could be a subject for further study.

### 4.4. Void Radius

In Figure 9, we show posteriors obtained with the $N_\nu + P^{vg} + P^{gg}$ likelihood, concentrating on void size. Cuts on effective radius are performed both in the VSF and $P^{vg}$ parts of the data vector. We observe that the posteriors are almost identical regardless of whether we cut at 30 (the baseline analysis) or $40 \ h^{-1} \ \text{Mpc}$. On the other hand, further increasing the minimum radius to $50 \ h^{-1} \ \text{Mpc}$ removes much of the effect of voids on the posterior. Figure 9 indicates that at least for the
<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Different void-only statistics combinations. Upper panel: while all posteriors are close to the prior, it appears that $P^{vg}$ carries somewhat more signal than $N_v$, and the combined signal is strongest. Lower panel: consistent with Figure 9, it appears that larger voids carry more neutrino mass signal, although the posteriors are quite close to the prior.</figcaption></figure>

present analysis voids with effective radii between 40 and 50 $h^{-1}$ Mpc are the most constraining. Smaller voids might be contaminated by spurious Poisson voids and perhaps also due to their shallower density profile being less affected by neutrinos. Larger voids presumably suffer from their low abundance.

In Figure 10 we show posteriors obtained from void statistics only. We show them mostly for completeness; in the present analysis these are entirely prior dominated. However, even in this plot we see the previously mentioned observation that larger voids appear to carry more signal.

## 4.5. Poissonian Void Size Function

As a final point of this section, we substantiate the previous claim that the VSF is very close to Poissonian distributed. While this seems to be a natural assumption, void exclusion makes it nontrivial. Indeed, previous works have assumed Poisson likelihoods (Sahlén 2019; Contarini et al. 2023); our simulations enable us to check this assumption.

In Figure 11, we show two checks performed with our fiducial mocks. The left panel shows the covariance matrix divided by the outer square of the Poissonian standard deviations; the result is close to the identity. The right panel shows a check using the variance-stabilizing Anscombe transform (Anscombe 1948). For each mock data vector $c^{(\alpha)}$ and bin $i$, we compute the transformed VSF count

$$ \tilde{c}_i^{(\alpha)} = 2 \left( \sqrt{c_i^{(\alpha)} + \frac{3}{8}} - \sqrt{\langle c_i \rangle + \frac{3}{8}} \right) + \frac{1}{4\sqrt{\langle c_i \rangle}} \quad (4) $$

In the limit of large counts the distribution of these transformed counts converges to the standard normal if the counts themselves are Poissonian. As can be seen, the agreement

with the standard normal is quite good indeed. These tests demonstrate that deviations from Poissonian distribution are small for the VSF, at least for the choice of binning considered here.

## 5. Conclusions

We have performed inference on galaxy clustering in the BOSS CMASS NGC sample, combining the VSF, the void–galaxy cross power spectrum, and the galaxy auto power spectrum. Our primary target was the neutrino mass sum, $\sum m_\nu$; thus, we imposed a tight prior on $\Lambda$CDM informed by primary CMB data.

We argued that analytic models for the considered void statistics are not mature enough and unsuitable for our specific problem, necessitating a simulation-based approach. To this end, we ran approximate gravity-only simulations and populated them with galaxies using an expressive HOD. Several factors motivated the use of ILI.

In our baseline analysis, we find $\sum m_\nu < 0.43$ eV from the galaxy auto power spectrum alone, and $\sum m_\nu < 0.35$ eV with the void statistics included (95% credible interval). We performed several tests to confirm the statistical and systematic validity of our likelihood approximation.

We performed a short investigation of the impact of voids on the neutrino mass posterior. It appears that the void statistics may be most effective in constraining $\sum m_\nu$ from below. This result would imply that future analyses aiming at measuring $\sum m_\nu$ may benefit from including void statistics.

Our results suggest that larger voids with effective radii $> 40 h^{-1}$ Mpc carry most of the signal despite their lower abundance. This has interesting implications for future analyses, since voids of this size should be detectable in photometric catalogs with relatively low redshift errors, such as the one expected for Rubin/LSST (Ivezić et al. 2019). Of course, spectroscopic surveys such as DESI (Dey et al. 2019), Euclid (Laureijs et al. 2011), SPHEREX (Doré et al. 2014), Subaru Prime Focus Spectrograph (Takada et al. 2014), and Roman (Spergel et al. 2015) will continue to be cornerstones of void science. The trade-off between volume, galaxy number density, and redshift precision warrants further investigation.

We also demonstrate that the VSF is very close to Poisson distributed, a feature that had been assumed in previous analyses but never explicitly confirmed.

Future work could improve upon our analysis in multiple ways. First, the cosmo-varied simulations should be run with different random seeds (we decided for a fixed seed in anticipation of an emulator-based analysis, which ultimately turned out to be very difficult). Second, it may be beneficial to normalize the void–galaxy cross power spectrum by void number. Although in principle this would contain the same information as our data vector once the VSF is included, the necessary transformation is nonlinear and thus potentially invisible to our data compression. Third, the HOD modeling could be improved. Some of our priors may not be optimal, and our $n(z)$ downsampling is simplistic. The CMASS sample’s completeness is quite well known and could be used to put a prior on the downsampling. Fourth, it turns out that the cosmological simulations did not dominate the computational cost. It may therefore be economical to increase the accuracy in FastPM or switch to a different solver.

Our results point toward a complicated picture with regard to the relationship between massive neutrinos and voids. Future
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Checks of the Poissonian nature of the VSF. In the left panel, we show a rescaled covariance matrix obtained from our fiducial simulations of the VSF part of the data vector. For an exactly Poissonian distribution, this matrix would be the identity. In the right panel, we show the distribution of Anscombe-transformed VSF counts. Overplotted is a standard normal. Both tests indicate a distribution that is very close to Poissonian.</figcaption></figure>

data sets, both spectroscopic as well as photometric, promise to bring tight cosmological constraints from void science, since it scales well with number.

## Acknowledgments

We thank the anonymous referee for the constructive report. We thank Sofia Contarini, Adrian Bayer, Jia Liu, Jo Dunkley, and Masahiro Takada for useful discussions. We thank Oliver Philcox for explaining the EFTofLSS likelihood. The work of L.T. is supported by the NSF grant AST 2108078. The authors are pleased to acknowledge that the work reported on in this paper was substantially performed using the Princeton Research Computing resources at Princeton University, which is a consortium of groups led by the Princeton Institute for Computational Science and Engineering (PICSciE) and Office of Information Technology's Research Computing.

## Appendix A
### Halo Occupation Distribution

In this appendix, we provide a more detailed description of the adopted HOD model.

First, for reference, the baseline five-parameter model only depends on halo mass M and has mean occupations

$$ \bar{N}_{\text{cen}} = \frac{1}{2} \left[ 1 + \text{erf} \left( \frac{\log M - \log M_{\text{min}}}{\sigma_{\log M}} \right) \right], \quad (\text{A1}) $$

for the central galaxy and

$$ \bar{N}_{\text{sat}} = \bar{N}_{\text{cen}} \left( \frac{M - M_0}{M_1} \right)^{\alpha}, \quad (\text{A2}) $$

for the satellites. A central galaxy is placed with probability $\bar{N}_{\text{cen}}$ and the number of satellites is drawn from a Poisson distribution with mean $\bar{N}_{\text{sat}}$. The central galaxy is placed at the halo's center and assigned the halo velocity. The satellites are distributed isotropically with a Navarro–Frenk–White (NFW) profile (Navarro et al. 1997) and the concentration model of Duffy et al. (2008), using the analytic solution for the inverse NFW CDF from Robotham & Howlett (2018). Satellite velocities are drawn from a distribution assuming virialization.

On top of this baseline model, we implement galaxy assembly bias using the decorated HOD (Hearin et al. 2016a, 2016b) with the ratio of kinetic to potential energy $r \equiv T/U$ as proxy for assembly history. In our preliminary tests $T/U$ outperformed halo concentration, possibly due to limited resolution within the FastPM halos. The decoration works by splitting halos into two groups according to $r$. In order to reduce the effect of any potential evolution of $r$ with halo mass, we do this splitting separately within 64 groups containing equal numbers of halos. The fraction $P_1$ of halos with lowest $r$ is assigned type 1 (2) for positive (negative) $a_{\text{bias}}$, while the rest is assigned type 2 (1). Then, the mean occupations are modified as

$$ \Delta \bar{N}_{\text{cen}} = |a_{\text{bias}}| \min \left[ 1 - \bar{N}_{\text{cen}}, \frac{1 - P_1}{P_1} \bar{N}_{\text{cen}} \right], \quad (\text{A3}) $$

$$ \Delta \bar{N}_{\text{sat}} = |a_{\text{bias}}| \frac{1 - P_1}{P_1} \bar{N}_{\text{sat}}, \quad (\text{A4}) $$

for type 1 and

$$ \Delta \bar{N}_{\text{cen}} = |a_{\text{bias}}| \max \left[ -\bar{N}_{\text{cen}}, \frac{1 - P_1}{P_1} (\bar{N}_{\text{cen}} - 1) \right], \quad (\text{A5}) $$

$$ \Delta \bar{N}_{\text{sat}} = |a_{\text{bias}}|(-\bar{N}_{\text{sat}}), \quad (\text{A6}) $$

for type 2.

Velocity bias for the centrals is implemented by adding $\eta_{\text{cen}}V_{\text{vir}}n$ where $n \sim N(0, 1)$. For the satellites, the velocity difference from the host halo is scaled by $\eta_{\text{sat}}$.

In terms of the redshift dependence for $M_{\text{min}}$, $M_1$ is approximated as linear in scale factor, such that

$$ \Delta \log M_i = \mu(M_i)(a - a_0), \quad (\text{A7}) $$

with $a_0 = 1/(1 + 0.53)$.

We adopt flat priors $12.5 < \log M_{\text{min}} < 13.2$, $0.1 < \sigma_{\log M} < 0.8$, $12.5 < \log M_{0.1} < 15.5$, $0.2 < \alpha < 1.5$, $-3 < P'_1 < 3$, $-1 < a_{\text{bias}} < 1$, $5 < \eta'_{\text{cen}} < 10$, $-1 < \eta'_{\text{sat}} < 1$, $-20 < \mu(M_{\text{min}}) < 20$, and $-40 < \mu(M_1) < 40$. Here, all masses are in units of $h^{-1} M_\odot$, and the primed parameters are defined as $2P_1 = (1 + \tanh P'_1)$, $\eta_{\text{sat}} = \exp(\eta'_{\text{sat}})$, $\eta_{\text{cen}} = \exp(-10+\eta'_{\text{cen}})$. The above intervals were found during preliminary inference runs. Note in particular the small values of $M_{\text{min}}$, compared to other analyses. This is partly explained by systematically lower halo masses in FastPM, and partly by the $n(z)$ downsampling described below. We picked the transformations given by the primed parameters based on the intuition that strictly (mathematically) bounded intervals often indicate that a uniform prior in a transformed quantity is a better choice.

## Appendix B
### Broadening of Posteriors

We have seen in our main posteriors that adding void statistics to the data vector sometimes broadens the posterior on $\sum m_\nu$. In this appendix, we attempt to better understand this observation, focusing on $k_{\text{max}} = 0.20$. For this, we will rely on inference on our fiducial mocks.

The first possible explanation could be a statistical fluctuation, and we cannot definitely exclude this hypothesis. One way, however, to test it is to look at *average* posteriors on our fiducial mocks. We perform inference on ~20 randomly chosen mocks and plot the CDFs of concatenated chains in
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Averaged posteriors on ~20 randomly chosen fiducial mocks. A widening of the Σ*m*<sub>v</sub> posterior upon including void statistics is also observed in simulations.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Schematic illustration of the proposed mechanism explaining the broadening of posteriors when voids are included. This figure is only meant as a guide and not as a literal depiction of posteriors. Correlations and nuisance parameters complicate the picture.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 14.** Posteriors on joint analyses of four randomly chosen fiducial mocks, averaged over ~30 groups. The solid and dashed lines correspond to likelihoods with two different sets of five nuisance parameters kept explicit. We see that the posteriors where void statistics are included have a slightly more pronounced bump at the true value Σ*m*<sub>v</sub> = 0.1 eV, consistent with the speculative picture in Figure 13.</figcaption></figure>

**Figure 12.** There, we observe that the expected, average behavior is for the posterior to broaden once void statistics are added to the data vector.

The second possible explanation could be that once void statistics are added our compression procedure becomes less

efficient. This could certainly be the case if at linear order the void statistics appear more constraining than they are globally, thus *P*<sub>00</sub><sup>gg</sup> would be unnecessarily downweighted. This hypothesis appears unlikely in light of the full posteriors presented in Figures 15 and 16. In these posteriors, we observe that for the parameters that are actually constrained (like *M*<sub>min</sub>) adding void statistics generically tightens the posteriors. It appears unlikely that Σ*m*<sub>v</sub> should be an exception.

Having found these two hypotheses unsatisfactory, we arrive at the third one: void statistics tend to constrain Σ*m*<sub>v</sub> from below. We illustrate this theory qualitatively in Figure 13, which should not be interpreted as a literal depiction. In fact, in Section 4.4 we show that void statistics alone yield posteriors close to the prior. Figure 13 provides merely an effective depiction.

We can investigate this hypothesis further by performing the following test. In order to increase the signal-to-noise ratio, we perform inference on four fiducial mocks at the same time, shown in Figure 14. For this, we use a different set of neural nets in which we leave five nuisance parameters explicit. The reason is that all implicitly marginalized nuisance parameters are effectively assumed to be different for each of the four mocks, an effect we would like to minimize. Of course, increasing the number of explicit nuisance parameters complicates the training and we have less confidence in the precise calibration of the posteriors. For this reason, our baseline results were obtained with only two explicit nuisance parameters. For reference, the real-data posteriors obtained with these alternative neural nets are shown in Figure 16. We perform this test with two different sets of nuisance parameters kept explicit in order to gauge robustness (corresponding to the solid and dashed lines in Figure 14). Similar to Figure 12, we average posteriors over ~30 randomly chosen groups of four mocks in order to decrease sample variance. We observe that, consistent with our theory, the posteriors that include void statistics show a more pronounced hint of a bump at the true Σ*m*<sub>v</sub>. In principle, one could increase the simulated volume further by combining more mocks, but our neural nets are not calibrated at the required level of precision and thus the resulting posteriors would not be robust.

In summary, the more mundane ideas to explain the observed broadening of posteriors appear questionable given the tests presented. On the other hand, the idea that void statistics are most effective at constraining Σ*m*<sub>v</sub> from below receives support from our experiments. A more in-depth examination of this issue would constitute a great starting point for future work.

## Appendix C: Corner Plots

This appendix collects posteriors in the full parameter spaces considered. Figure 15 shows the baseline parameter space with two explicit nuisance parameters. Figure 16 shows larger sections of parameter space (it should be mentioned, however, that the corresponding neural networks were trained without further hyperparameter optimization, implying a somewhat lower level of confidence in the validity of these posteriors). Figure 17 shows our EFTofLSS posteriors, demonstrating that the ΛCDM part of the parameter space is prior dominated.
<figure><img src="image_17.png" /><figcaption>**Figure 15.** Posteriors in the full parameter space considered, including the two HOD parameters we choose to keep explicit. The HOD posteriors are consistent between the different analysis choices.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>**Figure 17.** For reference, we show the full EFTofLSS posteriors.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>**Figure 16.** Posteriors with different sets of nuisance parameters kept explicit. Since no hyperparameter optimization was performed when training the corresponding networks and the larger parameter space makes training more difficult, we assess these posteriors as less robust than our baseline results. The neutrino mass posteriors, however, are quite consistent. The left panel corresponds to the solid lines in Figure 12 and the right panel to the dashed lines.</figcaption></figure>

<figure><img src="image_20.png" /></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 18.** Impact of relaxing the prior on $A_s$ to an essentially uninformative level. The posterior on $\sum m_\nu$ remains essentially unchanged. This test is performed with the EFTofLSS likelihood since the implicit-likelihood setup would require significant computational expense for a post hoc change of the prior.</figcaption></figure>

## Appendix D
### Impact of $A_s$ Prior

Given concerns about the consistency between early and late-time determinations of the clustering amplitude, we test whether our results are particularly sensitive to the choice of prior on $A_s$. Given the significant cost of changing the prior in our simulation-based analysis, we perform this test using the EFTofLSS likelihood. Figure 18 shows the change in $\sum m_\nu$ posterior upon relaxing the prior on $A_s$. We observe virtually no change. This result suggests that the neutrino mass constraint is driven by the BAO and not by the small-scale suppression. Although the HOD may constrain the small scales more than the EFT, we suspect similar results would be obtained.

## Appendix E
### Another Look at the Augmentations

At the level of the data vector, the augmentation procedure described in Section 3.3 appears to work relatively well, although the covariance matrix investigation in Figure 3 shows some problems. In order to further quantify how severe the approximation error due to the augmentations is, we focus on posterior-level summaries. In Figure 19, we histogram the standard deviation of the $\sum m_\nu$ posterior for inferences on fiducial and cosmo-varied simulations and compare them to the

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Comparison of standard deviations of neutrino mass posteriors, run on cosmo-varied and fiducial mocks.</figcaption></figure>

standard deviation obtained with real data. We point out that the fiducial simulations all have the same HOD realization, which is not very close to the best fit on real data, while the cosmo-varied simulations have randomly sampled HODs. Therefore, a perfect match between the histograms is likely not expected even if the augmentations reproduced the sample variance perfectly. Rather, Figure 19 serves to develop an upper bound on the approximation error due to the augmentations. We see that the real-data standard deviation is consistent with the distribution of the fiducials, whereas the distribution of the cosmo-varied simulations is somewhat offset.

## Appendix F
### Simulation Budget

One might worry that the 127 cosmo-varied simulations are not enough to properly sample the cosmological prior. We test this by discarding a third of the simulations and training on the rest. The resulting posterior, compared to our baseline result, is shown in Figure 20. Agreement between the two posteriors is almost perfect, demonstrating that our simulations cover the cosmological prior sufficiently well.
<figure><img src="image_23.png" /><figcaption>**Figure 20.** Comparison of posteriors obtained from neural networks trained with the full training set and with two-thirds of the training set. We observe very good agreement between these posteriors.</figcaption></figure>

## Appendix G
Simulation Data

About 50 TB of halo catalogs, light cones, void catalogs, and summary statistics have been saved (at 20 times between $z=0.44$ and $z=0.68$ in 127 different massive-neutrino cosmologies with various HODs and 69 different initial conditions with a fiducial model). We are currently finalizing how to make this data set publicly available.

## Appendix H
Code

In terms of new code, we have written C++ code to populate halo catalogs with galaxies and to generate light cones including survey realism.

We have also written a C implementation of the quasi-random sampling scheme for uniform and Gaussian priors.

This work necessitated several small modifications to public codes:

1. REPS: read files generated by the current CLASS version; write output in a user-defined directory.

2. FastPM: do not write neutrinos to disk.

3. bigfile: support for half-precision floats.

4. Rockstar: native reading of the bigfile snapshots generated by FastPM (using the file chunking to read in distributed fashion since Rockstar does not use MPI).

5. Rockstar/findParents: output to bigfile with lower priority fields in half precision.

6. cuboidremap: support for velocities.

7. sbi: custom splitting into training and validation data.

Since all these items are relatively obscure, we do not provide documentation. However, we are happy to share any of these with interested researchers. A repository with most of the code is available at https://github.com/leanderthiele/nuvoid_production.

## ORCID iDs

Leander Thiele https://orcid.org/0000-0003-2911-9163

Alice Pisani https://orcid.org/0000-0002-6146-4437

ChangHoon Hahn https://orcid.org/0000-0003-1197-0902

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

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

## References

Akiba, T., Sano, S., Yanase, T., Ohta, T., & Koyama, M. 2019, arXiv:1907.10902

Alcock, C., & Paczynski, B. 1979, Natur, 281, 358

Alsing, J., & Wandelt, B. 2019, MNRAS, 488, 5093

Anscombe, F. J. 1948, Biometrika, 35, 246

Aubert, M., Cousinou, M.-C., Escoffier, S., et al. 2022, MNRAS, 513, 186

Bahcall, J. N., & Davis, R., Jr 1976, Sci, 191, 264

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

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

Bayer, A. E., Villaescusa-Navarro, F., Massara, E., et al. 2021b, ApJ, 919, 24

Behroozi, P., Wechsler, R., & Wu, H.-Y., 2012 Rockstar: Phase-space halo finder, Astrophysics Source Code Library, ascl:1210.008

Behroozi, P. S., Wechsler, R. H., & Wu, H.-Y. 2013, ApJ, 762, 109

Berlind, A. A., & Weinberg, D. H. 2002, ApJ, 575, 587

Berlind, A. A., Weinberg, D. H., Benson, A. J., et al. 2003, ApJ, 593, 1

Blas, D., Lesgourgues, J., & Tram, T. 2011, JCAP, 2011, 034

Bond, J. R., Cole, S., Efstathiou, G., & Kaiser, N. 1991, ApJ, 379, 440

Carlson, J., & White, M. 2010, ApJS, 190, 311

Cautun, M., Paillas, E., Cai, Y.-C., et al. 2018, MNRAS, 476, 3195

Chen, S.-F., Vlah, Z., & White, M. 2022, JCAP, 2022, 008

Chudaykin, A., Ivanov, M. M., Philcox, O. H. E., & Simonović, M. 2020, PhRvD, 102, 063533

Colberg, J. M., Pearce, F., Foster, C., et al. 2008, MNRAS, 387, 933

Contarini, S., Marulli, F., Moscardini, L., et al. 2021, MNRAS, 504, 5021

Contarini, S., Pisani, A., Hamaus, N., et al. 2023, ApJ, 953, 46

Contarini, S., Pisani, A., Hamaus, N., et al. 2024, A&A, 682, A20

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

Contarini, S., Verza, G., Pisani, A., et al. 2022, A&A, 667, A162

Cooray, A., & Sheth, R. 2002, PhRvL, 372, 1

Cousinou, M. C., Pisani, A., Tilquin, A., et al. 2019, A&C, 27, 53

Cranmer, K., Brehmer, J., & Louppe, G. 2020, PNAS, 117, 30055

Cranmer, K., Pavez, J., & Louppe, G. 2015, arXiv:1506.02169

d'Amico, G., Gleyzes, J., Kokron, N., et al. 2020, JCAP, 2020, 005

Daya Bay Collaboration, An, F. P., Bai, J. Z., et al. 2012, PhRvL, 108 171803

de Lapparent, V., Geller, M. J., & Huchra, J. P. 1986, ApJL, 302, L1

Delaunoy, A., Hermans, J., Rozet, F., Wehenkel, A., & Louppe, G. 2022, arXiv:2208.13624

Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168

Doré, O., Bock, J., et al. 2014, arXiv:1412.4872

Duffy, A. R., Schaye, J., Kay, S. T., & Dalla Vecchia, C. 2008, MNRAS, 390, L64

Durkan, C., Murray, I., & Papamakarios, G. 2020, arXiv:2002.03712

Feldman, H. A., Kaiser, N., & Peacock, J. A. 1994, ApJ, 426, 23

Feng, Y., Chu, M.-Y., Seljak, U., & McDonald, P. 2016, MNRAS, 463, 2273

Foreman-Mackey, D., Conley, A., Meierjurgen Farr, W., et al., 2013a emcee: The MCMC Hammer, Astrophysics Source Code Library, ascl:1303.002

Foreman-Mackey, D., Hogg, D. W., Lang, D., & Goodman, J. 2013b, PASP, 125, 306

Gregory, S. A., & Thompson, L. A. 1978, ApJ, 222, 784

Guo, H., Zheng, Z., Zehavi, I., et al. 2015, MNRAS, 446, 578

Hahn, C., Eickenberg, M., Ho, S., et al. 2022, arXiv:2211.00723

Hahn, C., Eickenberg, M., Ho, S., et al. 2023, JCAP, 2023, 010

Hamaus, N., Cousinou, M.-C., Pisani, A., et al. 2017, JCAP, 2017, 014

Hamaus, N., Pisani, A., Choi, J.-A., et al. 2020, JCAP, 2020, 023

Hamaus, N., Pisani, A., Sutter, P. M., et al. 2016, PhRvL, 117, 091302

Hamaus, N., Sutter, P. M., & Wandelt, B. D. 2014a, arXiv:1409.7621

Hamaus, N., Sutter, P. M., & Wandelt, B. D. 2014b, PhRvL, 112, 251302

Hand, N., Feng, Y., Beutler, F., et al. 2018, AJ, 156, 160

Hand, N., Feng, Y., Beutler, F., et al., 2019 nbodykit: Massively Parallel, Large-scale Structure Toolkit, Astrophysics Source Code Library, ascl:1904.027

Hearin, A., Tollerud, E., Robitaille, T., et al., 2016a Halotools: Galaxy-Halo Connection Models, Astrophysics Source Code Library, ascl:1604.005

Hearin, A. P., Zentner, A. R., van den Bosch, F. C., Campbell, D., & Tollerud, E. 2016b, MNRAS, 460, 2552

Heavens, A. F., Jimenez, R., & Lahav, O. 2000, MNRAS, 317, 965

Hermans, J., Begy, V., & Louppe, G. 2019, arXiv:1903.04057

Hermans, J., Delaunoy, A., Rozet, F., et al. 2021, arXiv:2110.06581

Hotinli, S. C., Sabti, N., North, J., & Kamionkowski, M. 2023, PhRvD, 108, 103504

Howlett, C., Lewis, A., Hall, A., & Challinor, A. 2012, JCAP, 1204, 027

Hoyle, F., & Vogeley, M. S. 2004, ApJ, 607, 751

Icke, V. 1984, MNRAS, 206, 1P

Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020a, PhRvD, 101, 083504

Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020b, JCAP, 2020, 042

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

Jöveeer, M., Einasto, J., & Tago, E. 1978, MNRAS, 185, 357

K2K Collaboration, Ahn, M. H., Aliu, E., et al. 2006, PhRvD, 74, 072003
Ivezić, Ž., Kahn, S. M., Tyson, J. A., et al. 2019, ApJ, 873, 111

KamLAND Collaboration, Araki, T., Eguchi, K., et al. 2005, PhRvL, 94, 081801

KATRIN Collaboration, Aker, M., Balzer, M., et al. 2022, JPhG, 49, 100501

Kirshner, R. P., Oemler, A. J., Schechter, P. L., & Shectman, S. A. 1981, ApJL, 248, L57

Kobayashi, Y., Nishimichi, T., Takada, M., & Miyatake, H. 2022, PhRvD, 105, 083517

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

Kreisch, C. D., Pisani, A., Villaescusa-Navarro, F., et al. 2022, ApJ, 935, 100

Laureijs, R., Amiaux, J., Arduini, S., et al. 2011, arXiv:1110.3193

Lavaux, G., & Wandelt, B. D. 2012, ApJ, 754, 109

Lewis, A., Challinor, A., & Lasenby, A. 2000, ApJ, 538, 473

Lin, Y.-T., Mandelbaum, R., Huang, Y.-H., et al. 2016, ApJ, 819, 119

Mao, Q., Berlind, A. A., Scherrer, R. J., et al. 2017a, ApJ, 835, 160

Mao, Q., Berlind, A. A., Scherrer, R. J., et al. 2017b, ApJ, 835, 161

Massara, E., & Sheth, R. K. 2018, arXiv:1811.03132

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

Mikheyev, S. P., & Smirnov, A. Y. 1985, YaFiz, 42, 1441

Miller, B. K., Weniger, C., & Forré, P. 2022, arXiv:2210.06170

Moresco, M., et al. 2022, LRR, 25, 6

Nadathur, S. 2016, MNRAS, 461, 358

Nadathur, S., Woodfinden, A., Percival, W. J., et al. 2020, MNRAS, 499, 4140

Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493

Neyrinck, M. C. 2008, MNRAS, 386, 2101

Padilla, N. D., Ceccarelli, L., & Lambas, D. G. 2005, MNRAS, 363, 977

Pan, D. C., Vogeley, M. S., Hoyle, F., Choi, Y.-Y., & Park, C. 2012, MNRAS, 421, 926

Paranjape, A., Lam, T. Y., & Sheth, R. K. 2012, MNRAS, 420, 1429

Paranjape, A., Sheth, R. K., & Desjacques, V. 2013, MNRAS, 431, 1503

Paz, D., Lares, M., Ceccarelli, L., Padilla, N., & Lambas, D. G. 2013, MNRAS, 436, 3480

Pelliciari, D., Contarini, S., Marulli, F., et al. 2023, MNRAS, 522, 152

Philcox, O. H. E. 2021, PhRvD, 103, 103504

Philcox, O. H. E., & Ivanov, M. M. 2022, PhRvD, 105, 043517

Philcox, O. H. E., Ivanov, M. M., Simonović, M., & Zaldarriaga, M. 2020, JCAP, 2020, 032

Pisani, A., Lavaux, G., Sutter, P. M., & Wandelt, B. D. 2014, MNRAS, 443, 3238

Pisani, A., Massara, E., Spergel, D. N., et al. 2019, BAAS, 51, 40

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

Pisani, A., Sutter, P. M., & Wandelt, B. D. 2015b, arXiv:1506.07982

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

Press, W. H., & Schechter, P. 1974, ApJ, 187, 425

Pujol, A., & Gaztañaga, E. 2014, MNRAS, 442, 1930

Reid, B., Ho, S., Padmanabhan, N., et al. 2016, MNRAS, 455, 1553

Reid, B. A., Seo, H.-J., Leauthaud, A., Tinker, J. L., & White, M. 2014, MNRAS, 444, 476

Robotham, A. S. G., & Howlett, C. 2018, RNASS, 2, 55

Ronconi, T., Contarini, S., Marulli, F., Baldi, M., & Moscardini, L. 2019, MNRAS, 488, 5075

Ronconi, T., & Marulli, F. 2017, A&A, 607, A24

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

Sahlén, M., Zubeldía, Í., & Silk, J. 2016, ApJL, 820, L7

Schuster, N., Hamaus, N., Dolag, K., & Weller, J. 2023, JCAP, 2023, 031

Schuster, N., Hamaus, N., Pisani, A., et al. 2019, JCAP, 2019, 055

SDSS Collaboration, Dawson, K. S., Schlegel, D. J., et al. 2013, AJ, 145, 10

SDSS Collaboration, Eisenstein, D. J., Weinberg, D. H., et al. 2011, AJ, 142, 72

SDSS Collaboration, York, D. G., Adelman, J., et al. 2000, AJ, 120, 1579

Semenaite, A., Sánchez, A. G., Pezzotta, A., et al. 2023, MNRAS, 521, 5013

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

Sheth, R. K., & van de Weygaert, R. 2004, MNRAS, 350, 517

Simon, T., Zhang, P., Poulin, V., & Smith, T. L. 2023, PhRvD, 107, 123530

SNO Collaboration, Ahmad, Q. R., Allen, R. C., et al. 2002, PhRvL, 89, 011301

Spergel, D., Gehrels, N., Baltay, C., et al. 2015, arXiv:1503.03757

Stopyra, S., Peiris, H. V., & Pontzen, A. 2021, MNRAS, 500, 4173

Super-Kamiokande Collaboration, Fukuda, Y., Hayakawa, T., et al. 1998, PhRvL, 81, 1562

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. 2012a, ApJ, 761, 44

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

Sutter, P. M., Lavaux, G., Wandelt, B. D., et al. 2014a, MNRAS, 442, 3127

Sutter, P. M., Pisani, A., Wandelt, B. D., & Weinberg, D. H. 2014b, MNRAS, 443, 2983

Takada, M., Ellis, R. S., Chiba, M., et al. 2014, PASJ, 66, R1

Tejero-Cantero, A., Boelts, J., Deistler, M., et al. 2020, JOSS, 5, 2505

Tully, R. B., & Fisher, J. R. 1978, in IAU Symp. 79, Large Scale Structures in the Universe, ed. M. S. Longair & J. Einasto (Dordrecht: Reidel), 31

van de Weygaert, R., & Schaap, W. 2009, in Data Analysis in Cosmology, ed. V. J. Martínez et al., Vol. 665 (Berlin: Springer), 291

van den Bosch, F. C., Weinmann, S. M., Yang, X., et al. 2005, MNRAS, 361, 1203

Variu, A., Alam, S., Zhao, C., et al. 2024, MNRAS, 527, 11539

Verza, G., Carbone, C., Pisani, A., & Renzi, A. 2023, JCAP, 2023, 044

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

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

Wechsler, R. H., & Tinker, J. L. 2018, ARA&A, 56, 435

Wolfenstein, L. 1978, PhRvD, 17, 2369

Woodfinden, A., Percival, W. J., Nadathur, S., et al. 2023, MNRAS, 523 6360

Yoshikawa, K., Jing, Y. P., & Börner, G. 2003, ApJ, 590, 654

Zennaro, M., Bel, J., Villaescusa-Navarro, F., et al., 2016 REPS: REscaled Power Spectra for initial conditions with massive neutrinos, Astrophysics Source Code Library, ascl:1612.022

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

Zentner, A. R., Hearin, A. P., & van den Bosch, F. C. 2014, MNRAS, 443, 3044

Zhai, Z., Tinker, J. L., Banerjee, A., et al. 2023, ApJ, 948, 99

Zheng, Z., Berlind, A. A., Weinberg, D. H., et al. 2005, ApJ, 633, 791

Zhu, G., Zheng, Z., Lin, W. P., et al. 2006, ApJL, 639, L5