A robust measurement of the first higher-derivative bias of dark matter halos

TL;DR

This paper introduces a new simulation method to measure the first higher-derivative bias of dark matter halos, revealing a negative bias consistent with theoretical expectations and comparing well with other models.

Contribution

The authors develop a novel amplification technique for specific density modes in simulations, enabling precise measurement of higher-derivative halo bias.

Findings

The higher-derivative bias $b_{\nabla^2 \delta}$ is negative across all halo masses.

Measured bias roughly follows the $-R_L^2(M)$ relation, with $R_L$ as the Lagrangian radius.

Results agree with 1-loop power spectrum fits, previous studies, and peak theory predictions.

Abstract

We present a new simulation technique in which any chosen mode $\mathbf{k}$ of the density contrast field can be amplified by an amplitude $\Delta$. These amplified-mode simulations allow us to study the response of the halo density field to a long-wavelength mode other than the DC mode. In this sense they are a generalization of the separate-universe simulations to finite-wavelength modes. In particular, we use these simulations to obtain robust measurements of the first higher-derivative bias of dark matter halos $b_{\nabla^2 \delta}$. We find a negative bias at all mass considered, roughly following the $-R_L^2(M)$ relation, the Lagrangian radius of halos squared, as naively expected. We compare our results with those obtained from a fit to the 1-loop halo-matter power spectrum, as well as with the recent results from Abidi and Baldauf (2018), and to the prediction from the peak theory.

Figures

Equations

$${\delta }_h(\mathbf{x},\tau )=\sum _O{b}_O(\tau )O(\mathbf{x},\tau ),$$

$${\delta }_h(\mathbf{x},\tau )\supset b_1{\delta }_m(\mathbf{x},\tau ),$$

$$b_1(\tau ){\delta }_m(\mathbf{x},\tau )\to \int d^3\mathbf{y}F(\mathbf{y},\tau ){\delta }_m(\mathbf{x}+\mathbf{y},\tau ),$$

$$b_1(\tau ){\delta }_m(\mathbf{x},\tau )\to b_1(\tau ){\delta }_m(\mathbf{x},\tau )+b_{{\nabla }^2\delta }(\tau ){\nabla }^2{\delta }_m(\mathbf{x},\tau )+\cdots ,$$

$${\delta }_h(\mathbf{k})=[b_1-k^2{b}_{{\nabla }^2\delta }]{\delta }_m(\mathbf{k}),$$

$$\delta^{(1)}(\bm{x})=\delta^{(1)}(\bm{x})\big{|}_{\Delta=0}+\Delta{\rm Re}(\operatorname{e}^{i\bm{k}_{0}\cdot\bm{x}})=\delta^{(1)}(\bm{x})\big{|}_{\Delta=0}+\Delta\cos(\bm{k}_{0}\cdot\bm{x}).$$

$$\delta^{(1)}(\bm{k})=\delta^{(1)}(\bm{n}k_{F})=\delta^{(1)}(\bm{n}k_{F})\big{|}_{\Delta=0}+\frac{\Delta L^{3}}{2}(\delta^{K}_{\bm{n},\bm{m}}+\delta^{K}_{\bm{-n},\bm{m}}),$$

$${\delta }^{(1)}(\mathbf{x})=\Delta \cos (k_{0}x).$$

$${\delta }_h(\mathbf{x})=\frac{n_h(\mathbf{x})}{{\bar{n}}_h}-1=B_h\cos (k_{0}x),B_h=b(k_0)B_m=[b_1-b_{{\nabla }^2\delta }{k}_0^2+\cdots ]B_m,$$

$${\chi }^2=\sum _{x_i}\frac{1}{{\mathcal{N}}^2}{\left[n_h({\mathbf{x}}_i)-{\bar{n}}_h\left(1+B_h\cos (k_0{x}_i)\right)\right]}^2,$$

$$\frac{\partial }{\partial B_h}{\chi }^2(B_h)\stackrel{!}{=}0\iff \sum _{x_i}\cos (k_0{x}_i)\left[n_h({\mathbf{x}}_i)-{\bar{n}}_h\left(1+B_h\cos (k_0{x}_i)\right)\right]=0,$$

$${\hat{B}}_h=\frac{{\sum }_{x_i}\cos (k_0{x}_i)[n_h({\mathbf{x}}_i)-{\bar{n}}_h]}{{\bar{n}}_h{\sum }_{x_i}{\cos }^2(k_0{x}_i)}.$$

$$\sum _{x_i}{\cos }^2(k_0{x}_i)\to {\int }_0^L\frac{dx}{L}{\cos }^2(k_{0}x)=\frac{1}{2}.$$

$${\hat{B}}_h=\frac{2{\sum }_{\mathrm{halos}}\cos (k_0{x}_i)}{N_{\mathrm{halos}}},$$

$$\frac{{\hat{B}}_h^i(k_0,+\Delta ,M)-{\hat{B}}_h^i(k_0,-\Delta ,M)}{2\Delta }$$

$$\equiv b_1(M)-b_{{\nabla }^2\delta }^{\mathrm{eff}}(M)k_0^2,$$

$$\frac{P_{mm}(k)-P_{1-\mathrm{loop}}(k)}{P_{\text{L}}(k)}=-2C_{s,\mathrm{eff}}^2{k}^2,$$

$${\delta }^{(1)}(\mathbf{k})={\delta }_s^{(1)}(\mathbf{k})+\frac{1}{2}\Delta (2\pi {)}^3\left[e^{i\varphi }{\delta }_D(\mathbf{k}-{\mathbf{k}}_0)+e^{-i\varphi }{\delta }_D(\mathbf{k}+{\mathbf{k}}_0)\right],$$

$${\delta }_h({\mathbf{k}}_0)=$$

$$+\sum _{n=2}^{\infty }{\int }_{{\mathbf{p}}_1}\cdots {\int }_{{\mathbf{p}}_n}(2\pi {)}^3{\delta }_D({\mathbf{k}}_0-{\mathbf{p}}_{1\cdots n})F_n^{(h)}({\mathbf{p}}_1,\cdots ,{\mathbf{p}}_n){\delta }^{(1)}({\mathbf{p}}_1)\cdots {\delta }^{(1)}({\mathbf{p}}_n),$$

$$\widehat{\frac{{\rm d}\delta_{h}}{{\rm d}\Delta}}\equiv\frac{1}{2\Delta}\left[\delta_{h}(\bm{k}_{0})\Big{|}_{\Delta}-\delta_{h}(\bm{k}_{0})\Big{|}_{-\Delta}\right]\,.$$

$$\widehat{\frac{\mathrm{d}{\delta }_h}{\mathrm{d}\Delta }}=$$

$$+\frac{3}{2}{\int }_{\mathbf{p}}{P}_{\text{L}}(p)F_3^{(h)}(\mathbf{p},-\mathbf{p},{\mathbf{k}}_0)(2\pi {)}^3{\delta }_D(0)$$

$$+\frac{3}{8}(2\pi {)}^3{\delta }_D(0){\Delta }^2{F}_3^{(h)}({\mathbf{k}}_0,-{\mathbf{k}}_0,{\mathbf{k}}_0).$$

$$F_3^{(h)}(\mathbf{p},-\mathbf{p},{\mathbf{k}}_0)=$$

$$2\cdot \widehat{\frac{\mathrm{d}{\delta }_h}{\mathrm{d}\Delta }}({\mathbf{k}}_0)-3{\int }_{\mathbf{p}}{P}_{\text{L}}(p)F_3^{(h)}(\mathbf{p},-\mathbf{p},{\mathbf{k}}_0)=$$

$$\equiv$$

$$(3-3{)}^{\text{I}}\sim {\left({\int }_{\mathbf{p}}{F}_3^{(h)}(\mathbf{p},-\mathbf{p},\mathbf{k})P_{\text{L}}(p)\right)}^2\sim {\left(\frac{k}{k_{\mathrm{NL}}}\right)}^4,$$

$$(3-3{)}^{\text{II}}+(2-4)\sim {\int }_{\mathbf{p}}{F}_3^{(h)}(\mathbf{p},-\mathbf{p},\mathbf{k})P_{1-\mathrm{loop}}(p)\sim {\left(\frac{k}{k_{\mathrm{NL}}}\right)}^2,$$

$$(1-5)\sim {\int }_{\mathbf{p}}{\int }_{{\mathbf{p}}^{\prime }}{F}_5^{(h)}(\mathbf{p},-\mathbf{p},{\mathbf{p}}^{\prime },-{\mathbf{p}}^{\prime },\mathbf{k})P_{\text{L}}(p)P_{\text{L}}(p^{\prime })\sim {\left(\frac{k}{k_{\mathrm{NL}}}\right)}^2,$$

$$\frac{(k_0{R}_L{)}^4}{(k_0{R}_L{)}^2}=$$

$$P_{hm}^{1-\mathrm{loop}}(k)$$

$$+b_2{\int }_{\mathbf{p}}{F}_2(\mathbf{k}-\mathbf{p},\mathbf{p})P_{\text{L}}(p)P_{\text{L}}(|\mathbf{k}-\mathbf{p}|)$$

Full text

1 Introduction

The large-scale distribution of dark matter halos is one of the key ingredients of the theoretical description of large-scale structure (LSS). Since most observed tracers of LSS, such as galaxies, reside in halos, their statistics is determined by those of halos on large scales. In the context of perturbation theory, the statistics of halos are written in terms of bias parameters $b_{O}$ multiplying operators $O$ constructed out of the matter density field $\delta_{m}$ and the tidal field $K_{ij}$ (see [1] for a recent review)

[TABLE]

where $\delta_{h}$ is the fractional number density perturbation of a given halo sample.

Operators entering Eq. (1.1) can be divided into two categories: those which include exactly two net derivatives of the gravitational potential field $\Phi$, and higher-derivative operators which include four or more net derivatives of $\Phi$. Physically, these higher-derivative operators encapsulate the fact that halo formation involves the collapse of matter from a finite region in space, and thus, Eq. (1.1) cannot be completely local on all scales. Starting from a simple linear relation

[TABLE]

the way to incorporate the deviation from perfect locality of halo formation is to replace the local operators $\delta_{m}(\bm{x},\tau)$ with a functional [2, 3]

[TABLE]

where $F(\bm{y},\tau)$ is a kernel that is in general time dependent but has to be independent of $\bm{x}$ by homogeneity of the Universe. Performing a formal series expansion of $\delta_{m}$ around $\bm{x}$ leads to

[TABLE]

where statistical isotropy demands the absence of any preferred directions with which the derivative operators could be contracted. Hence the leading higher-derivative operator involves the Laplacian of $\delta_{m}(\bm{x},\tau)$, and the associated bias parameter has dimension $[{\rm length}]^{2}$. Its magnitude is expected to be of the order of $R_{*}^{2}$, where $R_{*}$ is the scale of the spatial support of the kernel $F(\bm{y},\tau)$, which we identify with the nonlocality scale of the tracer. For halos of mass $M$, this is given by their Lagrangian radius $R_{L}(M)$. In Fourier space, the term proportional to $b_{\nabla^{2}\delta}$ corresponds to a “scale-dependent bias” $-b_{\nabla^{2}\delta}k^{2}\delta$. However, let us emphasize that this is an expansion in powers of $k^{2}$ , rather than a general function $f(k)$, which is how the term “scale-dependent bias” is sometimes interpreted. To make this distinction clear, we will use the term higher-derivative bias throughout.

The peak model introduced by [4] allows for a theoretical prediction of $b_{\nabla^{2}\delta}$. Indeed, in this model, the real-space peak-matter 2-point correlation function can be derived in closed form in Lagrangian space, yielding analytic forms for the local Lagrangian bias $b^{L}_{1}$ and the peak bias $b_{01}^{L}$. We can then use a model for velocity bias (which arises from the correlation between linear velocities and density gradients, and reflects the fact that large-scale flows are more likely to be directed towards peaks than to be oriented randomly) to compute the Eulerian peak bias $b^{E}_{01}$ and, by taking into account the effect of the smoothing filter we can get an analytic prediction for the Eulerian higher-derivative bias.

On the measurements side, the first constraints on $b_{\nabla^{2}\delta}$ have been placed by studies testing the scale dependence of bias on large scales [5, 6]. More recent measurements include those of [7] who found a value of $3[R_{L}(M)]^{2}$ (but only ruled out $b_{\nabla^{2}\delta}=0$ at the $\sim 1.3\sigma$ level). On the other hand, [8] quote values for $b_{\nabla^{2}\delta}$ that are much smaller than $[R_{L}(M)]^{2}$. Very recently, [9] measured $b_{\nabla^{2}\delta}$ by fitting the 1-loop halo-matter power spectrum with one free parameter. They found a result consistent within errors with $-[R_{L}(M)]^{2}/\alpha$ with $\alpha\approx 2-3$. Thus, there is still large uncertainty on the magnitude of Eulerian higher-derivative bias for halos. This can be measured more easily in Lagrangian space, using either the halo-matter power spectrum [10, 11] or the projection method of [12, 13]. In particular, [10, 11] measured the so-called peak bias $b_{01}^{L}$, which contributes to $-b^{L}_{\nabla^{2}\delta}$ along with the leading contribution from the filtering kernel, and obtained $b_{01}^{L}\approx 2[R_{L}(M)]^{2}$ for halos with mass $M\geq 8\times 10^{12}h^{-1}M_{\odot}$ , with only a weak departure from the simple $[R_{L}(M)]^{2}$ scaling. However, as we explained above, it is not possible to relate Lagrangian higher-derivative biases to their Eulerian counterpart without using a model for the velocity bias.

In this paper, we propose a new technique to measure directly $b_{\nabla^{2}\delta}$ in Eulerian space using so-called amplified-mode simulations. The idea is to enhance a single mode $\bm{k}_{0}$ by adding a modulation $\Delta\cos(\bm{k}_{0}\cdot\bm{x})$ in the initial conditions for the density field of a gravitation-only N-body simulation, which, assuming linear growth, translates to $B_{m}\cos(\bm{k}_{0}\cdot\bm{x})$ at low redshift. This is a generalization of the separate universe simulations introduced in, e.g. [14], to non-DC modes with finite wavenumber. This enhancement amplifies the contribution of $b_{\nabla^{2}\delta}\nabla^{2}\delta_{m}$ in the bias expansion allowing for a clear detection of the linear higher-derivative bias parameter. More precisely, if the mode $\bm{k}_{0}$ is chosen to be small enough that linear theory still applies today, the same wavelength mode should be observed in the halo density field with a different amplitude, i.e. $\delta_{h}$ receives a contribution of the type $B_{h}\cos(\bm{k}_{0}\cdot\bm{x})$. Since the linear relation between $\delta_{m}$ and $\delta_{h}$ is given Fourier space by

[TABLE]

one can measure $b_{\nabla^{2}\delta}$ from a suite of amplified-mode simulations by measuring the ratio $B_{h}/\Delta$ for diverse values of $\bm{k}_{0}$ and fit a second order polynomial to this relation.

This paper is organized as follows: in section 2 we describe in more details the idea of amplified-mode simulations (section 2.1), how to estimate $b_{\nabla^{2}\delta}$ from them (section 2.2) and the higher-order corrections one needs to consider (sections 2.3–2.4). We present our set of simulations and shortly explain the halo finding procedure in section 3. Section 4 describes how to obtain the same results from the 1-loop power spectrum in perturbation theory (section 4.1) and reviews some aspects of the peak theory and how the higher-derivative bias can be computed using this model (section 4.2). Finally, section 5 presents and discusses our results, and we conclude in section 6. The appendices present some checks of our implementation of the simulations (Appendix A), detailed computation of various quantities in amplified-mode simulations (Appendix B), some considerations about measuring the effective sound speed for matter $C^{2}_{s,{\rm eff}}$ (Appendix C) and comparison between our amplified-mode technique and 1-loop power spectrum results (Appendix D). Throughout this paper we adopt the same cosmology as in [15], i.e. a flat $\Lambda{\rm CDM}$ cosmology with $\Omega_{m}=0.27$, $h=0.7$, $\Omega_{b}h^{2}=0.023$ and $\mathcal{A}_{s}=2.2\cdot 10^{-9}$.

2 Amplified-mode simulations

In this section, we introduce in more details the idea behind amplified-mode simulations. Here and in the following we drop the redshift argument from the equations since the results can be applied at any $z$.

2.1 Theoretical idea

The principle of amplified-mode simulations is to superimpose a plane wave of linear amplitude $\Delta$ and wavenumber $\bm{k}_{0}$ to the initial random density field \delta^{(1)}(\bm{x})\big{|}_{\Delta=0} coming from sampling the power spectrum in the absence of enhancement (we use the notation $\delta^{(1)}$ to denote the linear density field). Since the density field has to be real in configuration space it is given by

[TABLE]

Here, we have set the phase of the plane wave to zero. The value of the phase is arbitrary, given homogeneity and isotropy of the background. As mentioned above, the effect of the amplified mode is to amplify the contribution of $b_{\nabla^{2}\delta}\nabla^{2}\delta_{m}$ in the bias expansion, allowing for a clear detection of the higher-derivative bias factor. Since we are only interested in measuring the linear higher-derivative bias parameter, we choose the mode $\bm{k}_{0}$ to be on sufficiently large scales that linear theory still applies at redshift zero.

The implementation of the amplified mode in N-body simulations is straightforward since one only needs to modify the initial distribution of particles to incorporate the plane wave before running the simulation in a traditional way. Hence no modification needs to be done to the integration scheme. Using the fact that the density field is discrete and periodic in both configuration and Fourier space, the Fourier transform of Eq. (2.1) is

[TABLE]

where $k_{F}=(2\pi)/L$ is the modulus of the fundamental mode of the simulation box ($L$ is the one dimensional comoving box size), $\bm{n}$ a vector of integers, $\bm{m}$ another vector of integers such that $\bm{k}_{0}=\bm{m}k_{F}$ and $\delta^{K}_{\bm{n},\bm{m}}$ the Kronecker delta.111We see explicitly in this last expression that the reality condition on the configuration space density field involves the hermitianity of the Fourier space field with each of the two modes $\bm{m}$ and $-\bm{m}$ being amplified. Thus, after sampling the density field from the power spectrum in Fourier space, one simply needs to add a factor of $(\Delta L^{3})/2$ at the desired wavenumber and ensure hermitianity of the Fourier space density field. The positions and velocities of particles are then obtained from the 2LPT displacement field and the resulting distribution of particles given as an input to the cosmological simulation code. The integration scheme is then carried out using a standard N-body code without any further modifications to obtain the late-time particle distribution. In Appendix A we present some detailed tests to verify our implementation.

2.2 Estimating the higher-derivative bias

We present here an estimator for $B_{h}$ the amplified mode amplitude in the halo density field as well as the procedure to obtain $b_{\nabla^{2}\delta}$ from this estimator. We also discuss the higher-order corrections that we neglect in section 2.4. In this section, we absorb the small random contribution from sampling the power spectrum at $k_{0}$, P_{\text{L}}(k_{0})\big{|}_{\Delta=0}, in $\Delta$. Hence the Lagrangian real space density field is now given by

[TABLE]

The halo density field can be written in terms of the halo density $n_{h}$ as

[TABLE]

where, $B_{m}$ is the amplitude as measured in the late-time matter field and $\bar{n}_{h}$ is the mean halo density. We consider halos within a fixed mass range, but drop the mass argument since the results can be applied to any halo selection. Thus, we want to estimate $B_{h}$. Given the Gaussian nature of $\delta_{m}$ and $\delta_{h}$ at linear order, we can simply use a $\chi^{2}$ given by

[TABLE]

where $\bm{x}_{i}$ is the position of the $i^{th}$ halo, $x_{i}$ the $x$ component of $\bm{x}_{i}$ and $\mathcal{N}$ is a noise term which we assume to be constant in space (since we assume that the modulation of $n_{h}$ on the scale $k_{0}$ to be small). The least-squares estimator for $B_{h}$ is then

[TABLE]

which yields

[TABLE]

This can be implemented by simply summing over the halo positions $\bm{x}_{i}$, weighted by factors of $\cos(k_{0}x_{i})$. Further, if we assume that halos are approximately uniformly distributed (see the linear approximation above), the cosine-average over a constant vanishes, while the denominator yields

[TABLE]

We then obtain

[TABLE]

where $N_{\rm halos}$ is the total number of halos.

The same estimator can be applied to matter to obtain $\hat{B}_{m}$. One simply needs to replace $B_{h}\rightarrow B_{m}$, $B_{m}\rightarrow\Delta$ in the above derivation, and replace the halo density $n_{h}$ and total number of halos $N_{\rm halos}$ by the matter density $\rho_{m}$ and total number of particles $N$.

To get $b_{\nabla^{2}\delta}$ in practice we use the Fourier space relation Eq. (1.5) applied to our estimator. In order to partially cancel cosmic variance we actually compute the mean between the positive and negative amplitude results for each realization, i.e.

[TABLE]

where $i$ denotes the $i^{\rm th}$ realization, $C^{2}_{\rm s,eff}$ is the scaled effective sound speed of the dark matter fluid, and we use the results of [15] for the linear bias $b_{1}(M)$222It is crucial for the ratio in Eq. (2.10) to be computed for results from the same realization (i.e. the same seed for the random generation of the initial particle distribution) in order to cancel the effects of the random phase at $\bm{k}_{0}$.. We used the superscript “eff” to emphasize that we took the ratio with respect to the linear amplitude $\Delta$ and we neglected the $4^{\rm th}$ order term. We then simply fit a second-order polynomial in $k_{0}$ to this ratio to get an estimate for $b_{\nabla^{2}\delta}^{\rm eff}(M)$.

Finally, Eq. (2.10) involves the effective sound speed for matter $C^{2}_{\rm s,eff}$ that we will measure from the standard perturbation theory (SPT) 1-loop matter power spectrum prediction ($P_{\rm 1-loop}$) as

[TABLE]

where $P_{mm}$ is the matter power spectrum as measured in simulations and $P_{\text{L}}$ the linear one. We can then subtract the product $b_{1}C^{2}_{s,{\rm eff}}$ from $b_{\nabla^{2}\delta}^{\rm eff}$ to obtain $b_{\nabla^{2}\delta}(M)$.

2.3 Coupling between short and long wavelengths modes

In addition to what has been discussed in the previous section, there are nontrivial couplings between long and short wavelength modes that cannot be ignored. We present here the 1-loop calculation needed to compute them. We wish to stay concise and presents the details in Appendix B.

We want to compute the nonlinear matter and halo fields in the amplified-mode case. Throughout, we assume the infinite-volume limit. Then, the linear density field in the amplified-mode simulations in Fourier space is modified to

[TABLE]

where $\phi$ is the phase of the amplified-mode and we introduced the notation $\delta_{s}=\delta|_{\Delta=0}$ for shortness. The two Dirac delta functions ensure that the matter density field remains real, which requires $\delta^{(1)}(-\bm{k})=\delta^{(1)*}(\bm{k})$. In the following, we will set $\phi=0$ without loss of generality. We now consider the halo density field. In perturbation theory, the nonlinear halo density field at the mode $\bm{k}_{0}$ can be written as

[TABLE]

where $\int_{\bm{p}}\equiv\int d^{3}\bm{p}/(2\pi)^{3}$ and $F_{n}^{(h)}$ are the fully symmetrized kernels of the halo density field, which we will describe below. We now insert Eq. (2.12) into Eq. (2.13) and evaluate the result up to cubic order. The estimator applied to the halo density field is defined as the symmetric difference

[TABLE]

As we show in Appendix B, the final result for this quantity at cubic order and when averaging over many small modes (so that we can replace the small-scale modes with their ensemble average) is

[TABLE]

Now everything is multiplied by the same factor $(2\pi)^{3}\delta_{D}(\bm{0})$ (which simply gives $L_{\rm box}^{3}$ when restoring box normalization). We are interested in the terms in the first line. However, we see that there are further contributions from the second and third lines. The cubic kernel in the configuration $F_{3}^{(h)}(\bm{p},-\bm{p},\bm{k}_{0})$ is precisely what appears in the 1-loop halo power spectrum (e.g., [16]). There are two contributions: first, the cubic order matter kernel $b_{1}F_{3}$ multiplied by the linear bias. Second, there is a contribution from quadratic and cubic bias terms. We have

[TABLE]

In the limit $p\ll k_{0}$, $F_{3}^{(h)}(\bm{p},-\bm{p},\pm\bm{k}_{0})$ scales as $(k_{0}/p)^{2}$. The integrand in the second line of Eq. (B.10) peaks around $p\sim k_{\text{NL}}$, and hence we expect this contribution to be of order $(k_{0}/k_{\text{NL}})^{2}\Delta$, which is not negligible compared to $b_{\nabla^{2}\delta}k_{0}^{2}\Delta$ except possibly for the most massive halos. This kernel involves the bias combination $b_{K^{2}}+(2/5)b_{\rm td}$, which was recently measured by [17] using an optimal estimator for the trispectrum (see their figure 2 which shows $(5/2)b_{K^{2}}+b_{\rm td}$). We hence use their result multiplied by $2/5$ in Eq. (2.16). On the other hand, the last line in Eq. (2.15) is multiplied by an extra factor $\Delta^{2}$ which allows us to neglect it.

Hence the final expression allowing us to measure $b_{\nabla^{2}\delta}$ becomes

[TABLE]

2.4 Higher order corrections

In this section we discuss the various higher-order corrections that we do not take into account and we show that they can indeed be neglected.

First of all there is the obvious $\Delta^{2}$ contribution in the last line of Eq. (2.15). Since we choose for $\Delta$ a linearly extrapolated value of 0.05 at $z=0$ it is clear that this term can be neglected. As long as corrections of order $\Delta^{2}$ are negligible, the results for $b_{\nabla^{2}\delta}$ do not depend on the precise choice of $\Delta$.

Next, there are the 2-loop contributions that we neglected in our calculation at 1-loop order (terms up to cubic order in the linear density field) leading to Eq. (2.17). We note that showing these contributions to be negligible is also of importance for the estimator of $b_{\nabla^{2}\delta}$ using the 1-loop halo-matter power spectrum that we will present in section 4.1. There are several terms to consider and we will discuss them one after the other. They are shown in figure 1 of [18] and we will follow the same nomenclature as that reference. The first one is the $(3-3)^{\text{I}}$ contribution (i.e. two tadpoles). This is simply of the order of the $(1-3)$ contribution squared, i.e.

[TABLE]

where the last approximation is valid in the soft limit $\bm{k}\rightarrow 0$ (where $k_{\rm NL}$ is the nonlinear scale at a given redshift). In the same limit the (1-3) contribution scales as $(k/k_{\rm NL})^{2}$, and the $(3-3)^{\text{I}}$ one can hence be neglected. Then there are the $(3-3)^{\text{II}}$ and $(2-4)$ terms. These in fact simply sum up together to the integral of the 1-loop power spectrum

[TABLE]

in the soft limit again. The scaling of these terms is hence the same as the one of the 1-loop terms we consider. However, this contribution has to be renormalized with counterterms (in the same way as $C^{2}_{s,{\rm eff}}$ appears at 1-loop), and will not contribute to the 1-loop result in the end. Finally, there is the $(1-5)$ contribution

[TABLE]

in the double soft limit [19]. The scaling of this term is again the same as the one of the 1-loop terms. However, as in the case of the $(3-3)^{\text{II}}$ and $(2-4)$ terms, this contribution has to be renormalized with counterterms and will not contribute to the 1-loop result in the end. It can hence also be discarded. We conclude that all 2-loop corrections scale as $(k_{0}^{4}/k_{\rm NL})^{4}$, compared to the $k_{0}^{2}$ scaling of the desired higher-derivative contribution, provided that $k_{0}$ is sufficiently smaller than $k_{\rm NL}\sim 0.25\,h\,{\rm Mpc}^{-1}$ at $z=0$.

There is a further constraint on $k_{0}$ from the amplitude of the higher-derivative biases $\propto\partial^{4}\delta$ that we neglect. Estimating the latter to have a coefficient of order $R_{L}^{4}$ (recall that we expect $b_{\nabla^{2}\delta}\propto R_{L}^{2}$), we obtain

[TABLE]

Thus, by choosing a value of $k_{0}\ll 1/R_{L}$ and $k_{0}\ll k_{\rm NL}$, and making $\Delta$ sufficiently small to neglect the terms proportional to $\Delta^{2}$, the contribution $\propto b_{\nabla^{2}\delta}$ can be made to be the leading contribution.

Finally we note that technically, our measurement using the estimator $\hat{B}_{h}$ in Eq. (2.9) corresponds to a version of the “scatter-plot” technique to measure bias. That is, we calculate the weighted number of halos, corresponding to a plane-wave filter, for a range of values of the corresponding weighted matter density. As described in detail in Sec. 4.2 of [1], this technique exactly recovers the bias parameters relating the halo-matter moment (defined with the same filter) to the matter moments. Apart from the choice of filtering kernel, which is usually a spherical or cubic tophat filter but chosen to be plane-wave here, there is one further difference in our application of the technique: the value of the matter density is not random, but chosen deterministically as $\Delta$. This means that we cancel cosmic variance to leading order.

3 Simulations and halo finding

We present here the details of our set of simulations. We also provide a quick outline of the halo finding procedure.

We arbitrarily align the plane wave in the $x$ direction and choose $\Delta=0.05$ at redshift zero for the linear amplitude. We then run simulations where we amplify the modes $k_{0}=\{k_{F},2k_{F},3k_{F},4k_{F},5k_{F},8k_{F},10k_{F}\}$, and amplitude $\pm\Delta$ for each $k_{0}$ value. We choose a comoving box size $L=500\,h^{-1}{\rm Mpc}$ and number of particles $N=512^{3}$. These last two parameters are the same as for the “highres” set of simulations of [15] who computed the local bias parameters from separate universe simulations, and yield a mass resolution $m_{p}=7\cdot 10^{10}h^{-1}M_{\odot}$. Finally we ran 48 realizations of each simulation and initialized them with 2LPT at $z=49$. We refer to this set of simulations as amplified-mode simulations. In particular we refer to the fiducial set corresponding to no amplification as L500.

Furthermore, in order to cross-check our results with constraints from the power spectrum, we use another set of two simulations without amplified-mode, and with the same cosmological parameters but box size $L=2400\,h^{-1}{\rm Mpc}$ and $N=1536^{3}$ particles. This allows us to increase the signal-to-noise ratio on large-scales. The mass resolution in this set is $m_{p}=2.9\cdot 10^{11}h^{-1}M_{\odot}$ and we refer to it as L2400.

The halo finding procedure is same as the one used in [15]. Halos are identified at $z=0$, 0.5 and 1 using the spherical overdensity halo finder Amiga Halo Finder (AHF) [20, 21] with an overdensity threshold $200\rho_{m}$ for the halo definition ($\rho_{m}$ is the background density). We bin the mass range of halos in 11 tophat bins of width 0.2 in logarithmic scale centered from $\log M=12.55$ to $\log M=14.55$, where $\log$ is the base 10 logarithm. Hence the lowest mass bin is centered on halos with around 51 particles, with a lower limit around 40 particles. We refer the reader to [15] for more details and the justification of our choices.

Before moving on we shortly come back to the condition given by Eq. (2.21). For our simulation parameters, $10k_{F}=0.126\,h\,{\rm Mpc}^{-1}$ and $1/R_{L}$ is between 0.089 and 0.483 $h\,{\rm Mpc}^{-1}$ so that Eq. (2.21) is satisfied for $k\leq 8k_{F}(5k_{F})$ for objects of in bins $\log M=14.35(14.55)$ respectively and up to $10k_{F}$ for all less massive objects. This defines the range of $k_{0}$ values that we will use for the fit as a function of halo mass.

4 Other measurements and predictions

In this section we present how the same results can be obtained from the 1-loop power spectrum in SPT, as well as predictions from the peak model for $b_{\nabla^{2}\delta}$.

4.1 Power spectrum measurements

We start by describing how to measure the higher-derivative bias parameter from the 1-loop halo-matter power spectrum. This will provide a good cross-check of our results. This has already been done in [7, 9].

The one-loop halo-matter power spectrum is given by (see e.g. [1] and references therein)

[TABLE]

where we have neglected the stochastic contribution which is proportional to $k^{2}$ and expected to be smaller than the 1-loop order terms. Ref. [17] recently measured the bias parameters $b_{1}$, $b_{2}$, $b_{K^{2}}$, and $b_{\rm td}$ for the same cosmology as in this work, which allows us to fit the halo-matter power spectrum measured from simulations with a single free parameter to obtain a measurement of $b_{\nabla^{2}\delta}$.

In practice we use the L500 set of simulations to fit Eq. (4.1) up to $k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$ at all redshift. We choose this value for the same reason as for $C^{2}_{s,{\rm eff}}$ (see Appendix C) as well as to have a maximum $k$ roughly matching those of our amplified-mode simulations. We then follow the same procedure as outlined in [17] to obtain robust errorbars. Mainly, we first use a bootstrap technique to obtain errorbars on each data point as a function of the wavenumber $k$. We then use these errorbars to weight the points when fitting, and bootstrap the fit in order to obtain errorbars on the final result for $b_{\nabla^{2}\delta}$. We then use the two simulations of the L2400 set to obtain the final mean value of $b_{\nabla^{2}\delta}$ using errorbars scaled by the total effective volume (both on the data at each $k$ and on $b_{\nabla^{2}\delta}$ itself). This means that, knowing the $1\sigma$ error from the L500 set ([$\sigma(b_{\nabla^{2}\delta})]_{\rm L500}$), we infer the one for the L2400 set as

[TABLE]

where $V_{\rm L500}=48\cdot 500^{3}(\,h^{-1}{\rm Mpc})^{3}$ and $V_{\rm L500}=2\cdot 2400^{3}(\,h^{-1}{\rm Mpc})^{3}$. For more details and justification about this procedure, we refer the reader to section 3.1 of [17].

4.2 Prediction from peak theory

In this section, we introduce how the higher-derivative bias can be estimated from the peak model first introduced in [4]. Since the peak theory has already been extensively discussed in the literature, we refrain from giving a detailed description of this model here (we refer the reader to the original paper [4]). Notice that the apparition of a scale-dependent bias for peaks, as well as the concept of velocity bias (that we will introduce shortly) were first pointed out in [22, 23] and further studied in, e.g. [24, 11, 25]. We define the following spectral moments for a generic window function $W$

[TABLE]

where $P_{\text{L}}(k)$ is again the linear power spectrum, as well as the spectral shape parameter

[TABLE]

In the peak model, halos are in one-to-one correspondence with peaks of the Lagrangian density field. This assumption is expected to hold for halos with masses above a few $M_{*}$, where $M_{*}$ is the typical mass of halos that collapsed at redshift $z$. In Fourier space, the density of peaks $\delta_{\rm pk}$ is written in terms of the density field filtered on some scale $R$ as

[TABLE]

where $b^{E}_{10}=1+b^{L}_{10}=b^{E}_{1}$ is the local halo bias, $b^{E}_{01}$ contributes to $b_{\nabla^{2}\delta}$ and $\delta_{R}(\bm{k})=\delta(\bm{k})W(\bm{k})$. We choose for the filter the effective window function introduced in [26]

[TABLE]

with $W_{G}$ and $W_{TH}$ the gaussian and tophat filters respectively. Expanding $W$ in a Taylor series we get

[TABLE]

To get an expression for $b_{10}^{E}$ and $b_{01}^{E}$ we start from their Lagrangian counterpart given by

[TABLE]

where $\nu_{c}=\frac{\delta_{c}}{\sigma_{0}}$ with $\delta_{c}=1.686$ the critical threshold for collapse, $J=G_{1}(\gamma,\gamma\nu_{c})/G_{0}(\gamma,\gamma\nu_{c})$ is the mean peak curvature and

[TABLE]

with $f(x)$ the function defined in Eq. (A.15) of [4]. To obtain the Eulerian counterpart of Eqs. (4.8)–(4.9), we must take into account the so-called velocity bias, i.e. the fact that large-scale flows are more likely to be orientated towards peaks than in random directions. The velocity bias is defined in terms of the peak and linear matter displacement fields as [22, 23]

[TABLE]

where the displacement fields map the Lagrangian ($\bm{q}$) to Eulerian ($\bm{x}$) positions of matter particles and peaks of the density field

[TABLE]

with the subscript $i$ standing for peaks (pk) and matter (m), respectively. Integrating the continuity equation $\partial\delta_{\rm pk}(\bm{k})/\partial\tau=-\nabla\cdot\bm{v}_{\rm pk}(\bm{k})$ and evaluating the result in Fourier space then yields

[TABLE]

where $D(z)$ is the linear growth factor at redshift $z$ normalized so that $D(z_{0})=1$ for halos collapsing at redshift $z_{0}$. In these last two expressions it is important to understand that halo collapse is meant to happen at redshift $z_{0}$. Hence, to compute the Eulerian biases one should fix $z=z_{0}$ and compute the spectral moments Eq. (4.3) at this redshift (including the ones entering the Lagrangian biases), keeping $D=1$.

We can hence plug Eqs. (4.13)–(4.14) in Eq. (4.7) and obtain the higher-derivative bias parameter in the peak model

[TABLE]

with the dependence on the mass and the redshift being implicit.

5 Results and discussion

We now turn to presenting our results for the effective sound speed of matter, as well as for the leading higher-order derivative bias.

5.1 Effective matter sound speed

We show here results for the effective sound speed of matter using the 1-loop power spectrum for matter Eq. (2.11). To ensure that the 1-loop power spectrum accurately describes the matter power spectrum we perform the fit to $k_{\rm max}\ll k_{\rm NL}$, where $k_{\rm NL}$ is the nonlinear scale defined such that

[TABLE]

For our cosmology, $k_{\rm NL}=0.22,\,0.32$ and $0.47\,h/{\rm Mpc}$ at redshift $0,\,0.5$ and $1$ respectively. Hence we perform fits up to $k_{\rm max}=0.15\,h/{\rm Mpc}$ at all redshift (see Appendix C for a more detailed justification of this choice).

The relation corresponding to Eq. (2.11) at $z=0.0$ as a function of $k$ is presented in figure 1. The red points were obtained from the L2400 set of simulations where the 1-sigma errorbars have been rescaled from the ones obtained with the L500 set, as explained in section 4.1. Fitting a second order polynomial to this relation we determine $C^{2}_{\rm s,eff}=1.31\pm 0.06\,(\,h^{-1}{\rm Mpc})^{2}$ at $z=0$. We also get $C^{2}_{\rm s,eff}=0.65\pm 0.05\,(\,h^{-1}{\rm Mpc})^{2}$ and $C^{2}_{\rm s,eff}=0.33\pm 0.04\,(\,h^{-1}{\rm Mpc})^{2}$ at $z=0.5$ and $1$ respectively.

Several works presented results for $C^{2}_{s,{\rm eff}}$ measured in a similar way, including [27] who obtained $C^{2}_{s,{\rm eff}}=1.6\;({\rm Mpc}/h)^{2}$, and [8] who quote a value of $C^{2}_{s,{\rm eff}}=2.31\pm 0.02\;({\rm Mpc}/h)^{2}$. Ref. [18] also obtained results at various redshifts quoting, e.g. a value of $0.98\;({\rm Mpc}/h)^{2}$ at $z=0$. They furthermore provided a very detailed discussion about the effect of the 2-loop contributions that we neglect for this measurement. We elaborate on these results and show further tests about this quantity in Appendix C.

5.2 Halo higher-derivative bias

We now turn to results for $b_{\nabla^{2}\delta}(M)$. We first present the mean ratio in Eq. (2.17) (divided by $b_{1}$ to remove the mass dependent shift along the vertical axis) as a function of $k$ for three mass bins (color coded) in figure 2. The dashed lines on this figure show the best $k^{2}$ fits that are then used to determine $b_{\nabla^{2}\delta}$ together with the results from the previous section for $C^{2}_{\rm s,eff}$ as

[TABLE]

In order to respect the condition in Eq. (2.21) we use $k_{\rm max}=5\,(8)k_{F}$ as a maximum wavenumber for the fit at $\log M=14.35$ (14.55) respectively. For all lower mass bins we use the full $k$-range up to 10$k_{F}$. We checked the dependence of our results for $b_{\nabla^{2}\delta}$ on the maximum $k$ used for the fit, and found results consistent within $1\sigma$ errorbars, with of course a higher signal-to-noise ratio when including higher wavenumbers in the fit. This indicates that the 2-loop effects discussed in section 2.4 are sufficiently small to not significantly bias the measurement of $b_{\nabla^{2}\delta}$. In Appendix C, we confirm that $C^{2}_{\rm s,eff}$, which similarly receives 2-loop corrections, also does not show any significant scale dependence on the scales used for constraining $b_{\nabla^{2}\delta}$, $k<10k_{F}=0.126\,h\,{\rm Mpc}^{-1}$.

Figure 3 shows results for $b_{\nabla^{2}\delta}$ as a function of halo mass at redshift 0, 0.5 and 1 (blue circles, upper triangles and lower triangles respectively). We obtain a clear detection of a nonzero, negative bias at all mass, with the most precise constraints being at $z=0$. Interestingly, no strong evidence of a redshift dependence is seen in the relation $b_{\nabla^{2}\delta}(M)$. This is in strong contrast to other bias parameters, such as $b_{1}$, which are strongly redshift-dependent at fixed mass (they are approximately universal functions of $\sigma(M,z)$). For this reason, the fits described below are performed to the combined results at all redshifts. Notice that we do not have results for $\log M=14.55$ at $z=1$ since at this redshift the number of objects is too low to obtain a robust measurement.

We compare our results with the results from [9] who obtained their measurement from a fit to the 1-loop halo-matter power spectrum (as described in section 4.1) at $z=0$ only. We find consistent results within errors. The measurements reported here have smaller error bars and correspondingly higher signal to noise. We also obtained results from a fit to the 1-loop halo-matter power spectrum but with a very low signal-to-noise ratio. Hence we present them in figure 9 in Appendix D as a cross-check of our results. In addition we performed two fits of our combined results at all redshifts. Firstly, we fitted the relation $-R^{2}_{L}(M)/\alpha$ (the Lagrangian radius squared divided by some constant). Indeed, as explained in section 1, since $b_{\nabla^{2}\delta}$ is the first dimensionful bias parameter, we expect it to have a magnitude of the order of the only physical scale entering halo formation, i.e. $R^{2}_{L}$. We obtained a best fit value

[TABLE]

This relation is shown by the red dashed curve in figure 3. We see that this performs quite well, confirming the expected trend. However we obtained a reduced $\chi^{2}$ per degree of freedom $\chi^{2}=15$ which is quite poor; there are also deviations from a $R_{L}^{2}(M)$ scaling both at low and high masses. Hence we also performed a simple linear fit in $\log M$ shown by the black dotted curved. We obtained

[TABLE]

with a reduced $\chi^{2}=5.6$ hence performing much better (we checked that a quadratic fit does not improve the $\chi^{2}$). The physical interpretation of this logarithmic mass dependence however is unclear. Hence we caution against using this fitting relation outside the mass and redshift range probed here.

Finally, we compare our results with the peak prediction described in the previous section (magenta dashed-dotted line). This analytic prediction is the one that performs the worst quantitatively, especially at higher mass. However it is satisfying to see that this fairly simple model predicts the correct sign for this bias parameter, and a roughly correct overall amplitude. Moreover, it only predicts a weak redshift dependence of $b_{\nabla^{2}\delta}$ at fixed mass, which is confirmed by our results.

6 Conclusions

Amplified-mode simulations are a simple and efficient way to effectively modify the Laplacian of the matter density field $\delta_{m}$ on large scales by amplifying the real part of a Fourier mode at a given wavenumber $k_{0}$. This corresponds to adding a cosine wave in configuration space. We have performed various tests and computations in Appendix A and Appendix B to validate our implementation.

Using these simulations, and taking into account the nontrivial coupling between the amplified mode and short-wavelength modes, we have presented new measurements of the halo higher-derivative bias parameter $b_{\nabla^{2}\delta}$. This allowed us to obtain a clear detection of this parameter for all halo masses probed and up to $z=1$. Our results, which are the most robust to date, are in broad agreement with those of [9] with however a much higher signal-to-noise ratio.

We obtained a negative bias parameter at all mass with almost no redshift dependence in agreement with the analytic prediction from peak theory. Since $b_{\nabla^{2}\delta}$ has dimension of [length]2, it involves a physical length scale which, in addition to the nonlinear scale $k_{\rm NL}$, determines the range of scales on which rigorous perturbative approaches can be trusted. On physical grounds, we expect its magnitude to be of the order of the physical length scale in halo formation, the halo Lagrangian radius $R_{L}$. We indeed found $-R_{L}^{2}/2.5$ to be a reasonably good fit to our results. Our measurements thus confirm that nonlocal effects should only become important on scales $\lesssim R_{L}$. We also provided a simple empirical fit which performs better than the scaling with $R_{L}$ and that can be used quantitatively in future work in Eq. (5.4). It should only be trusted over the mass and redshift range probed here, however.

These measurements complete the ones of [15, 17], and we now have results for the complete set of bias parameters entering the 1-loop halo power spectrum and tree-level bispectrum for the cosmology adopted in this work at $z=0,$ 0.5 and 1. This enables future work such as studying the reach of perturbation theory by comparing, e.g. the halo-matter power spectrum as measured in simulation with the 1-loop SPT prediction with no free parameter. It will also allow us to study the stochasticity in halo formation in more detailed ways since we will be able to push the small scale-limit of the deterministic part (i.e. describe $\delta_{h}$ and its associated statistics accurately to higher $k$ values). Finally, it could also be helpful for future surveys. Indeed, in order to be able to use the observed distribution of discrete luminous tracers on smaller and smaller scales (thus increasing constraining power) to extract cosmological information, it is of crucial importance to use a robust bias model with as few free parameters as possible. The measurements of this work together with those of [15, 17], and the fitting functions provided, are precisely aiming towards this direction.

Acknowledgments

We thank Tobias Baldauf for useful discussions at the early stages of this project, Vincent Desjacques for his insight regarding the peak model, and Mehrdad Mirbabayi for his help regarding higher-order corrections. TL would like to thank KITP and MPA for their hospitality during some of this work. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1125915. FS acknowledges support from the Starting Grant (ERC-2015-STG 678652) “GrInflaGal” from the European Research Council.

Appendix A Simulation checks

This appendix presents a simple test to confirm that our implementation of the amplified-mode simulations is correct. We setup the initial condition of the simulation such that the final density field at redshift 0 is $\delta(\bm{k},z=0)=\Delta\,(L^{3}/2)\,\left[\delta_{D}(\bm{k}-\bm{k}_{0})+\delta_{D}(\bm{k}+\bm{k}_{0})\right]$, i.e. we set all modes of the density field to zero and only apply the mode amplification at $\bm{k}=\bm{k}_{0}$ (the factor $(L^{3}/2)$ comes from the fact that we want the amplitude to be $\Delta$ in real space and that we take discrete Fourier transforms, see Appendix B). We do this for a plane wave aligned with the $x$ axis, and various $\bm{k}_{0}$, initial redshift, amplitude $\Delta$ and box size $L$. We then output the density field in real space in order to check that we recover a cosine wave at the right wavenumber and with the right amplitude. This also allows us to verify that the redshift dependence of $\Delta$ is indeed given by $D(z)$ (the linear growth factor). We also output the results at redshift 0 both when running the full N-body code or when simply applying the 2LPT algorithm (since we work on scales that are still linear today we expect the two to be equal).

Results are presented in figures (4)–(5). In the first figure, we show the result of the 2LPT code used for the initial conditions at high redshift. The default setup is shown by the red crosses, i.e. output at redshift 49, box size $L=500\,h^{-1}{\rm Mpc}$, linear amplitude at $z=0$ $\Delta=0.05$ and fundamental mode $k_{F}$ being amplified. We then show results varying each of these parameters. The green dots show $\delta(x)$ when the mode $k_{0}=2k_{F}$ is amplified, the magenta squares are for an amplitude $\Delta=0.1$, blue triangles are at $z=99$, and the black circles are for a box size of $250\,h^{-1}{\rm Mpc}$.

The first point to note is that we do indeed retrieve a cosine wave at the required frequency, as shown by comparing the red and green curves. The amplitude of the wave is also correct since we obtain $\Delta(z=49)=0.0013$ and $\Delta(z=99)=0.00065$ when setting $\Delta(z=0)=0.05$ (red and blue curves respectively). This is expected as we work on linear scales and hence the amplitude should evolve with the growth factor $D(z)$. For our cosmology $D=0.76,\,0.02$ and $0.01$ at redshift 0, 49 and 99 respectively, and $D(z)/D(0)\cdot 0.05=0.0013(0.00065)$ at $z=49(99)$ respectively. We furthermore get $\Delta(z=49)=0.0026$ when setting $\Delta(z=0)=0.1$ (magenta curve) which is indeed twice the amplitude of the red curve. Finally, when dividing the box size by two and amplifying the fundamental mode we get the black circles which is also the correct result since the fundamental mode of this box corresponds to $2k_{F}$ for $L=500\,h^{-1}{\rm Mpc}$.

We now turn to figure 5, which shows the same as the red curve in figure 4 but at redshift 0. The red crosses show the result when the full N-body code is run. We see that we recover a cosine wave with the expected amplitude at redshift zero, $\Delta=0.05$, which is indicated by the black dashed line showing the linearly evolved density field. A small number of grid cells display outlying values, which are also present when generating the density field at redshift zero using the 2LPT code (blue points). The outliers predominantly happen in the low-density region and are presumably due to inaccuracies in the force calculation in a few specific locations. In addition, the simulation results show small oscillations at a high wavenumber close to the Nyquist frequency of the particle grid. Since they do not appear in the 2LPT results, these must be due to the N-body algorithm. We also see the same effect when outputting the results from Gadget-2 at $z=48$, close to the initial redshift. Notice that in this test setup, there are no small-scale initial density perturbations. In a cosmological N-body simulation, which has significant initial density perturbations on small scales, we expect this numerical artifact to become irrelevant.

Appendix B Density field in amplified-mode simulations at 1-loop order

In this appendix we present the detailed computation of the matter and halo density fields at 1-loop in amplified-mode simulations using SPT. We omit the redshift argument throughout since these computations are valid at all redshift.

We start with the matter density field $\delta_{m}$. This will be useful to compute the overdensity of halos later. We start again from the expression for the linear density field in the amplified-mode simulations Eq. (2.12) with $\phi=0$ without loss of generality. We then consider the nonlinear density field. In perturbation theory, it can be written as

[TABLE]

where $F_{n}$ are the fully symmetrized kernels of the matter density field. We now insert Eq. (2.12) into Eq. (B.1) and evaluate the result up to cubic order. This gives

[TABLE]

where we have used the symmetry of the kernels under permutation of their argument. We are interested in the response of the nonlinear matter field to the long-wavelength enhanced mode for which the estimator is defined as the symmetric difference

[TABLE]

By symmetry, only terms of Eq. (B.2) that are odd in $\Delta$ contribute to this. Since we typically choose $\Delta\sim 10^{-2}-10^{-1}$ we further neglect the terms that are cubic in $\Delta$ and we obtain

[TABLE]

We now average over many realizations. Terms with odd powers of $\delta_{s}^{(1)}$ average to zero, and are hence dropped as well, while terms with even powers of $\delta_{s}^{(1)}$ can be replaced by their ensemble average

[TABLE]

In the end, only terms that involve zero power of $\Delta$ and two powers of $\delta_{s}^{(1)}$ remain. For the first, there are terms in which the arguments of $\delta_{s}^{(1)}$ sum to zero, and those where they sum to $2\bm{k}_{0}$ when evaluating the expression at $\bm{k}_{0}$. The latter terms average to zero when considering many realizations since $\delta_{s}^{(1)}$ obeys Eq. (B.5) by definition, and we hence get

[TABLE]

Finally, evaluating this expression at $\bm{k}=\bm{k}_{0}$ leads to

[TABLE]

Now everything is multiplied by the same factor $(2\pi)^{3}\delta_{D}(\bm{0})$ (which simply gives $L_{\rm box}^{3}$ when restoring box normalization).

This last expression implies that it would be possible to measure $C^{2}_{s,{\rm eff}}$ from amplified-mode simulations by looking at the matter field response to the mode amplification. For simplicity, and in order to provide a more direct comparison with the literature on $C^{2}_{s,{\rm eff}}$, we measure this coefficient from the matter power spectrum. We however verified that Eq. (B.6) is respected in our simulations. This is shown in figure 6 at redshift 0 and for $k_{0}=k_{F}$. The blue crosses represent the l.h.s of Eq. (B.6) while the red line is the right one. We clearly see that the Dirac deltas on the r.h.s of the equation prevent modes at $k\neq k_{0}$ to respond to the amplification. More quantitatively, the value of the blue cross at $k=k_{F}$ is $0.9979\pm 5\cdot 10^{-4}$ while the red line is $0.9986\pm 2\cdot 10^{-4}$.

We now move to the halo density field and present details of the computation leading to Eq. (2.17) not shown in the main text. The computation follows the same line as for the matter case. We start from Eq. (2.14) and insert into it the nonlinear halo density field up to third order as given by Eq. (2.13). Following the same reasoning as to go from Eq. (B.3) to Eq. (B.6) yields the terms proportional to $\Delta$ (before dividing by $\Delta$ in the estimator). For the terms proportional to $\Delta^{3}$, we have

[TABLE]

The factor $\Delta^{3}/8$ comes from the last two terms in Eq. (2.12), while the square brackets contain the two terms in Eq. (2.14), taking into account that $(-\Delta)^{3}=-\Delta^{3}$. It is easy to see that only 3 combinations of $r,s,t$ contribute to teach term, so we obtain

[TABLE]

using the symmetries of the fully symmetrized $F_{3}$ kernel. Hence we get

[TABLE]

In the last line we have used the symmetry of the perturbation theory kernels under sign change of all momenta $F_{3}^{(h)}(-\bm{k}_{1},-\bm{k}_{2},-\bm{k}_{3})=F_{3}^{(h)}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$. Note that there are no odd contributions in $\Delta$ at second order, and that $\delta^{(1)}_{s}(\bm{p})\delta^{(1)}_{s}(-\bm{p})=|\delta_{s}^{(1)}(\bm{p})|^{2}$ is positive definite. When averaging over many small scale modes we can replace

[TABLE]

which finally leads to Eq. (2.15).

Appendix C On the effective sound speed of matter

This appendix presents further results and tests on $C^{2}_{s,{\rm eff}}$. In figure 7 we show the results for $C^{2}_{s,{\rm eff}}$ as a function of $k_{\rm max}$ used for the fit in Eq. (2.11). We also show the mean value inferred from results with $k_{\rm max}$ in the range $[0.063,0.15]\,\,h\,{\rm Mpc}^{-1}$. We decided to use the mean value in this range since at higher values of $k_{\rm max}$ the central value shifts systematically which is most likely due to higher-order effects (2-loop terms) that we do not take into account. As we discuss in the next paragraph, [18] already found these effects to be important. We hence obtain $C^{2}_{s,{\rm eff}}=1.31\pm 0.06\,(\,h^{-1}{\rm Mpc})^{2}$ at $z=0$.

We compared this to results the results from [27], [8], and [18]. These authors respectively quote values of $1.6$, $2.31\pm 0.02$, and $0.98\,(\,h^{-1}{\rm Mpc})^{2}$ at $z=0$. There are thus discrepancies between published results in the literature. In the case of [27] and [8], this could be due to the fact that they used a similar approach than us but with a single fit in the range $0.15-0.25(0.3)\,h\,{\rm Mpc}^{-1}$ respectively, where 2-loop terms induce a scale dependence in $C^{2}_{s,{\rm eff}}$. The most interesting result was found by [18] who found a result close to our one using the same technique as in this work with $k_{\rm max}=0.05\,h\,{\rm Mpc}^{-1}$. Their figure 6 is similar to our figure 7 at various redshifts and up to higher $k$. We found a very similar scale dependence of $C^{2}_{s,{\rm eff}}$ than them when we looked at similar $k_{\rm max}$. These authors showed that this is due to 2-loop terms that are not taken into account in our estimator Eq. (2.11), and they attribute to this the discrepancy with other works. On the other hand, as can be seen on figure 7, we do not see any strong dependence of $C^{2}_{s,{\rm eff}}$ with the maximum $k$ used for the fit up to $k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$ which is why we limit ourselves to this range. Notice that this is a bit larger than, but comparable to the maximum wavenumber used to measure $b_{\nabla^{2}\delta}$, which is $0.126\,h\,{\rm Mpc}^{-1}$ (or less, depending on mass). Some disagreement with values reported in the literature is also expected due to the fact that we use different cosmologies. However, the change in $C^{2}_{s,{\rm eff}}$ due to a change a cosmology should be mainly proportional to the amplitude of the linear power spectrum, i.e. $\sigma_{8}^{2}$, leading to expected differences only of order 10%.

Finally, we also look at the redshift dependence of $C^{2}_{s,{\rm eff}}$ in figure 8. In order to cancel divergences in perturbation-theory loop integrals, $C^{2}_{s,{\rm eff}}$ has to scale as $D^{4}(z)$, although the finite part which remains after the divergences are canceled and which is the parameter we are measuring could scale differently with redshift. In scale-free cosmologies (flat, matter-dominated cosmology with a power-law power spectrum), one expects the finite part to scale as $k_{\rm NL}^{-2}$ following dimensional reasoning. Indeed the $1/k_{\rm NL}^{2}$ relation reproduces the redshift dependence we find well, while the scaling $D^{4}$ is a bit too steep.

Appendix D Comparison with results from the 1-loop halo-matter power spectrum

Figure 9 presents a comparison between the measurements presented in section 5 obtained from amplified-mode simulations, and measurements obtained from a fit to the 1-loop halo-matter power spectrum with one free parameter, as presented in section 4.1. The results are at redshift 0. We see that we find a very good agreement between the two. However the signal-to-noise ratio is much higher for the results obtained from amplified-mode simulations which allow for a clear detection of $b_{\nabla^{2}\delta}$. Notice that we performed the fit of the power spectrum up to $k_{\rm max}=0.15\,h\,{\rm Mpc}^{-1}$, which roughly corresponds to $12k_{F}$, and we use the L2400 set of simulations. Errorbars are the 1$\sigma$ error obtained following the procedure outlined in section 4.1. We also emphasize that we use the CAMB linear power spectrum to compute the integrals multiplying $b_{2}$, $b_{K^{2}}$ and $\left(b_{K^{2}}+2/5\,b_{\rm td}\right)$, and use these as mean values (as well as the mean values presented in [17] for the bias parameters) when fitting Eq. (4.1). Cosmic variance is partially canceled by measuring $P_{\rm L}$ from the Zel’dovich density field at $z=99$ for each realization of L500, and using it in the first line of Eq. (4.1) (also in the computation of $P_{mm}^{1-{\rm loop}}$). Using $P_{\rm L}$ and the bias values as measured in each realization for each term in Eq. (4.1), and performing a fit realization by realization would allow to cancel more cosmic variance.