# Bye bye, local bias: the statistics of the halo field are not determined by the local mass density

Deaglan J Bartlett, Matthew Ho, Benjamin D Wandelt

► To cite this version:

Deaglan J Bartlett, Matthew Ho, Benjamin D Wandelt. Bye bye, local bias: the statistics of the halo field are not determined by the local mass density. The Astrophysical Journal Letters, 2024, 977 (2), pp.L44. 10.3847/2041-8213/ad97b9. hal-04576359

HAL Id: hal-04576359

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

Submitted on 17 Jan 2025

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

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

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

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

# Bye-bye, Local-in-matter-density Bias: The Statistics of the Halo Field Are Poorly Determined by the Local Mass Density

Deaglan J. Bartlett¹, Matthew Ho¹, and Benjamin D. Wandelt¹,²

¹CNRS & Sorbonne Université, Institut d'Astrophysique de Paris (IAP), UMR 7095, 98 bis bd Arago, F-75014 Paris, France; deaglan.bartlett@iap.fr

²Center for Computational Astrophysics, Flatiron Institute, 162 5th Avenue, New York, NY 10010, USA

Received 2024 May 31; revised 2024 November 25; accepted 2024 November 26; published 2024 December 16

## Abstract

Bias models relating the dark matter field to the spatial distribution of halos are widely used in current cosmological analyses. Many models predict halos purely from the local Eulerian matter density, yet bias models in perturbation theory require other local properties. We assess the validity of assuming that only the local dark matter density can be used to predict the number density of halos in a model-independent way and in the nonperturbative regime. Utilizing N-body simulations, we study the properties of the halo counts field after spatial voxels with near-equal dark matter density have been permuted. If local-in-matter-density (LIMD) biasing were valid, the statistical properties of the permuted and unpermuted fields would be indistinguishable since both represent equally fair draws of the stochastic biasing model. If the Lagrangian radius is greater than approximately half the voxel size and for halos less massive than ~$10^{15} h^{-1} M_{\odot}$, we find the permuted halo field has a scale-dependent bias with greater than 25% more power on scales relevant for current surveys. These bias models remove small-scale power by not modeling correlations between neighboring voxels, which substantially boosts large-scale power to conserve the field's total variance. This conclusion is robust to the choice of initial conditions and cosmology. Assuming LIMD halo biasing cannot, therefore, reproduce the distribution of halos across a large range of scales and halo masses, no matter how complex the model. One must either allow the biasing to be a function of other quantities and/or remove the assumption that neighboring voxels are statistically independent.

Unified Astronomy Thesaurus concepts: Large-scale structure of the universe (902); Dark matter distribution (356); Cosmology (343)

## 1. Introduction

Uncovering the statistical relation between the luminous tracers of large-scale structure and the underlying matter density field—the “bias”—is a crucial part of cosmological analyses since we observe the former but wish to understand the properties of the latter. It has long been known that clusters form at high-density peaks of the density field (N. Kaiser 1984), but an exact relation between these two quantities at the field level remains unknown.

On the largest scales, one expects that the bias is linear, namely, the overdensity of tracers is proportional to the overdensity of matter. This approximation breaks down as one moves from linear to quasi-linear scales, and more general functional relationships must be used. The simplest assumption one can make is that the tracer number density depends solely on the local Eulerian dark matter density—the so-called “local-in-matter-density” (LIMD) assumption (V. Desjacques et al. 2018), which, in a perturbative analysis, allows one to write the tracer number density as an expansion in the overdensity: the local bias expansion (J. N. Fry & E. Gaztanaga 1993). More generally, one can express the bias model as local in gravitational observables, namely, not just the density field but other terms related to the tidal field or higher derivative terms (M. Mirbabayi et al. 2015; V. Desjacques et al. 2018; T. Lazeyras & F. Schmidt 2018, 2019). One can enumerate such terms based on the order of perturbation theory such that,

if one goes to sufficiently high order and is in the large-scale, perturbative regime, then the result converges to the truth.

Relying on perturbation theory (F. Bernardeau et al. 2002) or effective field theory (EFT; J. J. M. Carrasco et al. 2012; L. Senatore & M. Zaldarriaga 2014; L. Senatore 2015; A. Perko et al. 2016) inherently limits our analysis to these quasi-linear scales. Although it is possible to study higher-order statistics (G. D'Amico et al. 2024) or perform field-level analyses in these frameworks (F. Schmidt et al. 2019; F. Elsner et al. 2020; A. Kostić et al. 2023; J. Stadler et al. 2023; N.-M. Nguyen et al. 2024), to extract information on nonlinear scales, one may need to go beyond such local bias expansions although some studies have shown that perturbative methods perform well when tested against N-body simulations (N. Roth & C. Porciani 2011; M. Schmittfull et al. 2019). Instead of relying on power series, many parametric forms for the biasing relation have been proposed in the literature based on the LIMD assumption (A. S. Szalay 1988; T. Matsubara 1995, 2011; N. Frusciante & R. K. Sheth 2012; M. C. Neyrinck 2014; M. Ata et al. 2015). As such, many reconstruction analyses typically assume that the number density of halos is solely a function of the smoothed density field (e.g., I. M. Schmoldt et al. 1999; P. Erdoğdu et al. 2004, 2006; F. S. Kitaura et al. 2009; J. Jasche et al. 2010, 2015; J. Jasche & B. D. Wandelt 2012, 2013; G. Lavaux & J. Jasche 2016; C. Modi et al. 2018; J. Jasche & G. Lavaux 2019; G. Lavaux et al. 2019; D. K. Ramanah et al. 2019). Alternatively, one can utilize machine learning methods to map the dark matter density field (using both local and nonlocal information) to the halo counts field (T. Charnock et al. 2020; see B. Dai & U. Seljak 2021 for a Lagrangian approach).

Unlike in the perturbative regime, we do not have theoretical guarantees that any of these bias models are accurate, and so this

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

Original content from this work may be used under the terms of the [Creative Commons Attribution 4.0 licence](https://creativecommons.org/licenses/by/4.0/). Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
must be tested empirically. M. Mirbabayi et al. (2015) argue that, at any order in perturbation theory, the bias can be expressed in terms of locally measurable observables and that this argument is closed under renormalization. However, it could be (as is often assumed; A. S. Szalay 1988; T. Matsubara 1995, 2011; N. Frusciante & R. K. Sheth 2012; M. C. Neyrinck 2014; M. Ata et al. 2015) that the only important observable is the local dark matter density field. Testing this assumption is the subject of this Letter. There is an important distinction, however, between testing the assumptions of the halo biasing model (i.e., which variables are needed to predict the tracer number density) and the particular model chosen (i.e., what is the functional form of the bias or the stochasticity in the tracer-matter relation). Previous tests of the latter seem to suggest that information beyond the local density is required (e.g., N.-M. Nguyen et al. 2021; F. Fang et al. 2024); however, this could just be due to insufficiently accurate LIMD biasing models. Initial studies of the former also seem to support this conclusion (A. Balaguera-Antolínez et al. 2019; A. Balaguera-Antolínez et al. 2020; M. Pellejero-Ibañez et al. 2020).

In this Letter, we study the former issue in the fully nonlinear regime.

We present a model-independent test of the assumption of LIMD halo biasing. Our test is independent of the functional form of the biasing model but simply tests whether halo biasing can be a stochastic function of only the local matter density. We find that LIMD halo biasing cannot reproduce the distribution of halos less massive than ~$10^{15} h^{-1} M_{\odot}$ when applied to gridded dark matter fields, when the Lagrangian radius (Equation (2)) of the halo is larger than approximately half the voxel size. Our permutation test naturally enforces conservation of all one-point moments and count statistics while demonstrating the insufficiency of the LIMD hypothesis to model halo clustering. We demonstrate that this conclusion is robust across initial conditions and cosmology by utilizing the QUIJOTE suite of N-body simulations. We conclude that no LIMD biasing model, however complex, should be used if one wishes to accurately predict the spatial distribution of halos given a matter density field.

Our test also provides a method for testing the convergence of halo catalogs from N-body simulations; we find that the minimum mass for which the catalogs are converged for QUIJOTE is ~$10^{13.1} h^{-1} M_{\odot}$.

In Section 2, we explicitly define LIMD biasing, introduce our model-independent test, and describe the simulations used to perform this test. Our results are presented in Section 3, where we demonstrate that our conclusions are highly robust to the one free parameter in our test, as well as cosmology and initial conditions, and investigate the adequacy of LIMD biasing as a function of scale and halo mass. Section 4 discusses these results and presents our conclusions.

## 2. Methods

### 2.1. Assessing the Assumption of Local Biasing

We define LIMD halo biasing to be the assumption that the number of halos in a given region of comoving space (voxel) is drawn from any probability distribution under the assumptions:

1. The parameters of this distribution in a given voxel are only functions of the dark matter density in that voxel.

2. Each voxel is statistically independent.

Note that we will not make any assumptions about the type of distribution (e.g., Poisson) nor the functional form of its parameters (e.g., the mean for a Poisson distribution) other than we expect that this function should be a continuous function of the local dark matter density. The statement about statistical independence is made within the context of conditional probability. Specifically, while the underlying quantities that influence the distribution in each voxel are indeed correlated among voxels (inherited from the correlation due to large-scale structure), the stochastic quantities drawn from this distribution at each voxel are treated as independent random variables. If the particular draw from the conditional distribution depended on other voxels in any other way than through this inherent large-scale structure correlation, then one would need information beyond the local density field to reconstruct the halo counts field, and thus, this would not be a model that is not local in matter density.

Under these assumptions, if we identify several voxels that have the same dark matter density and then permute the number of halos in these voxels, we should end up with an equally likely draw from the underlying distribution. If the statistical properties of the halo-count field are noticeably different after this permutation, then at least one of the assumptions of LIMD biasing cannot be true.

Since the dark matter density field takes continuous values, instead of performing permutations on voxels that have exactly equal dark matter density, we first sort the voxels into $n_{\text{bin}}$ bins of dark matter density and consider two voxels equivalent if they fall within the same bin, namely, if their density is approximately equal. In our fiducial analysis, we choose to use adaptive bin widths so that each bin has approximately the same number of voxels and use $n_{\text{bin}} = 1000$. We demonstrate that our results are highly robust to this choice in Section 3.

To assess the level of agreement between the permuted and unpermuted fields, we compute the power spectrum, $P_{\text{counts}}(k)$, of the counts field using the PYLIANS library (F. Villaescusa-Navarro 2018). We perform 100 permutations and report the mean and standard deviation at each wavenumber, then compare this to the true $P_{\text{counts}}(k)$.

### 2.2. Simulations

To test the assumption of LIMD biasing, we utilize the QUIJOTE suite of N-body simulations. These simulations were run in a periodic box of length $L = 1 h^{-1}$ Gpc using the GADGET-III code (V. Springel 2005). We focus on the high-resolution suite that contains $N_p^3 = 1024^3$ particles, as opposed to the standard resolution of $N_p^3 = 512^3$.

For our fiducial analysis, we consider simulation number 0 from the fiducial cosmology, which matches the Planck 2018 cosmological parameters (Planck Collaboration 2020): $\Omega_m = 0.3175$, $\Omega_b = 0.049$, $h = 0.6711$, $n_s = 0.9624$, and $\sigma_8 = 0.834$. In Section 3, we verify that our conclusions are robust if we chose different initial conditions (i.e., for different simulation numbers in the suite). We also study the sensitivity of our conclusions to the choice of cosmological parameters by performing our analysis with the (high-resolution) Latin hypercube suite of QUIJOTE simulations, consisting of 2000 dark-matter-only simulations with cosmological parameters arranged on a Latin hypercube in the range: $\Omega_m \in [0.1, 0.5]$, $\Omega_b \in [0.03, 0.07]$, $h \in [0.5, 0.9]$, $n_s \in [0.8, 1.2]$, and $\sigma_8 \in [0.6, 1.0]$.

We make all our comparisons at redshift zero, where we use the positions and masses of halos from the precomputed
Friends-of-Friends (FoF; M. Davis et al. 1985) halo catalogs, which were obtained using a linking length of $b=0.2$. For brevity, throughout this Letter, we only report results using FoF halos, but we have verified that our conclusions are robust to this choice by rerunning our analysis with the Rockstar halofinder (P. S. Behroozi et al. 2013) and found consistent conclusions.

Since one would expect halos of different masses to be distributed differently in space, we binned these halos into nine logarithmically spaced mass bins in the range $10^{13.1}-10^{15.8} h^{-1} M_{\odot}$ and performed our test separately for each bin. For halos within this mass range, we find that our results are consistent between the standard and high-resolution QUIJOTE simulations, whereas this is not true for lower masses. We therefore cannot be confident that the halo catalogs are converged for halos that are less massive than this, so we simply discard them in our analysis. Our test, therefore, additionally acts as a method of determining the minimum mass for which the halo catalogs from these simulations are converged.

We divide the simulation volume into $N^3$ cubic voxels and compute the density field in each voxel using a nearest grid point (NGP) estimator. This is the same estimator used to compute the gridded halo counts field, and thus, we mitigate effects that arise due to different discretization procedures. We have tested our method using a range of other density field estimators: cloud in cell, triangular-shaped cloud, and piecewise cubic spline kernel, as well as one inspired by smoothed particle hydrodynamics (J. J. Monaghan 1992; S. Colombi et al. 2007). For our fiducial analysis, we used $N=128$, corresponding to voxels with a side length of $\Delta x = 7.8 h^{-1} \text{ Mpc}$, but we show that our results are insensitive to this choice (and thus to the smoothing scale of the density field) in Section 3. At our fiducial grid resolution, we find that all density estimators give the same qualitative trends although the magnitude of the bias induced by LIMD models varies with this choice. For large voxel sizes, we find the performance is best for the NGP density estimator, but the degree of bias depends sensitively on the choice of mass assignment in the definition of the density field and is much worse for other choices. We thus choose the NGP estimator to ensure our results are conservative.

## 3. Results

### 3.1. LIMD Halo Biasing Cannot Reproduce the Halo Distribution

In Figure 1, we plot $P_{\text{counts}}(k)$ for our fiducial setup for halos in the mass range $10^{13.7}-10^{14.0} h^{-1} M_{\odot}$, both before and after permuting voxels with equal dark matter density. When compared to the truth, we see that the permuted fields have very different halo distributions than the original field on all scales. In particular, we observe a significant excess of power on scales $k \lesssim 0.3 h \text{ Mpc}^{-1}$ and too little power on smaller scales. We infer that this is due to the removal of small-scale power due to LIMD models ignoring local correlations, which then requires an overestimation of large-scale power, to respect the conservation of the field-variance under voxel permutations. This represents the main conclusion of this Letter: the assumption of LIMD halo biasing cannot produce the correct distribution of halos.

The vertical dashed lines in Figure 1 correspond to four characteristic wavenumbers that we use for the remainder of the

<figure><img src="image_4.png" /><figcaption>**Figure 1.** Power spectrum of the halo counts field before and after permuting voxels with approximately equal density, for halos of mass $10^{13.7}-10^{14.0} h^{-1} M_{\odot}$. The permuted fields have significantly greater power on large scales; thus, LIMD halo biasing cannot reproduce the distribution of halos. The mean and standard deviation are computed over 100 permutations. The vertical colored lines correspond to the wavenumbers used to characterize this discrepancy in later plots.</figcaption></figure>

analysis to quantify the level of discrepancy, corresponding to approximately 0.01, 0.03, 0.1, and $0.5 h \text{ Mpc}^{-1}$. For our fiducial analysis, we find that $P_{\text{counts}}(k)$ is a factor of $1.6 \pm 0.2$, $1.3 \pm 0.07$, $1.16 \pm 0.02$, and $0.97 \pm 0.01$ times the true $P_{\text{counts}}(k)$, respectively, for these wavenumbers. It is therefore clear that the discrepancy is significant for this simulation, with the true power spectrum at approximately the $7.6\sigma$ value for $k=0.1 h \text{ Mpc}^{-1}$.

In Figure 2, we investigate how this conclusion changes as a function of halo mass and voxel size. For the smaller voxel sizes ($\Delta x \lesssim 30 h^{-1} \text{ Mpc}$), one observes that the fractional difference is nonzero until $\sim 10^{15} h^{-1} M_{\odot}$, indicating that LIMD halo biasing is invalid for all masses below this point. For the largest halo masses and voxels, we find that the power spectrum of the permuted field is consistent with the truth. The number density of these objects is low—there are only 381 halos above $10^{15} h^{-1} M_{\odot}$ in the fiducial simulation—and these can only exist at the extreme peaks of the density field. Hence, given that these objects cannot exist in low-density environments, it is perhaps unsurprising that their number density can be well predicted by just the local density field. We therefore conclude that LIMD halo biasing is a reasonable assumption for halos above $\sim 10^{15} h^{-1} M_{\odot}$, but, for the majority of objects in the simulation, this framework cannot be used when the voxel size is below $\sim 30 h^{-1} \text{ Mpc}$.

In Figure 3, we investigate how the error on $P_{\text{counts}}(k)$ varies as a function of grid size for our fiducial mass bin. One observes that the level of disagreement is relatively unchanged on the largest scales as one changes the voxel size from $\sim 4 h^{-1} \text{ Mpc}$ to $\sim 30 h^{-1} \text{ Mpc}$, indicating that LIMD biasing is inappropriate for density fields smoothed on these scales. As one further increases the voxel size to $\sim 60 h^{-1} \text{ Mpc}$, the bias starts to diminish, with a bias of $\sim 1.4\sigma$ at $k = 0.03 h \text{ Mpc}^{-1}$. In this regime, the number of voxels is similar to the number of density bins (as indicated by the black line in Figure 3), but we have verified that similar behavior is seen if we reduce $n_{\text{bin}}$ to 100. However, even on scales $\sim 0.01 h \text{ Mpc}^{-1}$, we see a mean
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Fractional error on the power spectrum after permuting equal-density voxels as a function of halo mass and voxel size, $\Delta x$. Significant discrepancies are seen for a wide range of scales for halo masses below $10^{15} h^{-1} M_{\odot}$ and when $\Delta x \lesssim 30 h^{-1} Mpc$; hence, LIMD halo biasing is invalid in this regime. The points and error bars give the mean and standard deviation of the discrepancy across 100 permutations, respectively.</figcaption></figure>

fractional error of 25% for a voxel size of $\sim 62.5 h^{-1} Mpc$. Thus, the LIMD halo biasing is invalid for smoothed density fields with a smoothing scale smaller than $\sim 30 h^{-1} Mpc$ although it may be an acceptable approximation on larger scales.

To further investigate how the error on $P_{\text{counts}}(k)$ varies as a function of grid size and mass, we first define a way of summarizing the bias across the different wavenumbers. If the fractional error at each $k$ is $\epsilon(k)$ with standard deviation $\sigma(k)$, then we define

$$ \mathcal{B}^2 = \frac{\int \left( \frac{\epsilon(k)}{\sigma(k)} \right)^2 dk}{\int dk} - 1. \qquad (1) $$

If the fractional error had a constant mean of $n\sigma(k)$ for all wavenumbers and each wavenumber were independent, then

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Fractional error on the power spectrum for halos of mass $10^{13.7}-10^{14.0} h^{-1} M_{\odot}$ after permuting equal-density voxels as a function of voxel size. The black dashed line indicates the value of $N$ where we have as many voxels as density bins. For $N \ge 2^6$, we find significant disagreement on large scales, suggesting that one must smooth the density field over at least $30 h^{-1} Mpc$ for LIMD halo biasing to be valid. We do not plot values for wavenumbers beyond the Nyquist frequency; hence, some lines are truncated.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 4.** The significance of the deviation (Equation (1)) between the power spectra of the true and permuted fields as a function of halo mass, for different voxel sizes, $\Delta x$. The dashed line corresponds to $\mathcal{B} = 1$. Large values of $\mathcal{B}$ correspond to scenarios where LIMD is a poor approximation.</figcaption></figure>

this gives $\langle \mathcal{B}^2 \rangle = n^2$, and thus, $\mathcal{B}$ gives a dimensionless estimate of the bias, where $\mathcal{B} = 0$ corresponds to being unbiased.

We compute this quantity for the various mass bins and resolutions and plot the resulting $\mathcal{B}$ as a function of halo mass in Figure 4. As before, we find that the level of bias is small for $\Delta x \gtrsim 30 h^{-1} Mpc$, with $\mathcal{B} < 1$ for all halo masses. For the most massive halos (above $\sim 10^{15} h^{-1} M_{\odot}$), we again see no indication of bias for any grid resolution, presumably since these objects are rare and can only exist in the densest environments. As we move to smaller masses and $\Delta x < 30 h^{-1} Mpc$, we find that the bias increases, with $\mathcal{B} > 8$ for $\Delta x = 3.9 h^{-1} Mpc$.

Intriguingly, we find that the value of $\mathcal{B}$ does not monotonically increase as one decreases the halo mass, but for $\Delta x = 15.6 h^{-1} Mpc$, $\mathcal{B}$ reaches a maximum of $\sim 2$ at $\sim 10^{14.2} h^{-1} M_{\odot}$ before decreasing to below 1 at masses below $\sim 10^{13.4} h^{-1} M_{\odot}$. One can begin to see the same effect for $\Delta x = 7.8 h^{-1} Mpc$, with the maximum at $\sim 10^{13.5} h^{-1} M_{\odot}$. These maxima occur when the voxel size is approximately twice the Lagrangian radius of the halo, defined such that the halo mass, $M_h$, obeys

$$ M_h \equiv \frac{4\pi}{3} \rho_m R_L^3, \qquad (2) $$

where $\rho_m$ is the cosmic mean matter density. Therefore, one would expect to see a similar maximum for $\Delta x = 3.9 h^{-1} Mpc$
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Fractional error on the power spectrum of halo counts after permuting voxels within the same density bin as a function of the number of bins. We consider halos of mass $10^{13.7}-10^{14.0} h^{-1} M_{\odot}$ and choose a binning scheme such that each bin contains approximately an equal number of voxels. After an initial transient behavior, we find constant biases in the power spectrum across many orders of magnitude; hence, our conclusion is robust to the choice of $n_{\text{bin}}$.</figcaption></figure>

for a halo mass of ~$10^{12.4} h^{-1} M_{\odot}$. However, these halos are not resolved in the QUIJOTE simulations, so we cannot test this conjecture here. As one decreases the halo mass below this maximum, the voxel size becomes significantly larger than the Lagrangian radius of the halo, and thus LIMD biasing becomes an increasingly good approximation. If the Lagrangian patch is larger than the voxel size, then information from neighboring voxels is required to determine the number and positions of halos. Hence, we conclude that there exists an intermediate range of masses for a given voxel size where LIMD biasing breaks down, where the upper limit is due to the rarity of very massive halos, and the lower limit is set by the ability to resolve the Lagrangian radius of individual halos.

## 3.2. Sensitivity to Binning Scheme

If one did not bin the density sufficiently finely, then one would expect that this shuffling procedure would not reproduce the correct distribution of halos even if LIMD halo biasing were correct since voxels with very different densities would be considered equivalent, and therefore, permuting their halos would not represent any realistic biasing model. To test this, in Figure 5 we plot the fractional error on the power spectrum as a function of $n_{\text{bin}}$. When one has only one density bin, the permuted halo power spectrum is consistent with random noise. As we increase $n_{\text{bin}}$, one sees an initial transient behavior, where the fractional error becomes positive on large scales at approximately $n_{\text{bin}} = 10$ and then increases to a constant level by approximately $n_{\text{bin}} = 30$. Beyond this point, the fractional error on $P_{\text{counts}}(k)$ is independent across orders of magnitude of $n_{\text{bin}}$ suggesting that our conclusions are robust to this choice.

We have verified that our results are unchanged if we consider equally spaced bins in both density and the logarithm of density, instead of enforcing each bin to have approximately an equal number of voxels. The transient behavior disappears at $n_{\text{bin}} \approx 10$ when using logarithmic bins and at $n_{\text{bin}} \approx 300$ for the linear binning scheme. Above these values, we again find that the bias is insensitive to the choice of $n_{\text{bin}}$. We have also verified that our binning is sufficiently fine such that shuffling

the density field under this scheme yields a matter power spectrum that is consistent with the unshuffled field.

Of course, if one had as many bins as number of voxels ($n_{\text{bin}} = N^3$), then (assuming no voxels have identically equal densities) there would be only one voxel per bin, and thus, the permuted halo field is equal to the unpermuted field. For $N=128$, this corresponds to $n_{\text{bin}} \approx 10^{6.3}$. Since our conclusion holds at $n_{\text{bin}} = 10^5$, even with an average of just 21 voxels per bin we cannot reproduce the halo distribution. An LIMD halo biasing model that varies that rapidly with density seems implausible and highly fine-tuned. We have verified that using $n_{\text{bin}} = 100$ for a mock halo counts field generated from a Poisson distribution with a mean that depends on density according to a power law (with parameters optimized to fit the QUIJOTE simulations using a Poisson loss function in each voxel) yields unbiased $P_{\text{counts}}(k)$ upon shuffling, again indicating that our fiducial choice of $n_{\text{bin}} = 10^3$ is sufficient.

## 3.3. Sensitivity to Cosmology and Initial Conditions

To test the robustness of our conclusion to the initial conditions, we repeated our analysis across the 100 high-resolution QUIJOTE simulations at the fiducial cosmology. We find biases (mean divided by standard deviation across 100 permutations) of $1.4 \pm 0.9$, $4.5 \pm 1.0$, $6.6 \pm 0.9$, and $-4.7 \pm 1.0$ at the wavenumbers shown in Figure 1, indicating that our conclusion is insensitive to the initial white noise field.

Similarly, to investigate the sensitivity of the result to the cosmological parameters, we ran our analysis for all 2000 high-resolution simulations in the Latin hypercube suite. We find that, averaging over cosmological parameters, the biases for these wavenumbers are $1.7 \pm 1.4$, $4.2 \pm 1.5$, $6.4 \pm 1.2$, and $-4.5 \pm 1.1$ demonstrating that LIMD halo biasing does not hold for any of the cosmologies considered.

# 4. Discussion and Conclusion

In this Letter, we have shown that the distribution of halos cannot be modeled according to LIMD biasing: the assumption that the number of halos in a given voxel is drawn from a probability distribution that depends solely on the local dark matter density and where each voxel is statistically independent. At least one of these assumptions must be broken to accurately model the distribution of halos in the Universe when averaging over scales of ~$4-30 h^{-1}$ Mpc. Of course, one would expect that if one went to sufficiently small scales, halos would extend to multiple voxels, and thus, LIMD biasing must break once the voxels are made sufficiently small. This work quantifies this scale and its relation to the halo's Lagrangian radius (Equation (2)) using N-body simulations.

We have demonstrated this in a model-independent way, without considering a particular choice of probability distribution nor the dependence of its parameters on the local dark matter density. By permuting the number of halos found in voxels of approximately equal density, we showed that the permuted field has too much power on large scales, whereas this would be statistically indistinguishable from the unpermuted field if LIMD biasing were true. This conclusion is robust to the definition of “approximately equal density” as well as voxel size, cosmology, and initial conditions. We find that the most massive halos can be modeled under these assumptions, but this breaks down for halos less massive than ~$10^{15} h^{-1} M_{\odot}$.
The test described in this Letter is similar in approach to the first step of the Bias Assignment Method (BAM; A. Balaguera-Antolínez et al. 2019), which is used to fit nonparametric bias models. As here, A. Balaguera-Antolínez et al. (2019) find that information beyond the local density is required to recover the halo distribution but perform their analysis with a single-grid resolution of $\Delta x = 3 h^{-1}$ Mpc, which is slightly finer than what we consider. Our work extends this conclusion by systematically exploring the sensitivity to this smoothing scale and explicitly investigating the mass dependence of this conclusion in separate mass bins. BAM has recently been used to provide evidence for nonlocality in the bias model by studying the bispectrum (M. Pellejero-Ibañez et al. 2020), with up to $4.8\sigma$ significance (F.-S. Kitaura et al. 2022). The conclusions presented here are much stronger than this. Moreover, we demonstrate that the results are robust to cosmological parameters and provide more thorough tests of the robustness of our results to the binning scheme. Unlike A. Balaguera-Antolínez et al. (2019), we have demonstrated that our conclusions are robust to the resolution of the N-body simulation used in our test.

This conclusion is particularly important as the use of higher-order statistics and field-level approaches become increasingly popular in cosmological analyses since it is no longer sufficient to tune a bias model to only reproduce the power spectrum or other low-order statistics; one must be able to reproduce all moments of the field. In field-level studies, the initial conditions of the Universe are free parameters, so using an LIMD biasing model implies inferring an incorrect density field to obtain a predicted galaxy (or whichever tracer is used in the analysis) field that reasonably matches observations (N.-M. Nguyen et al. 2021). Although the inferred galaxy field at the redshift of the observations would be trustworthy, extrapolating to other redshifts could be problematic, and the underlying density field and its initial conditions would have artifacts. Given that it is often these initial conditions that one is interested in, it is imperative that future studies move beyond the LIMD halo biasing assumption to avoid these problems.

If the issue causing the deficiency of LIMD models was chiefly the dependence on local density (i.e., independent draws is a reasonable approximation), then one could apply a similar permutation technique to identify whether other halo biasing schemes are valid in a model-independent fashion. For example, by binning in multiple dimensions, one could assess whether the variables used in EFT based approaches are sufficient in the nonperturbative regime. Alternatively, one could search for transformations of the density field for which the permuted and unpermuted halo number counts are indistinguishable when binned in this variable and then construct bias models with these new parameters.

As discussed in Section 2, the degree to which LIMD biasing fails depends on the estimator used to obtain the density field. This high sensitivity to the estimator is problematic if one wishes to predict the halo distribution from an approximate gravity model. As such, using information beyond the local density field is essential to produce realistic halo catalogs for analysis pipelines for current and future surveys.

We conclude that the number of halos in a given region of space is a stochastic function of the local dark matter density field, but it cannot depend solely on this quantity. Models that rely on this assumption, no matter how complex, cannot

reproduce the distribution of halos and thus must be superseded by more accurate models. This Letter not only disproves the LIMD assumption of halo biasing for scales that are not significantly greater than the Lagrangian radius of a halo but provides a framework to aid the design of future approaches by testing their assumptions in a model-independent manner.

## Acknowledgments

We thank Andrés Balaguera, Yan-Chuan Cai, Neal Dalal, Simon Ding, Ludvig Doeser, Rozo Eduardo, Pedro Ferreira, Joachim Harnois-Deraps, Jens Jasche, Guilhem Lavaux, Titouan Lazeyras, Stuart McAlpine, Minh Nguyen, Cora Uhlemann, Francisco Villaescusa-Navarro, and the Implicit Likelihood Inference, Accelerated Forward Modeling and Robustness working groups of the Simons Collaboration on “Learning the Universe” for useful comments and suggestions.

D.J.B. and M.H. are supported by the Simons Collaboration on “Learning the Universe.” B.D.W. acknowledges the DIM-ORIGINES-2023, Infinity Next grant. The Flatiron Institute is supported by the Simons Foundation. This work was done within the Aquila Consortium (https://www.aquila-consortium.org/).

For the purposes of open access, the authors have applied a Creative Commons Attribution (CC BY) license to any Author Accepted Manuscript version arising.

*Software:* ASTROPY (Astropy Collaboration et al. 2013, 2018, 2022), MATPLOTLIB (J. D. Hunter 2007), NUMPY (C. R. Harris et al. 2020), PANDAS (pandas development team, T 2020; Wes McKinney 2010), PYLIANS (F. Villaescusa-Navarro 2018).

## ORCID iDs

Deaglan J. Bartlett https://orcid.org/0000-0001-9426-7723

Matthew Ho https://orcid.org/0000-0003-3207-8868

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

## References

Astropy Collaboration, Price-Whelan, A. M., Lim, P. L., et al. 2022, ApJ, 935, 167

Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al. 2018, AJ, 156, 123

Astropy Collaboration, Robitaille, T. P., Tollerud, J. E., et al. 2013, A&A, 558, A33

Ata, M., Kitaura, F.-S., & Müller, V. 2015, MNRAS, 446, 4250

Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibanez, M., Zhao, C., & Abel, T. 2019, MNRAS, 483, L58

Balaguera-Antolínez, A., Kitaura, F.-S., Pellejero-Ibañez, M., et al. 2020, MNRAS, 491, 2565

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

Bernardeau, F., Colombi, S., Gaztañaga, E., & Scoccimarro, R. 2002, PhR, 367, 1

Carrasco, J. J. M., Hertzberg, M. P., & Senatore, L. 2012, JHEP, 2012, 82

Charnock, T., Lavaux, G., Wandelt, B. D., et al. 2020, MNRAS, 494, 50

Colombi, S., Chodorowski, M. J., & Teyssier, R. 2007, MNRAS, 375, 348

Dai, B., & Seljak, U. 2021, PNAS, 118, e2020324118

D'Amico, G., Donath, Y., Lewandowski, M., Senatore, L., & Zhang, P. 2024, JCAP, 05, 059

Davis, M., Efstathiou, G., Frenk, C. S., & White, S. D. M. 1985, ApJ, 292, 371

Desjacques, V., Jeong, D., & Schmidt, F. 2018, PhR, 733, 1

Elsner, F., Schmidt, F., Jasche, J., Lavaux, G., & Nguyen, N.-M. 2020, JCAP, 2020, 029

Erdoğdu, P., Lahav, O., Huchra, J. P., et al. 2006, MNRAS, 373, 45

Erdoğdu, P., Lahav, O., Zaroubi, S., et al. 2004, MNRAS, 352, 939

Fang, F., Cai, Y.-C., Li, Z., et al. 2024, MNRAS, 530, 2355

Frusciante, N., & Sheth, R. K. 2012, JCAP, 2012, 016

Fry, J. N., & Gaztañaga, E. 1993, ApJ, 413, 447

Harris, C. R., Millman, K. J., van der Walt, S. J., et al. 2020, Natur, 585, 357
Hunter, J. D. 2007, CSE, 9, 90

Jasche, J., Kitaura, F. S., Wandelt, B. D., & Enßlin, T. A. 2010, MNRAS, 406, 60

Jasche, J., & Lavaux, G. 2019, A&A, 625, A64

Jasche, J., Leclercq, F., & Wandelt, B. D. 2015, JCAP, 2015, 036

Jasche, J., & Wandelt, B. D. 2012, MNRAS, 425, 1042

Jasche, J., & Wandelt, B. D. 2013, MNRAS, 432, 894

Kaiser, N. 1984, ApJL, 284, L9

Kitaura, F.-S., Balaguera-Antolínez, A., Sinigaglia, F., & Pellejero-Ibanez, M. 2022, MNRAS, 512, 2245

Kitaura, F. S., Jasche, J., Li, C., et al. 2009, MNRAS, 400, 183

Kostić, A., Nguyen, N.-M., Schmidt, F., & Reinecke, M. 2023, JCAP, 2023, 063

Lavaux, G., & Jasche, J. 2016, MNRAS, 455, 3169

Lavaux, G., Jasche, J., & Leclercq, F. 2019, arXiv:1909.06396

Lazeyras, T., & Schmidt, F. 2018, JCAP, 2018, 008

Lazeyras, T., & Schmidt, F. 2019, JCAP, 2019, 041

Matsubara, T. 1995, ApJS, 101, 1

Matsubara, T. 2011, PhRvD, 83, 083518

McKinney, W. 2010, Proc. 9th Python in Science Conf., ed. S. van der Walt & J. Millman, 56

Mirbabayi, M., Schmidt, F., & Zaldarriaga, M. 2015, JCAP, 2015, 030

Modi, C., Feng, Y., & Seljak, U. 2018, JCAP, 2018, 028

Monaghan, J. J. 1992, ARA&A, 30, 543

Neyrinck, M. C. 2014, in IAU Symp. 306, Statistical Challenges in 21st Century, ed. A. Heavens, J.-L. Starck, & A. Krone-Martins (Cambridge: Cambridge Univ. Press), 251

Nguyen, N.-M., Schmidt, F., Lavaux, G., & Jasche, J. 2021, JCAP, 2021, 058

Nguyen, N.-M., Schmidt, F., Tucci, B., Reinecke, M., & Kostić, A. 2024, arXiv:2403.03220

pandas development team, T 2020, pandas-dev/pandas: Pandas, v2.0.3, Zenodo, doi: 10.5281/zenodo.3509134

Pellejero-Ibañez, M., Balaguera-Antolínez, A., Kitaura, F.-S., et al. 2020, MNRAS, 493, 586

Perko, A., Senatore, L., Jennings, E., & Wechsler, R. H. 2016, arXiv:1610.09321

Planck Collaboration 2020, A&A, 641, A6

Ramanah, D. K., Lavaux, G., Jasche, J., & Wandelt, B. D. 2019, A&A, 621, A69

Roth, N., & Porciani, C. 2011, MNRAS, 415, 829

Schmidt, F., Elsner, F., Jasche, J., Nguyen, N. M., & Lavaux, G. 2019, JCAP, 2019, 042

Schmittfull, M., Simonović, M., Assassini, V., & Zaldarriaga, M. 2019, PhRvD, 100, 043514

Schmoldt, I. M., Saar, V., Saha, P., et al. 1999, AJ, 118, 1146

Senatore, L. 2015, JCAP, 2015, 007

Senatore, L., & Zaldarriaga, M. 2014, arXiv:1409.1225

Springel, V. 2005, MNRAS, 364, 1105

Stadler, J., Schmidt, F., & Reinecke, M. 2023, JCAP, 2023, 069

Szalay, A. S. 1988, ApJ, 333, 21

Villaescusa-Navarro, F., 2018 Pylians: Python Libraries for the Analysis of Numerical Simulations, Astrophysics Source Code Library, ascl:1811.008