Responses in Large-Scale Structure

TL;DR

This paper develops a rigorous framework for power-spectrum responses to long-wavelength perturbations in cosmology, enabling precise predictions of non-Gaussian covariance and extending previous results to nonlinear scales.

Contribution

It introduces a formalism for power-spectrum responses as resummed vertices, generalizes previous results, and provides leading-order expressions and nonlinear extrapolations validated by simulations.

Findings

Response coefficients accurately predict covariance in simulations

Model agrees with simulations for small soft modes and large hard modes

Formalism can be extended to galaxy tracers

Abstract

We introduce a rigorous definition of general power-spectrum responses as resummed vertices with two hard and $n$ soft momenta in cosmological perturbation theory. These responses measure the impact of long-wavelength perturbations on the local small-scale power spectrum. The kinematic structure of the responses (i.e., their angular dependence) can be decomposed unambiguously through a "bias" expansion of the local power spectrum, with a fixed number of physical response coefficients, which are only a function of the hard wavenumber $k$. Further, the responses up to $n$-th order completely describe the $(n+2)$-point function in the squeezed limit, i.e. with two hard and $n$ soft modes, which one can use to derive the response coefficients. This generalizes previous results, which relate the angle-averaged squeezed limit to isotropic response coefficients. We derive the complete expression of first- and second-order responses at leading order in perturbation theory, and present extrapolations to nonlinear scales based on simulation measurements of the isotropic response coefficients. As an application, we use these results to predict the non-Gaussian part of the angle-averaged matter power spectrum covariance ${\rm Cov}^{\rm NG}_{\ell = 0}(k_1,k_2)$, in the limit where one of the modes, say $k_2$, is much smaller than the other. Without any free parameters, our model results are in very good agreement with simulations for $k_2 \lesssim 0.06\ h/{\rm Mpc}$, and for any $k_1 \gtrsim 2 k_2$. The well-defined kinematic structure of the power spectrum response also permits a quick evaluation of the angular dependence of the covariance matrix. While we focus on the matter density field, the formalism presented here can be generalized to generic tracers such as galaxies.

Equations

$$⟨\delta (\mathbf{k})\delta ({\mathbf{k}}^{\prime })⟩=(2\pi {)}^3{\delta }_D(\mathbf{k}+{\mathbf{k}}^{\prime })⟨\delta (\mathbf{k})\delta ({\mathbf{k}}^{\prime }){⟩}^{\prime }=(2\pi {)}^3{\delta }_D(\mathbf{k}+{\mathbf{k}}^{\prime })P_m(k).$$

$$⟨\delta (\mathbf{k})\delta ({\mathbf{k}}^{\prime })\delta ({\mathbf{p}}_1)\delta ({\mathbf{p}}_2)\cdots \delta ({\mathbf{p}}_n){⟩}_c,\text{mbox}with{p}_i\ll k,k^{\prime }(i=1,2,\cdots ,n)\text{mbox}and{p}_{\mathrm{12..}n}\ll k,k^{\prime }.$$

$$\raisebox{0.0pt}{\parbox[c]{53.4798pt}{\includegraphics[scale={0.8}]{diag_Pm.pdf}}}\equiv P_{m}(k,t)\>(2\pi)^{3}\delta_{D}(\bm{k}+\bm{k}^{\prime})\,.$$

$$\lim_{p\to 0}\left(\raisebox{0.0pt}{\parbox[c]{53.4798pt}{\includegraphics[scale={0.8}]{diag_R1def.pdf}}}\right)\equiv\frac{1}{2}\mathcal{R}_{1}(k;\mu_{\bm{k},\bm{p}};t)P_{m}(k,t)\>(2\pi)^{3}\delta_{D}(\bm{k}+\bm{k}^{\prime}-\bm{p})\,.$$

$$\lim_{p\to 0}\left(\raisebox{0.0pt}{\parbox[c]{86.0013pt}{\includegraphics[scale={0.8}]{diag_sqB.pdf}}}+(\bm{k}\leftrightarrow\bm{k}^{\prime})\right)=\mathcal{R}_{1}(k;\mu_{\bm{k},\bm{p}})P_{m}(k)P_{\rm L}(p)\>(2\pi)^{3}\delta_{D}(\bm{k}+\bm{k}^{\prime}+\bm{p})\,,$$

$$\displaystyle B^{\rm tree}(\bm{k},\bm{k}^{\prime},\bm{p})=2\Big{[}F_{2}(\bm{k},\bm{p})P_{\rm L}(k)+(\bm{k}\leftrightarrow\bm{k}^{\prime})\Big{]}P_{\rm L}(p)+2F_{2}(\bm{k},\bm{k}^{\prime})P_{\rm L}(k)P_{\rm L}(k^{\prime}),$$

$${\mathcal{R}}_1^{\text{tree}}(k;\mu )$$

$$\underset{p\to 0}{\lim }B(\mathbf{k},{\mathbf{k}}^{\prime },\mathbf{p})=\underset{p\to 0}{\lim }⟨\delta (\mathbf{k})\delta ({\mathbf{k}}^{\prime })\delta (\mathbf{p}){⟩}_c^{\prime }={\mathcal{R}}_1(k,{\mu }_{\mathbf{k},\mathbf{p}})P_m(k)P_{\mathrm{L}}(p),$$

$$\displaystyle\lim_{\{p_{a}\}\to 0}\left(\raisebox{0.0pt}{\parbox[c]{86.0013pt}{\includegraphics[scale={0.8}]{diag_Rndef.pdf}}}\right)=\>$$

$$\times (2\pi {)}^3{\delta }_D(\mathbf{k}+{\mathbf{k}}^{\prime }-{\mathbf{p}}_{1\cdots n}),$$

$$\max \left\{\frac{p^2}{k^2},\frac{p^2}{k_{\text{NL}}^2}\right\}\text{mbox}wherep\equiv \max \{p_a{\}}_{a=1,\cdots n}.$$

$$\displaystyle\lim_{\{p_{a}\}\to 0}\left(\raisebox{0.0pt}{\parbox[c]{76.6062pt}{\includegraphics[scale={0.8}]{diag_sqnp2.pdf}}}+(\text{perm.})\right)$$

$$=n!{\mathcal{R}}_n(k;\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\})P_m(k)$$

$$\propto P_m(k)\underset{n\text{soft power spectra}}{\underbrace{P_{\mathrm{L}}(p_a)P_{\mathrm{L}}(p_b)P_{\mathrm{L}}(|{\mathbf{p}}_{ab}|)\cdots P_{\mathrm{L}}(|{\mathbf{p}}_{c\cdots d}|)}}.$$

$$\displaystyle\raisebox{0.0pt}{\parbox[c]{125.7498pt}{\includegraphics[scale={0.8}]{diag_sqnp2Rm.pdf}}}.$$

$$\displaystyle\raisebox{0.0pt}{\parbox[c]{125.7498pt}{\includegraphics[scale={0.8}]{diag_sqTR1.pdf}}}+(\mbox{perm.})$$

$$={\mathcal{R}}_1(k;{\mu }_{\mathbf{k},{\mathbf{p}}_{12}})P_m(k)\left[2F_2(-{\mathbf{p}}_{12},{\mathbf{p}}_2)P_{\mathrm{L}}(|{\mathbf{p}}_{12}|)P_{\mathrm{L}}(p_2)+({\mathbf{p}}_1\leftrightarrow {\mathbf{p}}_2)\right]$$

$$\times (2\pi {)}^3{\delta }_D(\mathbf{k}+{\mathbf{k}}^{\prime }+{\mathbf{p}}_{12}).$$

$$\propto P_{\mathrm{L}}(k)P_{\mathrm{L}}(k^{\prime })P_{\mathrm{L}}(p_2)\cdots P_{\mathrm{L}}(p_n),$$

$$\propto P_{\mathrm{L}}(k)P_{\mathrm{L}}(|\mathbf{k}+{\mathbf{p}}_a|)P_{\mathrm{L}}(|\mathbf{k}+{\mathbf{p}}_{ij}|)P_{\mathrm{L}}(p_3)\cdots P_{\mathrm{L}}(p_n),\text{mbox}andsoon,$$

$$⟨\delta (\mathbf{k},t)\delta ({\mathbf{k}}^{\prime },t)\delta ({\mathbf{p}}_1,t_1)\cdots \delta ({\mathbf{p}}_n,t_n){⟩}_{c,{\mathcal{R}}_n}^{\prime }$$

$$=n!{\mathcal{R}}_n\left(k,t,\{t_a\};\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\}\right)P_m(k,t)\prod _{a=1}^n{P}_{\mathrm{L}}(p_a,t_a),$$

$${\mathcal{R}}_n\left(k,t,\{t_a\};\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\}\right)={\mathcal{R}}_n\left(k,t;\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\}\right).$$

$$\displaystyle\mathcal{R}_{n}\big{(}k;\{\mu_{\bm{k},\bm{p}_{a}}\},\{\mu_{\bm{p}_{a},\bm{p}_{b}}\},\{p_{a}/p_{b}\}\big{)}=\frac{1}{n!\,P_{m}(k)}\frac{{\rm d}^{n}P_{m}(\bm{k}|\delta^{(1)}({\bm{p}_{1}}),\cdots,\delta^{(1)}({\bm{p}_{n}}))}{{\rm d}\delta^{(1)}({\bm{p}_{1}})\cdots{\rm d}\delta^{(1)}({\bm{p}_{n}})}\Big{|}_{\delta^{(1)}(\bm{p}_{a})=0}\,.$$

$$\frac{P_m(\mathbf{k}|\mathbf{x})}{P_m(k)}-1=\sum _O{R}_O(k)\left({\hat{k}}^i{\hat{k}}^j\cdots \right)O_{ij\cdots }(\mathbf{x}),$$

$$\max \left(\frac{p^2}{k^2},\frac{p^2}{k_{\text{NL}}^2}\right),$$

$$\left({\hat{k}}^i{\hat{k}}^j\cdots \right)O_{ij\cdots }^{(m)}(\mathbf{x})=\left[\prod _{a=1}^m{\int }_{{\mathbf{p}}_a}\right]e^{i\mathbf{x}\cdot {\mathbf{p}}_{12\cdots m}}{\mathcal{K}}_O^{(m)}(\{{\hat{\mathbf{p}}}_b\},\{p_b/p_c\}){\delta }^{(1)}({\mathbf{p}}_1)\cdots {\delta }^{(1)}({\mathbf{p}}_m),$$

$$\frac{P_m(\mathbf{k}|\mathbf{x})}{P_m(k)}-1\stackrel{\text{n-th order}}{=}\sum _O{R}_O(k)\left({\hat{k}}^i{\hat{k}}^j\cdots \right)O_{ij\cdots }^{(n)}(\mathbf{x}),$$

$$\displaystyle\mathcal{R}_{n}\big{(}k,t;\{\mu_{\bm{k},\bm{p}_{a}}\},\{\mu_{\bm{p}_{a},\bm{p}_{b}}\},\{p_{a}/p_{b}\}\big{)}=\>$$

$${\Pi }_{ij}^{[1]}(\mathbf{x},\tau )=\frac{2}{3{\Omega }_m{\mathcal{H}}^2}{\partial }_i{\partial }_j\Phi (\mathbf{x},\tau )=K_{ij}(\mathbf{x},\tau )+\frac{1}{3}{\delta }_{ij}\delta (\mathbf{x},\tau ),$$

$$K_{ij}\equiv {\Pi }_{ij}^{[1]}-\frac{1}{3}{\delta }_{ij}\mathrm{tr}[{\Pi }^{[1]}]=\left(\frac{{\partial }_i{\partial }_j}{{\nabla }^2}-\frac{1}{3}{\delta }_{ij}\right)\delta .$$

$${\Pi }_{ij}^{[n]}=\frac{1}{(n-1)!}\left[(\mathcal{H}f{)}^{-1}\frac{\text{D}}{\text{D}\tau }{\Pi }_{ij}^{[n-1]}-(n-1){\Pi }_{ij}^{[n-1]}\right],$$

$$\displaystyle{\rm 1^{st}\ order}\$$

Full text

1 Introduction

The large-scale distribution of matter in the Universe encodes a very rich set of observational imprints that can and have been used to test cosmological models. The standard way to describe its statistical properties is via $n$-point correlation functions of the matter density fluctuations field $\delta(\bm{x})$ [1]. The simplest such object is the $2$-point correlation function, $\xi(r)$, which measures the correlations of the density contrast in regions of the Universe separated by a distance $r$. Its Fourier counterpart is the matter power spectrum,

[TABLE]

Here and throughout, angle brackets denote an ensemble average and a prime on a correlator indicates that the overall momentum conserving $(2\pi)^{3}\delta_{D}(\bm{k}_{\rm tot})$ factor is dropped. For an initially Gaussian distributed density field, and during the linear regime of structure formation, the power spectrum encodes all the statistical information of the matter field. Higher order $n$-point functions become important during the later stages of structure formation, when mode-coupling between different Fourier modes develops specific non-Gaussian signatures; or when the primordial density fluctuations are themselves non-Gaussian, as predicted by a range of inflation models [2]. Accurate theoretical predictions of these higher-order correlations are thus necessary to properly exploit the statistical information of the observational data. Even studies that rely solely on comparing the predicted matter power spectrum with observations, such as gravitational lensing, require a good understanding of higher-order correlations because the mode coupling adds an important contribution to the matter power spectrum covariance. This is described by a specific configuration of the matter trispectrum, the Fourier transform of the connected $4$-point correlation function [3, 4].

The complexity of modeling $n$-point functions increases rapidly with $n$. This holds for perturbation theory methods (see Ref. [1] for a review) as well as measurements in N-body simulations. Perturbative approaches, including the effective field theory (EFT) approach [5, 6] (see Ref. [7] for a review), are limited to the quasi-linear regime, i.e., wavenumbers $k\lesssim k_{\text{NL}}$, where $k_{\text{NL}}=0.3\>h\,{\rm Mpc}^{-1}$ is the nonlinear scale at redshift $z=0$. Numerical simulations of structure formation are currently the only available tool to accurately predict matter fluctuations in the nonlinear regime, but these predictions do not come without downsides. Simulations require significant amounts of computational resources, which makes it harder to obtain quick predictions, needed for instance to cover multidimensional spaces of cosmological models. Higher-order $(n>2)$-point correlations require significantly larger volumes to obtain converged results, as they live in higher-dimensional parameter spaces. Furthermore, estimators of higher-order correlations become themselves more computationally demanding. Attempts at precision measurements of the bispectrum $B(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})$ (3-point correlation function) illustrate these challenges [8, 9, 10].

A simplification of the study of higher order correlation functions can be achieved by focusing on so-called squeezed limit configurations, i.e.,

[TABLE]

Here, a subscript c indicates that we are considering only the connected part of the correlator, and we adopt a shorthand notation for the sum of vectors: $\bm{k}_{12\cdots n}=\bm{k}_{1}+\bm{k}_{2}+\cdots\bm{k}_{n}$. We shall also denote magnitudes of vectors as $k=|\bm{k}|$. This squeezed-limit $(n+2)$-point function represents the coupling of $n$ long-wavelength (or soft) modes with two short-wavelength (or hard) modes. Symmetries of the large-scale structure provide strong constraints on this squeezed limit, a result that is known as “consistency relations” [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. Further, Ref. [22] linked a particular angle-average of these squeezed configurations of the matter $(n+2)$-point function to the $n$-th order response $R_{n}(k)$ of the local matter power spectrum to an initial density perturbation (more precisely, for the case of $n$ superimposed spherically symmetric soft modes). These responses, which are a subset of more general response coefficients defined below, can be measured accurately with separate universe simulations [23, 24, 25, 26], which incorporate spherically symmetric long-wavelength perturbations by simulating curved Friedmann-Robertson-Walker cosmologies (see also Refs. [27, 28, 29, 30, 31] for further applications of the separate universe approach).

One of the main goals of this paper is to show how the relation between squeezed-limit $(n+2)$-point correlation functions and responses can be generalized beyond the special case of $n$ spherically symmetric perturbations, to cover the full shape of these correlation functions at leading order in $p_{i}/k$. This is related to the multipoint-propagator formalism [32], although the latter works at the level of the density field, not the power spectrum considered here. For any given value of $n$, we will see that there is a well-defined, finite set of response coefficients $R_{O}(k)$, which includes the above-mentioned isotropic response coefficient $R_{n}(k)$. One can think of these coefficients as describing the response of the local nonlinear matter power spectrum $P_{m}(\bm{k};\bm{x},\tau)$ measured around position $\bm{x}$ at conformal time $\tau$ to the leading local gravitational observables. These observables include the density perturbation, tidal field, as well as convective time derivatives thereof. The $R_{O}(k)$ can therefore be regarded as the coefficients of a “bias expansion” of the local nonlinear matter power spectrum, in analogy with the expansion of galaxy bias (see Ref. [33] for a review).

The response coefficients $R_{O}(k)$ are physical observables that can be measured in simulations, and one can use physical considerations to evaluate their magnitude and scale dependence (see for example the detailed discussions in Refs. [28, 22]). We describe a procedure which uses perturbation theory results on large scales, together with the simulation measurements of the isotropic response coefficients presented in Ref. [22], to make physically well-motivated estimates for the $R_{O}(k)$. Note that these responses, which quantify the effect of long-wavelength perturbations on the nonlinear gravitational evolution of the small-scale power spectrum, are to be distinguished from the response of the power spectrum to changes in the initial power spectrum, which were measured in Refs. [34, 35].

The sequence of steps followed in this paper can be outlined as follows:

1. i.

First, we rigorously define general power spectrum responses, which we denote as $\mathcal{R}_{n}$ ($n=1,2,\cdots$) (Sec. 2). These are functions of several angles and momentum ratios and describe the general (i.e., not angle-averaged) squeezed limit of ($n+2$)-point functions involving two hard and $n$ soft modes. The small-scale modes are allowed to be fully nonlinear. A diagrammatic representation of these response-type interactions shows how these responses can be used in more general perturbative calculations. 2. ii.

We write the local power spectrum as a bias-like expansion of response coefficients $R_{O}(k)$ that multiply all leading gravitational observables on which the local power spectrum can depend (namely, the density and tidal fields and their time derivatives) at any given order in perturbation theory (Sec. 3). In this way, the $\mathcal{R}_{n}$ are decomposed into a finite set of functions of $k$ only, multiplied by kernels that are uniquely determined by perturbation theory. This decomposition enormously simplifies the description of the responses $\mathcal{R}_{n}$. 3. iii.

We then use these results to provide fully nonlinear predictions for $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$, which are exact on large scales, but use a physically motivated extrapolation of simulation results on small scales. First, we derive the large-scale predictions for the $R_{O}(k)$ relevant for $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ by matching to the tree-level matter bispectrum and trispectrum, respectively (Sec. 4). We then use the measured simulation results for the isotropic coefficients $R_{1}(k),\,R_{2}(k)$ presented in Ref. [22] to extrapolate all $R_{O}(k)$ to nonlinear scales, by employing a separation into “growth” and “dilation” effects (Sec. 5).

As an interesting first application of our framework, we use our description of $\mathcal{R}_{2}$ to predict the squeezed limit of the matter power spectrum covariance, as was already suggested by Ref. [36]. This corresponds to a special case of Eq. (1.2) with $n=2$, $\bm{k}^{\prime}=-\bm{k}$, $\bm{p}_{2}=-\bm{p}_{1}$, and $p_{1}=p_{2}\ll k_{\text{NL}}$ (Sec. 6). We shall see that this formalism, which combines analytical results with small-volume simulation measurements, allows one to match covariance matrix estimates based fully on numerical simulations to very good degree all the way up to $k=1\>h\,{\rm Mpc}^{-1}$ (the interested reader might have a quick look at Fig. 2 on p. 2). We also expand the squeezed covariance in a Legendre multipole expansion and briefly analyze its three non-vanishing moments (monopole, quadrupole and hexadecupole). In Sec. 7, we summarize our main conclusions and outline possible future applications of the framework. Finally, Appendix A lists the Feynman diagram conventions that we adopt in the main body of the paper; in Appendix B, we illustrate explicitly the equivalence between $\mathcal{R}_{2}$ and the $4$-point connected correlation function in the squeezed limit; and in Appendix C, we specify the distinction between Eulerian and Lagrangian definitions of power spectrum responses.

Throughout this paper, we always assume a flat $\Lambda{\rm CDM}$ cosmology with the following parameters (the same as those of the covariance matrix estimates from simulations in Ref. [37]): $h=0.72$, $\Omega_{m}h^{2}=0.1334$, $\Omega_{b}h^{2}=0.02258$, $n_{s}=0.963$, $\sigma_{8}(z=0)=0.801$.

2 Power spectrum response: definition and connection to squeezed $(n+2)$-point functions

In this section, we use diagrammatic representations of interactions (or mode-coupling) in cosmological perturbation theory to define the power spectrum responses, and show how they are directly related to certain squeezed limits of $(n+2)$-point correlation functions (see Appendix A for a summary of the Feynman rules). We begin with the simplest case of the first-order response, and then generalize to $n$-th order.

2.1 First-order response

Consider the nonlinear matter power spectrum, which is denoted as a “blob” propagator with two outgoing modes,

[TABLE]

Throughout, arrows denote the direction of momentum as well as time. Further, whenever it leads to no confusion, we suppress the time argument $t$ to shorten the notation. We will discuss aspects of the time dependence at the end of this section. The nonlinear power spectrum is a non-perturbative quantity, and the blob can be understood as resumming infinitely many perturbative contributions including counterterms. The linear power spectrum on the other hand will be denoted with a dot in our notation.

Now, consider the following, also non-perturbative, $3$-point interaction vertex, with two outgoing hard (high-momentum, or wavenumber) Fourier modes, and one ingoing soft (low-momentum) mode, all defined at a fixed time $t$:

[TABLE]

where $\mu_{\bm{k},\bm{p}}=\bm{k}\cdot\bm{p}/(kp)$. Here and throughout, thick and thin lines denote hard and soft lines, respectively. Further, the notation $\lim_{p\to 0}$ is not to be understood as mathematical limit, but signifies that only the leading contribution in this limit, i.e. the lowest power of $p$, is kept. The resummed vertex in Eq. (2.2) defines the first-order power spectrum response $\mathcal{R}_{1}$. The meaning of $\mathcal{R}_{1}$ can be elucidated as follows. In the limit $p/k\to 0$, we have $\bm{k}^{\prime}\approx-\bm{k}$. That is, up to corrections suppressed by $p/k$, the hard modes are in the same configuration as in the case of the nonlinear matter power spectrum Eq. (2.1). The comparison between Eqs. (2.2) and (2.1) then provides justification to call $\mathcal{R}_{1}(k;\mu_{\bm{k},\bm{p}})$ a power spectrum response, in that it describes the impact of the linearly evolved soft mode $\delta^{(1)}(\bm{p})$ on the nonlinear power spectrum. Note that we do not require that $\bm{k},\bm{k}^{\prime}$ be in the perturbative regime, i.e., they are allowed to be fully nonlinear. On the other hand, in addition to having to be much smaller than $k$, $p$ is also assumed to be in the perturbative regime, i.e. $p\ll k_{\text{NL}}$. One might wonder why we have not allowed for $\mathcal{R}_{1}$ to depend on $p$. The reason is that in the low-$p$ limit, a long-wavelength perturbation modulating the power spectrum appears as a uniform (spatially constant) contribution to the density and tidal field. Any dependence on the wavelength of the soft mode enters only at order $(p/k)^{2}$ or $(p/k_{\text{NL}})^{2}$ (whichever is larger). The same will correspondingly hold for higher-order responses. This will be explicitly justified in the next section.

The vertex corresponding to $\mathcal{R}_{1}$ has three lines, suggesting that its leading contribution will be at the three-point function (bispectrum) level. Indeed, we can obtain its contribution to the (equal-time) bispectrum by attaching a soft power spectrum as in

[TABLE]

where the black dot represents the linear power spectrum $P_{\rm L}(p)$ (see Appendix A). This constitutes the dominant contribution to the matter bispectrum, $B(\bm{k},\bm{k}^{\prime},\bm{p})$ in the squeezed limit $p/k\ll 1$. To demonstrate this, we can work at tree level in perturbation theory, in which case we have that

[TABLE]

where $F_{2}$ is the symmetrized second order perturbation theory kernel (see Eq. (3.18) below). The two terms in brackets in Eq. (2.4) correspond to the diagram of Eq. (2.3) at tree level.111At tree level in perturbation theory, the blob in Eq. (2.3) turns into an $F_{2}$ vertex with a propagator for momentum $\bm{k}$ (and similarly for $\bm{k}^{\prime}$ to account for the $\bm{k}\leftrightarrow\bm{k}^{\prime}$ permutation). The second term in Eq. (2.4), $\propto P_{\rm L}(k)P_{\rm L}(k^{\prime})$, corresponds to a diagram where two hard ingoing modes form a soft outgoing mode, and it is not classifiable as a power spectrum response. However, in the squeezed limit that we are considering, this interaction is suppressed by $F_{2}(\bm{k},\bm{k}^{\prime})\propto\left(p/k\right)^{2}$. This is required by mass and momentum conservation: by definition, small-scale perturbations can only redistribute the mass within a region that is much smaller than the scale $1/p$ of the long mode. This means that fluid momentum generated by the coupling of small-scale modes is of order $i\bm{p}/k^{2}$, which leads to a contribution to the large-scale density suppressed by $(p/k)^{2}$ (see Appendix B of Ref. [38] for a more detailed discussion). Beyond this, there is a further suppression by $P_{\rm L}(k)/P_{\rm L}(p)$ which, given $p\ll k$, is much less than 1 for $p\gtrsim 10^{-2}\,h\,{\rm Mpc}^{-1}$ due to the shape of the power spectrum.

As shown in more detail in Sec. 4, the tree-level prediction for $\mathcal{R}_{1}$ can be read off from equating the first two contributions in Eq. (2.4) to Eq. (2.3):

[TABLE]

Note that indeed, $\mathcal{R}_{1}^{\text{tree}}$ is independent of $p$, as it should be. Moreover, the dependence on $\mu$ simply consists of a monopole and a quadrupole. We will see below that the fully nonlinear response $\mathcal{R}_{1}$ retains a very similar structure, and is completely described by two functions of $k$, which reduce to Eq. (2.5) on large scales. The physical interpretation of the angular structure (i.e., the $\mu$ dependence) is addressed in the next section. Beyond tree level, we can write for the squeezed-limit bispectrum,

[TABLE]

where corrections away from the limit $p\to 0$ are suppressed by $(p/k)^{2}$ and $(p/k_{\text{NL}})^{2}$.

2.2 Generalization to $n$-th order responses

Having gained intuition with the simpler first-order case, we now provide the definition of the general $n$-th order response $\mathcal{R}_{n}$. It is defined analogously to $\mathcal{R}_{1}$ in Eq. (2.2) as

[TABLE]

which corresponds to the modulation of the nonlinear matter power spectrum $P_{m}(k)$ by $n$ linearly evolved long-wavelength modes $\bm{p}_{1},\,\cdots,\,\bm{p}_{n}$. Here and throughout, the notation $\lim_{\{p_{a}\}\to 0}$ implies that only the lowest powers in all of the $p_{a}$ are kept. Specifically, the leading correction to Eq. (2.7) for finite soft momenta is suppressed by

[TABLE]

As we will see below in concrete examples, all singular terms in the $p_{a}\to 0$ limit cancel in Eq. (2.7), so that the limit can be taken in any order. Note that, in addition to $k$, $\mathcal{R}_{n}$ depends on $n(n+1)/2$ cosines of wavenumbers, as well as on $n(n-1)/2$ relative magnitudes of soft momenta, adding up to a total of $n^{2}+1$ arguments including $k$. Recall that $\mathcal{R}_{n}$ does not depend on the overall scale of the $p_{a}$ in the limit $p_{a}/k\to 0$. Analogously to the relations in Eqs. (2.3) and (2.6) between the first-order response and the $3$-point correlation function, the $n$-th order response contributes to the connected $(n+2)$-point function $\langle\delta(\bm{k})\delta(\bm{k}^{\prime})\delta(\bm{p}_{1})\cdots\delta(\bm{p}_{n})\rangle_{c}$ in the kinematic regime where the modes $\bm{p}_{1},\bm{p}_{2},\cdots\bm{p}_{n}$ are soft (the condition that $|\bm{p}_{1\cdots n}|\ll k$ must also hold), through the following diagram:

[TABLE]

The factor $n!$ arises from the permutations among the $\bm{p}_{a}$. The subscript $\mathcal{R}_{n}$ on the correlator indicates that we are considering only the contributions that involve $\mathcal{R}_{n}$, and the squeezed limit is implicitly assumed. In addition to Eq. (2.9), there are two other types of terms that contribute to the $(n+2)$-point function in this limit:

Contributions of the form

[TABLE]

These diagrams can be broken down into vertices involving only soft lines (which can be treated in standard perturbation theory) that interact with the nonlinear power spectrum in a response vertex of lower-order $\mathcal{R}_{m}$ ($m<n$ is the number of soft lines that attach to the $\mathcal{R}_{m}$ response vertex, and these must include at least one whose momentum is the sum of two or more soft momenta). A schematic example of such a term is

[TABLE]

Explicitly, at $n=2$, there is a single such contribution given by (see Appendix B)

[TABLE] 2. 2.

Contributions of the form

[TABLE]

where we have only written tree-level contributions. These terms are characterized by two or more hard power spectra (and correspondingly $n-1$ or less soft power spectra). They necessarily result from two or more ingoing hard modes combining to soft modes, and are the generalization of the squeezed-bispectrum contribution $\propto P_{\rm L}(k)P_{\rm L}(k^{\prime})$ in Eq. (2.4). Therefore, they do not correspond to response vertices. However, as discussed above, in the squeezed regime, they are highly suppressed both by mass-momentum conservation and by the shape of the matter power spectrum. If $n_{h}$ is the number of hard power spectra in a given contribution, they are suppressed by $(p/k)^{2(n_{h}-1)}[P_{\rm L}(k)/P_{\rm L}(p)]^{n_{h}-1}$, where again $p\equiv\max\{p_{1},p_{2},\cdots,p_{n}\}$. On the other hand, loop corrections are instead suppressed by factors of $(p/k_{\text{NL}})^{3+n_{P}}$, where $n_{P}=d\ln P_{\rm L}/d\ln k$ is the linear power spectrum index, evaluated at a scale of order $k_{\text{NL}}$ (to be distinguished from the primordial spectral index $n_{s}$).

To summarize, the leading contributions to the squeezed-limit $(n+2)$-point function with two hard and $n$ soft modes are completely described by the responses $\mathcal{R}_{m}$ with $m\leq n$. This contains the contribution of Eq. (2.9), proportional to $\mathcal{R}_{n}$, as well as those of point 1 above, which are given in terms of $\mathcal{R}_{m}$, $1\leq m<n$. We stress that this holds for fully nonlinear hard modes. The only restriction is that $p\ll{\rm min}\{k,k^{\prime},k_{\text{NL}}\}$. In Appendix B, we illustrate this in detail by listing all the terms for $n=2$ at tree level.

In the discussion presented thus far, the responses $\mathcal{R}_{n}$ are vertex interactions that depend, in addition to $k$, also on the orientations and relative magnitudes of the soft modes, for a total of $n^{2}+1$ arguments. One may therefore argue that, at least at first sight, the problem of calculating the squeezed limit of correlation functions has not been significantly reduced. Crucially however, as we will see in the next section, there is a well-defined decomposition of the $\mathcal{R}_{n}$ into a finite number of response coefficients $R_{O}(k)$ which are only a function of $k$, multiplied by specific kinematic shapes, i.e. functions of the angles and relative magnitudes of the soft modes. For $\mathcal{R}_{1}$, there are two such coefficients and kinematic shapes, while for $\mathcal{R}_{2}$, there are six. Thus, instead of being a free function of five arguments, we will see that $\mathcal{R}_{2}$ is completely determined by six functions of $k$ only, which represents a significant reduction in functional freedom. These discussions are the subject of the next section.

Before continuing however, for completeness, we restore the time dependences in the relation between squeezed $(n+2)$-point functions and $\mathcal{R}_{n}$. Specifically, the responses describe the contributions to the unequal time $(n+2)$-point functions in the limit $\{p_{a}\}\to 0$ as

[TABLE]

where, we recall, the subscript $\mathcal{R}_{n}$ indicates that only those contributions that are captured by $\mathcal{R}_{n}$ are considered, and the squeezed limit is implied. At tree level, i.e. to zeroth order in $k/k_{\text{NL}}$, the $\mathcal{R}_{n}$ are directly related to the perturbation theory kernels, and hence independent of time in an Einstein-de Sitter (EdS) universe. This also holds in $\Lambda$CDM to better than percent-level accuracy. Beyond tree level, the time independence no longer holds. In the EdS limit however, this time dependence only appears in the response coefficients $R_{O}(k)$, which are multiplied by time-independent kernels. Thus, to percent-level accuracy, we can write

[TABLE]

We stress that the response approach does not rely on the EdS approximation, and can be analogously performed using the exact $\Lambda$CDM expansion history as well.

It is also worth noting that additional terms enhanced by $k/p_{a}$ appear in the squeezed limit of unequal-time $(n+2)$-point functions, which are induced by the displacement of the small-scale modes by the large-scale modes [39] (note that the terms derived in Ref. [39] for the bispectrum appear only if the hard modes are evaluated at different times, whereas the response vertices always describe two hard modes at equal time). It can be shown that these contributions are also described by responses, and are of the type of Eq. (2.11), i.e., they involve $\mathcal{R}_{m}$ with $m<n$.

3 Decomposition into response coefficients

We now show how the responses $\mathcal{R}_{n}$ can be decomposed into a set of $k$-dependent coefficients multiplying distinct kinematic shapes. For this, we first derive how the $\mathcal{R}_{n}$ measure the response, or modulation, of the matter power spectrum in the presence of $n$ long-wavelength, linearly evolved density fluctuations. This can be formulated precisely by defining the response as the functional derivative of the local nonlinear power spectrum with respect to the amplitude of the long-wavelength density perturbations, $\delta^{(1)}({\bm{p}_{a}})$ [22, 36]

[TABLE]

Here, we implicitly require that $p_{a}\ll k$, and $p_{a}\ll k_{\text{NL}}$, since the notion of a local power spectrum is only meaningful in this limit. The meaning of $P_{m}$ in Eq. (3) can be rigorously understood as follows. Let us consider a region of size $L$, with $1/k\ll L\ll 1/p$ and $L\gg 1/k_{\text{NL}}$, and measure the Fourier modes $\delta(\bm{k})$ within this volume. Since the large-scale modes can be approximated as constant over this region, it corresponds to a homogeneous (though not necessarily isotropic) patch of space, so that the local, anisotropic power spectrum of small-scale modes in the region of size $L$, $\langle\delta(\bm{k})\delta(\bm{k}^{\prime})\rangle_{L}$, is meaningful.

In order to derive the functional derivative in Eq. (3), we treat the local, small-scale power spectrum as a particular case of a biased tracer (see Sec. 2 of Ref. [33] for a detailed introduction). For this, it is convenient to work in a mixed real/Fourier-space representation. While the small-scale modes are described in Fourier space (within the region of size $L$ as described above), we treat the large-scale modes in real space. In the large-scale bias expansion, we have to include all local gravitational observables that can be constructed out of the long-wavelength modes. As shown in Refs. [40, 41], these consist of all nonlinear local combinations of the tidal tensor $\partial_{i}\partial_{j}\Phi(\bm{x}_{\rm fl}[\tau],\tau)$, as well as its convective time derivatives. Here, $\bm{x}_{\rm fl}(\tau)$ denotes the fluid trajectory leading to a given Eulerian position $\bm{x}$ at time $\tau$, while convective time derivatives are defined as $D/D\tau=\partial/\partial\tau+v^{i}\partial/\partial x^{i}$, where $v^{i}=dx_{\rm fl}^{i}/d\tau$ is the velocity of the matter fluid. Reference [41] derived a compact set of complete, linearly independent terms in the bias expansion at any order, on which we will rely in our considerations below.

To start with, we write the local, anisotropic power spectrum within a region of size $L$ around position $\bm{x}$ as a sum of operators $(\hat{k}^{i}\hat{k}^{j}\cdots)O_{ij\cdots}(\bm{x})$ describing the local gravitational observables corresponding to the long-wavelength modes, multiplied by “bias,” or response coefficients $R_{O}(k)$:

[TABLE]

where we have continued to suppress the time argument (recall that $R_{O}\equiv R_{O}(k,t)$ in general). The notation $\left(\hat{k}^{i}\hat{k}^{j}\cdots\right)$ stands for a product of unit vectors that are contracted with $O_{ij\cdots}$, and which effectively take into account the preferred directions induced in the local patch by the long-wavelength modes. Let us discuss a few noteworthy aspects of the expansion of Eq. (3.2). As before, we define $p\equiv\max\{p_{1},p_{2},\cdots,p_{n}\}$.

The large-scale modes provide preferred directions, so that the growth of small-scale perturbations within the region considered depends on the angle of the small-scale modes with these directions. We have therefore allowed the local power spectrum to be direction-dependent. Note however that $P_{m}(-\bm{k}|\bm{x})=P_{m}(\bm{k}|\bm{x})$ still has to hold, since the density field is real. Contributions of the form $\hat{k}^{i}O_{i}$, or in fact any term involving odd powers of $\hat{\bm{k}}$, are therefore forbidden. Given these constraints and the results of the bias expansion of general tracers, there is a well-defined set of operators appearing in the expansion of Eq. (3.2) at any given order in perturbations. 2. 2.

By definition, the operators $O_{ij\cdots}$ describe the modulation of the local power spectrum within the region of size $L\ll 1/p$. They should therefore be interpreted as coarse-grained on the scale $L$, although we do not emphasize this in the notation. Moreover, in analogy to the bias expansion for galaxy or halo overdensities [42, 43], the operators $O_{ij\cdots}$ are to be considered as renormalized. This implies in particular that $\langle O_{ij\cdots}\rangle=0$. 3. 3.

In real space and evaluated along a fixed fluid trajectory, each of the leading operators $O_{ij\cdots}$ is given by a combination of terms which each involve exactly two spatial derivatives acting on the gravitational potential $\Phi$ [40, 41]. Note that $O_{ij\cdots}$ does in general contain terms that are nonlocal in $\partial_{i}\partial_{j}\Phi$, such as, at second order, $(\partial_{i}\partial_{j}/\nabla^{2})(\nabla^{2}\Phi)^{2}$. When evaluating these operators in perturbation theory at a fixed Eulerian position $(\bm{x},t)$, we also obtain terms that correspond to the displacement of the fluid trajectory from the Lagrangian position. These terms are of the form $s^{k}\partial_{k}O_{ij\cdots}$, where $\bm{s}(\bm{q},t)=\bm{x}_{\rm fl}(\bm{q},t)-\bm{q}$ is the Lagrangian displacement and $\bm{q}$ is the Lagrangian coordinate of the fluid trajectory. 4. 4.

Beyond the leading local gravitational observables mentioned above, there are also contributions to the bias expansion that have more than two spatial derivatives of $\Phi$ in Eq. (3.2). Each additional derivative yields a power of $p_{a}\leq p$. Since odd powers of $k$ are ruled out by the constraint that the density field is real, the leading such higher-derivative term is of order $p^{2}$. The coefficients of these higher-derivative terms contain a spatial scale which describes the size of the region “probed” by the tracer. In our case, the tracer is the matter power spectrum $P_{m}(\bm{k}|\bm{x})$, which means that the spatial scale has to be at least $1/k$. Furthermore, even very small-scale modes $k\gg k_{\text{NL}}$ probe scales that are of order the nonlinear scale $2\pi/k_{\text{NL}}\sim 20\,h^{-1}\,{\rm Mpc}$, which corresponds to the typical distance traveled by dark matter particles over the course of structure formation. In summary, the leading higher-derivative contributions to Eq. (3.2) are expected to scale as

[TABLE]

and are therefore small. In fact, these higher-derivative terms are of the same order as the terms in squeezed-limit $(n+2)$-point functions that are not captured by the response definition of Eq. (2.2) and Eq. (2.7). We therefore do not include them throughout. In principle, it is possible to do so, thereby allowing the response approach to recover as well subleading contributions to squeezed-limit $(n+2)$-point functions, but we do not do this here. 5. 5.

At any given order $m$ in standard perturbation theory, the operators $O_{ij\cdots}$ contracted with $\hat{k}^{i}\hat{k}^{j}\cdots$ can be written as

[TABLE]

where $\mathcal{K}^{(m)}_{O}$ is a kernel that is specific to the $m$-th order contribution in perturbation theory to the operator $O$. Similarly to the perturbation theory kernels $F_{n}$, we assume without loss of generality that $\mathcal{K}^{(m)}_{O}$ is fully symmetrized in its arguments. The Fourier-space kernel only involves the unit vectors $\hat{\bm{p}}_{a}$ as well as relative magnitudes $p_{b}/p_{c}$ of the soft modes. Even though, as mentioned above, nonlocal terms are present in the $O_{ij\cdots}$ in real space, the Fourier representation is an algebraic function of these quantities. We will encounter concrete examples illustrating Eq. (3.4) and the properties of $\mathcal{K}_{O}$ below. In Eq. (3.4), as throughout, we have implicitly assumed that all quantities are evaluated at a fixed time $t$. In general, $\mathcal{K}_{O}$ is a function of $t$ as well, although it is independent of time in the EdS approximation, which we will assume throughout. If better than percent-level precision in $\Lambda$CDM is desired, then this time dependence can be straightforwardly incorporated.

We are now in a position to decompose the responses. At fixed order $n$ in perturbation theory, Eq. (3.2) is given by

[TABLE]

where the sum runs over all operators which start at $n$-th or lower order in perturbation theory (we will see concrete examples below). Crucially, the coefficients $R_{O}(k)$ are independent of the order $n$ in perturbation theory that is employed to describe the soft modes $\bm{p}_{a}$ that make up the operators $O_{ij\cdots}^{(n)}$. After inserting Eq. (3.4) into Eq. (3.5), and the resulting equation into Eq. (3), we can write the response $\mathcal{R}_{n}$ as a sum over kernels multiplied by coefficients that only depend on $k$:

[TABLE]

where the sum runs over the same operators as in Eq. (3.5). This exercise therefore organizes the cumbersome multi-dimensional dependencies of $\mathcal{R}_{n}$ into a simpler linear combination of the response coefficients $R_{O}(k,t)$, which depend only on the scale $k$ and time $t$, and that multiply the kernels $\mathcal{K}^{(n)}_{O}$, which unambiguously fix the possible angular and configuration dependencies. We stress that while the kernels $\mathcal{K}^{(n)}_{O}$ are independent of time in the commonly used EdS-approximated standard perturbation theory approach, as mentioned above, the $R_{O}(k,t)$ always depend on time, as they capture fully nonlinear structure formation.

Before continuing, we note that the $R_{O}$ can be seen as “Eulerian” response coefficients, as we have expanded the matter power spectrum in terms of operators at the final time. In Appendix C, we introduce a slightly different definition of “Lagrangian” response coefficients $R_{O}^{L}$, which are a generalization of those introduced in Ref. [22], and analogous to the Lagrangian bias coefficients of dark matter halos. These correspond to a different basis in the same vector space, i.e., at any given order $n$, the Lagrangian $R_{O}^{L}$ can be expressed in terms of the Eulerian $R_{O}$, and vice versa. Crucially however, the expression for the full response $\mathcal{R}_{n}$ is independent of the choice of Eulerian or Lagrangian, or any other equivalent basis.

3.1 General power spectrum expansion

As mentioned above, the expansion of Eq. (3.2) can be constructed out of $\partial_{i}\partial_{j}\Phi(\bm{x},\tau)$ and its convective time derivatives. We now provide a recipe to write down this expansion by following Ref. [41], who defined222Note that the prefactor $2/(3\Omega_{m}\mathcal{H}^{2})$ is absorbed into the definition of $\Phi$ there.

[TABLE]

where $\Omega_{m}(\tau)$ is the ratio of matter density to critical density and $\mathcal{H}=a^{-1}{\rm d}a/{\rm d}\tau$ is the conformal Hubble rate. $\Pi^{[1]}$ contains the density perturbation and the tidal field via $\delta=\operatorname{tr}[\Pi^{[1]}]$ and

[TABLE]

The superscript $[1]$, which is to be distinguished from $(1)$, refers to the fact that $\Pi^{[1]}$ starts at first order in perturbation theory, but contains higher order terms as well. We then recursively define higher-order tensors $\Pi^{[n]}$ by taking convective time derivatives as (see Ref. [41] for more details)

[TABLE]

where $f=d\ln D/d\ln a$ is the linear growth rate and $D(\tau)$ is the linear growth factor. As shown in Ref. [41] for the case of the galaxy number density, the complete bias expansion consists of all scalar combinations of the $\Pi^{[n]}_{ij}$ that are relevant at any given order. However, the operators $\operatorname{tr}[\Pi^{[n]}]$ with $n>1$ are degenerate with lower-order operators, since they are completely determined by the equations of motion for matter. For this reason, they do not need to be included in the bias expansion.

We now wish to generalize this expansion to the local matter power spectrum. Thus, we include all terms that have an even number of indices to contract with $\hat{k}^{i}\hat{k}^{j}\cdots$, i.e. 2-tensors, 4-tensors, and so on, again excluding $\operatorname{tr}[\Pi^{[n]}]$ with $n>1$. Up to second order, we have

[TABLE]

In the expansion of Eq. (3.2), the ${\rm 1^{st}\ order}$ terms above are to be contracted with $\hat{k}^{i}\hat{k}^{j}$, while the ${\rm 2^{nd}\ order}$ terms are contracted with $\hat{k}^{i}\hat{k}^{j}\hat{k}^{k}\hat{k}^{l}$. For clarity of presentation, we have left the ${\rm 2^{nd}\ order}$ terms written above unsymmetrized. Further, note that the term $\delta_{ij}\delta_{kl}\operatorname{tr}[\Pi^{[2]}]$ is absent following the discussion above. Although we only list the corresponding terms up to second order, the expansion in Eqs. (3.10)–(3.11) can be continued to any desired order (accompanied by a proliferation of the number of terms: at third order for example, there are 14 terms).

A very similar expansion can be performed in Lagrangian space, using the Lagrangian deformation tensor

[TABLE]

In this case, since convective time derivatives reduce to ordinary time derivatives in Lagrangian space, the basis can be simply constructed out of the $m$-th order contributions $M_{ij}^{(m)}$ to $M_{ij}$. This is in contrast to the $\Pi^{[m]}_{ij}$ introduced in the Eulerian expansion above: due to the fact that convective time derivatives are nontrivial in Eulerian coordinates, the tensor $\Pi^{[m]}_{ij}$ is not simply the $m$-th order contribution to $\Pi^{[1]}_{ij}$. As we have noted already above, the choice of Lagrangian vs. Eulerian expansions amounts to a basis choice in the same vector space, so that the final result for $\mathcal{R}_{n}$ is independent of the choice of expansion. We will used the Eulerian expansion Eqs. (3.10)–(3.11) in the following, as they are more directly related to standard perturbation theory results.

3.2 First-order expansion

Let us begin by deriving the decomposition of the linear response $\mathcal{R}_{1}$. From Eq. (3.10), we see that there are only two relevant operators, which we parametrize as $O_{ij}=\{\delta\,\delta_{ij},\,K_{ij}\}$. Equation (3.5) then becomes

[TABLE]

where $R_{1}\equiv R_{\delta}$. Here and in the following, any long-wavelength operator is evaluated at linear order in perturbation theory, unless explicitly stated otherwise, i.e., we drop the superscript (1) in the following. Performing a Fourier transform in the long-wavelength mode, and using that in Fourier space,

[TABLE]

we obtain the following decomposition of $\mathcal{R}_{1}$:

[TABLE]

where $\mu=\bm{p}\cdot\bm{k}/(pk)$. The coefficients $R_{1}(k)$ and $R_{K}(k)$ can be derived at tree level by plugging the above definition of $\mathcal{R}_{1}$ into Eq. (2.6). We shall do this explicitly in the next section.

3.3 Second-order expansion

Including all second-order operators from Eq. (3.11) in Eq. (3.5) (slightly rewritten to single out the density and tidal fields) and performing the Fourier transform on the two long-wavelength modes, we obtain

[TABLE]

where $R_{2}\equiv 2R_{\delta^{2}}$, and

[TABLE]

is directly related to the trace-free part of $\Pi^{[2]}_{ij}$ (see Appendix C of Ref. [33]). Note that we only need to consider the trace-free part of $\Pi^{[2]}$, as the trace is equivalent to other terms already taken into account.

As noted in the discussion above, the first two terms on the right-hand side of Eq. (3.3) are required by the physical definition of the first-order response, since two long-wavelength modes which combine to a second-order density field and tidal field are to be multiplied by the corresponding first-order response coefficients. This is why the above second-order expansion contains eight terms, as opposed to only six as one could naively expect from Eq. (3.11). This is analogous to the term $b_{1}\delta^{(2)}$ appearing in the expansion of biased tracers at second order (see Sec. 2.2 in Ref. [33]). The explicit kinematic dependence encoded in Eq. (3.3) follows from standard perturbation theory, using

[TABLE]

where $\mu_{1}=\hat{\bm{k}}\cdot\hat{\bm{p}}_{1}$, $\mu_{2}=\hat{\bm{k}}\cdot\hat{\bm{p}}_{2}$ and $\mu_{12}=\hat{\bm{p}}_{1}\cdot\hat{\bm{p}}_{2}$. The kinematic shape of $\hat{\varPi}_{ij}$ and $K^{(2)}_{ij}$ contains the quantity $(\hat{\bm{k}}\cdot\bm{p}_{12})^{2}/p_{12}^{2}$, which is a consequence of the fact that these terms are nonlocally related to $\partial_{i}\partial_{j}\Phi$ in real space. Defining $f_{12}=p_{2}/p_{1}$, the said term can be written as

[TABLE]

where in the second equality we have specialized to the case $p_{1}=p_{2}$. Interestingly, the rational function in the angles reduces to a polynomial shape in this particular limit. Combining Eqs. (3.18), (3.19), (3.20), (3.3) and (3), we can then write

[TABLE]

We stress again that this expression for the full second-order response $\mathcal{R}_{2}$ is valid for fully nonlinear $k$.

A similar expression has been derived in Ref. [36] using the squeezed limit of the tree-level four-point function. In particular, when restricting to $f_{12}=1$, Eq. (3.22) matches Eq. (3.16) of Ref. [36], but with a different kinematic basis. That is, their coefficients $A_{n}$ are linear combinations of our $R_{O}$. Note however that the kinematic basis in Eq. (3.22) is derived from operators that correspond to physical, local observables. In Sec. 5, we shall see that this plays an important role in guiding the extrapolation of the response coefficients to the nonlinear regime.

4 Tree-level response coefficients

In this section, we put the results of the last two sections together to explicitly derive the shape of the $R_{O}(k)$ response coefficients in the first-order and second-order response, by matching to the tree-level bispectrum and trispectrum, respectively.

4.1 First order

The first-order response coefficients $R_{1}(k)$ and $R_{K}(k)$ can be read off from the equation that is obtained by plugging Eq. (3.15) into Eq. (2.6) at tree level:

[TABLE]

where we have expanded the power spectrum as P_{\rm L}(|\bm{k}+\bm{p}|)\approx P_{\rm L}(k)\big{(}1+\mu pP_{\rm L}^{\prime}(k)/P_{\rm L}(k)\big{)}, and here and throughout, a prime on a power spectrum denotes a derivative with respect to $k$. Note that the terms $\propto k/p$ in $F_{2}$ that diverge in the limit $p\to 0$ cancel exactly. This is expected, since these terms correspond to an overall displacement due to the long mode, which is not locally observable and hence does not physically modulate the local power spectrum. The functions $R_{1}^{\rm tree}(k)$ and $R_{K}^{\rm tree}(k)$ are obtained straightforwardly by matching to the constant and the $\mu^{2}$ terms in Eq. (4.1):

[TABLE]

In this section, the superscript tree serves to emphasize that the result is valid only at tree level. These expressions match the equivalent ones derived in Refs. [22] (in the case of $R_{1}$) and Refs. [36, 44]. Note also that, as expected, $R_{1}(k)=\int_{-1}^{1}\mathcal{R}_{1}(k,\mu){\rm d}\mu$, and that this tree-level result for $R_{1}$ is recovered by the simulation measurements of Ref. [22] on large scales.

4.2 Second order

The derivation of the second-order response coefficients follows the exact same steps, but applied to the $n=2$ case of the tree-level limit of Eq. (2.7), which describes the 4-point function or trispectrum. Explicitly, we plug the $\mathcal{R}_{2}$ expansion of Eq. (3.22) into the following equation:

[TABLE]

where (see Appendix B, and Appendix A. 2. of Ref. [22] for the angle-averaged case)

[TABLE]

Diagrammatically, the tree-level $\mathcal{R}_{2}^{\rm tree}$ interaction vertex can be represented as the following sum of $F_{3}$ and $(F_{2})^{2}$ kernels (clarifying the meaning of the open circle in Eq. (4.2)):

[TABLE]

where $p_{1}$,$p_{2}$ are to be understood as much softer than $k$,$k^{\prime}$. Equation (4.4) is obtained from Eq. (4.5) by the Feynman rules, after attaching propagators to the external momenta $\bm{p}_{1}$, $\bm{p}_{2}$.

After Taylor-expanding $P_{\rm L}(|\bm{k}+\bm{p}|)$ and $P_{\rm L}(|\bm{k}+\bm{p}_{12}|)$ in Eq. (4.4) to second order in the amplitude of the soft modes, so that the IR divergences cancel out (cf. Eqs. (B.13) in Appendix B), and specializing to the case $f_{12}=1$, we obtain

[TABLE]

The general result for $f_{12}\neq 1$ is much more lengthy and can be found in Eq. (B.9); however, the case $f_{12}=1$ is sufficient to unambigously determine all second-order response coefficients.

The second-order response coefficients can now be derived by equating Eq. (3.22), specializing to $f_{12}=1$, to Eq. (4.2) through Eq. (4.2). After matching the relevant kinematic shapes (analogously to the first-order responses above; performing this matching for a general $f_{12}$ leads to the same result), we find

[TABLE]

As expected, we recover the first-order response coefficients via the distinct shapes of the second-order density and tidal fields, which serves as a cross-check of the calculation. Recall that, as discussed in Sec. 3.3, the presence of these first-order coefficients in the second-order expansion is required by the fact that the response coefficients are independent of the perturbative order used to describe the long-wavelength modes.

5 Response coefficients in the nonlinear regime

The response coefficients $R_{O}(k)$ are physical quantities that can be measured with appropriately setup N-body simulations, thereby offering a means to accurately determine them at any $k$. In the case of the isotropic response coefficients $R_{n}(k)=n!R_{\delta^{n}}(k)$, which correspond to the response with respect to $\delta^{n}$, the task can be carried out effectively using the separate universe simulation formalism [23, 24, 25, 26], in which the effect of the long-wavelength density perturbation is included by following structure formation in a spatially curved Friedmann-Robertson-Walker universe. The latter can in turn be neatly incorporated into N-body simulations by an appropriate modification of the cosmological parameters. The coefficients $R_{n}$ can be seen as the matter power spectrum analogue of the so-called “local-in-matter-density” (LIMD; previously known simply as “local”) bias parameters of large-scale structure tracers, which can similarly be measured through separate universe simulations [29, 45, 46]. In Ref. [22], the authors used separate universe simulations to measure the first three isotropic response coefficients $R_{1}(k)$, $R_{2}(k)$ and $R_{3}(k)$ in the nonlinear regime.

The remaining response coefficients, one at first order and five at second order, which disappear under a full angle average but are relevant in general, cannot be measured using the standard separate-universe setup. Instead, simulations including large-scale tidal fields, which break isotropy, are required. Setting up such simulations, while clearly well motivated, is beyond the scope of this paper. Instead, we shall use the tree-level expressions of the previous section together with the measured isotropic response coefficients from Ref. [22] to obtain physically motivated extrapolations of all $R_{O}(k)$, as described in this section.

In the separate universe picture, which is exact for $\Lambda$CDM [47], the isotropic responses studied in Ref. [22] can be broken down into three physically distinct contributions [48, 14, 49, 22, 20]

The reference density contribution, which accounts for the fact that the local matter power spectrum is defined with respect to the local mean density, while the global power spectrum is defined with respect to the fiducial background. This in practice corresponds to an overall rescaling of the amplitude of the power spectrum, and is uniquely determined at all orders by the perturbative evolution of the large-scale density field. 2. 2.

The dilation effect contribution, which arises because the long-wavelength modes perturb the local scale factor, which therefore rescales the size of the small-scale modes. This corresponds to a rescaling of the argument of the matter power spectrum. Note that if the long-wavelength perturbations are anisotropic, then the perturbations to the scale factor will also be anisotropic. This contribution is completely determined by the perturbative evolution of the large-scale modes, as well as the shape of the nonlinear matter power spectrum. 3. 3.

The growth effect, which quantifies the changes in structure formation induced by the long-wavelength perturbations, and corresponds to the actual physical coupling between the long- and short-wavelength modes. Reference [22] defined the growth-only response coefficients $G_{n}(k)$, which measure the impact of this effect alone.

Effects 1 and 2 above can be calculated exactly at any order given the nonlinear matter power spectrum in the fiducial cosmology, i.e., no dedicated separate universe simulations are necessary. Separate universe simulations are however required to measure effect 3 on nonlinear scales. Next, we will attempt to single out the contribution from the growth effect in all coefficients $R_{O}(k)$, and extrapolate this contribution to the nonlinear regime of structure formation using the measured growth-only isotropic responses in the separate universe simulations of Ref. [22].

5.1 Extrapolating the first order response coefficients

We start with the first order response coefficients, $R_{1}$ and $R_{K}$. Following Ref. [22], at tree level, $R_{1}$ can be decomposed as

[TABLE]

where the first, second and third terms on the second line are the reference density, dilation and growth effects, respectively. At tree level, $G^{\rm tree}_{1}=26/21$. By replacing $G_{1}$ in this expression with its full nonlinear shape, as well as replacing the linear matter power spectrum with the nonlinear one, we obtain an expression for $R_{1}(k)$ which, because of the exact validity of the separate universe, perfectly matches that from the separate universe simulations [22]. Here and below, we use the superscript tree to indicate expressions valid at tree level only, and implicitly assume the general nonlinear result otherwise. The same steps can be followed to extrapolate $R_{K}$ to nonlinear scales as well. In particular, the tidal response at tree level is

[TABLE]

There is no reference density effect here, and the first term corresponds to the physical growth effect of a long-wavelength tidal field [50], while the second term is the dilation effect. Thus, in order to obtain the nonlinear $R_{K}$, we can adopt the following form:

[TABLE]

The coefficient $\alpha=12/13$ is determined by the requirement that $R_{K}$ recovers the tree-level limit as $k\to 0$, so that $\alpha G_{1}^{\rm tree}=8/7$. Note also that we have replaced $P_{\rm L}$ by $P_{\rm m}$. Hence, our extrapolation for $R_{K}(k)$ becomes

[TABLE]

This assumes that the scale dependence of the physical modulation, i.e. the growth effect, of the local power spectrum by a long-wavelength tidal field is the same as that of a long-wavelength density perturbation. Beyond the large-scale limit, where the tree-level result is exact, this expression is not expected to be perfect. In the halo model description of the nonlinear matter density (see Ref. [51] for a review), the scale dependence of the growth response is determined by the sensitivity of the abundance of massive halos to long-wavelength perturbations, which leads to the increased response on scales of $k\sim 0.7\,h\,{\rm Mpc}^{-1}$ (see e.g. Fig 2 of Ref. [22]). On the other hand, the structure of halos, which determines the power spectrum in the deeply nonlinear “1-halo” regime, is largely insensitive to long-wavelength perturbations, which suppresses the response at $k\gtrsim 1\,h\,{\rm Mpc}^{-1}$. These considerations, which are described in more detail in Ref. [22], are expected to hold for anisotropic large-scale perturbations as well. Thus, our extrapolation should serve as a reasonable first approximation that has no free parameters. The tree-level and nonlinear shapes of $R_{1}(k)$ and $R_{K}(k)$ are shown in Fig. 1.

5.2 Extrapolating the second order response coefficients

The same considerations as above apply, but become more involved in the case of the second-order anisotropic response coefficients, $R_{K\delta}$, $R_{K^{2}}$, $R_{K.K}$, $R_{KK}$ and $R_{\hat{\varPi}}$. We can split these coefficients into those that do not depend on the second derivative of the matter power spectrum and those that do. For those that do not, which are $R_{K^{2}}$, $R_{K.K}$ and $R_{\hat{\varPi}}$, we choose to extrapolate them in the exact same way as $R_{K}$ above, i.e., by using $G_{1}(k)$ to rescale the constant term.333We have checked that using $G_{2}(k)$ instead gives the same result up to $10-20\%$ for $k\lesssim 1\>h\,{\rm Mpc}^{-1}$. We have also checked that assuming a reference density term to extrapolate these coefficients makes practically no difference up to $k\approx 0.3\>h\,{\rm Mpc}^{-1}$, and never induces differences larger than $20\%$ for higher values of $k$. Those that do depend on the second derivative, which are $R_{K\delta}$ and $R_{KK}$, require more deliberation. We can use clues provided by the decomposition of $R_{2}$ into reference density, dilation, and growth contributions derived in Ref. [22] (using their Eulerian definitions, not Lagrangian ones, cf. Appendix C):

[TABLE]

The term $G_{1}^{\prime}(k)$ vanishes at tree level and only accounts for less than $\approx 10\%$ of the total value of $R_{2}(k)$ on small scales. Here, we ignore this term and take the following ansatz for the extrapolation of $R_{K\delta}$ and $R_{KK}$:

[TABLE]

where the coefficients $\beta$ and $\gamma$ are determined by the conditions that $\beta\left[\frac{8}{21}G_{1}^{\rm tree}+G_{2}^{\rm tree}\right]$ and $\gamma\left(-\frac{2}{9}-\frac{2}{3}G_{1}^{\rm tree}\right)$ match the corresponding coefficients $\propto 1$ and $\propto P_{\rm L}^{\prime}(k)/P_{\rm L}(k)$ of the tree-level expressions in Eq. (4.2). Note that $\beta$ and $\gamma$ take on different values in $R_{K\delta}$ and $R_{KK}$, and $G_{2}^{\rm tree}=3002/1323$. We emphasize again that there are no free parameters in the resulting extrapolated $R_{O}(k)$.

Putting this all together, the nonlinear extrapolations of the response coefficients are given by

[TABLE]

where $G_{1}(k)$ and $G_{2}(k)$, as well as $R_{1}(k)$ and $R_{2}(k)$, are taken directly from the measurements of Ref. [22]. We note that these extrapolations assume similar scale dependencies of the response coefficients for different operators, which is an approximation that should be assessed using simulations implementing the different long-wavelength configurations.

The tree-level and nonlinear extrapolation of the response coefficients are shown in Fig. 1 at $z=0$. We have used CAMB [52] to evaluate $P_{\rm L}(k)$, the Coyote emulator [53] to evaluate $P_{m}(k)$, and compute the derivatives via finite-differencing. By construction, the curves on the left and right panels agree on large scales. When compared to the tree-level result, the nonlinear expressions exhibit the expected suppression of the amplitude of the BAO oscillations, and the broad-band suppression on small scales as one begins to probe the interior of virialized structures. Further, one can identify the following hierarchy between the response coefficients:

[TABLE]

This hierarchy is determined by the tree-level result, and it is preserved by our nonlinear extrapolation. This again should also be confirmed by full simulation measurements of the anisotropic response coefficients.

6 Application: squeezed matter power spectrum covariance

In this section, we apply the response formalism described in the previous sections to predict the non-Gaussian part of the matter power spectrum covariance in the squeezed limit. We focus on the frequently considered case of the equal-time power spectrum covariance here, although this can be straightforwardly generalized to the covariance of power spectra at different times (relevant for line-of-sight integrated observables such as lensing). The power spectrum covariance is defined by

[TABLE]

where the hat on top of the power spectra indicates that these are power spectrum estimators. Specifically, the estimators average over a set of Fourier modes within bins centered on $\bm{k}_{1},\,\bm{k}_{2}$ (we use the central wavenumbers $\bm{k}_{1},\bm{k}_{2}$ to label the bins, for clarity). The covariance can be decomposed into three contributions:

[TABLE]

described in the following.

The Gaussian contribution ${\rm Cov}^{\rm G}$. This represents the trivial disconnected diagonal contribution to Eq. (6.1). If we consider the covariance of the angle-averaged power spectrum within a wavenumber bin of width $\Delta k$, the Gaussian contribution is given by

[TABLE]

where $N_{k}=VV_{k}/(2\pi)^{3}$ is the number of modes that are averaged over, with $V_{k}=4\pi k^{2}\Delta k$ and $V$ is the survey volume. 2. 2.

The non-Gaussian (trispectrum) contribution, ${\rm Cov}^{\rm NG}$, which (assuming Gaussian initial conditions) is induced by mode-coupling due to nonlinear gravitational evolution. Specifically, this involves the “parallelogram” configuration of the trispectrum, $T(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{2},-\bm{k}_{2})$, and will be described in more detail below. 3. 3.

The super-sample covariance contribution ${\rm Cov}^{\rm SS}$ [54, 55, 27, 28, 49]. This second non-Gaussian contribution accounts for the interaction of modes observable within the survey with modes whose wavelength is larger than the scale of the survey [27]. This contribution is in fact completely captured by the first-order response $\mathcal{R}_{1}$ [28, 44]. For the case of angle-averaged spectra and isotropic survey window functions one obtains

[TABLE]

where $\tilde{W}(p)$ is the Fourier transform of the survey window function. This contribution formally arises from the convolution of the trispectrum with the window function, and scales differently with survey volume than the other two contributions. This is why it is sensible to treat it separately from the parallelogram trispectrum that contributes to ${\rm Cov}^{\rm NG}$. Note also that the super-sample contribution is not included when estimating the power spectrum covariance from an ensemble of standard N-body simulations, which do not include fluctuations on scales larger than the simulation box.

We note in passing that when estimating the covariance from N-body simulations, there is also a purely numerical contribution from the particle shot noise, which we have not included here.

Both the first and third contributions in Eq. (6.2) are well understood. The second contribution, ${\rm Cov}^{\rm NG}$, is however considerably more challenging to predict. For the case of angle-averaged power spectra, this contribution is given as a bin-average over the trispectrum through

[TABLE]

where $\mu_{12}={\bm{k}}_{1}\cdot{\bm{k}}_{2}/(k_{1}k_{2})$ and the simplifying approximation follows from assuming very narrow bins ($\Delta k/k_{i}\ll 1$) around $\bm{k}_{1},\,\bm{k}_{2}$ such that the trispectrum can be treated as constant over the bins (see e.g. Ref. [56] for an explicit demonstration showing that this is valid). Equation (6) provides the precise relation between the covariance of the angle-averaged power spectrum and the matter trispectrum. Moreover, the trispectrum configuration appearing in Eq. (6) is a function of $k_{1}$, $k_{2}$, and the angle $\mu_{12}$ between them. We can thus generalize Eq. (6) and define the non-Gaussian power spectrum covariance before angle-averaging, $\operatorname{Cov}^{\rm NG}(k_{1},k_{2},\mu_{12})=T(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{2},-\bm{k}_{2})/V$ (distinguished here from ${\rm Cov}^{\rm NG}(k_{1},k_{2})$ by the presence of $\mu_{12}$ in the argument), which we can decompose in multipoles as

[TABLE]

and $\mathcal{P}_{\ell}(x)$ is the Legendre polynomial of order $\ell$. Note that due to the fact that $P_{m}(-\bm{k})=P_{m}(\bm{k})$ has to hold, the odd multipole moments vanish identically. Note also that the $\ell=0$ case (called the monopole) corresponds exactly to ${\rm Cov}^{\rm NG}(k_{1},k_{2})$ in Eq. (6). This is the most commonly considered case of the power spectrum covariance, i.e., that of the angle-averaged power spectrum.

Here, we shall be interested in the squeezed limit of ${\rm Cov}^{\rm NG}(k_{1},k_{2},\mu_{12})$, where $k_{1}\ll\min\{k_{2},\,k_{\text{NL}}\}$, i.e., the non-Gaussian covariance of the small-scale with the large-scale power spectrum. Naturally, by symmetry, the same results trivially hold in the limit $k_{2}\ll\min\{k_{1},\,k_{\text{NL}}\}$. Using Eq. (6) as well as Eq. (2.9) for $n=2$, one can write

[TABLE]

where $k_{\text{hard}}=\max\{k_{1},k_{2}\}$, $k_{\text{soft}}=\min\{k_{1},k_{2}\}$. Note that the response-type contribution $\propto\mathcal{R}_{1}$ of Eq. (2.12) does not contribute in this particular configuration because $p_{12}=0$. The leading corrections come from beyond-squeezed-limit as well as from nonlinear terms, as discussed in Sec. 2 and Sec. 3. Some of the nonlinear terms involve loop interactions that are not captured by the $\mathcal{R}_{2}$ interaction vertex alone, and which become important if $k_{\text{soft}}$ is not sufficiently smaller than $k_{\text{NL}}$. Here, we work only at tree level in $k_{\text{soft}}$, whose leading contribution is the first term in the second line of Eq. (6.7). For completeness, we note also that Eq. (6.7) can naturally be generalized to the unequal-time covariance as

[TABLE]

where $t_{\rm hard}$ and $t_{\rm soft}$ correspond to the times of the hard and soft momenta, respectively. Using Eq. (3.22), we can write the configuration of $\mathcal{R}_{2}$ that enters in Eq. (6.7) as

[TABLE]

where $\mathcal{P}_{2}(\mu_{12})=(3\mu_{12}^{2}-1)/2$ is the second Legendre polynomial. Recall that Eqs. (6.7)–(6) are exact in the limit where the soft mode is sufficiently small, with any restrictions placed on the hard mode coming only from the maximum wavenumbers probed by the simulations which provide the response coefficient measurements.

Naturally, the multipole expansion of Eq. (6.6) remains valid in the squeezed limit, and so we can define

[TABLE]

Given the kinematic dependences appearing in Eq. (6), only $\ell=0,2,4$ are non-vanishing. These multipoles are given as linear combinations of $R_{O}(k)$ as we discuss next.

Taking $\ell=0$ in Eq. (6), we arrive at the following expression for the monopole squeezed-limit covariance,

[TABLE]

When inserting the tree-level expressions for the $R_{O}(k)$ derived in Sec. 4 into Eq. (6.11), we recover exactly the tree-level squeezed covariance derived in Ref. [36].

In Fig. 2, we show Eq. (6.11) (red line), as well as the angle-averaged covariance matrix estimated using N-body simulations in Ref. [37] (black; see e.g. Refs.[57, 58, 59, 28] for other estimations of the covariance from simulations). The result corresponds to $z=0$ and we take $V=656.25\>h^{-3}{\rm Mpc}^{3}$, which is the volume of the simulation boxes. The covariance matrix is shown as a function of $k_{1}$ for the fixed $k_{2}$ values indicated above each panel. Note that in the figure we are not restricting to squeezed configurations, i.e., we also show the result for $k_{1}\approx k_{2}$, which is a regime in which the result is not expected to be accurate. The grey shaded area shows the region where $1/2\leq k_{1}/k_{2}\leq 2$, which serves to mark roughly these non-squeezed configurations. Moreover, note also that in each panel, the meaning of $k_{\text{hard}}$ and $k_{\text{soft}}$ switches on either side of the $k_{1}=k_{2}$ equality. Specifically, when $k_{1}<k_{2}$ (left of the center of the grey area) $k_{\text{hard}}=k_{2}$, whereas $k_{\text{hard}}=k_{1}$ when $k_{1}>k_{2}$ (right of the center of grey area).

Figure 2 shows that Eq. (6.11) provides a good description of both the amplitude and shape of the simulation results for all configurations shown. For $k_{\text{soft}}=k_{2}\lesssim 0.06\,h\,{\rm Mpc}^{-1}$ and $k_{\text{soft}}<k_{\text{hard}}/2$ (i.e., outside of the shaded area), the agreement is remarkable. In this squeezed regime, we expect all other contributions to the total covariance to be negligible, such that comparisons to the result of Eq. (6.11) may be a useful diagnostic of any systematics in covariance estimates based fully on N-body simulations. For higher values of $k_{2}$, shown in the lower panels, one can observe that our result slightly underpredicts the simulation results at high $k_{1}$. This is due to nonlinear corrections which become noticeable as $k_{\text{soft}}/k_{\text{NL}}$ becomes sizeable. By including the 1-loop contribution to the covariance, we expect to obtain a significantly improved match.

Note also that the response coefficients were measured in Ref. [22] for a slightly different cosmology than the one used for the covariance measurements of Ref. [37]. In particular, Ref. [22] used: $h=0.7$, $\Omega_{m}h^{2}=0.1323$, $\Omega_{b}h^{2}=0.023$, $n_{s}=0.95$, $\sigma_{8}(z=0)=0.8$. We have verified that the results obtained using the cosmology of Ref. [22] to compute power spectra are indistinguishable from those using the cosmology of Ref. [37]. Relative to the power spectrum, we expect the responses to be less sensitive to changes in cosmological parameters, hence we expect that this slight inconsistency in the parameters is entirely negligible for our comparisons here.

Figure 3 shows the relative magnitude of the response coefficient terms that contribute to the monopole squeezed-limit covariance in Eq. (6.11). The figure shows that, although the isotropic second-order response coefficient $R_{2}$ provides the dominant contribution (blue line), the other response coefficients still contribute appreciably. This highlights the importance of the anisotropic responses, even when considering the angle-averaged covariance. Recall that we have relied on extrapolations for the anisotropic responses (cf. Sec. 5), as these have so far not been measured directly in N-body simulations. The fact that their contribution is non-negligible, together with the good agreement between model and simulations displayed in Fig. 2, then suggests that our extrapolation steps are, at the very least, not grossly wrong.

One might also wonder how well the prediction performs in Fig. 2 when using the tree-level (cf. Eqs. (4.2)) rather than fully nonlinear response coefficients (cf. Eqs. (5.2)). We find that the lack of both BAO damping and suppression at high $k$ of the tree-level coefficients leads to a significantly worse description of the covariance, as expected.

Finally, we give the quadrupole ($\ell=2$) and hexadecupole ($\ell=4$) of the squeezed-limit covariance, obtained from Eq. (6):

[TABLE]

As expected, the isotropic responses $R_{n}(k)$ do not contribute for $\ell>0$. Figure 4 compares the $l=0,2,4$ cases. We see that, for the squeezed covariance, the quadrupole is predicted to be larger than the monopole. These results reiterate the power and usefulness of the response approach, which immediately allows us to derive all multipoles of the squeezed matter power spectrum covariance. We leave the comparison of the higher multipoles of the power spectrum covariance with simulation results (e.g. Ref. [58]) to future work.

7 Summary and conclusions

We have presented a rigorous definition of power spectrum responses $\mathcal{R}_{n}$, which can be used to directly describe the squeezed limit of $(n+2)$-point functions with two hard and $n$ soft modes. More specifically, the dominant tree-level contributions to the $(n+2)$-point matter correlation functions in this limit are those that are linear in the hard power spectrum (cf. Sec. 2.2 and Appendix B), and these are precisely the terms which can be represented through power spectrum responses. The response approach can be employed whenever long-short mode coupling is relevant in perturbation theory. The response functions can be measured accurately in the nonlinear regime with only a few N-body simulations, and as a result, by replacing the relevant interactions in perturbation theory with the corresponding response vertices (cf. Eq. (2.7)), one effectively resums infinitely many higher loop contributions (of a specific form), dramatically extending the range of scales where perturbative approaches can be employed and are predictive.

The $\mathcal{R}_{n}$ can be decomposed into a finite set of response coefficients $R_{O}$ (which at fixed time only depend on the hard mode $k$) that multiply well-defined kernels that specify the kinematic dependences and are fully determined by perturbation theory. This organization of the structure of the responses simplifies their description significantly (cf. Sec. 3). The determination of the kernels follows from a complete enumeration of local gravitational observables, which was recently developed in the context of general perturbative bias expansions [40, 41]. The shape of the $R_{O}$ can be determined at tree level by matching to the corresponding squeezed matter correlation functions (cf. Sec. 4). On nonlinear scales, their shape can be determined with the aid of dedicated small-volume simulations, as was shown recently for the isotropic response coefficients in Ref. [22].

As an application of this framework, we have considered the non-Gaussian (connected) part of the matter power spectrum covariance. This is a crucial component in the cosmological interpretation of, for example, the two-point function of gravitational lensing shear [60, 55, 61, 62], and several steps in modeling it have been made recently (see e.g. Ref. [63] for a perturbation theory calculation at 1-loop order, and Refs. [64, 56] for a more phenomenological approach). Here, we have considered the squeezed limit of the power spectrum covariance ${\rm Cov}^{\rm NG}(k_{\text{soft}},k_{\text{hard}})$ at tree level in the soft mode, i.e. for $k_{\text{soft}}\ll k_{\text{NL}}$, $k_{\text{soft}}\ll k_{\text{hard}}$, but any $k_{\text{hard}}$. This covariance corresponds to a configuration of the matter trispectrum whose relevant contributions are completely determined by $\mathcal{R}_{2}$ (cf. Sec. 6). The response formalism therefore allows us to readily evaluate the covariance of the angle-averaged power spectrum, as well as its anisotropic part. The interpretation of this squeezed limit of the covariance as a specific response had been already pointed out by Ref. [36]. The specific decomposition of $\mathcal{R}_{2}$ into a set of physical response coefficients that we performed here, however, has the advantage of permitting the use of physical considerations to estimate the nonlinear behavior of the response. This includes guiding the design of the simulation setups that are needed to measure the anisotropic response coefficients, as well as guiding the extrapolation of the anisotropic $R_{O}$ to nonlinear scales using the isotropic responses measured in simulations (cf. Sec. 5).

We have compared the resulting prediction of the angle-averaged ${\rm Cov}^{\rm NG}(k_{\text{soft}},k_{\text{hard}})$ in the squeezed limit to fully numerical estimates of the covariance, and found very good agreement when the soft mode is sufficiently linear, $k_{\text{soft}}\lesssim 0.06\ \>h\,{\rm Mpc}^{-1}$ (cf. Sec. 6). For higher values of $k_{\text{soft}}$, terms beyond tree level become relevant. We stress that our predictions here have relied on an extrapolation, rather than measurement, of the anisotropic response coefficients in the nonlinear regime. Indeed, this is the only approximation made in our prediction for the squeezed-limit covariance that is relevant when comparing with gravity-only simulations. A simulation measurement of these will require nontrivial technical work, since the anisotropy is incompatible with standard cubic simulation boxes. Nevertheless, the actual measurement of the coefficients will involve substantially less computational resources compared to full numerical measurements of the covariance (e.g., compare the total simulation volume used in Ref. [22] to that in Ref. [37]). Given such measurements, the prediction of the squeezed-limit covariance becomes exact. The response formalism applied to the covariance remains true even beyond the gravity-only case we considered so far. A measurement of the responses, and hence squeezed covariance, including baryonic effects is entirely feasible due to the greatly reduced computational demands. For the same reasons, estimating the covariance for a range of cosmological models would also be within reach. The issue of how responses can be useful for the covariance beyond the squeezed limit is subject of ongoing work.

We finish with a brief list on possible further applications and extensions of the response approach:

• •

Covariance of power spectra in redshift space, relevant for the application to spectroscopic tracers, and at different times (cf. Eq. (6)), relevant for applications to weak lensing measurements.

• •

Response of the matter bispectrum: this can be defined in close analogy with the power spectrum response defined in Sec. 2, and corresponds to resummed vertices with three hard outgoing lines $\bm{k}_{1},\bm{k}_{2},\bm{k}_{3}$ and $n$ ingoing soft lines $\bm{p}_{a}$. The main difference to the power spectrum case is that the kinematic structure becomes more complicated.

• •

Response of the galaxy power spectrum: there is in principle no obstacle to applying this formalism to biased tracers instead of matter, for example to predict the galaxy power spectrum covariance [65, 66, 67]. However, because of the absence of mass/momentum conservation of galaxies, the coupling of two hard modes $k$ to a soft galaxy density perturbation $p$ is not necessarily suppressed by $(p/k)^{2}$. This leads to additional stochastic contributions which need to be included to provide an accurate description of squeezed-limit $n$-point functions.

The squeezed-limit covariance results presented here therefore form only a small subset of a wide range of applications of the response approach, which remains to be explored in the future.

Acknowledgments

We thank Linda Blot for providing the numerical measurements of the power spectrum covariance, and Daniele Bertolini, Joachim Harnois-Déraps, Wayne Hu, and Mikhail Solon for useful discussions.

FS acknowledges support from the Marie Curie Career Integration Grant (FP7-PEOPLE-2013-CIG) “FundPhysicsAndLSS,” and Starting Grant (ERC-2015-STG 678652) “GrInflaGal” from the European Research Council.

Appendix A Feynman rules for cosmological perturbation theory

We now spell out the Feynman rules we employ in calculating perturbative predictions for cosmological correlations. Different conventions are used in the literature, and ours is based on the one used by Ref. [38]. The rules are as follows:

An $n$-point correlation function is represented by a collection of diagrams with $n$ outgoing external legs. 2. 2.

Interaction vertices have $m\geq 2$ ingoing lines $\bm{p}_{1},\cdots,\bm{p}_{m}$ coupling to a single outgoing line $\bm{p}$. Each such vertex is assigned a factor

[TABLE]

In this paper, we assign a positive (negative) sign to outgoing (ingoing) momenta. Each ingoing line has to be directly connected to a propagator (linear power spectrum).444This is because diagrams that involve interaction vertices directly connected to each other are absorbed into higher-order interaction vertices. 3. 3.

Propagators are represented in our notation as vertices with 2 outgoing lines of equal momentum $k$ as

, and they are assigned a factor $P_{\rm L}(k)$. To ease the notation, we often skip labeling these two outgoing lines (which line is which can always be inferred from momentum conservation). 4. 4.

All momenta that are not fixed in terms of momentum constraints are integrated over as

[TABLE]

A diagram without any loop integral is said to be a tree-level diagram. 5. 5.

Each diagram is multiplied by the symmetry factor, which accounts for the number of all nonequivalent labelings of external lines and degenerate configurations of the diagram.

In order to include response-type interactions, we augment these rules by one additional rule:

A second class of interaction vertices is allowed, which have 2 (instead of 1) outgoing lines with momenta $\bm{k},\bm{k}^{\prime}$, and $n\geq 1$ incoming lines with momenta $\bm{p}_{a}$. These vertices are only predictable in the limit where $p=\max\{p_{a}\}\ll\min\{k,k_{\text{NL}}\}$, and $\Sigma_{a}p_{a}\ll\min\{k,k_{\text{NL}}\}$, but no restriction is placed on the magnitude of the outgoing momenta, which can be fully nonlinear. In this regime, each such vertex is assigned a factor (cf. Eq. (2.7))

[TABLE]

We have not written the arguments of $\mathcal{R}_{n}$ explicitly here, which are described in detail in Sec. 2. The factor $1/2$ cancels the trivial permutation $\bm{k}\leftrightarrow\bm{k}^{\prime}$, which is always present.

Appendix B Tree-level matter trispectrum

In this appendix, we explicitly list all the terms that contribute to the $4$-point matter correlation function at tree level. This serves as a concrete example of the discussion in Sec. 2.2 around which types of contributions to the $(n+2)$-point connected correlation function can be represented as power spectrum responses. Explicitly,

[TABLE]

At tree level, we have that

[TABLE]

where

[TABLE]

The terms $T_{\rm a}$, $T_{\rm c}$ and $T_{\rm d}$ belong to type 2 discussed in Sec. 2.2, and $T_{\rm b}$ is an example of a term that belongs to type 1. As we noted there, this term can be described with the lower order power spectrum response $\mathcal{R}_{1}$. This can be straightforwardly verified by inserting (cf. Eqs. (2.4) and (2.6))

[TABLE]

into Eq. (2.12) at tree level to get Eq. (B). Finally, the term $T_{\mathcal{R}_{2}}$ corresponds to the trispectrum terms that are captured by $\mathcal{R}_{2}$ in the squeezed limit, i.e.,

[TABLE]

with

[TABLE]

and where we have used the expansions

[TABLE]

with $\varepsilon_{1}=p_{1}/k$ and $\varepsilon_{2}=p_{2}/k$. When restricting Eq. (B.9) to the case $f_{12}=1$, the result simplifies considerably and we obtain Eq. (4.2).

Appendix C Lagrangian response coefficients

An equivalent expansion to Eq. (3.2) can be performed to yield what we call here the “Lagrangian” response coefficients $R_{O}^{L}(k)$. These are in some sense the analog of Lagrangian bias parameters for the bias expansion of the matter power spectrum. Consider an operator $O^{[m]}$ whose lowest-order contribution is at $m$-th order. We will drop the indices and contraction with $\hat{\bm{k}}$ here, as it is of no importance for this discussion. At order $m+1$ in perturbation theory, $O^{(m+1)}$ is given by

[TABLE]

where the sum runs over the same operators that appear in our basis at order $m+1$ (i.e. these operators start at order $m+1$), with coefficients $Z_{O,O^{\prime}}$ which in the EdS approximation are constants. As an example, consider $O^{[2]}=\delta^{2}$. Then, we have $(\delta^{2})^{(3)}$ on the left-hand side, and the sum on the right-hand side runs over cubic operators such as $\delta^{3}$. The last term is the displacement contribution which, at this order, translates the operator $O^{[m]}$ from the Eulerian to the Lagrangian, or fluid coordinate. A similar division of terms into those which involve the operators in our basis, and those corresponding to the displacement from Eulerian to Lagrangian coordinates, can be made analogously at any order greater than $m+1$. We thus write

[TABLE]

Note that $O^{(n)}_{\rm disp}$ only appears for $n\geq m+1$. The displacement contribution can be rigorously defined as containing those terms which are not invariant under a time-dependent, uniform coordinate transformation, $\bm{x}\to\bm{x}+\bm{\xi}(\tau)$. We can now define the Lagrangian response coefficients, by writing the power spectrum expansion at $n$-th order (cf. Eq. (3.5)) as

[TABLE]

The crucial difference to Eq. (3.5), and the Eulerian coefficients, is that the $R_{O}^{L}$ for operators $O^{[m]},\,m<n$, that start at lower order than $n$, only multiply the displacement part of these lower-order operators. The significance of this expansion is that, when taking the angle-averaged squeezed limit of the equal-time $(n+2)$-point function, the displacement terms cancel by symmetry, as they involve single powers of $\hat{\bm{p}}_{a}\cdot\hat{\bm{p}}_{b}$.

Specifically, as shown in Ref. [22], after performing an angle average over all soft modes in the squeezed-limit $(n+2)$-point function, a single isotropic response coefficient $R_{n}^{L}(k)$ remains:

[TABLE]

where the second equality is obtained using Eq. (2.9). The second relation has in fact been derived in Ref. [22], who however did not emphasize that only the $\mathcal{R}_{n}$-type contributions as written in Eq. (2.9) are to be included in the $(n+2)$-point connected correlator (that is, not the contributions which involve lower $\mathcal{R}_{m}$ and PT couplings between soft modes). Nevertheless, Eq. (C) establishes a well-defined relation between the angle-averaged squeezed $(n+2)$-point function and the $n$-th order isotropic Lagrangian response coefficient $R_{n}^{L}(k)$.

Crucially, both $R_{O}^{E}$ and $R_{O}^{L}$ are physical coefficients, which moreover can be related to each other unambiguously in perturbation theory. For this, we separate Eq. (3.5) into operators which start at $n$-th order, and those which start at lower order. For the latter, we then insert Eq. (C.2). In this matching, we only need to consider the boost-invariant (non-displacement) terms. We then obtain

[TABLE]

and hence

[TABLE]

where the sum runs over all lower-order operators in whose $n$-th order expression in perturbation theory the $n$-th order operator $O$ appears.

For completeness, we now give the explicit expressions for the Lagrangian response coefficients in Eq. (C.3) at second order. From Eq. (C.6), it follows that the first-order responses $R_{1}$ and $R_{K}$ are identical in both Eulerian and Lagrangian expansions. The only difference to Eq. (3.2) is that the contributions from nonlinear evolution of the first-order operators disappear. Instead, we only explicitly include terms that displace the first-order terms from the Eulerian (with respect to which the long-wavelength fields $\delta(\bm{p}_{1})$ and $\delta(\bm{p}_{2})$ are defined) to the Lagrangian position, specifically $-s_{(1)}^{k}\partial_{k}K_{ij}^{(1)}$ and $-s_{(1)}^{k}\partial_{k}\delta^{(1)}$. Then, Eq. (3.3) is replaced with

[TABLE]

The explicit kinematic dependence of the Lagrangian response Eq. (C) becomes, for $p_{1}=p_{2}$,

[TABLE]

By matching to the tree-level trispectrum as in the case of the Eulerian expansion, the Lagrangian response coefficients are obtained as

[TABLE]

where we have dropped the superscript “tree” for clarity, and did not repeat the trivially identical results for $R_{1}$ and $R_{K}$ (note however that $R^{L}_{KK}=R^{E}_{KK}$, as well). The Eulerian and Lagrangian coefficients can therefore be related as

[TABLE]

This relation in fact is not restricted to the tree-level responses, but also holds for the fully nonlinear response coefficients, as derived above. In fact, the first relation was already obtained by Ref. [22]. Using that

[TABLE]

where all quantities on the right-hand sides are evaluated at linear order, it is straightforward to verify that Eq. (C.6) indeed yields all relations in Eq. (C).

Figure 5 compares the second-order Eulerian and Lagrangian response coefficients that are different. We see that, in most cases, the Eulerian coefficients are smaller than the Lagrangian ones. This is because the Lagrangian coefficients contain some of the coupling between long-wavelength modes that is not included in the Eulerian coefficients, as has been already noticed in the case of $R_{2}$ and $R_{3}$ in Ref. [22].