TY - JOUR
AU1 - Mead, A J
AU2 - Verde, L
AB - ABSTRACT We derive a simple prescription for including beyond-linear halo bias within the standard, analytical halo-model power spectrum calculation. This results in a corrective term that is added to the usual two-halo term. We measure this correction using data from N-body simulations and demonstrate that it can boost power in the two-halo term by a factor of ∼2 at scales |$k\sim 0.7\, h\mathrm{Mpc}^{-1}$|⁠, with the exact magnitude of the boost determined by the specific pair of fields in the two-point function. How this translates to the full power spectrum depends on the relative strength of the one-halo term, which can mask the importance of this correction to a greater or lesser degree, again depending on the fields. Generally, we find that our correction is more important for signals that arise from lower mass haloes. When comparing our calculation to simulated data, we find that the underprediction of power in the transition region between the two- and one-halo terms, which typically plagues halo-model calculations, is almost completely eliminated when including the full non-linear halo bias. We show improved results for the autospectra and cross-spectra of galaxies, haloes, and matter. In the specific case of matter–matter or matter–halo power, we note that a large fraction of the improvement comes from the non-linear biasing between low- and high-mass haloes. We envisage our model being useful in the analytical modelling of cross-correlation signals. Our non-linear bias halo-model code is available at https://github.com/alexander-mead/BNL. large-scale structure of Universe, cosmology: theory 1 INTRODUCTION The halo model (reviewed by Cooray & Sheth 2002) is widely used in the interpretation of data from cosmological large-scale structure surveys. The model is not derived from first principles, and is instead phenomenological in that it is a description of the properties of the universe as seen in N-body simulations. In the model, all matter is taken to reside in haloes that trace the large-scale matter fluctuations in a biased way, with this ‘halo bias’ usually taken to be linear with respect to the underlying linear matter field. The model also makes a number of other simplifying assumptions; usually that haloes are spherical, devoid of substructure, and that the halo mass determines all of the halo properties with no scatter. In addition, choices must be made for the mass function, biasing recipe, and profile for haloes. Fields, be they sourced by point tracers or else via some emissivity, can then be assumed to occupy haloes in different ways depending on the halo properties; the model will then make predictions for the n-point correlation functions. The halo model has been successful in explaining the broad shape of the galaxy–galaxy correlation function (e.g. Peacock & Smith 2000; Seljak 2000), as well as making inroads in the understanding of the connection between galaxy formation and the haloes in which the galaxies reside (e.g. Mandelbaum et al. 2005; Cacciato et al. 2012) and is also widely used to describe other cross-correlations (e.g. Addison, Dunkley & Spergel 2012; Hill & Spergel 2014; Ma et al. 2015; Padmanabhan, Refregier & Amara 2017; Hill et al. 2018; Tanimura et al. 2019; Wolz et al. 2019; Koukoufilippas et al. 2020). It has been 20 yr since the modern formulation of the halo model came to prominence, and since then the quality of data to which the halo model is exposed have increased dramatically. With this in mind, it is reasonable to revisit some of the foundational assumptions. The two-point model breaks the clustering signal down in to a sum of two parts: a ‘two-halo’ term that originates from the clustering between different pairs of haloes, and a ‘one-halo’ term that originates from clustering within the same halo. For all cosmological observables, the two-halo component of the clustering dominates the signal at large scales, while the one-halo component dominates at small scales. One perennial problem with the model has been the ‘transition’ region, where both two- and one-halo terms have similar magnitude and both contribute to the predicted signal. In general, the model underpredicts the strength of clustering in this region, for example Fedeli et al. (2014) and Mead et al. (2015) showed that this underprediction can be as much as 30 per cent for the matter–matter power spectrum, with the exact amount depending on redshift and cosmology. Tinker et al. (2005) also note similar problems in halo-model descriptions of the transition region of the galaxy–galaxy correlation function. This region of the spectrum is also called the ‘quasi-linear’ regime because the evolution of perturbations at these scales is not exactly governed by linear perturbation theory, but there is hope that it can be understood via higher order perturbative schemes. The fact that the transition and quasi-linear regions coincide is no coincidence as it is inevitable that halo formation is linked to scales where perturbative descriptions break down. In the special case of the matter–matter power spectrum, the inaccuracies of the halo model are ‘remedied’ by devising fitting functions (e.g. halofit of Smith et al. 2003; Takahashi et al. 2012) or else by adding phenomenological parameters (e.g. hmcode Mead et al. 2015, 2016, 2021) and fitting these to power spectra measured from high-resolution cosmological N-body simulations. In this way, models of the power spectrum that are accurate at the 5 per cent level for z < 2 and |$k\lt 10\, h\mathrm{Mpc}^{-1}$| have been developed. Alternatively, Valageas & Nishimichi (2011), Mohammed & Seljak (2014), Seljak & Vlah (2015), and Philcox, Spergel & Villaescusa-Navarro (2020) have all improved predictions at quasi-linear scales by incorporating perturbation theory for the matter field, with per cent level accuracy for |$k\lt 1\, h\mathrm{Mpc}^{-1}$| reported for some models. However, it could be argued that predictions of the matter–matter power spectrum are one of the least useful applications of the halo model, given that this spectrum can be accurately measured from N-body simulations (at least, if one ignores the inconvenient issue of baryonic feedback). Indeed, the most accurate predictors for the matter–matter spectrum as a function of cosmology for |$k\lt 10\, h\mathrm{Mpc}^{-1}$| come from ‘emulation’ schemes (e.g. Lawrence et al. 2010; Agarwal et al. 2012, 2014; Lawrence et al. 2017; Knabenhans et al. 2019) in which measured power spectra from simulations that span a range of cosmologies are interpolated between. The models discussed for the matter–matter power in the previous paragraph are difficult to generalize to spectra other than matter–matter. For example, if one is interested in the halo-model prediction for anything connected with galaxies, including the connection between matter and galaxies, then one is reduced to using the standard halo model, which comes with the standard problems. In this paper, we are interested in improving the halo-model predictions around the transition region for any field pair, and we do not restrict our focus to the matter–matter spectrum. We focus our attention on the non-linear portion of the halo bias, a treatment of which is almost always absent from standard halo-model calculations. We aim to do this while doing minimal damage to the existing halo-model apparatus because this is widely used by the community in its standard form. In Section 2, we present the standard derivation of the halo-model power spectrum and we demonstrate how to include non-linear halo bias in the calculation in a clean, isolated way. In Section 3, we show how we calculate the new non-linear halo bias term, which involves measurements of this new term from N-body simulations. In Section 4, we present the results of including this new term in halo-model calculations for matter–matter, halo–matter, and halo–halo power spectra, and we show an application to power spectra involving galaxies. Finally, we conclude in Section 5 and propose some ideas for future work. Appendix A discusses the technical details of how we incorporate the new non-linear bias term in our numerical calculations. Appendix B presents results at redshifts other than those presented in the main body of the paper. Appendix C presents a pedagogical example mock universe where linear biasing is exactly respected, and we demonstrate that this can be modelled essentially perfectly by the standard halo-model calculation. 2 HALO MODEL 2.1 Standard derivation We shall first present the standard derivation of the halo model power spectrum (e.g. Cooray & Sheth 2002). This subsection can be skipped by those readers familiar with the halo model, but is included for completeness. We carry out this derivation in a comoving, periodic volume V because we believe this is the most transparent. Note that this means that the Fourier modes are discrete. We eventually take V → ∞ to retrieve the continuum limit. To keep the notation simple, we suppress any time dependence in function arguments, but the reader should recall that most functions can be time as well as (Fourier) space dependent. We consider the power spectrum between a pair of three-dimensional cosmological ‘fields’, which could be identical. Examples of such fields would be ‘matter’, ‘halo’, or ‘galaxy’ overdensities that vary from place-to-place in the Universe, and that could require smoothing to be reasonably defined. These fields could be further restricted to be haloes in a specific mass range or galaxies of a specific type. They could also be fields that are generated by some emissive process, such as infrared light or electron pressure. We define fields |$\theta _u(\boldsymbol {x})$|⁠, where |$\boldsymbol {x}$| is comoving position and the label u stands for the field. We will be interested in power spectra between pairs of such fields, Puv(k), which has units of the product of the units of fields u and v and an additional factor of volume. We start from the assumption that all fields reside in haloes, which allows us to write the total field at some position as a sum of contributions, θu, i, from haloes with positions |$\boldsymbol {x}_i$| $$\begin{eqnarray*} \theta _u(\boldsymbol {x})=\sum _i \theta _{u,i}(\boldsymbol {x}-\boldsymbol {x}_i), \end{eqnarray*}$$(1) if we allow our definition of ‘halo’ to be sufficiently general then this is always true: we can always break down the total field into a sum of contributions. For example, θu, i could be the ‘halo profile’ for matter, but for other fields it could be an emission profile. In what follows we refer to θu, i as the halo profile, but one should keep in mind that this is a more general quantity. We can write the Fourier transform, in terms of comoving wavenumber |$\boldsymbol {k}$|⁠, as $$\begin{eqnarray*} \theta _{u,\boldsymbol {k}}=\frac{1}{V}\int {\theta _u(\boldsymbol {x})}\mathrm{ e}^{-i\boldsymbol {\boldsymbol {k}}\cdot \boldsymbol {\boldsymbol {x}}}\, \mathrm{d}^3{\boldsymbol {x}}= \frac{1}{V}\sum _i\int {\theta _{u,i}(\boldsymbol {x}-\boldsymbol {x}_i)}e^{-i\boldsymbol {\boldsymbol {k}}\cdot \boldsymbol {\boldsymbol {x}}}\, \mathrm{d}^3{\boldsymbol {x}}, \nonumber\\ \end{eqnarray*}$$(2) where the factor of 1/V ensures that the units of |$\theta _{u,\mathbf {k}}$| are the same of those of |$\theta _u(\boldsymbol {x})$|⁠. If we make the variable change |$\boldsymbol {y}=\boldsymbol {x}-\boldsymbol {x}_i$|⁠, then we can write $$\begin{eqnarray*} \theta _{u,\mathbf {k}}=\frac{1}{V}\sum _i\mathrm{e}^{-i\boldsymbol {k}\cdot \boldsymbol {x_i}}\int {\theta _{u,i}(\boldsymbol {y})}e^{-i\boldsymbol {\boldsymbol {k}}\cdot \boldsymbol {\boldsymbol {y}}}\, \mathrm{d}^3{\boldsymbol {y}}= \frac{1}{V}\sum _i\mathrm{e}^{-i\boldsymbol {k}\cdot \mathbf {x_i}}W_{u,i}(k), \nonumber\\ \end{eqnarray*}$$(3) where we have defined $$\begin{eqnarray*} W_{u,i}(k)=\int _0^\infty \frac{\sin (kr)}{kr}\theta _{u,i}(r)4\pi r^2\, \mathrm{d}r, \end{eqnarray*}$$(4) as a Fourier transform of the halo profile, which has units of field u multiplied by volume. We have made the additional assumption that the halo profile is spherically symmetric, which allows us to write Wu, i in terms of |$k=|\boldsymbol {k}|$| and |$r=|\boldsymbol {r}|$| and also ensures that Wu, i is real. Let us now construct the power spectrum between two fields (which could be identical), θu and θv: $$\begin{eqnarray*} P_{uv}(k)=\left\langle \theta ^{*}_{u,\mathbf {k}}\theta ^{}_{v,\boldsymbol {k}}\right\rangle =\left\langle \frac{1}{V^2}\sum _i\sum _j \mathrm{e}^{-i\boldsymbol {k}\cdot \boldsymbol {(x_i-x_j)}}W_{u,i}(k)W_{v,j}(k)\right\rangle.\nonumber\\ \end{eqnarray*}$$(5) We can break this equation up into two pieces, the i = j piece, which we call the one-halo term and the i ≠ j piece, which we call the two-halo term. The one-halo term is $$\begin{eqnarray*} P^\mathrm{1H}_{uv}(k)=\left\langle \frac{1}{V^2}\sum _i W_{u,i}(k)W_{v,i}(k)\right\rangle , \end{eqnarray*}$$(6) where the expectation value is taken over all modes in the volume and we have made the standard assumptions about isotropy and homogeneity that ensure that we can write expressions as a function of k only. We now assume a continuum of haloes, labelled with a mass M, and that the halo mass is the sole determinant of the halo properties. We further assume that these haloes are distributed according to a mass-distribution function n(M) (sometimes denoted dn/dM in the literature), where |$n(M)\, \mathrm{d}M$| describes the number density of haloes in the range M to M + dM. We can then write the sum as an integral by taking V → ∞ using $$\begin{eqnarray*} \frac{1}{V}\sum _i \rightarrow \int n(M)\, \mathrm{d}M, \end{eqnarray*}$$(7) to get $$\begin{eqnarray*} P^\mathrm{1H}_{uv}(k)=\int _0^\infty W_{u}(M,k)W_{v}(M,k) n(M)\, \mathrm{d}M. \end{eqnarray*}$$(8) We can apply the same reasoning to the two-halo term, where we convert both sums to integrals $$\begin{eqnarray*} P^\mathrm{2H}_{uv}(k)&=&\int _0^\infty \int _0^\infty P_\mathrm{hh}(M_1,M_2,k)\nonumber\\ &&\times\, W_{u}(M_1,k)W_{v}(M_2,k)n(M_1)n(M_2)\, \mathrm{d}M_1\mathrm{d}M_2, \end{eqnarray*}$$(9) and we recognize the expectation of the complex exponential to be the power spectrum of the halo centres: Phh(M1, M2, k). Note that the functions Wu and Wv are common between equations (8) and (9). At this stage of the derivation, it is common to make the approximation that haloes are linearly biased tracers of the underlying linear matter field $$\begin{eqnarray*} P_\mathrm{hh}(M_1,M_2,k)\simeq b(M_1)b(M_2)P^\mathrm{lin}_\mathrm{mm}(k), \end{eqnarray*}$$(10) where b(M) is the linear halo bias of haloes with mass M and |$P^\mathrm{lin}_\mathrm{mm}(k)$| is the linear-theory matter power spectrum. This approximation has the virtue of being correct at very large scales and also is mathematically convenient because it allows us to split the double integral in equation (9) into a product of two one-dimensional integrals of similar form. The final result is then $$\begin{eqnarray*} P^{2\mathrm{H}}_{uv}(k)=P^\mathrm{lin}_\mathrm{mm}(k) \prod _{n=u,v}\left[\int _0^\infty W_n(M,k)b(M)n(M)\, \mathrm{d}M\right]. \end{eqnarray*}$$(11) The adopted halo mass function and linear halo bias must satisfy the following properties for any power spectrum involving the matter to have the correct large-scale limit1: $$\begin{eqnarray*} \int _0^\infty Mn(M)\, \mathrm{d}M=\bar{\rho }, \end{eqnarray*}$$(12) $$\begin{eqnarray*} \int _0^\infty Mb(M)n(M)\, \mathrm{d}M=\bar{\rho }, \end{eqnarray*}$$(13) where |$\bar{\rho }$| is the mean comoving cosmological matter density. In words, these equations enforce that all matter is contained in haloes and that, on average, matter is unbiased with respect to itself. As an example of how the two-halo term works in practice: In the special case of the matter–matter overdensity power spectrum, we can write |$\theta _\mathrm{m}(M,r) = \rho _\mathrm{m}(M,r)/\bar{\rho }$| where ρm(M, r) is the halo matter density profile. We then note that |$W_\mathrm{m}(M,k\rightarrow 0)=M/\bar{\rho }$|⁠, and it is usual to factorize this normalization, such that |$W_\mathrm{m}(M,k)=MU_\mathrm{m}(M,k)/\bar{\rho }$|⁠, with Um(M, k → 0) = 1. We can then write equation (11) as $$\begin{eqnarray*} P^{2\mathrm{H}}_{uv}(k)=P^\mathrm{lin}_\mathrm{mm}(k) \left[\frac{1}{\bar{\rho }}\int _0^\infty M U_\mathrm{m}(M,k)b(M)n(M)\, \mathrm{d}M\right]^2, \end{eqnarray*}$$(14) and we see that |$P^{2\mathrm{H}}_{uv}(k\rightarrow 0)=P^\mathrm{lin}_\mathrm{mm}(k\rightarrow 0)$| automatically as the term in the square brackets equals unity in this limit. For spectra other than matter–matter this is no longer true, and in general the large-scale limit of the two-halo term will be equal to the linear spectrum multiplied by amplitude factors that account for the field content and the bias that arises from how these fields populate haloes. 2.2 Including non-linear halo bias It is worth examining some of the approximations that lead to the standard halo-model equations (8) and (11): It has been assumed that halo profiles are perfectly spherical with no substructure, that there is no scatter in profile properties at fixed host halo mass, and that halo properties depend only on halo mass and not on other properties, for example, halo location. These approximations will break down, and the errors in the eventual power spectrum that they contribute will vary with the fields that are being considered (e.g. Smith & Watts 2005; Giocoli et al. 2010; Smith & Markovic 2011; Chen & Afshordi 2020; Voivodic, Rubira & Lima 2020). It is also usually assumed that haloes trace the underlying linear matter distribution with a linear halo bias. In this paper, we focus on including scale-dependent, non-linear halo bias within the halo-model calculation. A previous attempt to include non-linear halo bias has been made by Smith, Scoccimarro & Sheth (2007) where the combination of standard perturbation theory (SPT) and an Eulerian bias expansion were used to model the matter–matter, matter–halo, and halo–halo power spectra. This approach was demonstrated to be successful at very large scales, where perturbation theory is a good description of clustering. A similar model is presented by McDonald (2006) where the renormalization of coefficients in the bias expansion was considered for the first time. Both of these models rely on perturbation theory, and they fail on smaller scales where perturbation theory breaks down and where much of the constraining power of a contemporary cosmological survey lies. Alternatively, Fedeli et al. (2014) investigate a phenomenological non-linear biasing model where b(M) → b(M, k) in equation (11) and the k dependence is fitted to N-body data. In this paper, we follow a different approach, and use measurements of the non-linear halo bias from simulations to push more deeply in to the non-linear regime. Our study is of particular relevance to galaxy–galaxy lensing, which is the study of the two-point function between galaxy and matter overdensities, where the matter clustering is accessed via weak gravitational lensing. The halo model is very widely used in the interpretation of data from such observations, and is used to understand the fundamental link between haloes and galaxy formation therein (e.g. Mandelbaum et al. 2005; Cacciato et al. 2009). Most halo models of galaxy–galaxy lensing assume a linear halo bias (e.g. Cacciato et al. 2012; Dvornik et al. 2018), but the halo model of the galaxy–galaxy lensing signal developed by van den Bosch et al. (2013) includes some non-linear effects of halo biasing in two ways. First, the idea that haloes cannot overlap (so-called halo exclusion) is considered by forcing the halo–halo correlation function to −1 on scales below the sum of the virial radii of the haloes contributing to the two-point function. Secondly, scale-dependent halo bias is included using a fitting function for the ‘radial-dependent’ bias taken from Tinker et al. (2005). This has the advantage that the Tinker et al. (2005) result is returned automatically on quasi-liner scales, but the disadvantage that the applicability of results is limited to a specific galaxy population. In the van den Bosch et al. (2013) model, and in some other galaxy–galaxy lensing prescriptions, the halo bias is defined to be relative to the non-linear matter field. In our work, we avoid this, since in principle a good halo model should be able to predict the non-linear power spectra of matter–matter, matter–galaxies, and galaxies–galaxies all from the same set of founding assumptions. Different assumptions regarding non-linear halo bias impact on halo-model predictions in the transition region between the two- and one-halo terms, an accurate modelling of which is becoming important given the quality of contemporary data; traditional halo models generally underpredict either power or correlation in the transition region. In Hayashi & White (2008), the halo–matter cross-correlation is modelled by adopting the maximum of either the two- or one-halo term in the transition region. In Garcia et al. (2020), it is demonstrated that the transition region can be well modelled if halo boundaries are tailored to the cross-correlation and if a halo exclusion term is included. In the Koukoufilippas et al. (2020) measurement of the the thermal Sunyaev–Zel’dovich (tSZ)–galaxy cross-correlation the transition region was scaled by the ratio of halofit to the standard halo model prediction for the matter–matter power spectrum. In Hang et al. (2021), the galaxy auto correlation is assumed to be a constant bias multiplying a linear part, plus a different constant bias multiplying a non-linear part, which is taken from halofit. None of these solutions are particularly appealing from a theoretical perspective. We note that it is not obvious that the same correction that works in halofit for matter–matter would apply more generally, and indeed later we show that this is not the case. In the standard halo-model derivation, the calculation of the two-halo term is made tractable by making the approximation that the haloes are linearly biased tracers of the underlying linear matter field (equation 10). At large scales, the shape of the two-halo term of any power spectrum computed in this way will be exactly that of the linear spectrum; scale-dependent deviations from this arise due to the factors of Wn(M, k) in equation (11). This means that the shape of the standard two-halo term is only different from linear theory on scales corresponding to the virial radii of the most massive haloes where Wn(M, k) starts to deviate from constant. We note that it seems unlikely that a linear halo bias is a good description of the clustering relation between haloes and matter on scales comparable to the sizes of individual haloes. Given that we know that in reality halo bias is not linear, let us instead not make this approximation, but write the halo power spectrum in the following way: $$\begin{eqnarray*} P_\mathrm{hh}(M_1,M_2,k)=b(M_1)b(M_2)P^\mathrm{lin}_\mathrm{mm}(k)[1+\beta ^\mathrm{NL}(M_1,M_2,k)], \end{eqnarray*}$$(15) where the function βNL captures all the things missing from the standard linear-bias–linear-field model. We know some things about βNL, specifically the large-scale limit: βNL(M1, M2, k → 0) = 0, and also that it must obey symmetry with respect to mass arguments: βNL(M1, M2, k) = βNL(M2, M1, k). The new idea presented in this paper is to include βNL within semi-analytical calculations using the halo model. If we substitute equation (15) into equation (9) we have $$\begin{eqnarray*} P^{2\mathrm{H}}_{uv}(k)&=&P^\mathrm{lin}_\mathrm{mm}(k)\prod _{n=u,v}\left[\int _0^\infty W_n(M,k)b(M)n(M)\mathrm{d}M\right] \nonumber\\ &&+ \, P^\mathrm{lin}_\mathrm{mm}(k)I^\mathrm{NL}_{uv}(k). \end{eqnarray*}$$(16) The first term is standard, and is identical to equation (11) and the new content is captured by |$I^\mathrm{NL}_{uv}(k)$| in the second term: $$\begin{eqnarray*} I^{\mathrm{NL}}_{uv}(k)&=&\int _0^\infty \int _0^\infty \beta ^\mathrm{NL}(M_1,M_2,k) W_u(M_1,k)\nonumber\\ &&\times\, W_v(M_2,k)b(M_1)b(M_2)n(M_1)n(M_2)\, \mathrm{d}M_1\mathrm{d}M_2. \end{eqnarray*}$$(17) It is worth considering what this new content represents. The standard approximation is that haloes are linearly biased tracers of an underlying linear matter field. Since we know the linear matter field to be a Gaussian random field this implies that the halo fields themselves must be Gaussian random. The new content is all departures from this simple picture, including enhanced halo clustering, filamentary, sheet, and void structure; all of which can be realized by moving haloes from their linear Gaussian locations. This is often termed ‘scale-dependent’ bias. There is no reason to assume that the function βNL should be particularly simple, indeed, given the complexity of the new content we might expect it to be a complicated function. We also have no reason to believe that it should be independent of either redshift or cosmological parameters. Since we know perturbation theory to be a good description of the Universe for |$k\lesssim 0.2\, h\mathrm{Mpc}^{-1}$|⁠, the βNL function must also contain these perturbative corrections if it is to be a good, general model of clustering. 2.3 Relation to other non-linear bias definitions Non-linear halo bias is often discussed in terms of the scale-dependent function, bhh(M1, M2, k), defined by $$\begin{eqnarray*} P_\mathrm{hh}(M_1,M_2,k)=b^2_\mathrm{hh}(M_1,M_2,k)P^\mathrm{lin}_\mathrm{mm}(k), \end{eqnarray*}$$(18) where $$\begin{eqnarray*} b^2_\mathrm{hh}(M_1,M_2,k\rightarrow 0)=b(M_1)b(M_2). \end{eqnarray*}$$(19) This is can be computed using our formalism (via equation 15) as $$\begin{eqnarray*} b^2_\mathrm{hh}(M_1,M_2,k)=b(M_1)b(M_2)[1+\beta ^\mathrm{NL}(M_1,M_2,k)], \end{eqnarray*}$$(20) if M1 = M2 there is an additional one-halo contribution (in this case called the shot noise) of |$1/\bar{n}_\mathrm{h}$| where |$\bar{n}_\mathrm{h}$| is the mean number density of the halo sample. Some authors define a halo bias via the cross-spectrum between the halo–number density field and the matter overdensity field: $$\begin{eqnarray*} P_\mathrm{hm}(M,k)=b_\mathrm{hm}(M,k)P^\mathrm{lin}_\mathrm{mm}(k). \end{eqnarray*}$$(21) This halo bias will be generally different from that measured through the halo–number density auto spectrum, although they coincide at large scales. This ‘cross-bias’ can be written in our notation (setting u = h and v = m) as $$\begin{eqnarray*} b_\mathrm{hm}(M,k) &=& W_\mathrm{m}(M,k)/P^\mathrm{lin}_\mathrm{mm}(k)+b(M)\int _0^\infty \mathrm{d}M_2\nonumber\\ &&\times\,[1+\beta ^\mathrm{NL}(M,M_2,k)]W_\mathrm{m}(M_2,k)b(M_2)n(M_2). \end{eqnarray*}$$(22) The first term on the right-hand side of equation (22) arises from the one-halo term (equation 8), while the second is from the two-halo term (equation 16). The equations in this subsection can be derived by taking |$W_\mathrm{h}(M^{\prime },k)=\delta _\mathrm{D}(M-M^{\prime })/\bar{n}_\mathrm{h}$| as the window function for the halo therefore in both equations (20) and (22) we have assumed that we are considering thin halo mass bins. Some authors define the halo biases in equation (18) or (21) with respect to the non-linear matter–matter power spectrum, rather than the linear matter spectrum. In this paper, we always define it with respect to the linear spectrum and this distinction is important in our work because: First, it is halo bias defined in this way that enters standard halo-model calculations, and secondly a general model should be able to predict the non-linear matter–matter spectrum, and since all matter is contained in haloes this itself must come from the haloes. If one wanted to work with halo bias defined relative to the underlying non-linear matter–matter spectrum, then this itself is computable from our model by taking matter overdensity profiles in equations (16) and (17). 3 MEASURING NON-LINEAR HALO BIAS 3.1 Measurement It would be ideal if we could calculate βNL from first principles. However, while perturbation theory may be able to give us some insight at large scales, we anticipate βNL to have structure on scales that are out of reach of even state-of-the-art perturbation theories. Therefore, we decide to measure the required function from N-body simulations. If we simply rewrite equation (15) as $$\begin{eqnarray*} 1+\beta ^\mathrm{NL}(M_1,M_2,k)=\frac{P_\mathrm{hh}(M_1,M_2,k)}{b(M_1)b(M_2)P^\mathrm{lin}_\mathrm{mm}(k)}, \end{eqnarray*}$$(23) we can see that a sensible way to measure βNL is to measure Phh(M1, M2, k) and then to divide it by the linear power spectrum multiplied by the product of linear halo bias factors for M1 and M2. Recall that Phh(M1, M2, k) is the cross-power spectrum measured between haloes of mass M1 with those of mass M2. In real measurements, we have to bin haloes in mass to measure this function. If we consider the autospectrum (M1 = M2), then we will also have to subtract (halo) shot noise from the measured halo power spectra because, in our approach, this is the one-halo contribution to the halo autospectrum, and we are interested only in the two-halo contribution. This shot noise is not a problem if we consider the cross-spectrum between haloes in two different mass bins because it is automatically zero when the haloes in each leg of the cross-spectrum are different. The shot-noise contribution to a given autohalo-power-spectrum measurement is given by $$\begin{eqnarray*} P^\mathrm{SN}_\mathrm{hh}=\frac{1}{\bar{n}_\mathrm{h}}, \end{eqnarray*}$$(24) where |$\bar{n}_\mathrm{h}$| is the mean halo number density measured for the mass bin in question. Some authors consider a non-Poissonian shot noise, not given by equation (24) and may also consider halo exclusion (the spatially exclusivity of haloes) to affect the shot-noise term. In our picture, halo exclusion will enter in the non-linear halo bias portion of the two-halo term, since it pertains to the way that haloes trace the underlying matter field, rather than the structure within the haloes themselves. In fact, regardless of the position one takes on a shot-noise correction for haloes, subtracting power according to equation (24) is correct in our work given that this is exactly what we take for the one-halo term when we evaluate halo–halo autopower spectra using our model (equation 8). Subtracting equation (24) can therefore be seen as a method for isolating the two-halo term from the halo–halo auto power measurement. 3.2 Simulations We measure βNL(M1, M2, k) using data taken from the multidark simulation data base2 (Klypin, Trujillo-Gomez & Primack 2011; Prada et al. 2012; Riebe et al. 2013). We use data from the original multidark simulation, which simulated N = 20483 particles in a |$L=1000\, h^{-1}\mathrm{Mpc}$| cube. We consider the combination of a high number of particles in the large volume of multidark advantageous for our measurement at both large and small scales. We utilize haloes that have been identified via the rockstar phase-space algorithm of Behroozi, Wechsler & Wu (2013)3 and defined using the ‘virial’ criterion, with a spherical-overdensity threshold given by |$\simeq 360\bar{\rho }$| at z = 0, which comes from spherical-collapse calculations (e.g. Bryan & Norman 1998) for the simulated WMAP 5 level accurate for standard (ΛCDM) cosmology: Ωm = 0.27, Ωb = 0.0469, Ωv = 1 − Ωm, h = 0.7, n = 0.95, and σ8 = 0.82. We prefer to use the virial halo definition of haloes because of results presented in Mead et al. (2021), where it was demonstrated that halo-model calculations are more robust with respect to cosmological dependence when haloes are defined using the virial condition (as opposed to |$200\bar{\rho }$| or 200ρc). However, in principle one could work with any desired halo definition as long as consistency is maintained.4 To measure βNL (equation 23) we calculate halo–halo cross-power between eight mass bins, which leads to 36 unique cross-combinations. The choice of eight mass bins represents a reasonable compromise between having enough bins so that we gain in halo mass resolution, while also allowing each bin to contain enough haloes such that the measurement is not too noisy.5 The lower limit of our lowest mass bin corresponds to haloes with 50 particles, which in multidark is |$\simeq 10^{11.6}\, h^{-1}\mathrm{M_\odot }$|⁠. In this paper, we only care about the halo position and mass being accurately measured, and 50 particles can be considered a minimum for this (Knebe et al. 2011). We compared the mass function of our halo sample with theoretical expectations from Tinker et al. (2010) and find good agreement, even at our 50 particle lower limit.6 The limits of our mass bins are defined to be equally spaced in peak height, ν, defined via $$\begin{eqnarray*} \nu =\frac{\delta _\mathrm{c}}{\sigma (M)}, \end{eqnarray*}$$(25) between the ν value corresponding to the lowest mass haloes (ν ≃ 0.74 at z = 0) and that corresponding to the most massive halo (ν ≃ 3.97 at z = 0). We prefer to work with ν bins, rather than M or log (M), because many quantities of cosmological interest, such as the halo mass function, bias and concentration–mass relation, have been shown to be more independent of further cosmology dependence when expressed in terms of ν. While we do not investigate the cosmology dependence of βNL in this paper, we feel that expressing it in terms of ν may be useful in the future. In equation (25), we take the (cosmology-dependent) critical threshold for collapse, δc from Nakamura & Suto (1997), although this cosmology dependence has a negligible impact (see Mead et al. 2021) on the eventual results in this paper. The standard deviation in the linear matter field, σ(M), is calculated analytically in the usual way using a top-hat filter. 3.3 Constructing βNL We compute Phh(M1, M2, k) via fast Fourier transform with 5123 cells and subtract halo shot noise for auto spectra as per equation (24) in order to isolate the two-halo component. We only keep wavenumbers up to half of the Nyquist frequency (for our number of mesh cells |$k_\mathrm{Ny}/2\simeq 0.8\, h\mathrm{Mpc}^{-1}$|⁠) in order to avoid the effects of aliasing. Following equation (23), we construct βNL(M1, M2, k) by dividing the measured halo–halo power spectra by the linear power spectrum. For |$k\lt 0.08\, h\mathrm{Mpc}^{-1}$|⁠, we take the linear matter–matter spectrum to be that measured from the simulation itself.7 At large scales, this allows us to cancel some cosmic variance that would otherwise dominate our measurement. We are unable to use the simulation measurement for |$k\gt 0.08\, h\mathrm{Mpc}^{-1}$| because non-linear contributions to the measured matter–matter field become important, so we use a smooth theory calculation from camb (Lewis, Challinor & Lasenby 2000). In practice, this is not a problem at these smaller scales because the measurement is not noise dominated, and in any case the noise would not benefit from cosmic-variance cancellation since it is not Gaussian. We then also divide by the product of the linear halo bias factors, which we take from fitting a linear bias model to the |$k\lt 0.08\, h\mathrm{Mpc}^{-1}$| portion of the measurement of |$P_\mathrm{hh}(M_1, M_2, k)/P^\mathrm{lin}_\mathrm{mm}(k)$|⁠. We also investigated taking the linear bias from the Tinker et al. (2010) fitting function, and found similar results.8 However, we noticed discrepancies if we take the linear bias from a less accurate source, for example, the peak-background split model of Tinker et al. (2010). In Fig. 1, we show the βNL function measured from multidark for four different mass bins at z = 0. The measurement is noisy at large scales, but appears to asymptote to zero. Structure in βNL becomes visible for |$k\gt 0.08\, h\mathrm{Mpc}^{-1}$| and we observe a prominent, positive detection of the function for |$k\gt 0.1\, h\mathrm{Mpc}^{-1}$| with an amplitude that is dependent on the mass bin. At smaller scales still, the function seems to decay, particularly for the higher mass bins, which we suspect may be because the spatially exclusivity of haloes ensures that the correlation function must be −1 at scales smaller than the sum of the virial radii of the two-halo populations being correlated. This halo-exclusion condition for small scales in the correlation function will translate in to a small-scale depletion in power for wavenumbers greater than that corresponding to the sum of the virial radii. The two highest mass bins shown in Fig. 1 correspond to rv ≃ 0.94 and |$1.26\, h^{-1}\mathrm{Mpc}$|⁠, respectively, which correspond to the dips visible in the right-hand panel to within factors of |$\sim \pi$|⁠. Figure 1. Open in new tabDownload slide Non-linear halo bias function βNL(M1, M2, k) measured from the multidark simulation at z = 0. We show this function measured for autocombination and cross-combination of four halo-mass bins centred at 1012.2, 1013.2, 1013.8, and |$10^{14.3}\, h^{-1}\mathrm{M_\odot }$|⁠. The left-hand panel shows the four autospectra, while the other panels show the six possible cross-spectra, with each point/error set coloured according to the bin colours denoted in the left-hand panel. The dashed vertical line |$0.08\, h\mathrm{Mpc}^{-1}$| indicates our split between those modes we take to be linear and those we take to be non-linear. Error bars, shown only for |$k\gt 0.08\, h\mathrm{Mpc}^{-1}$|⁠, show error-on-the-mean power in each k bin in the halo–halo power spectra measured from the simulation, the errors are not shown for |$k\lt 0.08\, h\mathrm{Mpc}^{-1}$| where we enjoy some cosmic variance cancellation because here we divide by the measured linear spectrum. At low k, we see that βNL tends to zero, as per our expectation. We see non-zero structure start to emerge in βNL for |$k\gtrsim 0.08\, h\mathrm{Mpc}^{-1}$|⁠. For |$k\sim 0.8\, h\mathrm{Mpc}^{-1}$|⁠, we see the function approaches unity for some halo-mass bins, which indicates that this will provide order-unity corrections to analytical calculations around these scales. That the function turns over and starts to decrease at small scales for spectra that involve the higher halo masses is probably because of halo exclusion, which limits how close two haloes can physically be, thus killing the power below scales that correspond to the sum of the two virial radii. Figure 1. Open in new tabDownload slide Non-linear halo bias function βNL(M1, M2, k) measured from the multidark simulation at z = 0. We show this function measured for autocombination and cross-combination of four halo-mass bins centred at 1012.2, 1013.2, 1013.8, and |$10^{14.3}\, h^{-1}\mathrm{M_\odot }$|⁠. The left-hand panel shows the four autospectra, while the other panels show the six possible cross-spectra, with each point/error set coloured according to the bin colours denoted in the left-hand panel. The dashed vertical line |$0.08\, h\mathrm{Mpc}^{-1}$| indicates our split between those modes we take to be linear and those we take to be non-linear. Error bars, shown only for |$k\gt 0.08\, h\mathrm{Mpc}^{-1}$|⁠, show error-on-the-mean power in each k bin in the halo–halo power spectra measured from the simulation, the errors are not shown for |$k\lt 0.08\, h\mathrm{Mpc}^{-1}$| where we enjoy some cosmic variance cancellation because here we divide by the measured linear spectrum. At low k, we see that βNL tends to zero, as per our expectation. We see non-zero structure start to emerge in βNL for |$k\gtrsim 0.08\, h\mathrm{Mpc}^{-1}$|⁠. For |$k\sim 0.8\, h\mathrm{Mpc}^{-1}$|⁠, we see the function approaches unity for some halo-mass bins, which indicates that this will provide order-unity corrections to analytical calculations around these scales. That the function turns over and starts to decrease at small scales for spectra that involve the higher halo masses is probably because of halo exclusion, which limits how close two haloes can physically be, thus killing the power below scales that correspond to the sum of the two virial radii. In the top row of Fig. 2, we show βNL at z = 0 as a function of two mass variables as measured from multidark at k = 0.1, 0.3, and |$0.8\, h\mathrm{Mpc}^{-1}$| separately. We parametrize the function in terms of the ‘peak height’, ν, rather than M (equation 25) because this pertains directly to our numerical implementation (Section 4). For reference, ν = 0.5, 1, 2, and 3 correspond to ≃1010.4, 1012.5, 1014.2, and |$10^{14.9}\, h^{-1}\mathrm{M_\odot }$| haloes for this cosmology at z = 0. We see that, in general, the function increases in amplitude as k increases, as can also be seen in Fig. 1. We also see the expected reflection symmetry about the ν1 = ν2 line. We measured this function from multidark for ν ≳ 0.75, and for lower values of ν we linearly extrapolate to get an estimate of βNL (we discuss the validity of this later). We are not able to measure βNL for ν ≲ 0.75 because this corresponds to haloes below our 50 particle limit. Note that we only have 25 bins in k and 8 bins in ν for the measurement, and the rest of βNL is calculated either via linear interpolation or extrapolation in three dimensions. The locations of high signal in Fig. 2 indicate cross-spectra that have particularly strong non-linear halo biases. We note particular intensity in the cross-between very low ν ∼ 0.4 and high ν ∼ 2.5 haloes. This plausibly originates from low-mass haloes falling on to higher mass haloes. We also see a negative signal between ν ∼ 2 and ν ∼ 3 haloes, which could be a result of high-mass haloes having formed through major mergers leaving a deficit of the haloes that needed to have merged. Figure 2. Open in new tabDownload slide The upper set of three panels show the function βNL(ν1, ν2, k) at k = 0.2 (left), 0.4, and |$0.8\, h\mathrm{Mpc}^{-1}$| (right), all at z = 0. The lower three panels show the integrand for the function |$I^\mathrm{NL}_\mathrm{mm}(k)$|⁠, defined in equation (17). This integrand is βNL(ν1, ν2, k)Um(ν1, k)Um(ν2, k)b(ν1)b(ν2)g(ν1)g(ν2), and it is shown here for the case of the matter–matter power spectrum; for most of the masses and scales shown here Um(ν, k) ≃ 1. βNL is measured directly from the multidark simulation in ν bins and we linearly interpolate between these bins to acquire a continuous function, which leaves some square residual patterns, most visible in the top left-hand panel. The lowest mass haloes that we measure correspond to ν ≃ 0.75, which corresponds to the dashed vertical and horizontal lines in each panel. To get values for βNL outside this limit, we linearly extrapolate from our measurements at higher halo masses. For reference, for this cosmology at this redshift ν = 0.5, 1, 2, and 3 correspond to halo masses: |$\log _{10}(M/\, h^{-1}\mathrm{M_\odot })\simeq 10.4$|⁠, 12.5, 14.2, and 14.9. Figure 2. Open in new tabDownload slide The upper set of three panels show the function βNL(ν1, ν2, k) at k = 0.2 (left), 0.4, and |$0.8\, h\mathrm{Mpc}^{-1}$| (right), all at z = 0. The lower three panels show the integrand for the function |$I^\mathrm{NL}_\mathrm{mm}(k)$|⁠, defined in equation (17). This integrand is βNL(ν1, ν2, k)Um(ν1, k)Um(ν2, k)b(ν1)b(ν2)g(ν1)g(ν2), and it is shown here for the case of the matter–matter power spectrum; for most of the masses and scales shown here Um(ν, k) ≃ 1. βNL is measured directly from the multidark simulation in ν bins and we linearly interpolate between these bins to acquire a continuous function, which leaves some square residual patterns, most visible in the top left-hand panel. The lowest mass haloes that we measure correspond to ν ≃ 0.75, which corresponds to the dashed vertical and horizontal lines in each panel. To get values for βNL outside this limit, we linearly extrapolate from our measurements at higher halo masses. For reference, for this cosmology at this redshift ν = 0.5, 1, 2, and 3 correspond to halo masses: |$\log _{10}(M/\, h^{-1}\mathrm{M_\odot })\simeq 10.4$|⁠, 12.5, 14.2, and 14.9. 4 RESULTS 4.1 Calculation detail The main results of this paper are halo-model calculations for the matter–matter, matter–halo, and halo–halo power spectra when the two-halo term is calculated using equation (16). From equation (17), we see that this requires us to evaluate a double integral over the βNL(M1, M2, k) function weighted by halo bias, mass function, and profile functions. Difficulty arises in the numerical integration when one of the fields is ‘matter’ because of substantial contributions from low-mass haloes, well below the mass typically resolved by N-body simulations, we discuss this in detail in Appendix A. In practice, we evaluate this integral by converting M to the dimensionless ‘peak height’, defined in equation (25). We then use the mass functions g(ν), normalized such that the integral over all ν gives unity, which is related to n(M) via $$\begin{eqnarray*} g(\nu)\, \mathrm{d}\nu =\frac{M}{\bar{\rho }}n(M)\, \mathrm{d}M. \end{eqnarray*}$$(26) For our halo-model calculations, we take the form of g(ν) and b(ν) from Tinker et al. (2010) and we take the overdensity threshold, Δv(z), to be that used for halo identification in the multidark simulations (taken from Bryan & Norman 1998: |$\simeq 360\bar{\rho }$| at z = 0, asymptoting to |$\simeq 178\bar{\rho }$| at high z). The mass function and halo bias are defined in such a way that equations (12) and (13) are automatically satisfied. For any power spectrum involving the matter field we must also choose a halo profile, and we adopt the profile of Navarro, Frenk & White (1997): $$\begin{eqnarray*} \rho (M,r)\propto \frac{1}{r/r_\mathrm{s}(1+r/r_\mathrm{s})^2}, \end{eqnarray*}$$(27) which is truncated at the virial radius defined such that this encloses an average density of Δv(z) times the mean density: $$\begin{eqnarray*} M=4\pi r_\mathrm{v}^3\Delta _\mathrm{v}(z)\bar{\rho }. \end{eqnarray*}$$(28) The virial radius is related to the halo-scale radius, rs via the mass-dependent concentration parameter c = rv/rs, which we take from Duffy et al. (2008) and use the appropriate redshift-dependent relation for their full sample of haloes identified using a virial criterion. We are primarily interested in intermediate ‘quasi-linear’ scales and neither the adopted halo profile nor specific concentration-mass relation are important for our results,9 which only start to have a significant impact for |$k\gtrsim 1\, h\mathrm{Mpc}^{-1}$|⁠. The lower row of Fig. 2 shows the |$I^\mathrm{NL}_\mathrm{mm}$| integrand, defined in equation (17), for the special case of the matter–matter power spectrum. For the scales and mass ranges shown, |$W_\mathrm{m}(M, k)\simeq M/\bar{\rho }$| and therefore the most significant change when going from βNL to the integrand for |$I^\mathrm{NL}_\mathrm{mm}$| is the suppression at high ν caused by the halo-mass function. Fig. 2 therefore shows the halo mass ranges that give additional contributions in our two-halo term from the non-linear halo bias. When evaluating the non-linear bias correction in equation (17), we force the correction to be zero for |$k\lt 0.08\, h\mathrm{Mpc}^{-1}$|⁠. This is consistent with our earlier choices because we measure our linear halo bias using these scales, and so we are making an implicit assumption that scale-dependent halo bias should be zero here. This choice makes only a small difference to our results because the correction, even when evaluated, is very small for |$k\lt 0.08\, h\mathrm{Mpc}^{-1}$|⁠. However, because βNL is noise dominated at these scales, if we do not force the correction to zero we transfer this noise into our halo-model calculation. In Fig. 3, we compare the halo-model predictions for the matter–matter power spectrum decomposed in to the different two- and one-halo terms. Scale-dependence compared to linear theory in the standard two-halo term can only ever be a suppression, which arises only due to the scale dependence in the halo window profiles [Wn(M, k) in equation 11]. This small suppression can be seen at the smallest scales shown in Fig. 3. In contrast, the inclusion of non-linear halo bias invokes a strong scale dependence in the two-halo term, which boosts the power compared to the standard, starting at |$k\sim 0.1\, h\mathrm{Mpc}^{-1}$|⁠, and this boost can reach a factor of ∼2 at |$k\sim 0.7\, h\mathrm{Mpc}^{-1}$|⁠. However, such a strong effect is not seen in full in the total power prediction because the maximum of this occurs on scales where the one-halo term (which is identical in each model) starts to dominate the power spectrum prediction. Despite this, the total power in our new model is still boosted by ∼50 per cent in the transition region between the two terms compared to the standard model. That power in our two-halo term starts to decay for |$k\gtrsim 0.7\, h\mathrm{Mpc}^{-1}$| is a combination of halo exclusion effects that are included in our measurement of βNL, window profiles truncation from equation (17), and the fact that we only measure βNL for |$k\lesssim 0.8\, h\mathrm{Mpc}^{-1}$| and extrapolate it beyond this. In any case, at such small scales the one-halo term dominates the overall power spectrum. Figure 3. Open in new tabDownload slide Halo model matter–matter power spectrum predictions (dashed) decomposed in to two-halo (solid) and one-halo terms (solid black) for the standard halo model (red) and our improved non-linear halo bias model (blue) at z = 0. The one-halo term is identical for each model. The linear power (dashed black) is also shown. The top panel shows the power spectra, while the lower panel shows the ratio of the predictions to linear theory. We see that the power in our new two-halo term is approximately double that of the standard at |$k\simeq 0.7\, h\mathrm{Mpc}^{-1}$|⁠. Figure 3. Open in new tabDownload slide Halo model matter–matter power spectrum predictions (dashed) decomposed in to two-halo (solid) and one-halo terms (solid black) for the standard halo model (red) and our improved non-linear halo bias model (blue) at z = 0. The one-halo term is identical for each model. The linear power (dashed black) is also shown. The top panel shows the power spectra, while the lower panel shows the ratio of the predictions to linear theory. We see that the power in our new two-halo term is approximately double that of the standard at |$k\simeq 0.7\, h\mathrm{Mpc}^{-1}$|⁠. 4.2 Halo-model power spectra Fig. 4 is the main result of this paper, and shows halo-model calculations with and without including βNL. We show power for autocombination and cross-combination for three halo mass bins: 1012.5–1013.0, 1013.0–1013.5, 1013.5–|$10^{14.0}\, h^{-1}\mathrm{M_\odot }$| and the matter field, all taken from multidark. The diagonal panels show the autospectra, while the off-diagonal panels show the cross-spectra. The error bars show error-on-the-mean power in each k bin. Figure 4. Open in new tabDownload slide The upper triangular set shows cross-power spectra computed between different cosmological fields using a standard halo-model calculation (red) and our improved method (blue and green) compared to measurements from simulations (the black points with errors). The cosmological fields we show are matter overdensity and three sequential halo-mass bins of halo overdensity. The two- (long-dashed coloured) and one-halo (short-dashed black) terms are also shown, but note that the one-halo term is identical for both models. The linear theory matter power spectrum is also shown (solid black). The effect of the improvement to the calculation that is discussed in this paper is prominent in the transition region between the two- and the one-halo terms. This can be better appreciated in the residual panels in the lower triangular set where the model is divided by the simulated data. For each spectrum, we see a clear improvement when using the new halo-model calculation with non-linear halo bias included. The error bars are errors-on-the-mean taken from power spectrum measurements from the simulations. Figure 4. Open in new tabDownload slide The upper triangular set shows cross-power spectra computed between different cosmological fields using a standard halo-model calculation (red) and our improved method (blue and green) compared to measurements from simulations (the black points with errors). The cosmological fields we show are matter overdensity and three sequential halo-mass bins of halo overdensity. The two- (long-dashed coloured) and one-halo (short-dashed black) terms are also shown, but note that the one-halo term is identical for both models. The linear theory matter power spectrum is also shown (solid black). The effect of the improvement to the calculation that is discussed in this paper is prominent in the transition region between the two- and the one-halo terms. This can be better appreciated in the residual panels in the lower triangular set where the model is divided by the simulated data. For each spectrum, we see a clear improvement when using the new halo-model calculation with non-linear halo bias included. The error bars are errors-on-the-mean taken from power spectrum measurements from the simulations. The top triangle set show a comparison of the raw power spectrum measurement from simulations to the model predictions. This set of plots is not very useful for investigating the fine details of the performance of our method but does demonstrate that we are generating predictions that are broadly realistic when compared to the simulated data. Note that in each panel the one-halo term is the same for both the standard and non-linear bias model calculations. The one-halo term is only present in some panels: those where there is an overlap between the haloes in the two fields that make up the power spectrum. The lower triangle set show the ratio of our model predictions to the simulation measurements, with the error bar translated from the simulations into this new space. The standard halo-model prediction is shown together with that from two versions of our calculation that include βNL. First, we discuss the halo–halo autospectra: in this case we see a small improvement when non-linear halo bias is included in the calculation, but an improvement that corrects the ∼10 per cent low residual that would otherwise be present. If we instead consider the cross-spectra between the halo-mass bins, we see a more dramatic improvement, where underpredictions in power of ∼40 per cent are almost perfectly cured when one includes the non-linear halo bias. This dramatic difference between the autospectra and cross-spectra is due to the halo shot noise, which appears in the autospectra only and is accounted for via the one-halo term of the standard halo model calculation. The fact that this is absent in the cross-spectra between different halo-mass bins allows us to see the true deficit in two-halo power when compared to the linear-bias–linear-power assumption in the standard halo model, which is obscured in the autospectra by the relatively powerful one-halo term at smaller scales. In some ways, the success of our model for the halo autospectra and cross-spectra is not much of a success, given that we essentially use the difference between the standard halo model prediction and the simulations to inform the correction in the non-linear halo bias model. What saves these plots from complete triviality is the fact that the mass bins used are different from those used in our measurement of βNL and also that the model predictions come via the full halo-model apparatus, including the choice of Tinker et al. (2010) for the mass function and bias. If there were any serious discrepancies between these choices and reality, these would manifest in the Figure. In this sense, the halo–halo panels of Fig. 4 provide a useful sanity check and inform us that our halo-model implementation is performing as expected. That the most serious discrepancies occur in the auto spectra is a result of our theoretical mass function not agreeing perfectly with the halo population seen in the simulations, which leads to a slightly different amplitude of the one-halo term in each case. This discrepancy indicates that more accurate halo mass function predictions10 would be useful in further application of this work. The more interesting panels of Fig. 4 are those that show the matter–matter and matter–halo power. These demonstrate the utility of our method, given that the non-linear halo bias correction to the matter originates from an integral over all halo masses. We show two versions of this calculation, one where we restrict the limits of the integral in equation (17) to be only over halo masses that we have actually measured from multidark. The upper limit of multidark is ν ≃ 4, which is effectively infinite from the point of view of the calculation as the results are unchanged if we extrapolate βNL above this or fix it to zero. The lower limit of multidark is ν ≃ 0.75 and our results do change depending on if we either extrapolate below this limit or fix βNL to zero, and we show the impact of these two choices in Fig. 4. Note that this choice has no impact on the halo–halo spectra that we show since all of these correspond to the correction evaluated only in sub-squares of βNL and INL shown in Fig. 2 that have been measured well. The halo–matter spectra instead correspond to slices and the matter–matter spectra corresponds to the whole plane. In all cases where matter forms part of a two-point function, we see a dramatic improvement in the accuracy of the calculations when the non-linear halo bias is included, particularly when we allow ourselves to extrapolate the measured βNL function to all values of ν. In this case, the problem of an underprediction in power in the transition region that plagues the standard halo-model calculation is dramatically ameliorated – this is the main result of this paper. In all cases, there are only very small corrections predicted by the model for |$0.08\lt k/\, h\mathrm{Mpc}^{-1}\lesssim 0.1$|⁠, which is a result of linear theory and constant bias being a very good approximation for these scales. There is perhaps some small improvement, even for |$k\sim 0.1\, h\mathrm{Mpc}^{-1}$|⁠, which could arise from βNL including some ‘perturbative’-type corrections, e.g. pre-virialization, but the data are quite noisy and we are unwilling to draw any firm conclusions about this. For |$k\gt 0.1\, h\mathrm{Mpc}^{-1}$|⁠, the correction has a significant impact on our predictions, and this demonstrates that the lack of a proper incorporation of non-linear halo bias in the standard halo model is mostly responsible for the underprediction of power in the transition region. From Fig. 4, we note that the required correction in any power spectra that involves haloes is larger for the lower mass halo bins (both in halo–halo and halo–matter spectra). The relatively good performance of halo models that pertain to higher mass haloes must therefore be due to two effects: First, that non-linear halo biasing seems to be intrinsically less important for high-mass haloes, and secondly that the one-halo term is relatively larger due to the increased shot-noise amplitude arising from the low number density of rare objects, and this one-halo term then obscures the two-halo term on scales where non-linear halo biasing is important. This suggests that traditional halo-model approaches will be more successful when describing the power spectra of fields dominated by higher mass haloes, and this is broadly the trend seen in the literature. For example, the power spectrum of the tSZ effect has a well-understood halo-model interpretation (e.g. Komatsu & Seljak 2002; Refregier & Teyssier 2002; Battaglia et al. 2012; Horowitz & Seljak 2017). In the case of tSZ, the dominance of the one-halo term is even more extreme due to the extra mass weighting that the electron pressure brings for high-mass haloes (e.g. Mead et al. 2020). On the other hand, we would expect a poorer performance of the standard halo model when considering fields that arise from lower mass haloes. For example, Addison, Dunkley & Bond (2013) note that including non-linear halo bias is necessary when modelling the autospectrum of the cosmic infrared background (CIB), which would make sense given the |$\sim 10^{12.5}\, h^{-1}\mathrm{M_\odot }$| peak halo mass for star formation efficiency (e.g. Viero et al. 2013; Planck Collaboration XVIII 2014; Maniyar, Béthermin & Lagache 2021). Neutral hydrogen (H i) is known to trace relatively low-mass haloes due to its destruction by energetic feedback processes in massive galaxies, which would suggest that traditional halo model approaches (e.g. Padmanabhan et al. 2017) may be less successful for modelling the H i distribution. 4.3 Extrapolation We note that a significant fraction of the improvement that we see in the transition region for spectra that involve matter arises only when we allow ourselves to extrapolate our measured βNL to lower halo masses: ν ≲ 0.75. The exact fraction of the improvement that depends on these low halo masses varies, but in all of our matter–halo power spectra it is at least half and in our matter–matter it is a factor of ∼5. Taken at face value, this would suggest that it is the non-linear bias of low-mass haloes relative to high-mass haloes, and the non-linear bias of low-mass haloes with themselves, that are primarily responsible for the power deficit in the transition region. This agrees somewhat with conclusions drawn by van Daalen & Schaye (2015) who demonstrated that significant power arises when matter that is ‘just outside’ the halo virial radius is accounted for. It it plausible that this matter is comprised of infalling, low-mass haloes, and that it is the non-linear biasing (in our language) of these objects that is adding considerable power. Of course, we should also be wary of trusting our extrapolation, particularly given that we are inferring the non-linear bias of roughly half of the mass from the other half, and there may well be complexities regarding the distribution of low-mass haloes that are not captured by linear extrapolation. To partially address this we included data from the bolshoi simulation in our βNL measurement.11bolshoi is similar to multidark, but is a smaller |$250\, h^{-1}\mathrm{Mpc}$| cube, so the mass resolution is better by a factor of 64, which allows us to measure the non-linear bias of haloes down to ν ∼ 0.5, albeit with increased noise at any given scale due to the smaller box size. If we do this we get broadly consistent results with those shown in Fig. 4 with the amount of extrapolation lessened, which is shown in Fig. 5. The difference between the extrapolation in the two cases is around 5 per cent in power. The degraded performance that we see for |$k\sim 0.1\, h\mathrm{Mpc}^{-1}$| is due to noise from (the small volume of) bolshoi transferring in to our measurements. This gives us hope that our extrapolation is moderately robust. In any case, even if one is reluctant to trust our extrapolation, the fact that we see any improvement at all in the transition region tells us that non-linear halo bias is at least part of the puzzle surrounding missing power in this region, and does not disprove the hypothesis that it may be entirely responsible. Figure 5. Open in new tabDownload slide Residual power spectra from halo models compared to measurements from simulations at z = 0 for the matter–matter power spectrum (top left) and matter–halo power spectra for three different halo-mass bins (other panels). In each case, we show the standard halo model (red) and then the non-linear halo bias model with βNL taken from multidark (green); this extrapolated to low halo mass (blue); the combination of multidark and bolshoi (yellow); this extrapolated to low halo mass (purple). The difference between green and yellow therefore shows the difference when including the ∼0.5 < ν < 0.75 haloes from bolshoi. The difference between purple and blue gives some indication of the robustness of our extrapolation. These panels correspond to the left-hand column of the lower triangle of Fig. 4. Figure 5. Open in new tabDownload slide Residual power spectra from halo models compared to measurements from simulations at z = 0 for the matter–matter power spectrum (top left) and matter–halo power spectra for three different halo-mass bins (other panels). In each case, we show the standard halo model (red) and then the non-linear halo bias model with βNL taken from multidark (green); this extrapolated to low halo mass (blue); the combination of multidark and bolshoi (yellow); this extrapolated to low halo mass (purple). The difference between green and yellow therefore shows the difference when including the ∼0.5 < ν < 0.75 haloes from bolshoi. The difference between purple and blue gives some indication of the robustness of our extrapolation. These panels correspond to the left-hand column of the lower triangle of Fig. 4. 4.4 Relation to perturbation theory The model presented in this paper does not explicitly include beyond-linear perturbation theory, and instead takes a pragmatic approach by measuring the required non-linear halo bias correction directly from N-body simulation data. We include here a discussion of how our results relate to beyond-linear perturbation theory (see e.g. Bernardeau et al. 2002 and references therein), which is known to be an excellent description of the matter clustering on large scales. For example, one-loop SPT is sub-per cent ΛCDM spectra for |$k\lesssim 0.08\, h\mathrm{Mpc}^{-1}$|⁠, whereas linear theory provides a ∼2 per cent over prediction of power on these ‘pre-virializataion’ scales. Higher order perturbative schemes and effective field theories have the potential to push this accuracy to smaller scales (for recent incarnations, see e.g. Seljak & Vlah 2015; Foreman, Perrier & Senatore 2016; Philcox et al. 2020). Taking a schematic approach, assume that one believes that a good large-scale model of the quasi-linear matter–matter power is given by $$\begin{eqnarray*} P^\mathrm{QL}_\mathrm{mm}(k)\simeq P^\mathrm{lin}_\mathrm{mm}(k)+P^\mathrm{cor}_\mathrm{mm}(k), \end{eqnarray*}$$(29) where the ‘cor’ term on the right-hand side contains all corrections to linear theory. One could explicitly include this in our approach in two ways: First, the large-scale limit of βNL could be changed to enforce |$I^\mathrm{NL}_\mathrm{mm}(k\rightarrow 0)=P^\mathrm{cor}_\mathrm{mm}(k)/P^\mathrm{lin}_\mathrm{mm}(k)$| (in the previous discussion this limit was assumed to be zero). Secondly, one could replace the factors of |$P^\mathrm{lin}_\mathrm{mm}(k)$| that appear on the right-hand side of equation (16) with |$P^\mathrm{QL}_\mathrm{mm}(k)$|⁠, such that all of the halo biasing is defined relative to a spectrum that one believes to be more correct at smaller scales. This would entail a redefinition of βNL, again with the linear power replaced by the quasi-linear power (in equation 23), and would need the limit |$I^\mathrm{NL}_\mathrm{mm}(k\rightarrow 0)=0$| to be enforced on scales where |$P^\mathrm{QL}_\mathrm{mm}(k)$| was thought to be a good description of the power. This redefinition may reduce the scale-dependence of our βNL function. We attempted to include the above by taking |$P_\mathrm{mm}^\mathrm{cor}(k)$| to be given by one-loop SPT and redefining our biasing relative to this. However, the problem we encountered is that while SPT provides a slightly improved prediction of power for |$k\lt 0.1\, h\mathrm{Mpc}^{-1}$| when it does go wrong, it goes wrong quite spectacularly – grossly overpredicting the non-linear matter–matter power for |$k\gtrsim 0.2\, h\mathrm{Mpc}^{-1}$|⁠. Given this, we found using linear theory to be more convenient because when it does go wrong it does not go very wrong. In addition, the SPT correction on scales where we might notice it is much smaller than the statistical errors from the limited simulation volume, so would not be noticeable in any of our plots. In the perturbative bias models of McDonald (2006) or Smith et al. (2007), one-loop SPT and an Eulerian bias expansion are taken to be good descriptions of the matter and halo clustering on large scales. This has the advantage that the integral constraints on the Eulerian bias coefficients ensure that the SPT result is returned when integrating over all halo masses (⁠|$I^\mathrm{NL}_\mathrm{mm}(k)=0$| in our language). The disadvantage is that this model fails when perturbation theory fails, and this failure can be quite extreme compared to the failures endured by linear theory (with linear bias) at the same scale. Still, it may be possible to incorporate the successes of the Smith et al. (2007) approach with our model by making a split in method based on wavenumber. One may also consider using the result from a fitting function (halofit, hmcode, or an emulator) for the full non-linear matter–matter power in place of the linear spectrum in equations (16) and (23). We do not do this as we wish for our model to be a good description of the power of any combination of fields, including matter–matter. Including the full non-linear matter–matter in the two-halo term therefore leads to a recursion that we would prefer to avoid. 4.5 Galaxy clustering example Finally, in Fig. 6, we demonstrate the importance of our non-linear halo bias correction when comparing analytical halo-model predictions to mock galaxy catalogues that we generate from the multidark halo catalogues using a simple halo-occupation distribution (HOD) prescription. Our simple HOD catalogues are generated12 by taking 1012.5 and |$10^{14}\, h^{-1}\mathrm{M_\odot }$| as the minimum and maximum halo masses that can host galaxies. Haloes of |$10^{12.5}\, h^{-1}\mathrm{M_\odot }$| are given a single central galaxy, while haloes above this mass are assigned a number of additional satellite galaxies that depends linearly on the halo mass. For example, a |$4\times 10^{12.5}\, h^{-1}\mathrm{M_\odot }$| halo would be assigned one central and three satellite galaxies. These satellite galaxies are then set to follow an isothermal radial distribution around the halo centre,13 extending out to the halo virial radius. We make a realization of this galaxy distribution using the multidark halo catalogues and we measure the three-dimensional galaxy–galaxy, galaxy–matter, and central–satellite power spectra. We compare these to predictions from the halo model described in this paper, both with and without the non-linear bias correction in Fig. 6. Once again, it can be seen that the inclusion of the non-linear halo bias makes a significant improvement to the halo-model predictions around the quasi-linear transition region. We extrapolate βNL to low halo masses as described in Section 4.3, but this only affects the galaxy–matter power spectrum because all galaxies live in haloes that have been well measured in βNL. These results demonstrate the power of the approach advocated in this paper, in that once the βNL correction has been measured and characterized, which only needs to be done once using the halo catalogue, any type of HOD prescription can be applied and the novel halo model can be expected to make reasonable predictions. To contrast with this, in Fig. 6 we also show a halo model where we replace the linear power in the standard two-halo term (equation 11) with the non-linear matter–matter power (which we take from halofit of Takahashi et al. 2012), thus assuming that the haloes in which the galaxies reside are linearly biased with respect to the underlying non-linear matter distribution. This scheme has no physical motivation whatsoever, but is occasionally considered in the literature. We see that this leads to dramatic over-predictions of power for |$k\gt 0.1\, h\mathrm{Mpc}^{-1}$|⁠, which is perhaps not surprising given that this model lacks a physical basis and that two one-halo contributions now enter the calculation. Figure 6. Open in new tabDownload slide Residual power spectra from the standard halo model (red) and our new model that includes non-linear halo bias (dark blue) compared to power spectra from HOD realizations in simulations at z = 0. We also show power from a halo model where the linear power has been replaced by the non-linear matter–matter power (light blue), thus assuming a linear halo bias with respect to the underlying non-linear distribution. The generation of the HOD galaxies is described in the text. We show the galaxy–galaxy (top), galaxy–matter, and the central–satellite (bottom) power spectra. The error bars show error-on-the-mean power in each k bin. Figure 6. Open in new tabDownload slide Residual power spectra from the standard halo model (red) and our new model that includes non-linear halo bias (dark blue) compared to power spectra from HOD realizations in simulations at z = 0. We also show power from a halo model where the linear power has been replaced by the non-linear matter–matter power (light blue), thus assuming a linear halo bias with respect to the underlying non-linear distribution. The generation of the HOD galaxies is described in the text. We show the galaxy–galaxy (top), galaxy–matter, and the central–satellite (bottom) power spectra. The error bars show error-on-the-mean power in each k bin. 5 SUMMARY We have proposed an extension to the analytical halo-model formalism that allows for beyond-linear halo bias to be included within an otherwise standard calculation. We have expressed our correction in such a way that the a single new term is added to the (otherwise standard) two-halo term (equations 16 and 17), so that it can be easily integrated within existing halo-model implementations. The new correction requires the evaluation of a double integral over halo masses of a new function, which we call βNL, that accounts for all aspects of the two-point function of haloes that are not accounted for by the simple ‘linear bias with respect to linear matter field’ model. If one prefers to think of halo bias as being defined with respect to the underlying non-linear field then our βNL correction can be thought of as an effective bias that incorporates both standard halo biasing and the non-linear evolution of the underlying matter. We have demonstrated that our new corrective term generally boosts power in the transition region between the two- and one-halo terms, which is known to be poorly modelled by traditional halo-model approaches. Our results suggest that most of the power deficit in this region arises from the lack of consideration of beyond-linear halo bias, and that this can be a ∼50 per cent under prediction of power in some cases. We show that when including our correction we get essentially perfect auto spectra for galaxies, and much improved power spectra for matter–galaxy cross-spectra. For the matter–matter power spectrum we also get much improved predictions, but we noted a sensitivity in our results to the extrapolation to low-mass haloes that we are required to make of βNL. The advantage of our new approach over traditional approaches is that we capture a genuine, physical effect that is missing from the classic calculation. We also demonstrate that this effect is important, in that it provides corrections of order tens of per cent at scales that are relevant to contemporary surveys. Finally, we show that once the correction has been measured from haloes, it then applies to any two-point function one could calculate via the halo model – it does not need to be recalibrated each time. Of course, the disadvantage is that we must measure a new ingredient from simulations and that currently there exists no fitting form for βNL. Recently, Nishimichi et al. (2019) have developed the dark quest emulator, which provides the halo–halo and halo–matter correlation functions with the aim of providing an accurate basis for HOD modelling for galaxy–galaxy and galaxy–matter clustering studies. This has strong similarities to the approach we advocate in this paper, in that the non-linear biasing of haloes is automatically accounted for, and these terms are extracted from N-body simulations. However, our method is different in some important ways: We retain exactly the standard halo model apparatus on the largest scales, where all power spectra are described by the standard ‘linear power multiplied by a linear biasing factor’. This allows our non-linear biasing term to appear as a simple addition to an otherwise standard two-halo term, essentially concentrating all new content in the single new βNL term. dark quest, also recovers this large-scale limit, but does so less transparently. We feel that this gives our approach a pedagogical advantage. We also demonstrated that good results can be obtained for spectra that involve the matter field using information from a measured correction that comes only from the halo field, which is not considered by Nishimichi et al. (2019). We note that this was not guaranteed to work, for example the deficit in power in the transition region could have arisen from ‘one-halo’ effects, such as dispersion in halo profiles at fixed mass, or large-scale correlations in shapes. The fact that we see such dramatic improvement tells us that non-linear halo bias is the most important missing ingredient at quasi-linear scales. We have not investigated the cosmology dependence of the non-linear halo bias. We suggest that writing the correction in the form given in equation (23) may remove some of the cosmology dependence because it is a ratio of power spectra (see Mead 2017), and we also suggest that expressing the correction in terms of peak height, ν, rather than halo mass, may remove additional cosmology dependence. However, we noted that even when expressed in this way our correction retained some redshift dependence (Appendix B), which hints at a more complicated relationship with the underlying cosmology. The cosmology dependence could also be investigated using the dark quest emulator of Nishimichi et al. (2019) and this could even be used to build an emulator for βNL, but dark quest is not currently public. One could also use a set of simulations such as the quijote of Villaescusa-Navarro et al. (2020) to build an emulator for βNL from scratch. We strongly recommend building such an emulator as a way forward for cosmology, which will allow the halo model to continue to be useful as we enter the epoch of precision measurements. Note carefully that this emulator only needs to be build using information about halo clustering, and then the method outlined in this paper can be used to apply the correction to any halo model two-point function. One may also consider the application of our work to halo-model descriptions of higher order statistics of cosmological fields. In this case, we note that βNL originates from a two-point function and therefore that beyond-two-point information about halo clustering is absent. Therefore, to successfully describe the n-point correlation of a field would require a measurement of the n-point version of βNL. This is in contrast to the standard halo model, where in principle the n-point functions all follow from the class halo-model ansatz. However, at large scales it is probable that a perturbative model for non-linear bias could be used, and from this one could consistently build-up the n-point equivalents of βNL at large scales. This could then be augmented with simulation measurements at smaller scales. Our work was motivated by a consideration of the way halo models are traditionally used to infer properties of halo occupation via cross-correlation, be it from tracers or via some emissivity (e.g. galaxy or CMB lensing; tSZ; galaxy clustering; CIB), but we note that the error bars on current measurements are quite large (e.g. Hill & Spergel 2014; Hojjati et al. 2017; Dvornik et al. 2018; Tanimura et al. 2019). This contrasts with more precise work using spectroscopic galaxy surveys, particularly in redshift space, where careful scale cuts can be made. In order to assess the importance of our correction in cases with more exacting accuracy requirements would require us to measure βNL more accurately, for example, using an ensemble of simulations. These simulations could be quite low resolution since all that is required for βNL is the halo clustering; our method requires no knowledge of the internal structure of haloes. Once βNL has been measured from the halo population, it can then be used in any subsequent halo-model calculation. In principle, one could attempt to calculate βNL from perturbation theory, but we instead measure the function from the multidark simulations. We note that βNL provides additional power between |$0.1\lesssim k/\, h\mathrm{Mpc}^{-1}\lesssim 1$|⁠, thus casting doubt on whether it could ever be described fully perturbatively. However, we suggest that in future a perturbative description at low k could be combined with measurements from simulations at smaller scales. In our calculations, we evaluate βNL from an interpolator that is constructed from our measurements, but in principle one could also create a fitting function or an emulator for βNL, in the same spirit that McClintock et al. (2019a) have provided for the linear halo bias and Valcin et al. (2019) have provided for the large-scale non-linear halo bias. We end this paper with a warning: We have demonstrated that the standard halo model systematically underestimates power in the transition region between the two- and one-halo terms, and that there is compelling evidence that all the missing power in this region originates from the beyond-linear biasing of haloes. The amount of power underestimated depends on the particular two-point combination under consideration, but seems to be higher for spectra that include lower mass haloes, and those for which the one-halo term is subdued (e.g. cross-spectra between fields that arise from different haloes), because a powerful one-halo term has the potential to obscure the effect of non-linear halo biasing. In the specific three-dimensional power spectra investigated in this paper, the amount of missing power varied between 10 and 50 per cent. Despite this, the standard halo model is often used to draw conclusions from cosmological data sets about the connection between an observable (e.g. galaxies, tSZ, CIB) and the host haloes and even about cosmological parameters. These halo models can often be quite complicated, with a plethora of parameters that govern the distribution (be it of galaxies or some emissivity) within the halo. It is entirely possible that conclusions drawn from such models will be significantly biased if they rely on data that is modelled by the transition region of a standard halo model calculation. This will be particularly true if there are ‘physical’ terms in these models that can add power in the transition region, which is often the case when dealing with two-dimensional models of projected three-dimensional fields, because one does not then enjoy the k localization of physical effects as in three dimensions. We therefore recommend that βNL (or something like it) be included in future data analyses that use the halo model. At least, if it is not used, we then recommend excising data that are modelled by the transition region between the two- and one-halo terms to avoid biasing results. Clearly, the potential constraining power of data that is depleted by this excision will vary between two-point functions, and we suggest that it will be greater for data sets that probe lower mass haloes. Finally, if the accuracy of our βNL measurement is not sufficient for precision analyses then one could at least use our results to gauge the importance of the quasi-linear regime to a specific measurement. ACKNOWLEDGEMENTS AJM has received funding from the Horizon 2020 research and innovation programme of the European Union under the Marie Skłodowska-Curie Actions grant number 702971 and support from the European Research Council under grant number 647112. LV acknowledges support by European Union’s Horizon 2020 research and innovation programme ERC (BePreSySe, grant number 725327) and Spanish MINECO under project PGC2018-098866-B-I00, FEDER, UE. AJM acknowledges useful conversations with Samuel Brieden, Catherine Heymans, Alex Hall, John Peacock, Tilman Tröster, Oliver Philcox, and David Alonso. The CosmoSim data base used in this paper is a service by the Leibniz-Institute for Astrophysics Potsdam (AIP). The MultiDark data base was developed in cooperation with the Spanish MultiDark Consolider Project CSD2009-00064. The authors gratefully acknowledge the Gauss Centre for Supercomputing e.V. (www.gauss-centre.eu) and the Partnership for Advanced Supercomputing in Europe (PRACE, http://www.prace-ri.eu) for funding the MultiDark simulation project by providing computing time on the GCS Supercomputer SuperMUC at Leibniz Supercomputing Centre (LRZ, http://www.lrz.de). The multidark data base used in this paper and the web application providing online access to it were constructed as part of the activities of the German Astrophysical Virtual Observatory as result of a collaboration between the Leibniz-Institute for Astrophysics Potsdam (AIP) and the Spanish MultiDark Consolider Project CSD2009-00064. The Bolshoi and MultiDark simulations were run on the NASA’s Pleiades supercomputer at the NASA Ames Research Center. The MultiDark-Planck (MDPL) and the BigMD simulation suite have been performed in the Supermuc supercomputer at LRZ using time granted by PRACE. DATA AVAILABILITY STATEMENT The code and some data used in this work are available at https://github.com/alexander-mead/BNL. The remaining data underlying this article will be shared on request to the corresponding author. Footnotes 1 Achieving these limits is difficult numerically because of the large amount of mass contained in low-mass haloes according to most popular mass functions. Special care must be taken with the two-halo integral in the case of power spectra that involve the matter field. See Appendix A. 2 https://www.cosmosim.org 3 We obtain very similar results if we utilize haloes identified via the bound-density maximum (BDM) algorithm of Klypin & Holtzman (1997). 4 This would necessitate remeasuring βNL from a different halo catalogue with the new mass definition, or else translating between halo mass definitions (e.g. via the NFW profile). 5 We obtain very similar results if we increase the number of mass bins from 8 to 12. 6 We also use data from the bolshoi simulation, which is the same as multidark in all respects apart from being smaller, |$250\, h^{-1}\mathrm{Mpc}$|⁠, and thus able to resolve lower mass haloes. 7 In fact, we take the z ≃ 3 density field grown to the desired redshift because the initial conditions of these simulations no longer exist. 8 In principle, one could also use a linear bias emulator (e.g. Valcin et al. 2019; McClintock et al. 2019a). 9 We get essentially identical results for everything presented in this paper if we replace the NFW profile with an isothermal, ∝1/r2, profile. 10 In future, one could use the mass-function emulators of either McClintock et al. (2019b) or Bocquet et al. (2020). 11 Note that here we have used the BDM halo catalogues from multidark and bolshoi because bolshoi has no public rockstar catalogue. 12 We could have used a more complicated HOD, such as Zheng et al. (2005), but choose to keep the example simple for the sake of pedagogy. 13 Our model is quite insensitive to the exact radial distribution of satellite galaxies on the scales we show in Fig. 6, instead being primarily sensitive to the occupation number of each halo. 14 This condition is often, but not always, imposed on functional forms that are fitted to halo mass functions measured from N-body simulations. 15 Note that this means that haloes can overlap. REFERENCES Addison G. E. , Dunkley J., Spergel D. N., 2012 , MNRAS , 427 , 1741 10.1111/j.1365-2966.2012.21664.x Crossref Addison G. E. , Dunkley J., Bond J. R., 2013 , MNRAS , 436 , 1896 10.1093/mnras/stt1703 Crossref Agarwal S. , Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2012 , MNRAS , 424 , 1409 Agarwal S. , Abdalla F. B., Feldman H. A., Lahav O., Thomas S. A., 2014 , MNRAS , 439 , 2102 10.1093/mnras/stu090 Crossref Battaglia N. , Bond J. R., Pfrommer C., Sievers J. L., 2012 , ApJ , 758 , 75 10.1088/0004-637X/758/2/74 Crossref Behroozi P. S. , Wechsler R. H., Wu H.-Y., 2013 , ApJ , 762 , 109 10.1088/2041-8205/762/2/L31 Crossref Bernardeau F. , Colombi S., Gaztañaga E., Scoccimarro R., 2002 , Phys. Rep. , 367 , 1 10.1016/S0370-1573(02)00135-7 Crossref Bocquet S. , Heitmann K., Habib S., Lawrence E., Uram T., Frontiere N., Pope A., Finkel H., 2020 , ApJ , 901 , 5 Bryan G. L. , Norman M. L., 1998 , ApJ , 495 , 80 : 10.1086/305262 Crossref Cacciato M. , van den Bosch F. C., More S., Li R., Mo H. J., Yang X., 2009 , MNRAS , 394 , 929 10.1111/j.1365-2966.2008.14362.x Crossref Cacciato M. , Lahav O., van den Bosch F. C., Hoekstra H., Dekel A., 2012 , MNRAS , 426 , 566 : 10.1111/j.1365-2966.2012.21762.x Crossref Chen A. Y. , Afshordi N., 2020 , Phys. Rev. D , 101 , 103522 10.1103/PhysRevD.101.103522 Crossref Cooray A. , Sheth R., 2002 , Phys. Rep. , 372 , 1 10.1016/S0370-1573(02)00276-4 Crossref Duffy A. R. , Schaye J., Kay S. T., Dalla Vecchia C., 2008 , MNRAS , 390 , L64 10.1111/j.1745-3933.2008.00537.x Crossref Dvornik A. et al. , 2018 , MNRAS , 479 , 1240 10.1093/mnras/sty1502 Crossref Fedeli C. , Semboloni E., Velliscig M., Van Daalen M., Schaye J., Hoekstra H., 2014 , J. Cosmol. Astropart. Phys. , 8 , 28 10.1088/1475-7516/2014/08/028 Crossref Foreman S. , Perrier H., Senatore L., 2016 , J. Cosmol. Astropart. Phys. , 2016 , 027 Garcia R. , Rozo E., Becker M. R., More S., 2020 , preprint (arXiv:2006.12751) Giocoli C. , Bartelmann M., Sheth R. K., Cacciato M., 2010 , MNRAS , 408 , 300 10.1111/j.1365-2966.2010.17108.x Crossref Hang Q. , Alam S., Peacock J. A., Cai Y.-C., 2021 , MNRAS , 501 , 1481 Hayashi E. , White S. D. M., 2008 , MNRAS , 388 , 2 10.1111/j.1365-2966.2008.13371.x Crossref Hill J. C. , Spergel D. N., 2014 , J. Cosmol. Astropart. Phys. , 2 , 030 10.1088/1475-7516/2014/02/030 Crossref Hill J. C. , Baxter E. J., Lidz A., Greco J. P., Jain B., 2018 , Phys. Rev. D , 97 , 083501 Hojjati A. et al. , 2017 , MNRAS , 471 , 1565 10.1093/mnras/stx1659 Crossref Horowitz B. , Seljak U., 2017 , MNRAS , 469 , 394 Klypin A. , Holtzman J., 1997 , preprint (arXiv:e-print) Klypin A. A. , Trujillo-Gomez S., Primack J., 2011 , ApJ , 740 , 102 10.1088/0004-637X/740/2/102 Crossref Knabenhans M. et al. , 2019 , MNRAS , 484 , 5509 10.1093/mnras/stz197 Crossref Knebe A. et al. , 2011 , MNRAS , 415 , 2293 10.1111/j.1365-2966.2011.18858.x Crossref Komatsu E. , Seljak U., 2002 , MNRAS , 336 , 1256 10.1046/j.1365-8711.2002.05889.x Crossref Koukoufilippas N. , Alonso D., Bilicki M., Peacock J. A., 2020 , MNRAS , 491 , 5464 Lawrence E. , Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010 , ApJ , 713 , 1322 Lawrence E. et al. , 2017 , ApJ , 847 , 50 Lewis A. , Challinor A., Lasenby A., 2000 , ApJ , 538 , 473 10.1086/309179 Crossref Ma Y.-Z. , Van Waerbeke L., Hinshaw G., Hojjati A., Scott D., Zuntz J., 2015 , J. Cosmol. Astropart. Phys. , 9 , 046 Mandelbaum R. , Tasitsiomi A., Seljak U., Kravtsov A. V., Wechsler R. H., 2005 , MNRAS , 362 , 1451 Maniyar A. , Béthermin M., Lagache G., 2021 , A&A , 645 , A40 McClintock T. et al. , 2019a , preprint (arXiv:1907.13167) McClintock T. et al. , 2019b , ApJ , 872 , 53 McDonald P. , 2006 , Phys. Rev. D , 74 , 103512 10.1103/PhysRevD.74.129901 Crossref Mead A. J. , 2017 , MNRAS , 464 , 1282 10.1093/mnras/stw2312 Crossref Mead A. J. , Peacock J. A., Heymans C., Joudaki S., Heavens A. F., 2015 , MNRAS , 454 , 1958 10.1093/mnras/stv2036 Crossref Mead A. J. , Heymans C., Lombriser L., Peacock J. A., Steele O. I., Winther H. A., 2016 , MNRAS , 459 , 1468 10.1093/mnras/stw681 Crossref Mead A. J. , Tröster T., Heymans C., Van Waerbeke L., McCarthy I. G., 2020 , A&A , 641 , A130 10.1051/0004-6361/202038308 Crossref Mead A. J. , Brieden S., Tröster T., Heymans C., 2021 , MNRAS , 502 , 1401 Mohammed I. , Seljak U., 2014 , MNRAS , 445 , 3382 10.1093/mnras/stu1972 Crossref Nakamura T. T. , Suto Y., 1997 , Prog. Theor. Phys. , 97 , 49 10.1143/PTP.97.49 Crossref Navarro J. F. , Frenk C. S., White S. D. M., 1997 , ApJ , 490 , 493 10.1086/304888 Crossref Nishimichi T. et al. , 2019 , ApJ , 884 , 29 Padmanabhan H. , Refregier A., Amara A., 2017 , MNRAS , 469 , 2323 10.1093/mnras/stx979 Crossref Peacock J. A. , Smith R. E., 2000 , MNRAS , 318 , 1144 10.1046/j.1365-8711.2000.03779.x Crossref Philcox O. H. E. , Spergel D. N., Villaescusa-Navarro F., 2020 , Phys. Rev. D , 101 , 123520 10.1103/PhysRevD.101.123520 Crossref Planck Collaboration XVIII , 2014 , A&A , 571 , A30 Prada F. , Klypin A. A., Cuesta A. J., Betancort-Rijo J. E., Primack J., 2012 , MNRAS , 423 , 3018 10.1111/j.1365-2966.2012.21007.x Crossref Refregier A. , Teyssier R., 2002 , Phys. Rev. D , 66 , 043002 10.1103/PhysRevD.66.043002 Crossref Riebe K. et al. , 2013 , Astron. Nachr. , 334 , 691 10.1002/asna.201211900 Crossref Schmidt F. , 2016 , Phys. Rev. D , 93 , 063512 10.1103/PhysRevD.93.063512 Crossref Seljak U. , 2000 , MNRAS , 318 , 203 10.1046/j.1365-8711.2000.03715.x Crossref Seljak U. , Vlah Z., 2015 , Phys. Rev. D , 91 , 123516 10.1103/PhysRevD.91.123516 Crossref Sheth R. K. , Tormen G., 1999 , MNRAS , 308 , 119 10.1046/j.1365-8711.1999.02692.x Crossref Smith R. E. , Markovic K., 2011 , Phys. Rev. D , 84 , 063507 10.1103/PhysRevD.84.063507 Crossref Smith R. E. , Watts P. I. R., 2005 , MNRAS , 360 , 203 10.1111/j.1365-2966.2005.09053.x Crossref Smith R. E. et al. , 2003 , MNRAS , 341 , 1311 10.1046/j.1365-8711.2003.06503.x Crossref Smith R. E. , Scoccimarro R., Sheth R. K., 2007 , Phys. Rev. D , 75 , 063512 10.1103/PhysRevD.75.063512 Crossref Takahashi R. , Sato M., Nishimichi T., Taruya A., Oguri M., 2012 , ApJ , 761 , 152 10.1088/0004-637X/761/2/152 Crossref Tanimura H. , Aghanim N., Douspis M., Beelen A., Bonjean V., 2019 , A&A , 625 , A67 10.1051/0004-6361/201833413 Crossref Tinker J. L. , Weinberg D. H., Zheng Z., Zehavi I., 2005 , ApJ , 631 , 41 Tinker J. L. , Robertson B. E., Kravtsov A. V., Klypin A., Warren M. S., Yepes G., Gottlöber S., 2010 , ApJ , 724 , 878 10.1088/0004-637X/724/2/878 Crossref Valageas P. , Nishimichi T., 2011 , A&A , 527 , A87 10.1051/0004-6361/201015685 Crossref Valcin D. , Villaescusa-Navarro F., Verde L., Raccanelli A., 2019 , J. Cosmol. Astropart. Phys. , 2019 , 057 van Daalen M. P. , Schaye J., 2015 , MNRAS , 452 , 2247 10.1093/mnras/stv1456 Crossref van den Bosch F. C. , More S., Cacciato M., Mo H., Yang X., 2013 , MNRAS , 430 , 725 10.1093/mnras/sts006 Crossref Viero M. P. et al. , 2013 , ApJ , 779 , 32 Villaescusa-Navarro F. et al. , 2020 , ApJS , 250 , 2 10.3847/1538-4365/ab9d82 Crossref Voivodic R. , Rubira H., Lima M., 2020 , J. Cosmol. Astropart. Phys. , 2020 , 033 Wolz L. , Murray S. G., Blake C., Wyithe J. S., 2019 , MNRAS , 484 , 1007 10.1093/mnras/sty3142 Crossref Zheng Z. et al. , 2005 , ApJ , 633 , 791 10.1086/466510 Crossref APPENDIX A: NUMERICAL CALCULATIONS In this paper, we have demonstrated that the accuracy of analytical halo-model calculations can be substantially improved by including the non-linear bias of haloes in the two-halo term. Unfortunately, this increases the computational complexity of a calculation and also provides some additional numerical hurdles to overcome in comparison with the standard calculation. In this Appendix, we present some of the technical details of how we evaluate our new non-linear bias term. We first remind the reader of a well-known numerical problem (and a solution) that is encountered when trying to evaluate the standard two-halo term for matter–field combinations: $$\begin{eqnarray*} P^{2\mathrm{H}}_{uv}(k)=P^\mathrm{lin}_\mathrm{mm}(k) \prod _{n=u,v}\left[\int _0^\infty W_n(M,k)b(M)n(M)\, \mathrm{d}M\right]. \end{eqnarray*}$$(A1) In principle, one ought to integrate over all halo masses but usually the integration range is restricted to be finite so that the properties of very high- and low-mass haloes do not need to be specified. The mass function purports to describe all haloes, however, most popular mass functions (e.g. Sheth & Tormen 1999; Tinker et al. 2010) propose that a large fraction of mass is locked within very low-mass haloes, if taken literally. The reality is that very little is known about the ultralow-mass halo population at z = 0, and the inferences about this population that are made by standard mass functions are a result of the condition that all mass reside within haloes, $$\begin{eqnarray*} \int _0^\infty Mn(M)\, \mathrm{d}M=\bar{\rho }, \end{eqnarray*}$$(A2) being imposed on the fitting function by hand.14 For a matter field, we know that in the k → 0 limit the integral in equation (A1) must equal unity, equivalently: $$\begin{eqnarray*} \int _0^\infty Mb(M)n(M)\, \mathrm{d}M=\bar{\rho }. \end{eqnarray*}$$(A3) but this requires an unreasonably low lower limit on the integral. As discussed in Schmidt (2016) and Mead et al. (2020) one solution is to enforce equation (A3), which can be achieved in practice by assuming that all mass in haloes below the lower integration limit, Mmin, is contained in haloes of mass exactly Mmin: $$\begin{eqnarray*} n(M) \rightarrow n(M)\Theta (M-M_\mathrm{min})+\frac{A(M_\mathrm{min})\delta _\mathrm{D}(M-M_\mathrm{min})}{b(M_\mathrm{min})M_\mathrm{min}/\bar{\rho }}, \end{eqnarray*}$$(A4) where $$\begin{eqnarray*} A(M_\mathrm{min})=1-\frac{1}{\bar{\rho }}\int _{M_\mathrm{min}}^\infty Mb(M)n(M)\, \mathrm{d}M, \end{eqnarray*}$$(A5) which can be easily evaluated and should be some number between zero and unity, with ∼0.5 being typical for standard mass functions and standard integration limits. This procedure is valid because the linear halo bias tends to be constant for low-mass haloes and is sufficient as long as the calculation is not sensitive to the profiles (k dependence) of haloes with M < Mmin. Note that this problem generally does not affect the one-halo term (equation 8) because the extra factor of W within the integral ensures that the this term is dominated by comparatively high-mass haloes that will be included within the range of haloes that is typically integrated over. Our improved halo-model calculation requires us to evaluate the double integral $$\begin{eqnarray*} I^{\mathrm{NL}}_{uv}(k) &=& \int _0^\infty \int _0^\infty \beta ^\mathrm{NL}(M_1,M_2,k)W_u(M_1,k)\nonumber\\ &&\times \, W_v(M_2,k)b(M_1)b(M_2)n(M_1)n(M_2)\, \mathrm{d}M_1\mathrm{d}M_2. \end{eqnarray*}$$(A6) If either Wu or Wv pertain to matter we run into the same problems as for the standard two-halo term because we would usually evaluate this term over a finite range in mass. Once again, the upper limit is usually not a concern, so here we will treat it as effectively infinite (⁠|$M=10^{16}\, h^{-1}\mathrm{M_\odot }$| is a sensible upper limit for a standard ΛCDM model at z = 0) and it is the lower limit that is problematic. However, we can use the same trick and replace both n(M) that appear in equation (A6) with that in equation (A4). This splits equation (A6) into four terms: $$\begin{eqnarray*} I^{\mathrm{11}}_{uv}(k)= A^2(M_\mathrm{min}) \frac{W_u(M_\mathrm{min},k)}{M_\mathrm{min}/\bar{\rho }} \frac{W_v(M_\mathrm{min},k)}{M_\mathrm{min}/\bar{\rho }}, \end{eqnarray*}$$(A7) $$\begin{eqnarray*} I^{\mathrm{12}}_{uv}(k) &=& A(M_\mathrm{min})\frac{W_u(M_\mathrm{min},k)}{M_\mathrm{min}/\bar{\rho }}\int _{M_\mathrm{min}}^\infty \beta ^\mathrm{NL}(M_\mathrm{min}, M_2, k)\nonumber\\ &&\times \, W_v(M_2, k)b(M_2)n(M_2)\, \mathrm{d}M_2, \end{eqnarray*}$$(A8) $$\begin{eqnarray*} I^{\mathrm{21}}_{uv}(k) &=& A(M_\mathrm{min})\frac{W_v(M_\mathrm{min},k)}{M_\mathrm{min}/\bar{\rho }}\int _{M_\mathrm{min}}^\infty \beta ^\mathrm{NL}(M_1, M_\mathrm{min}, k)\nonumber\\ &&\times \, W_u(M_1, k)b(M_1)n(M_1)\, \mathrm{d}M_1, \end{eqnarray*}$$(A9) $$\begin{eqnarray*} I^{\mathrm{22}}_{uv}(k) &=& \int _{M_\mathrm{min}}^\infty \int _{M_\mathrm{min}}^\infty \beta ^\mathrm{NL}(M_1,M_2,k)W_u(M_1,k)\nonumber\\ &&\times \, W_v(M_2,k)b(M_1)b(M_2)n(M_1)n(M_2)\, \mathrm{d}M_1\mathrm{d}M_2, \end{eqnarray*}$$(A10) which can all be evaluated over M ∈ [Mmin, ∞]. The final result is then given by |$I^{\mathrm{NL}}_{uv}=I^{\mathrm{11}}_{uv}+I^{\mathrm{12}}_{uv}+I^{\mathrm{21}}_{uv}+I^{\mathrm{22}}_{uv}$|⁠. The terms other than |$I^{\mathrm{22}}_{uv}(k)$| only ever have an impact if signal arises from haloes with M < Mmin – they are identically zero if W(Mmin, k) = 0 (e.g. for galaxies when Mmin is below the occupation threshold). The additional terms are typically only evaluated when ‘matter’ appears in a two-point function because significant matter exists in low-mass haloes. We checked that our calculation is only minimally sensitive to Mmin with results only changing by ∼2 per cent if we change our fiducial lower limit from 108 to either 106 or |$10^{10}\, h^{-1}\mathrm{M_\odot }$|⁠. APPENDIX B: ADDITIONAL REDSHIFTS In Fig. B1, we show the ratio of halo model predictions for matter–matter, matter–halo, and halo–halo power spectra measured from the multidark simulations. This is the same information as in the lower triangle of Fig. 4, but for z = 0.53 and z = 1 (the upper and lower triangles of Fig. B1, respectively). We see very similar results to those presented in Fig. 4 with the accuracy being nearly perfect for the halo–halo cross-power (as it should be). The halo autopower show slightly larger discrepancies, but these must arise from the one-halo term, which is telling us that the theoretical mass function must not be a perfect description of the simulations here. For the halo–matter power we see very encouraging results, especially when we allow ourselves to extrapolate the non-linear bias correction to low halo masses. The matter–matter power is also encouraging, although we see a slight tendency for the extrapolated non-linear bias model to overpredict the power in the quasi-linear regime, particularly at z = 1, which may be indicative of a short coming of our extrapolation. Figure B1. Open in new tabDownload slide As the lower triangle set of Fig. 4 but for z = 0.53 (the upper triangle) and z = 1 (the lower triangle). In each case the βNL function has been calculated from multidark data at the redshift in question. The red points with the error bars show the standard halo model prediction, while other coloured points show the new model prediction, either extrapolated (blue) or not (green). Figure B1. Open in new tabDownload slide As the lower triangle set of Fig. 4 but for z = 0.53 (the upper triangle) and z = 1 (the lower triangle). In each case the βNL function has been calculated from multidark data at the redshift in question. The red points with the error bars show the standard halo model prediction, while other coloured points show the new model prediction, either extrapolated (blue) or not (green). For both redshifts shown in Fig. B1, βNL has been re-measured from the simulation data at the redshift in question. We also looked at using the βNL function measured at z = 0 in making predictions at other z, but we found that this did not work as well (although the results were still okay), hinting that the correction has a significant redshift dependence. This in turn suggests that βNL would have cosmology dependence. This cosmology and redshift dependence is despite the fact that we have chosen to parametrize the correction as a function of ν (rather than M, or log M) and that we have defined βNL as in equation (23) where we might have hoped that the division by Plin would null the cosmology and redshift dependence. We leave any further investigation of this to the future. APPENDIX C: EXACT LINEAR HALO BIASING In this Appendix, we present a contrived scenario in which the standard halo-model calculation (equations 11 and 8; neglecting non-linear halo bias) is exactly correct. We consider this example to have great pedagogical value, despite its absurdity. We create a mock universe by first creating 643 ‘haloes’, which we distribute uniform randomly15 within a |$100\, h^{-1}\mathrm{Mpc}$| cube. We consider these haloes to make up all of the mass in our ‘universe’, and we give them each an isothermal density profile corresponding to Δv = 3, which we fill with 512 particles per halo, such that each halo is of equal mass. We then make a Gaussian realization of a linear displacement field, corresponding to z = 19 for our WMAP 5 cosmology, and displace the haloes accordingly, thus giving them some large-scale clustering. In our example, because we work at z = 19 in a |$100\, h^{-1}\mathrm{Mpc}$| cube with only 643 haloes, the displacements of each halo are very small compared to the mean interhalo separation. This makes the Zel’dovich approximation almost perfect, and generates a genuinely linearly biased sample of haloes (albeit with b = 1). We then measure the power spectrum of haloes–haloes and particles–particles and compare these to analytical halo-model calculations. In the case of particles–particles, we subtract the particle shot noise, but we do not do this for haloes–haloes since, in this case, the (halo) shot-noise has physical meaning – the haloes are the field. The upper panel of Fig. C1 shows the power spectra averaged over 10 realizations. On large scales, the halo and particle power spectra agree perfectly, but at smaller scales we see the particle power spectrum turn over relative to the halo spectrum, which we can attribute to the physical extent of the halo profiles. The advantage of making each halo so bloated (Δv = 3, compared to the 200 standard value) is that we can actually see the turn over in the one-halo term caused by the halo profiles, which would otherwise be at too small a scale. Intriguingly, we see no power deficit in the two- to one-halo transition region in our mock universe. This is confirmed in the lower panel, where we show the ratio of particle–particle model power to simulations, with the error bars transferred from the simulation measurements. We see that, in this case, the standard halo-model calculation is accurate at the few per cent level for all wavenumbers shown. Note that there is certainly no ∼30 per cent under prediction, which is what is typically seen in measurements from N-body simulations. Figure C1. Open in new tabDownload slide The points with the error bars in the upper panel show the average particle–particle (green) and halo–halo (purple) power spectra measured from 25 realizations of our mock universe, which is described in the text. The error bars show error-on-the-mean power in each k bin. The black lines shows a halo model calculation for the particle–particle power (solid), and this broken down into two-halo (dashed) and one-halo (dotted) terms. The lower panel shows the ratio of the model to the simulation for the particle–particle data, with the error bars transferred from the top panel. Here, we see that the standard halo model provides an almost perfect description of the data from the (absurd) mock universes. Figure C1. Open in new tabDownload slide The points with the error bars in the upper panel show the average particle–particle (green) and halo–halo (purple) power spectra measured from 25 realizations of our mock universe, which is described in the text. The error bars show error-on-the-mean power in each k bin. The black lines shows a halo model calculation for the particle–particle power (solid), and this broken down into two-halo (dashed) and one-halo (dotted) terms. The lower panel shows the ratio of the model to the simulation for the particle–particle data, with the error bars transferred from the top panel. Here, we see that the standard halo model provides an almost perfect description of the data from the (absurd) mock universes. The results presented in this Appendix demonstrate that the analytical halo-model calculation is exactly correct when compared to our mock universe, where all of the approximations that go into deriving it are satisfied. Haloes are linearly biased, exactly spherical with known mass and virial radius (which in our example are the same for all haloes, but we could have drawn these from a mass distribution function) with smooth profiles. This illustrates that there is nothing strange going on when the calculation breaks down when compared to realistic simulations (or the realistic Universe) – it is simply a reflection of the assumptions behind the calculation not being valid. As demonstrated in this paper, the major missing ingredient at intermediate scales between the standard two- and one-halo terms is the non-linearity of the halo bias. © 2021 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited. © 2021 The Author(s) Published by Oxford University Press on behalf of Royal Astronomical Society.
TI - Including beyond-linear halo bias in halo models
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stab748
DA - 2021-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/including-beyond-linear-halo-bias-in-halo-models-wiPVihaqFB
SP - 3095
EP - 3111
VL - 503
IS - 2
DP - DeepDyve
ER -