Including beyond-linear halo bias in halo models

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 $\sim2$ at scales $k\sim0.7\,h 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 under-prediction 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 auto 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


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., Seljak 2000;Peacock & Smith 2000), as well as making inroads in the understanding of the E-mail: alexander.j.mead@googlemail.com 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 et al. 2012;Hill & Spergel 2014;Ma et al. 2015;Padmanabhan et al. 2017;Hill et al. 2018;Wolz et al. 2019;Tanimura et al. 2019;Koukoufilippas et al. 2020).
It has been 20 years 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 onehalo 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., HM-CODE Mead et al. 2015, 2016, 2020b 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 < 10 hMpc −1 have been developed. Alternatively, Valageas & Nishimichi (2011), Mohammed & Seljak (2014), Seljak & Vlah (2015) and Philcox et al. (2020) have all improved predictions at quasi-linear scales by incorporating perturbation theory for the matter field, with per-cent level accuracy for k < 1 hMpc −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 mattermatter spectrum as a function of cosmology for k < 10 hMpc −1 come from 'emulation' schemes (e.g., Lawrence et al. 2010;Agarwal et al. 2012Agarwal et al. , 2014Lawrence 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 generalise 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 halomodel power spectrum and we demonstrate how to include nonlinear 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.

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 threedimensional cosmological 'fields', that 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 θ u (x), where x is comoving position and the label u stands for the field. We will be interested in power spectra between pairs of such fields, P uv (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 x i 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 k, as where the factor of 1/V ensures that the units of θ u,k are the same of those of θ u (x). If we make the variable change y = x − x i then we can write where we have defined 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 W u,i in terms of k = |k| and r = |r| and also ensures that W u,i is real. Let us now construct the power spectrum between two fields (which could be identical), θ u and θ v : 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 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) dM 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 to get We can apply the same reasoning to the two-halo term, where we convert both sums to integrals and we recognise the expectation of the complex exponential to be the power spectrum of the halo centres: P hh (M 1 , M 2 , k). Note that the functions W u and W v 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 where b(M) is the linear halo bias of haloes with mass M and P lin 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 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 limit 1 : 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 whereρ 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 θ m (M, r) = ρ m (M, r)/ρ where ρ m (M, r) is the halo matter density profile. We then note that W m (M, k → 0) = M/ρ, and it is usual to factorise this normalisation, such that W m (M, k) = MU m (M, k)/ρ, with U m (M, k → 0) = 1. We can then write equation (11) as and we see that P 2H uv (k → 0) = P lin mm (k → 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.

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 2019;Voivodic et al. 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 scaledependent, 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 renormalisation 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. Second, 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 quasiliner 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. 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 mattermatter power spectrum, but we note that it is not obvious that the same correction that works in HALOFIT for matter-matter applies more generally.
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 W n (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 W n (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 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 (M 1 , M 2 , k → 0) = 0, and also that it must obey symmetry with respect to mass arguments: β NL (M 1 , M 2 , k) = β NL (M 2 , M 1 , 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 The first term is standard, and is identical to equation (11) and the new content is captured by I NL uv (k) in the second term 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 realised 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 < ∼ 0.2 hMpc −1 the β NL function must also contain these perturbative corrections if it is to be a good, general model of clustering.

Relation to other non-linear bias definitions
Non-linear halo bias is often discussed in terms of the scaledependent function, b hh (M 1 , M 2 , k), defined by This is can be computed using our formalism (via equation 15) as if M 1 = M 2 there is an additional one-halo contribution (in this case called the shot noise) of 1/n h wheren 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: 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 nota- 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 h (M , k) = δ D (M − M )/n 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: Firstly, 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 mattermatter 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).

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 Nbody simulations. If we simply rewrite equation (15) as we can see that a sensible way to measure β NL is to measure P hh (M 1 , M 2 , k) and then to divide it by the linear power spectrum multiplied by the product of linear halo bias factors for M 1 and M 2 . Recall that P hh (M 1 , M 2 , k) is the cross-power spectrum measured between haloes of mass M 1 with those of mass M 2 . In real measurements we have to bin haloes in mass to measure this function.
If we consider the auto-spectrum (M 1 = M 2 ) 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 auto spectrum, 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 wheren 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 auto power spectra using our model (equation 8). Subtracting equation (24) can therefore be seen as a method for isolating the twohalo term from the halo-halo auto power measurement.

Simulations
We measure β NL (M 1 , M 2 , k) using data taken from the MUL-TIDARK simulation database 2 (Klypin et al. 2011;Prada et al. 2012;Riebe et al. 2013). We use data from the original MUL-TIDARK simulation, which simulated N = 2048 3 particles in a L = 1000 h −1 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 utilise haloes that have been identified via the ROCKSTAR phasespace algorithm of Behroozi et al. (2013) 3 and defined using the 'virial' criterion, with a spherical-overdensity threshold given by 360ρ at z = 0, which comes from spherical-collapse calculations (e.g., Bryan & Norman 1998) for the simulated WMAP 5 ΛCDM cosmology: Ω m = 0.27; Ω b = 0.0469; Ω v = 1 − Ω m ; h = 0.7; n = 0.95; σ 8 = 0.82. We prefer to use the virial halo definition of haloes because of results presented in Mead et al. (2020b), 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ρ 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 8 mass bins, which leads to 36 unique cross-combinations. The choice of 8 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 10 11.6 h −1 M . 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 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. 2020b) 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.

Constructing β NL
We compute P hh (M 1 , M 2 , k) via fast Fourier transform with 512 3 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 Ny /2 0.8 hMpc −1 ) in order to avoid the effects of aliasing. Following equation (23), we construct β NL (M 1 , M 2 , k) by dividing the measured halo-halo power spectra by the linear power spectrum. For k < 0.08 hMpc −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 > 0.08 hMpc −1 because non-linear contributions to the measured matter-matter field become important, so we use a smooth theory calculation from CAMB (Lewis et al. 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 < 0.08 hMpc −1 portion of the measurement of P hh (M 1 , M 2 , k)/P lin 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 MULTI-DARK for 4 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 > 0.08 hMpc −1 and we observe a prominent, positive detection of the function for k > 0.1 hMpc −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 r v 0.94 and 1.26 h −1 Mpc respectively, which correspond to the dips visible in the right-hand panel to within factors of ∼ π.
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 hMpc −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 10 10.4 , 10 12.5 , 10 14.2 , and 10 14.9 h −1 M 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 MUL-TIDARK 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 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). 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.

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 (M 1 , M 2 , 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(ν), normalised such that the integral over all ν gives unity, which is related to n(M) via 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: 360ρ at z = 0, asymptoting to 178ρ 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), which is truncated at the virial radius defined such that this encloses an average density of ∆ v (z) times the mean density: The virial radius is related to the halo-scale radius, r s via the massdependent concentration parameter c = r v /r s , 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 > ∼ 1 hMpc −1 . The lower row of Fig. 2 shows the I NL mm integrand, defined in equation (17), for the special case of the matter-matter power spectrum. For the scales and mass ranges shown, W m (M, k) M/ρ and .08 hMpc −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 > 0.08 hMpc −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 < 0.08 hMpc −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 > ∼ 0.08 hMpc −1 . For k ∼ 0.8 hMpc −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.  (17). This integrand is β NL (ν 1 , ν 2 , k)U m (ν 1 , k)U m (ν 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 U m (ν, 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 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 M ) 10.4, 12.5, 14.2, and 14.9. therefore the most significant change when going from β NL to the integrand for I NL mm is the suppression at high ν caused by the halomass function. Fig. 2 therefore shows the halo-mass ranges that give additional contributions in our two-halo term from the nonlinear halo bias.
When evaluating the non-linear bias correction in equation (17) we force the correction to be zero for k < 0.08 hMpc −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 < 0.08 hMpc −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 twoand 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 (W n (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 ∼ 0.1 hMpc −1 , and this boost can reach a factor of ∼ 2 at k ∼ 0.7 hMpc −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 > ∼ 0.7 hMpc −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 < ∼ 0.8 hMpc −1 and extrapolate it beyond this. In any case, at such small scales the one-halo term dominates the overall power spectrum. Fig. 4 is the main result of this paper, and shows halo-model calculations with and without including β NL . We show power for auto and cross combinations for three halo mass bins: 10 12.5 -10 13.0 , 10 13.0 -10 13.5 , 10 13.5 -10 14.0 h −1 M and the matter field, all taken from MULTIDARK. Diagonal panels show the auto spectra while off-diagonal panels show the cross spectra. Error bars show erroron-the-mean power in each k bin.

Halo-model power spectra
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 auto spectra: 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 under predictions in power of ∼ 40 per cent are almost perfectly cured when one includes the non-linear halo bias. This dramatic difference between the auto-and cross-spectra is due to the halo shot noise, which appears in the auto spectra 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 auto spectra by the relatively powerful one-halo term at smaller scales. In some ways, the success of our model for the halo auto-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 discrep-  . 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 (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. Error bars are errors-on-the-mean taken from power spectrum measurements from the simulations.
ancy indicates that more accurate halo mass function predictions 10 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 I NL 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 under-prediction 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 < k/ hMpc −1 < ∼ 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 ∼ 0.1 hMpc −1 , which could arise from β NL including some 'perturbative' type corrections, e.g., pre virialisation, but the data are quite noisy and we are unwilling to draw any firm conclusions about this. For k > 0.1 hMpc −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 under-prediction 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 second 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 twohalo 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. 2020a). On the other hand, we would expect a poorer performance of the 10 In future one could use the mass-function emulators of either McClintock et al. (2019b) or Bocquet et al. (2020).  Figure 5. 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 standard halo model when considering fields that arise from lowermass haloes. For example, Addison et al. (2013) note that including non-linear halo bias is necessary when modelling the auto spectrum of the cosmic infrared background (CIB), which would make sense given the ∼ 10 12.5 h −1 M peak halo mass for star formation efficiency (e.g., Viero et al. 2013;Planck Collaboration 2014;Maniyar et al. 2020). Neutral hydrogen (HI) 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 HI distribution.

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 matterhalo 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, lowmass 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 11 . BOLSHOI is similar to MULTIDARK, but is a smaller 250 h −1 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, this 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 ∼ 0.1 hMpc −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.

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 standard perturbation theory (SPT) is sub-per-cent level accurate for standard ΛCDM spectra for k < ∼ 0.08 hMpc −1 whereas linear theory provides a ∼ 2 per cent over prediction of power on these 'pre-virialisataion' 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., Foreman et al. 2016;Seljak & Vlah 2015;Philcox et al. 2020). Taking a schematic approach, assume that one believes that a good large-scale model of the quasi-linear mattermatter power is given by 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 NL mm (k → 0) = P cor mm (k)/P lin mm (k) (in the previous discussion this limit was assumed to be zero). Second, one could replace the factors of P lin mm (k) that appear on the right-hand side of equation (16) with P QL 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 NL mm (k → 0) = 0 to be enforced on scales where P QL 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 cor mm (k) to be given by one-loop SPT and redefining our biasing relative to this. 11 Note that here we have used the BDM halo catalogues from MULTIDARK and BOLSHOI because BOLSHOI has no public ROCKSTAR catalogue.  Figure 6. 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 realisations 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. Error bars show error-on-the-mean power in each k bin.
However, the problem we encountered is that while SPT provides a slightly improved prediction of power for k < 0.1 hMpc −1 when it does go wrong, it goes wrong quite spectacularly -grossly overpredicting the non-linear matter-matter power for k > ∼ 0.2 hMpc −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 NL 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 mattermatter 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.

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 MULTI-DARK halo catalogues using a simple halo-occupation distribution (HOD) prescription. Our simple HOD catalogues are generated 12 by taking 10 12.5 and 10 14 h −1 M as the minimum and maximum halo masses that can host galaxies. Haloes of 10 12.5 h −1 M 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 × 10 12.5 h −1 M 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 realisation 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 nonlinear 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 characterised, 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 > 0.1 hMpc −1 , which is perhaps not surprising given that this model lacks a physical basis and that two one-halo contributions now enter the calculation.

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 corrective 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 twoand 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 .
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 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 an important way: we retain the standard halo model apparatus on very large scales, where all power spectra are well described by the standard linear power multiplied by a linear biasing factor. This allows our non-linear biasing term to work as a simple addition to an otherwise standard two-halo term, essentially concentrating all new content in to a single new term. 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 cosmology dependence. 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 the emulator is not currently public. One could also use a set of simulations such as the QUIJOTE of Villaescusa-  to build an emulator for β NL from scratch. 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 may require the n-point version of β NL .
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, by 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 < ∼ k/ hMpc −1 < ∼ 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 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 twopoint 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 localisation 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.

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.
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 < M min . 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 If either W u or W v 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 M 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) in to 4 terms: which can all be evaluated over M ∈ [M min , ∞]. The final result is then given by I NL uv = I 11 uv + I 12 uv + I 21 uv + I 22 uv . The terms other than I 22 uv (k) only ever have an impact if signal arises from haloes with M < M min -they are identically zero if W (M min , k) = 0 (e.g., for galaxies when M min 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 M min with results only changing by ∼ 2 per cent if we change our fiducial lower limit from 10 8 to either 10 6 or 10 10 h −1 M .

APPENDIX B: ADDITIONAL REDSHIFTS
In Fig. B1 we show the ratio of halo model predictions for mattermatter, 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 (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 auto power 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.
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 parameterise 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 P lin 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 nonlinear 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 64 3 'haloes', which we distribute uniform randomly 15 within a 100 h −1 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 realisation 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 Mpc cube with only 64 3 haloes, the displacements of each halo are very small compared to the mean inter-halo 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 haloeshaloes 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 realisations. 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  Figure B1. As the lower-triangle set of Fig. 4 but for z = 0.53 (upper triangle) and z = 1 (lower triangle). In each case the β NL function has been calculated from MULTIDARK data at the redshift in question. Red points with error bars show the standard halo model prediction, while other coloured points show the new model prediction, either extrapolated (blue) or not (green).  Figure C1. Points with error bars in the upper panel show the average particle-particle (green) and halo-halo (purple) power spectra measured from 25 realisations of our mock universe, which is described in the text.
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 particleparticle 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.
standard value) is that we can actually see the turn over in the onehalo 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.
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.