A Robust BAO Extractor

TL;DR

This paper introduces a fast, analytical method to extract and model the nonlinear evolution of Baryonic Acoustic Oscillations (BAO) from power spectra, applicable to simulations and real data without fitting parameters.

Contribution

A novel BAO extractor operator that accurately models nonlinear effects across all scales without requiring fitted parameters from simulations.

Findings

UV effects are irrelevant for the oscillating part's evolution.

The method accurately reproduces BAO evolution down to redshift zero.

The approach eliminates the need for fitting functions or nuisance parameters.

Abstract

We define a procedure to extract the oscillating part of a given nonlinear Power Spectrum, and derive an equation describing its evolution including the leading effects at all scales. The intermediate scales are taken into account by standard perturbation theory, the long range (IR) displacements are included by using consistency relations, and the effect of small (UV) scales is included via effective coefficients computed in simulations. We show that the UV effects are irrelevant in the evolution of the oscillating part, while they play a crucial role in reproducing the smooth component. Our "extractor" operator can be applied to simulations and real data in order to extract the Baryonic Acoustic Oscillations (BAO) without any fitting function and nuisance parameter. We conclude that the nonlinear evolution of BAO can be accurately reproduced at all scales down to $z=0$ by our fast analytical method, without any need of extra parameters fitted from simulations.

Equations

$$({\delta }_{ab}{\partial }_{\eta }+{\Omega }_{ab}){\phi }_b^R(\mathbf{k},\eta )=I_{\mathbf{k},{\mathbf{q}}_1,{\mathbf{q}}_2}{e}^{\eta }{\gamma }_{abc}({\mathbf{q}}_1,{\mathbf{q}}_2){\phi }_b^R({\mathbf{q}}_1,\eta ){\phi }_c^R({\mathbf{q}}_2,\eta )-h_a^R(\mathbf{k},\eta ),$$

$${\phi }_1^R\left(\mathbf{k},\eta \right)\equiv {\mathrm{e}}^{-\eta }{\bar{\delta }}_R\left(\mathbf{k},\eta \right),{\phi }_2^R\left(\mathbf{k},\eta \right)\equiv {\mathrm{e}}^{-\eta }\frac{-{\bar{\theta }}^R\left(\mathbf{k},\eta \right)}{\mathcal{H}f},$$

$${\bf\Omega}=\left(\begin{array}[]{cc}1&-1\\ -\frac{3}{2}\frac{\Omega_{m}}{f^{2}}&\frac{3}{2}\frac{\Omega_{m}}{f^{2}}\end{array}\right).$$

$${\gamma }_{121}\left(\mathbf{p},\mathbf{q}\right)=\frac{\left(\mathbf{p}+\mathbf{q}\right)\cdot \mathbf{p}}{2p^2},{\gamma }_{112}\left(\mathbf{q},\mathbf{p}\right)={\gamma }_{121}\left(\mathbf{p},\mathbf{q}\right),{\gamma }_{222}=\frac{{\left(\mathbf{p}+\mathbf{q}\right)}^2\mathbf{p}\cdot \mathbf{q}}{2p^2{q}^2}.$$

$$e^{\eta }{\gamma }_{abc}({\mathbf{q}}_1,{\mathbf{q}}_2){\phi }_b^R({\mathbf{q}}_1,\eta ){\phi }_c^R({\mathbf{q}}_2,\eta )=e^{\eta }{\delta }_{a2}\tilde{\gamma }({\mathbf{q}}_1,{\mathbf{q}}_2){\phi }_2^R({\mathbf{q}}_1,\eta ){\phi }_2^R({\mathbf{q}}_2,\eta )$$

$$-i{k}^j\frac{{\bar{v}}_R^j({\mathbf{q}}_1)}{\mathcal{H}f}{\phi }_a^R({\mathbf{q}}_2,\eta ),$$

$$\tilde{\gamma }({\mathbf{q}}_1,{\mathbf{q}}_2)=\frac{({\mathbf{q}}_1\cdot {\mathbf{q}}_2{)}^2}{q_1^2{q}_2^2}-1.$$

$$e^{\eta }\frac{\mathbf{k}\cdot {\mathbf{q}}_1}{q_1^2}{\phi }_2^R({\mathbf{q}}_1,\eta ){\phi }_a^R({\mathbf{q}}_2,\eta ),$$

$$h_a^R(\mathbf{k},\eta )\equiv -i\frac{k^i{J}_R^i(\mathbf{k},\eta )}{{\mathcal{H}}^2{f}^2}{e}^{-\eta }{\delta }_{a2},$$

$$J_R^i(\mathbf{x},\eta )=J_{1,R}^i(\mathbf{x},\eta )+J_{\sigma ,R}^i(\mathbf{x},\eta ),$$

$$J_{1,R}^i(\mathbf{x},\eta )=\frac{1}{[1+\delta {]}_R(\mathbf{x},\eta )}[(1+\delta ){\nabla }^i\varphi {]}_R(\mathbf{x},\eta )-{\nabla }^i[\varphi {]}_R(\mathbf{x},\eta ),$$

$$J_{\sigma ,R}^i(\mathbf{x},\eta )=\frac{1}{[1+\delta {]}_R(\mathbf{x},\eta )}\frac{\partial }{\partial x^k}\left[[1+\delta {]}_R(\mathbf{x},\eta ){\sigma }_R^{ik}(\mathbf{x},\eta )\right].$$

$$P_{ab}^R(k)=⟨{\phi }_a^R(\mathbf{k}){\phi }_b^R(-\mathbf{k}){⟩}^{\prime },$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\partial_{\eta}P^{R}_{ab}(k)=\Bigg{[}-\Omega_{ac}P^{R}_{cb}(k)+e^{\eta}I_{{\bf{k}};{\bf{p}}_{1},{\bf{p}}_{2}}\gamma_{acd}({\bf{p}}_{1},{\bf{p}}_{2})B^{R}_{bcd}(k,p_{1},p_{2})-\langle h^{R}_{a}({\bf{k}})\varphi^{R}_{b}(-{\bf{k}})\rangle^{\prime}$$

$$\displaystyle+(a\leftrightarrow b)\Bigg{]},$$

$$B_{abc}^R(q_1,q_2,q_3)=⟨{\phi }_a^R({\mathbf{q}}_1){\phi }_b^R({\mathbf{q}}_2){\phi }_c^R({\mathbf{q}}_3){⟩}^{\prime }.$$

$$\Delta P_{ab}^R(k)\equiv P_{ab}^R(k)-P_{ab}^{R,ss}(k),$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\partial_{\eta}\Delta P^{R}_{ab}(k)=\Bigg{[}-\Omega_{ac}\Delta P^{R}_{cb}(k)+e^{\eta}I_{{\bf{k}};{\bf{p}}_{1},{\bf{p}}_{2}}\gamma_{acd}(k,p_{1},p_{2})\Delta B^{R}_{bcd}(k,p_{1},p_{2})$$

$$\displaystyle\qquad-\Delta\langle h_{a}^{R}({\bf{k}})\varphi^{R}_{b}(-{\bf{k}})\rangle^{\prime}+(a\leftrightarrow b)\Bigg{]},$$

$${\partial }_{\eta }\Delta B_{abc}^R(k,q,p)=$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\Bigg{[}-\Omega_{ad}\Delta B^{R}_{dbc}(k,q,p)-\Delta\langle h_{a}^{R}({\bf{k}})\varphi^{R}_{b}({\bf{q}})\varphi^{R}_{c}({\bf{p}})\rangle^{\prime}$$

$$+2e^{\eta }{\gamma }_{aef}(k,q,p)\left(\Delta P_{eb}^R(q)P_{fc}^{R,ss}(p)+P_{eb}^{R,ss}(q)\Delta P_{fc}^{R,ss}(p)+\Delta P_{eb}^R(q)\Delta P_{fc}^R(p)\right)$$

$$+e^{\eta }{I}_{\mathbf{k};{\mathbf{p}}_1,{\mathbf{p}}_2}{\gamma }_{ade}(k,p_1,p_2)\Delta T_{debc}^R({\mathbf{p}}_1,{\mathbf{p}}_2,\mathbf{q},\mathbf{p})$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!+{\mathrm{cyclic\;\;permutations\;\;of}}\;\;(a,{\bf{k}}),\;(b,{\bf{q}}),\;(c,{\bf{p}})\Bigg{]}\,,$$

$$P_{ab}^R(k)\simeq P_{ab}^{R,1-loop}(k)+\Delta P_{ab}^{R,1-loop}(k),$$

$$⟨h_a^R(\mathbf{k}){\phi }_b^R(-\mathbf{k}){⟩}^{\prime }={\alpha }^R(\eta )\frac{k^2}{k_m^2}{P}_{1b}^R(k;\eta ){\delta }_{a2},$$

$$\Delta \alpha (\eta )={\alpha }^R(\eta )-{\alpha }^{R,ss}(\eta ),$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\partial_{\eta}\Delta P_{ab}(k;\eta)=\Bigg{[}-\Omega_{ac}\,\Delta P_{cb}(k;\eta)$$

$$\displaystyle\qquad-\Delta\alpha(\eta)\frac{k^{2}}{k_{m}^{2}}\left[P^{1-loop}_{1b}(k;\eta)+\Delta P_{1b}(k;\eta)\right]\,\delta_{a2}+(a\leftrightarrow b)\Bigg{]}\,.$$

$$\displaystyle\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!e^{\eta}\;\int\frac{d^{3}q}{(2\pi)^{3}}\;\frac{{\bf{k}}\cdot{\bf{q}}}{q^{2}}\Bigg{[}\langle\varphi_{2}({\bf{q}})\varphi_{a}({\bf{k}}-{\bf{q}})\varphi_{b}(-{\bf{k}})\rangle^{\prime}+\langle\varphi_{a}({\bf{k}})\varphi_{b}(-{\bf{k}}+{\bf{q}})\varphi_{2}(-{\bf{q}})\rangle^{\prime}\Bigg{]}\,.$$

$$⟨{\phi }_2(\mathbf{q}){\phi }_a(\mathbf{k}-\mathbf{q}){\phi }_b(-\mathbf{k}){⟩}^{\prime }\simeq -e^{\eta }\frac{\mathbf{k}\cdot \mathbf{q}}{q^2}{P}^0(q)\left(P_{ab}(k)-P_{ab}(|\mathbf{k}-\mathbf{q}|)\right)+O\left({\left(\frac{q}{k}\right)}^0\right),$$

$$-2e^{2\eta }{\int }^{\Lambda (k)}\frac{d^{3}q}{(2\pi {)}^3}{\left(\frac{\mathbf{k}\cdot \mathbf{q}}{q^2}\right)}^2{P}^0(q)\left(P_{ab}(k)-P_{ab}(|\mathbf{k}-\mathbf{q}|)\right)$$

$$\simeq -2e^{2\eta }\frac{k^2}{(2\pi {)}^2}{\int }^{\Lambda (k)}dq{P}^0(q){\int }_{-1}^{1}dx{x}^2\left(P_{ab}(k)-P_{ab}(k-qx)\right)$$

$$=-2e^{2\eta }\frac{k^2}{(2\pi {)}^2}{\int }^{\Lambda (k)}dq{P}^0(q){\bar{P}}_{ab}^1(k;q)F^1(q{r}_{bao})$$

$${\bar{P}}_{ab}^n(k;q)\equiv \frac{{\int }_{-1}^{1}dx{x}^{2n}\left(1-\frac{P_{ab}(k-qx)}{P_{ab}(k)}\right)}{F^n(q{r}_{bao})}{P}_{ab}(k),$$

$$F^n(q{r}_{bao})\equiv {\int }_{-1}^{1}dx{x}^{2n}\left(1-\cos (q{r}_{bao}x)\right),$$

Full text

1 Introduction

The evolution of the large scale structure of the Universe (LSS) beyond the linear regime must be accurately described in order to test cosmological models with present and forthcoming observations. Besides N-body simulations, analytical tools have been investigated and developed in the last years (for a review, see [1]). It is well known that standard perturbation theory (SPT) [2] fails to achieve percent accuracy at low redshifts and at BAO scales, and smaller. The failure is due to the inaccurate description of the couplings of the intermediate scales with very large (IR) and very small (UV) ones.

At large scales, the dominant effect is given by long range displacements, of typical rms $\sim 6~{}\mathrm{Mpc/h}$, induced by velocity fields coherent on scales as large as $\sim 110~{}\mathrm{Mpc/h}$. This effect can be well described in Lagrangian Perturbation Theory, already at the level of the Zel’dovich approximation, but is badly reproduced by lowest order SPT. Methods to overcome this problem have been envisaged, and after some debate regarding the Galilean invariance of the proposed schemes, a consensus has emerged on the way to resum the leading IR effects at all SPT orders [3, 4, 5, 6].

The effect of very small scales on intermediate ones can be studied by looking at the “response function” , which quantifies the sensitivity of the nonlinear Power Spectrum (PS) at a given scale $k$ on a slight modification of the linear one at a different scale $q$. N-body simulations show a progressive screening of the effect of UV scales on intermediate ones, a behavior which is completely missed by SPT [7]. The UV failure of SPT is expected, as nonlinear effects become more and more relevant at small scales and, at a more fundamental level, it is based on the single-stream approximation, in which velocity dispersion and all higher order momenta of the Dark Matter (DM) particle distribution function are artificially set to zero.

The shortcomings of SPT in describing the effect of small scales can be overcome by methods such as the Coarse Grained Perturbation Theory (CGPT) [8, 9] or the Effective Field Theory of the Large Scale Structure EFToLSS [10, 11] (see also [12]), in which the “wrong” UV behavior of SPT is replaced by counterterms measured from N-body simulations. These procedures have been shown to increase the maximum wavenumber $k_{max}$ up to which they match the results from simulations from $k_{max}^{SPT}\sim 0.05\;\mathrm{h/Mpc}$ to, typically, $k_{max}\sim 0.25-0.3\;\mathrm{h/Mpc}$ at $z=0$. Moreover, in [9] it was shown that the cosmology dependence of these counterterms can be efficiently reproduced by SPT, so that one does not need to run a simulation for each different cosmological model.

In this paper we present a scheme in which the three physical effects discussed above are clearly identified and combined. The mode-coupling between nearby intermediate scales is taken into account by SPT, the effect of UV scales is treated via effective counterterms and the long range displacements are taken into account by IR resummations based on nonlinear consistency relations. Our scheme improves the Time Renormalization Group (TRG) of ref. [13] and the time-evolution equations of refs. [14, 15] by consistently including these effects while conserving the flexibility of these approaches in describing, for instance, models in which the scale factor is scale dependent, as in massive neutrino cosmologies [16].

By clearly identifying the different physical effects we are able to assess the dependence of a given observable on each of them. Therefore, we focus on BAO, and derive an evolution equation for the oscillating component of the PS alone. By using exact consistency relations first derived in [17, 18], we include the effect of long range displacements at all SPT orders. Moreover, the structure of the equation itself suggest a procedure to “extract” the oscillating part from a given nonlinear PS, as the part of the PS which evolves independently from the smooth one. This “extractor operator” can be also applied to the PS obtained in N-body simulations and in real data, to directly compare theory and observations. This is to be contrasted with the usual approach, in which to extract the BAO scale one usually fits the data with a function containing a number of nuisance parameters modelling the nonlinear evolution of the broadband spectrum and the damping of the BAO fluctuations.

Comparing our extracted BAO with simulations we can clearly show that the BAO evolution is largely insensitive to the details of UV physics: our IR resummed equation for the oscillating part accurately describes BAO at all scales and down to $z=0$, and the results do not change appreciably even by setting to zero the UV counterterms (in EFToLSS language, the “sound speed” and other coefficients). The latter are, on the other hand, crucial in order to reproduce the broadband spectrum, and we detail the procedure to be followed to systematically include them at higher orders.

The paper is organised as follows: in sect. 2 we introduce the evolution equations derived in refs. [8, 9] and introduce the formalism and notations, in sect. 3 we discuss the inclusion of UV effects and show the results for the broadband PS at the lowest nontrivial order in our approximations. In sect. 4 we focus on the BAO and derive an evolution equation for the oscillating part of the PS including large scale flows by means of consistency relations. In sect. 5 we define our “BAO extractor” and apply it both to our results and to the PS from N-body simulations, discussing the comparison between the two. Finally, in sect. 6 we summarise our results and present our outlook on future work. In A we describe the procedure for extracting the UV counterterms.

2 Evolution equation for the PS

Neglecting vorticity, the first two moments of the coarse-grained Vlasov equations [8, 9] give the continuity and Euler equations, which, in Fourier space, can be cast in compact form as

[TABLE]

where $\eta=\log D(\tau)$ and $I_{k;p_{1},\cdots,p_{n}}\equiv\int\frac{d^{3}p_{1}}{(2\pi)^{3}}\cdots\frac{d^{3}p_{n}}{(2\pi)^{3}}(2\pi)^{3}\delta_{D}({\bf{k}}-\sum_{i=1}^{n}{\bf{p}}_{i})$.

We have introduced the fields

[TABLE]

where we have also used the linear growth function $f\left(\eta\right)\equiv\frac{1}{\cal H}\frac{d\eta}{d\tau}$. Repeated indices are summed over $b,c=1,2$. The overbars indicate fields obtained by Fourier transforming $[\delta]_{R}({\bf{x}},\eta)$ and the divergence of $[(1+\delta)v^{i}]_{R}({\bf{x}},\eta)/[1+\delta]_{R}({\bf{x}},\eta)$, respectively, where $[\cdots]_{R}$ indicates filtering up to some scale $R$ (see sect. 3 and refs. [8, 9] for details).

The left hand side of this relation is the linearized part of the equation for the two dynamical modes, and it is characterized by the matrix

[TABLE]

The first term at the right hand side of eq. (1) encodes the mode-coupling between the fields. The only nonvanishing components of the vertex functions are

[TABLE]

The mode-coupling term can also be written as

[TABLE]

with

[TABLE]

When $q_{1}\ll k\simeq q_{2}$ the second term at the RHS of (5) singles out the leading contributions to the mode-coupling induced by long (and time-dependent) velocity modes, whose form is dictated by the Galilean invariance of the system [19, 17, 20]. Expressing the velocity through its divergence, it gives the (formally) IR divergent term

[TABLE]

where we have used the definition (2). In sect. 4 we will discuss the resummation of these effects at all orders. On the other hand, the residual vertex function, eq. (6), amounts to $\cos(\theta_{12})^{2}-1$, where $\theta_{12}$ is the angle between ${\bf{q}}_{1}$ and ${\bf{q}}_{2}$, and therefore it is never divergent, and moreover vanishes for ${\bf{k}}={\bf{q}}_{1}+{\bf{q}}_{2}\to 0$.

The last term on the right hand side of eq. (1) is the contribution from the short modes that have been integrated out in the coarse-graining procedure [8, 9],

[TABLE]

where the UV source $J_{R}^{i}({\bf{k}},\eta)$ is the Fourier transformed of

[TABLE]

which depend on the gravitational potential $\phi$ and the velocity dispersion $\sigma_{R}^{ik}$,

[TABLE]

Applying the equation of motion (1) to the (equal-time) PS,

[TABLE]

– where the prime indicates that we have divided by $(2\pi)^{3}$ times the overall momentum delta function – gives

[TABLE]

where we have omitted the $\eta$-dependence, and where the bispectrum is given by

[TABLE]

Before proceeding, we emphasize that the only approximation in the equation above is in the way we deal with the vorticity of the coarse-grained velocity field. This has two components: a microscopic one, related to UV scales smaller than $R$, and one induced by the coarse-graining procedure itself. While the first one is completely included in the source terms $h_{a}^{R}$, we deal with the second one at a perturbative level. While in this section we have set the second vorticity component to zero from the beginning, one can show, using the methods of [9], that including it perturbatively would give exactly the same equations as those considered in the next sections. The effect of vorticity on the PS was investigated in [21], where it was found to be negligible at all scales and redshifts of interest.

No other approximation has been imposed so far. In particular, we are not assuming the single stream approximation, as it is usually done in SPT and other semi-analytic methods. Eq. (1), and its PS counterpart, eq. (12), contain all the relevant physics: the effect of the UV scales on the intermediate ones, through the source $h_{a}^{R}$, the mode-coupling between the intermediate scales, through the vertex functions, and the IR displacements in the terms containing the vertex (7) for $q_{1}\ll k$. In the following, we will discuss how to deal with all these effects.

3 UV effects

The source term $h_{a}^{R}$ is responsible for all deviations from the single stream approximation and all the nonlinear effects occurring at small scales. It is therefore cleaner to consider the subtracted PS,

[TABLE]

where $P^{R,ss}_{ab}(k)$ is the PS computed in the single stream approximation. It solves an equation analogous to (12), in which all the quantities, including the $\langle h_{a}^{R}\varphi_{b}^{R}\rangle^{\prime}$ correlator are obtained in the single stream approximation, in practice, by considering SPT or other approximation schemes at some finite order. This correlator vanishes in the $R\to 0$ limit while its value at nonvanishing $R$ takes into account all nonlinear effects due to modes q~{}\mbox{\raisebox{-2.58334pt}{\stackrel{{\scriptstyle>}}{{\sim}}}}~{}1/R in the single stream approximation, to be subtracted from the fully nonperturbative correlator measured in N-body simulations. Working at finite $R$ is essential in practice, in order to extract the source terms from the simulation, but the final results should be independent on the value chosen for $R$. We discuss this point below.

The evolution equation for the subtracted PS therefore reads

[TABLE]

where $\Delta\langle h_{a}^{R}({\bf{k}})\varphi^{R}_{b}(-{\bf{k}})\rangle^{\prime}\equiv\langle h_{a}^{R}({\bf{k}})\varphi^{R}_{b}(-{\bf{k}})\rangle^{\prime}-\langle h_{a}^{R,ss}({\bf{k}})\varphi^{R,ss}_{b}(-{\bf{k}})\rangle^{\prime}$. The subtracted bispectrum, $\Delta B^{R}_{abc}(k,p_{1},p_{2})\equiv B^{R}_{abc}(k,p_{1},p_{2})-B^{R,ss}_{abc}(k,p_{1},p_{2})$, in turn, solves the equation

[TABLE]

where $\Delta\langle h_{a}^{R}\varphi^{R}_{b}\varphi^{R}_{c}\rangle^{\prime}\equiv\langle h_{a}^{R}\varphi^{R}_{b}\varphi^{R}_{c}\rangle^{\prime}-\langle h_{a}^{R,ss}\varphi^{R,ss}_{b}\varphi^{R,ss}_{c}\rangle^{\prime}$ and we have defined the deviation from the single stream approximation of the trispectrum, $\Delta T^{R}_{debc}\equiv T^{R}_{debc}-T^{R,ss}_{debc}$, which, in turn, solves an evolution equation which can be straightforwardly derived.

The system of coupled evolution equations solved by the subtracted correlation functions must be truncated at some order. However, unlike the original TRG proposal [13], where the evolution equations for the unsubtracted functions were considered, there is a clear hierarchy in these equations which leads to a natural criterium for the truncation. Indeed, we first notice that the sources of these equations are given by the $\Delta\langle h_{a}^{R}\varphi^{R}_{b_{1}}\varphi^{R}_{b_{2}}\cdots\varphi^{R}_{b_{n}}\rangle^{\prime}$ correlators, since, if they all vanish, the single stream approximation is exact. Moreover, the role of these differences between correlators is to replace the “wrong” behavior of the UV modes in the single stream approximation with the “correct” ones, encoded in the correlators measured, for instance in N-body simulations. However, while the two-point correlator $\langle h_{a}^{R,ss}\varphi^{R,ss}_{b_{1}}\rangle^{\prime}$ starts contributing to the UV loops for $P^{R,ss}_{ab}(k)$ at 1-loop order, $\langle h_{a}^{R,ss}\varphi^{R,ss}_{b_{1}}\varphi^{R,ss}_{b_{2}}\rangle^{\prime}$ does it only from 2-loop, since it corrects the single-stream bispectrum at 1-loop order, and so on. In summary, if we want to correct the UV-loops behaviour of the $l-$loop order PS, then we need consider only the correlators up to $\Delta\langle h_{a}^{R}\varphi^{R}_{b_{1}}\varphi^{R}_{b_{2}}\cdots\varphi^{R}_{b_{l}}\rangle^{\prime}$, and, correspondingly, only the first $l$ equations of the system.

At the lowest order, we will therefore have

[TABLE]

where $P^{R,1-loop}_{ab}(k)$ is the PS computed at 1-loop SPT with the linear PS filtered at the scale $R$, while $\Delta P^{R,-1loop}_{ab}(k)$ is obtained from eq. (15) with $\Delta B^{R}_{bcd}(k,p_{1},p_{2})=0$ and $\Delta\langle h_{a}^{R}({\bf{k}})\varphi^{R}_{b}(-{\bf{k}})\rangle^{\prime}$ subtracted at 1-loop.

In order to go to next order, one has to take into account the 2-loop PS and include eq. (16), with $\Delta T^{R}=0$ and the three point correlator $\langle h_{a}^{R}\varphi^{R}_{b_{1}}\varphi^{R}_{b_{2}}\rangle^{\prime}$ subtracted at 1-loop order.

As we anticipated above, the total PS, eq. (17), should not depend on the coarse graining scale $R$ apart from an overall dependence on the filter function. For instance, the density-density PS should depend on $R$ only through a $W[kR]^{2}$ factor, where $W[kR]$ is the filter function in Fourier space. The same holds for the PS computed in the single stream approximation at a finite loop order (see [9] for a 1-loop check). Therefore, one has to check that no spurious $R$-dependence is induced by the difference between the source correlator measured in simulations and that computed in SPT. Indeed, each of them, taken separately, has a strong $R$-dependence, as shown in the left panel of fig. 1, where we use the parameterization

[TABLE]

and we plot the quantities $\alpha^{R}(\eta)$, $\alpha^{R,ss}(\eta)$, and the difference between the two (for $b=1$). The simulations we use were presented already in [9], see A for details. They are based on the TreePM code GADGET-II [22]. They follow the evolution, until z = 0, of $n_{part}=512^{3}$ CDM particles within a periodic box of $L_{box}=512h^{-1}\mathrm{Mpc}$ comoving. The initial conditions were generated at z = 99 by displacing the positions of the CDM particles, that were initially set in a regular cubic grid, using the Zel’dovich approximation. We assumed the “REF” cosmological parameters of [9], namely, $\Omega_{m}=0.271$, $\Omega_{b}=0.045$, $\Omega_{\Lambda}=0.729$, $h=0.703$, $n_{s}=0.966$, and $A_{s}=2.42\times 10^{-9}$.

As we see, for the scales of interest, the $R$-dependence of the correlator measured in N-body simulations is mostly cancelled by the perturbative one already at 1-loop, so that a $R$-independent function,

[TABLE]

can be defined. The residual $R$-dependence, is given by two contributions: nonlinear effects from scales q~{}\mbox{\raisebox{-2.58334pt}{\stackrel{{\scriptstyle>}}{{\sim}}}}~{}1/R not captured by the 1-loop subtraction, and non-perturbative (that is, beyond single-stream) effects from scales q~{}\mbox{\raisebox{-2.58334pt}{\stackrel{{\scriptstyle<}}{{\sim}}}}~{}1/R. The magnitude of these effects decreases with the external momentum $k$, and, as we see from this plot, they are clearly subdominant for $k$ in the BAO range of scales. The complementary information to the left panel of Fig. 1 is given in the right panel of the figure, where the scale-dependence of the ratio between the correlator and the PS is given, and the $k^{2}$ dependence clearly emerges.

Notice that the $\Delta\alpha(\eta)$ function bears some analogies with the sound speed coefficient in the EFToLSS [11], however, with some key differences, that we now outline. First of all, the stress tensor in [11], whose divergence gives basically our source terms, eq.  (10), is expanded in terms of the linear filtered fields, and therefore the coefficients of this expansion are obtained from the cross-correlator of the sources with the linear fields. Here, on the other hand, the relevant cross-correlators involve the sources and the nonlinear fields $\varphi_{a}^{R}$, and therefore include effects, like short scale displacements and source-source correlators ( $\langle h^{R}_{a}h^{R}_{b}\rangle^{\prime}$, where the second $h^{R}_{b}$ is contained in the nonlinear evolution of the $\varphi_{b}^{R}$ field) which are not included at 1-loop order in EFToLSS. Moreover, we will directly measure these sources, and the coefficient $\Delta\alpha(\eta)$, from N-body simulation and then include the result into our evolution equations for the PS. While this procedure can be followed also in the EFToLSS, more often, in practical applications, one first derives an expression for the PS containing the “sound speed” and other counterterms as parameters to be fitted from the PS measured in simulations. In doing so, the physical meaning of these counterterms is less transparent, and the amount of “overfitting”, in order to get the PS right, is difficult to estimate.

In summary, we can choose $R$ in the plateau region, or equivalently, take the formal $R\to 0$ limit, and consider the evolution equation

[TABLE]

In fig. 2 we show the ratios between eq. (17) (orange line) and the PS computed with the Coyote interpolator of N-body simulations [23]. The agreement clearly improves over the 1-loop SPT result (blue line) showing that the UV correction represented by the $h_{a}^{R}$ sources plays a decisive role, already at the lowest order considered here, namely, correcting the UV of the 1-loop PS. At $k=0.1\,h\,{\rm Mpc}^{-1}$, the source correlator modifies the PS by $\simeq-0.6\%$ at $z=1$, by $\simeq-1.1\%$ at $z=0.5$, and by $\simeq-1.8\%$ at $z=0$. The agreement in the PS shape degrades at low redshifts, where higher loop orders should be taken into account along the lines discussed below eq. (16). However, in the following, we focus our attention on the residual BAO oscillations exhibited by the orange lines, which indicate that only improving the UV effects does not account for the BAO damping well enough. In the next section we will discuss how to deal with this issue.

4 IR resummation and BAO wiggles

The damping of the BAO wiggles in the PS, and of the corresponding peak in the correlation function is mainly caused by random long range displacements [24]. It is well known that such effects are badly reproduced in SPT at any finite order, while they are much better taken into account by the Zel’dovich approximation, which provides a resummation at all SPT orders (see for instance, [25]). In our approach, the effect of these long range displacements on an intermediate scale $k$ are encoded in the last term at the RHS of eq. (5) when the momentum of the velocity field is $q_{1}\ll k$. We will now discuss how they can be naturally resummed.

Indeed, if one applies again the equation of motion, eq. (1), to the correlator $\langle\varphi_{a}({\bf{k}},\eta)\varphi_{b}(-{\bf{k}},\eta)\rangle^{\prime}$, to get the evolution equation for the PS, then, the last term at the RHS of eq. (5) gives

[TABLE]

The consistency relations first derived in [17, 18] for the bispectrum give

[TABLE]

in the $q\ll k$ configuration. The consistency relation is depicted diagrammatically in fig. 3. The second term in (21) gives the same contribution.

Notice that, while $P^{0}(q)$ in eq. (22) is taken to be the linear PS, the other two are fully nonlinear. Moreover, the RHS vanishes at the leading order in $q/k$, that is, if one sets the argument of the second PS inside parentheses to $k$. This is in agreement with the form of the consistency relations derived in [17], which vanish if, as in this case, all the fields are taken at equal times. However, as we now show, the different arguments of the two PS is crucial when they have an oscillatory component.

Inserting (22) in (21) gives

[TABLE]

where we have inserted a UV cutoff \Lambda(k)~{}\mbox{\raisebox{-2.58334pt}{\stackrel{{\scriptstyle<}}{{\sim}}}}~{}k, in order to enforce the validity range of the consistency relation (22). We have defined

[TABLE]

for q~{}\mbox{\raisebox{-2.58334pt}{\stackrel{{\scriptstyle<}}{{\sim}}}}~{}k, with

[TABLE]

so that, $F^{1}(q\;r_{bao})=2\,(1-j_{0}(q\,r_{bao})+2j_{2}(q\,r_{bao}))/3$, where the $j_{n}(x)$ are the spherical Bessel functions.

If the nonlinear PS has an oscillatory component,

[TABLE]

with $A_{ab}(k)$ a smooth modulating function which damps the oscillations beyond the Silk scale, then eq. (24) returns the oscillatory component itself plus a smooth contribution,

[TABLE]

and other terms proportional to derivatives of $A_{ab}(k)$, which are suppressed since the oscillatory part is proportional to $\Omega_{b}/\Omega_{m}$. Therefore, eq. (23) gives

[TABLE]

with

[TABLE]

The exact value of the cut-off $\Lambda(k)$ (slightly) affects the amplitude of the evolved BAO’s but not their scale, as one can see by comparing the three red lines (dashed, solid, and dotted) in Fig. 7, in which we plot the results for the oscillating part of the PS obtained by multiplying the integrand in (29) by $\exp(-q/(c\,k)^{2})$ with $c=1/2,1,\infty$, respectively, and integrating in $q$ from [math] to $\infty$. It should be noted that the integrand of eq. (23), as far as the oscillatory part is concerned, is naturally cut-off at the Silk scale. Therefore, even removing the UV cutoff altogether (setting $c=\infty$), scales q~{}\mbox{\raisebox{-2.58334pt}{\stackrel{{\scriptstyle>}}{{\sim}}}}~{}k_{Silk}\simeq 0.12\,{\mathrm{h/Mpc}} do not contribute to the resummation.

If we now consider the equation for the PS derived in the previous section, and we add to it the IR resummation term, eq. (28), we have completed our goal: we have an evolution equation in which IR, intermediate, and UV scales are taken into account.

We note that the last line of eq. (23) has been obtained by multiplying and dividing the previous line by the function $F^{1}\left(q\,r_{bao}\right)$ specified in eq. (25). For PS of the form (26), the final line of eq. (23) simplifies further into (27), where the scale $r_{bao}$, which is the comoving sound horizon at recombination, appears. For a given cosmology, we can compute the $r_{bao}$ using eq. (6) of [26]. However, notice that the BAO extraction procedure defined in the next section is quite insensitive to the input $r_{bao}$ value, see discussion after eq. (36).

We can consider also the same equation for a “smooth” cosmology, in which the initial PS has no BAO feature. They will not be generated by the evolution equation itself. If we subtract this equation from the one for the real cosmology, we get an evolution equation for the oscillatory component of the PS, in which the $O\left(P^{nw^{\prime\prime}}_{ab}/r_{bao}^{2}\right)$ terms in (28) cancel out, which, splitting the oscillatory part as

[TABLE]

gives the two equations

[TABLE]

with initial conditions $P^{w,ss}_{ab}(k;\eta_{in})=P^{1-loop,w}_{ab}(k;\eta_{in})$, $\Delta P^{w}_{ab}(k;\eta_{in})=0$. If $\eta_{in}$ is taken early enough, we can approximate $P^{1-loop,w}_{ab}(k;\eta_{in})\simeq P^{0,w}(k)u_{a}u_{b}$, and the first equation has the analytical solution

[TABLE]

The total PS is obtained by adding the smooth component,

[TABLE]

with

[TABLE]

where $\Delta P^{nw}_{ab}(k;\eta)$ solves an equation analogous to eq. (20),

[TABLE]

5 Extracting the BAO

The solution for the total PS in which the effect of the IR displacement on the oscillating component has been taken into account by eqs. (31) are given by the green lines in fig. 2 (for which we took $\Lambda=\infty$ in eq. (29)). The residual BAO oscillations are greatly reduced with respect to the results obtained in sect. 3, while the performance on the overall PS shape (the “broadband”) is basically the same. This fact suggests that the evolutions of the two components of the PS, namely the smooth broadband shape and the oscillatory one, are governed by different physical effects: the mode-coupling with intermediate and UV scales for the former and long range displacements for the latter, with just a moderate amount of interference between them. To explore this possibility, we first look for a procedure to extract the oscillatory PS from a given PS, linear or nonlinear. Of course, this procedure has some degree of arbitrariness, as any small smooth function vanishing outside the BAO range of scales can be assigned either to the oscillatory or to the smooth part. However, the evolution equation itself suggests an optimal way to define such a procedure, as it extracts the oscillatory component which evolves (mostly) independently from the smooth one. Indeed, as we have seen, eq. (23) led us to define the wiggly component as in eq. (27), using eq. (24) (for $n=1$). Therefore we will consider the operation,

[TABLE]

where the parameter $\bar{\alpha}$, that we will call “range”, determines the width of the momentum interval, centered on $k$ and expressed in units of the BAO scale, over which we take the integral. The parameter $n$, on the other hand, is set to 1 in the evolution equation, and has no effect when the operation $R[P]$ is applied on theoretical PS’s, but it might be useful on real data, when different bins are measured with different errors. Moreover, notice that the quantity $r_{bao}$ does not affect the scale of the oscillations extracted from the PS via eq. (36), as it is shown explicitly in Fig. 4. In other terms, the $R[P](k;\bar{\alpha},n)$ procedure extracts the scale of the true oscillations contained in $P(k)$, that is the $r_{bao}$ of eq. (26), while the $r_{bao}$ parameter in eq. (36), which could also differ from the former, only affects the amplitude of the oscillations.

As indicated in (27), the ratio $R(k;\bar{\alpha},n)$ is “contaminated” by smooth terms coming from the non-oscillatory component of the PS, of order $\sim P^{\prime\prime}/(Pr^{2}_{bao})\sim 1/(kr_{bao})^{2}$. This is clearly seen in fig. 5, where we show the action of the operation $R[P^{0}]$ on the linear PS (black dashed lines), and on a smooth interpolation of it, in which BAO oscillations are absent $R[P^{0,nw}]$ (red dotted lines)111We thank Tobias Baldauf for sharing with us a code to extract the smooth PS.. The difference between the two (black solid lines) oscillates around the x-axis. In the following plots, we will always subtract the same $R[P^{0,nw}]$ from all the different $R[P]$, in order to visualize oscillatory behaviors along the horizontal axis and, at the same time, keep the extraction procedure as simple as possible.

From the expression above, we can define a family of estimators for the ratio above. Assuming the data are binned, we can write

[TABLE]

where the value of the maximum $|l|$ in the sum, $L(\bar{\alpha})$, is chosen such that

[TABLE]

Assuming that the errors on the PS at different bins are uncorrelated, the error on $\hat{R}[P](k_{m};\bar{\alpha},n)$ is given by

[TABLE]

In fig. 6 we show the comparison between our results and the N-body simulations of ref. [27], which are optimised for the PS on the BAO scales. The UV cutoff in eq. (29) has been implemented by multiplying the integrand by the gaussian exponential $\exp(-(q/k)^{2})$. To show the behavior of the discrete $R[P]$ (37), we have binned the N-body data in $\Delta k$ intervals, and we have assumed a relative error $\Delta P_{m}/P_{m}=0.01$ for the PS in each bin (the N-body curves shown in the figure are obtained from an interpolation of the binned data obtained from (37) and (39)).

In fig. 7 we show the sensitivity of our results to this choice, by removing the UV cutoff entirely (that is, by taking the gaussian exponential to 1). In the same plot, we show also the effect of neglecting the UV contribution, that is, of setting $\Delta P^{w}_{ab}(k)=\Delta P^{nw}_{ab}(k)=0$ in eqs. (30), (34). This corresponds to using the 1-loop SPT results for the broadband part of the PS, which, as shown in fig. 2, reproduces the broadband part quite badly. Nevertheless, the effect on the oscillating part is given by the difference between the red and the black lines, which is barely recognisable on this plot scale.

6 Conclusions

The results of this paper somehow pair with those of ref. [25] on the damping of the BAO peak in the correlation function. In that case, it was shown that the Zel’dovich approximation (and slights improvements of it) can account for the widening of the BAO peak which, rather than being a disturbance to be marginalised over, is a well controlled physical phenomenon which can be used to extract cosmological information. Thanks to our procedure to extract the oscillating part of the PS, one can discuss the same issue in Fourier space and show that the damping of the oscillatory component of the PS depends, in practice, only on the $\Xi$ function defined in eq. (29), and therefore carries information on the linear PS and the linear growth factor. These IR effects are orthogonal to the well known limitations of SPT in the UV, which are quantitatively taken into account by effective “sound speed”, and other counterterms, which we have shown to be basically irrelevant on the BAO evolution. As a result, the BAO oscillations are robustly reproduced at all scales and at all redshifts by means of the semi-analytical computation presented here, which requires only the 1-dimensional momentum integral defining the $\Xi$ function.

When extracting BAO from data one usually fits the total PS with a model containing a certain number of nuisance parameters (see, for instance, [28]). The BAO extractor introduced in (36), (37), contains no free parameter, ($\bar{\alpha}$ and $n$ define an extraction procedure, they are not parameters to be fitted to compare theory and simulations or data). Moreover, since the ratio $R$ depends on ratios of PS measured at nearby scales, it should be insensitive to most of the systematics (including bias): this will be explored in a forthcoming publication, where also the extension to redshift space will be investigated. The procedure can be equally applied to the reconstructed PS, obtained after the long range displacements have been undone [29, 30, 31, 32]: therefore it is not an alternative to reconstruction, but rather, it provides a parameter independent procedure to extract BAO information from reconstructed data.

As for the broadband part of the PS, we showed that a 1-loop SPT computation supplemented with just one UV counterterm gives results in agreement with N-body simulations up to $k_{max}\sim 0.4\;\mathrm{h/Mpc}$ for $z\geq 0.5$, rapidly degrading at lower redshifts. We have discussed how to systematically improve our approximation, by including higher order SPT corrections and more correlators between the UV sources and the density and velocity fields. The use of time-evolution equations considered here is particularly fit to deal with models beyond $\Lambda$CDM in which the boradband part of the PS carries a distinctive signature, like cosmologies with massive neutrinos or based on modified GR, as the scale-dependence of the growth factor can be directly implemented in the linear evolution matrix in eq. (3).

Acknowledgments

We thank Francisco Villaescusa-Navarro and Matteo Viel for discussions and for providing us with the N-body simulations from which the UV source terms have been extracted. We also thank M. Sato and T. Matsubara for providing us with the N-body simulations used in the comparisons in figs. 6. M. Pietroni thanks T. Baldauf for discussions. E. Noda and M. Pietroni acknowledge support from the European Union Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreements Invisible- sPlus RISE No. 690575, Elusives ITN No. 674896 and Invisibles ITN No. 289442. The work of M. Peloso was supported in part by DOE grant de-sc0011842 at the University of Minnesota.

Appendix A Computation of $\Delta\alpha\left(\eta\right)$

In this Appendix we briefly discuss how we computed the coefficient $\frac{\Delta\alpha\left(\eta\right)}{k_{m}^{2}}$ used in the TRG equations (31) and (35). We recall from equations (18) and (19) that it is obtained from the correlator between the source $h_{a}^{R}$ and the course grained fields $\varphi^{R}$ measured from N-body simulations, minus the corresponding quantity evaluated in single stream approximation.

We use the N body simulations of the ‘reference’ cosmology of [9] (see Section 4 of that work for details). We use the simulations to evaluate $\frac{\langle h_{2}^{R}\left({\bf{k}},\eta\right)\varphi^{R}_{1}\left(-{\bf{k}},\eta\right)\rangle^{\prime}}{\langle\varphi^{R}_{1}\left({\bf{k}},\eta\right)\varphi^{R}_{1}\left(-{\bf{k}},\eta\right)\rangle^{\prime}}$ at the redshifts $z=0,\,0.25,\,0.5,\,1,\,1.5,\,2,\,3,\,5$. We fix the filter scale at $R=2\,{\rm Mpc}\,h^{-1}$, which is the smallest value of $R$ at which the resulting $\Delta\alpha$ is in the plateau region visible in the left panel of Figure 2 (we verified that using larger values of $R$ does not change the final values of $\Delta\alpha$ significantly). We then evaluate the coefficient $\frac{\alpha^{R}}{k_{m}^{2}}$ at $k=0.1\,h/{\rm Mpc}$ (as visible from the right panel of Figure 2, the result does not significantly change if we use a diffierent but comparable scale).

From this quantity we subtract the analogous ratio $\frac{\langle h_{2}^{R}\left({\bf{k}},\eta\right)\varphi^{R}_{1}\left(-{\bf{k}},\eta\right)\rangle^{\prime}}{\langle\varphi^{R}_{1}\left({\bf{k}},\eta\right)\varphi^{R}_{1}\left(-{\bf{k}},\eta\right)\rangle^{\prime}}$, with the numerator evaluated in single stream approximation through a one loop computation, and with the Coyote Emulator [23] PS in the denominator. The one loop computation is performed following Appendix A of [9]. 222Specifically, we add up “22” and “13” contributions. The “22” contributions are obtained from the terms ${\cal C}^{h0}$ (only half of it) and ${\cal C}^{hh}$ written in eqs. (A.13) and (A.14) of [9]; the “13” conttibutions are obtained from eqs. (A.20)-(A.21)-(A.22)-(A.23) of [9]. The correlators written in [9] are given in terms of the quantity denoted as $\bar{\varphi}^{(h)}$ in that work. To obtain the correlators in terms of $h_{a}$ we formally performed the replacement $g_{ac}\left(\eta-s\right)\rightarrow\delta_{ac}\delta^{(D)}\left(\eta-s\right)$ in the expressions of [9].

From this subtraction, we obtain the coefficient $\frac{\Delta\alpha}{k_{m}^{2}}$ at the various redshifts listed above. In the ETRG system a continuous function $\alpha\left(\eta\right)$ is needed. We find that the time dependence is well reproduced by the three-parameter fit:

[TABLE]

where $\eta_{i}=\eta\left(z_{i}\right)$. This parametrization ensures that $\Delta\alpha\left(\eta\right)$ is continuous, and that it grows faster than the 1-loop single stream term at $z>z_{2}$ (as suggested by the fact that this subtracted quantity contains terms that are of higher order in perturbation theory. This behavior is in qualitative agreement with our highest redshift data). We perform a fit of the data at $z=0,\,0.25,\,0.5,\,1$ in a $\left\{\log\eta,\,\log\Delta\alpha\right\}$ plane, and obtain the best fit values ${\cal C}_{0}=8.11\,h^{-2}\,{\rm Mpc}^{2}$ and $n_{1}=1.48$. We then fit the data at $z=1,\,1.5,\,2$ to obtain the best fit value $n_{2}=2.09$ (see Figure 8). The parametrization (40), with these values for the three fitting parameters ${\cal C}_{0},\,n_{1},\,n_{2}$, is used in the TRG evolutions of the main text.