Response Approach to the Matter Power Spectrum Covariance

TL;DR

This paper introduces a response-based method to accurately compute the matter power spectrum covariance matrix, extending previous approaches by including loop contributions and validated against simulations.

Contribution

The paper presents a novel response approach that calculates the covariance matrix beyond perturbation theory without fitting parameters, applicable at all orders.

Findings

Captures most of the covariance up to k ~ 1 h/Mpc

Accurate when the softer mode is linear, k < 0.08 h/Mpc

Extends response methods to include loop contributions

Abstract

We present a calculation of the matter power spectrum covariance matrix ${\rm Cov}(\bf{k}_1,\bf{k}_2)$ that uses power spectrum responses to accurately describe the coupling between large- and small-scale modes beyond the perturbative regime. These response functions can be measured with (small-volume) N-body simulations, which is why the response contributions to the covariance remain valid and predictive at all orders in perturbation theory. A novel and key step presented here is the use of responses to compute loop contributions with soft loop momenta, which extends the application of the response approach beyond that of previously considered squeezed $n$-point functions. The calculation presented here does not involve any fitting parameters. When including response-type terms up to 1-loop order in perturbation theory, we find that our calculation captures the bulk of the total covariance as estimated from simulations up to values of $k_1, k_2 \sim 1\ h/{\rm Mpc}$. Moreover, the prediction is guaranteed to be accurate whenever the softer mode is sufficiently linear, ${\rm min}\{k_1,k_2\} \lesssim 0.08\ h/{\rm Mpc}$. We identify and discuss straightforward improvements in the context of the response approach, which are expected to further increase the accuracy of the calculation presented here.

Figures

Equations

$$\displaystyle\operatorname{Cov}({\bm{k}_{1},\bm{k}_{2}})\equiv\big{<}\hat{P}_{m}(\bm{k}_{1})\hat{P}_{m}(\bm{k}_{2})\big{>}-\big{<}\hat{P}_{m}(\bm{k}_{1})\big{>}\big{<}\hat{P}_{m}(\bm{k}_{2})\big{>},$$

$$k_{\mathrm{soft}}\equiv \min \{k_1,k_2\},k_{\mathrm{hard}}\equiv \max \{k_1,k_2\}.$$

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

$$=\frac{1}{2}{\mathcal{R}}_n(k;\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\})P_m(k)(2\pi {)}^3{\delta }_D(\mathbf{k}+{\mathbf{k}}^{\prime }-{\mathbf{p}}_{1\cdots n}),$$

$$\displaystyle\lim_{\{p_{a}\}\to 0}\left(\raisebox{0.0pt}{\parbox[c]{78.0516pt}{\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)$$

$${\mathcal{R}}_n(k;\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\})=\sum _O{R}_O(k){\mathcal{K}}_O(\{{\mu }_{\mathbf{k},{\mathbf{p}}_a}\},\{{\mu }_{{\mathbf{p}}_a,{\mathbf{p}}_b}\},\{p_a/p_b\}).$$

$${\mathcal{R}}_2(k;\mu ,-\mu ,-1,1)$$

$$\mathrm{Cov}({\mathbf{k}}_1,{\mathbf{k}}_2)$$

$$⟨\delta ({\mathbf{k}}_1)\delta ({\mathbf{k}}_2)⟩$$

$$⟨\delta ({\mathbf{k}}_a)\delta ({\mathbf{k}}_b)\delta ({\mathbf{k}}_c)\delta ({\mathbf{k}}_d){⟩}_c$$

$${\hat{P}}_m(k_1)=V_f{\int }_{V_s(k_1)}\frac{{\mathrm{d}}^3\mathbf{k}}{V_s(k_1)}\delta (\mathbf{k})\delta (-\mathbf{k}),$$

$${\mathrm{Cov}}^{\mathrm{G}}(k_i,k_j)=\frac{2}{N_k}{P}_m(k_i{)}^2{\delta }_{ij},$$

$$\displaystyle\operatorname{Cov}^{\rm NG}(k_{1},k_{2},\mu_{12})=\operatorname{Cov}^{\rm NG}_{\rm tree}(k_{1},k_{2},\mu_{12})+\operatorname{Cov}^{\rm NG}_{\rm 1loop}(k_{1},k_{2},\mu_{12})+\big{(}{\rm higher\ loops}\big{)}.$$

$$\mathrm{Cov}(k_1,k_2,{\mu }_{12})$$

$${\mathrm{Cov}}^{\ell }(k_1,k_2)$$

$$V{\mathrm{Cov}}_{\mathrm{SPT}-\mathrm{tree}}^{\mathrm{NG}}(k_1,k_2,{\mu }_{12})$$

$$\displaystyle\raisebox{0.0pt}{\parbox[c]{57.0933pt}{\includegraphics[scale={0.8}]{diag_R2tree.pdf}}}=\>$$

$${\mathrm{Cov}}_{\mathcal{R}\text{-tree}}^{\mathrm{NG}}(k_1,k_2,{\mu }_{12})$$

$${\mathrm{Cov}}_{\text{st-tree}}^{\mathrm{NG}}(k_1,k_2,{\mu }_{12})=\{\begin{array}{ll}{\mathrm{Cov}}_{\text{SPT-tree}}^{\mathrm{NG}}(k_{\text{hard}},k_{\text{soft}},{\mu }_{12})& ,\mathrm{if}{k}_{\text{soft}}>f_{\text{sq}}{k}_{\text{hard}}\\ {\mathrm{Cov}}_{\mathcal{R}\text{-tree}}^{\mathrm{NG}}(k_{\text{hard}},k_{\text{soft}},{\mu }_{12})& ,\mathrm{otherwise}\end{array}.$$

$${\mathrm{Cov}}_{1loop,P_{\mathrm{L}}(p)}^{\mathrm{NG}}(k_1,k_2,{\mu }_{12})\stackrel{p\ll k_1,k_2}{\propto }\int {\mathrm{d}}^{3}p\left[\dots \right]P_{\mathrm{L}}(p)P_{\mathrm{L}}(k_i)P_{\mathrm{L}}(k_j)P_{\mathrm{L}}(k_k),$$

$${\mathrm{Cov}}_{1loop,[P_{\mathrm{L}}(p){]}^2}^{\mathrm{NG}}(k_1,k_2,{\mu }_{12})\stackrel{p\ll k_1,k_2}{\propto }\int {\mathrm{d}}^{3}p\left[\dots \right][P_{\mathrm{L}}(p){]}^2{P}_{\mathrm{L}}(k_i)P_{\mathrm{L}}(k_j).$$

$$=$$

$$+$$

$$\displaystyle\raisebox{0.0pt}{\parbox[c]{101.178pt}{\includegraphics[scale={0.8}]{diag_cov1l.pdf}}}+(\text{perm.})=$$

$$=2\int \frac{{\mathrm{d}}^{3}p}{(2\pi {)}^3}[P_{\mathrm{L}}(p){]}^2{\mathcal{R}}_2(k_1,{\mu }_1,-{\mu }_1,-1,1)P_m(k_1){\mathcal{R}}_2(k_2,{\mu }_2,-{\mu }_2,-1,1)P_m(k_2)$$

$$\displaystyle=\frac{2P_{m}(k_{1})P_{m}(k_{2})}{(2\pi)^{3}}\Bigg{[}\int_{0}^{p_{\text{max}}}p^{2}[P_{\rm L}(p)]^{2}{\rm d}p\Bigg{]}$$

$$\times {\int }_{-1}^1\mathrm{d}{\mu }_1{\int }_0^{2\pi }\mathrm{d}\phi {\mathcal{R}}_2(k_1,{\mu }_1,-{\mu }_1,-1,1){\mathcal{R}}_2(k_2,{\mu }_2,-{\mu }_2,-1,1),$$

$$\displaystyle\operatorname{Cov}^{{\rm NG}}_{\mathcal{R}\text{-1loop}}(k_{1},k_{2},\mu_{12})=V^{-1}\frac{2P_{m}(k_{1})P_{m}(k_{2})}{(2\pi)^{2}}\Bigg{[}\int_{0}^{p_{\text{max}}}p^{2}[P_{\rm L}(p)]^{2}{\rm d}p\Bigg{]}$$

$$\displaystyle\hskip 56.9055pt\times\Bigg{[}2\mathcal{A}_{1}\mathcal{A}_{2}+\frac{2}{5}\mathcal{B}_{1}\mathcal{B}_{2}\mathcal{P}_{2}(\mu_{12})+\frac{2}{5}\big{(}\mathcal{A}_{1}\mathcal{C}_{2}+\mathcal{A}_{2}\mathcal{C}_{1}\big{)}+\frac{4}{35}\big{(}\mathcal{B}_{1}\mathcal{C}_{2}+\mathcal{B}_{2}\mathcal{C}_{1}\big{)}\mathcal{P}_{2}(\mu_{12})$$

$$\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +\frac{2}{35}\big{(}1+2\mathcal{P}_{2}(\mu_{12})^{2}\big{)}\mathcal{C}_{1}\mathcal{C}_{2}\Bigg{]},$$

$$p_{\text{max}}=\min \{f_{\text{sq}}{k}_1,f_{\text{sq}}{k}_2,k_{\text{NL}}\}.$$

$$\frac{P_{\mathrm{L}}(k_{\text{soft}}){\int }_0^{p_{\text{max}}}{p}^2{P}_{\mathrm{L}}(p)\mathrm{d}p}{{\int }_0^{p_{\text{max}}}{p}^2[P_{\mathrm{L}}(p){]}^2\mathrm{d}p},$$

$$\text{mbox}Eq.(\text{ref}eq:ratioest)\approx \{0.27,0.11,0.02\}\text{mbox}for{k}_{\text{soft}}=\{0.1,0.3,1\}h{\mathrm{Mpc}}^{-1},$$

$$\displaystyle\raisebox{0.0pt}{\parbox[c]{101.178pt}{\includegraphics[scale={0.8}]{diag_cov2l.pdf}}}+(\text{perm.})=$$

Full text

1 Introduction

The unprecedented statistical precision that upcoming large-scale structure surveys are expected to attain requires cosmologists to develop equally precise methods to predict the various observables. The simplest and most widely applied way to describe the statistical information encoded in the large-scale structure is via the $2$-point correlation function [2]. This includes the $2$-point galaxy correlation function, or $2$-point correlations of cosmic shear maps in the case of gravitational lensing, as well as their cross-correlation. The starting point to predicting both these observables is the $2$-point correlation function of matter $\xi_{m}$, or its Fourier transform, the power spectrum $P_{m}$. The matter power spectrum is very well understood in the context of gravity-only N-body simulations (that is, neglecting baryonic effects on the total matter distribution). The simulation requirements for a given pre-specified precision have been studied in Ref. [3], and simulations have also allowed for the calibration of semi-analytical models such as Halofit [4] and construction of efficient interpolations such as that of the Coyote project [5]. Baryonic effects are known to have an impact on the small-scale matter power spectrum [6, 7, 8] and work on modeling these effects has also been carried out recently [9, 10, 11, 12].

An accurate model of the matter power spectrum alone is however insufficient to properly exploit upcoming surveys, especially when inferring cosmological parameter values from the data, for which one also needs accurate determinations of the covariance matrix of the power spectrum,

[TABLE]

in order to quantify the statistical error of the measurements. In the equation above, angle brackets denote ensemble averages and $\hat{P}_{m}(\bm{k})$ is an estimate of the matter power spectrum within some wavemode bin centered at $\bm{k}$. The power spectrum covariance, hereafter simply referred to as matter covariance, therefore measures the correlation between the power at wavemodes $\bm{k}_{1}$ and $\bm{k}_{2}$. For Gaussian initial conditions, different Fourier modes evolve independently in the linear stages of structure formation. In this regime, only the diagonal terms are non-vanishing and they are trivially related to the matter power spectrum itself. At later stages, nonlinear structure formation effectively couples different Fourier modes, which leads to important off-diagonal ($k_{1}\neq k_{2}$) terms in $\operatorname{Cov}({\bm{k}_{1},\bm{k}_{2}})$ through a special configuration of the matter trispectrum (the Fourier transform of the $4$-point correlation function), which we will describe in more detail below. The matter trispectrum is a cumbersome quantity to predict, which is why our current knowledge of the covariance matrix is far poorer than that of the power spectrum. Given that inaccurate estimates of the covariance matrix can result in wrong interpretations of the data (see e.g. Refs. [13, 14, 15, 16, 17]), this naturally motivates research in obtaining accurate theoretical predictions of the covariance, including its dependence on cosmological parameters [18, 19, 20, 21], as well as on baryonic effects. This is important for current data sets, but even more so for future large-volume surveys such as Euclid [22], LSST [23] and DESI [24].

One frequently employed tool to estimate $\operatorname{Cov}({\bm{k}_{1},\bm{k}_{2}})$ are direct estimates of the covariance via Eq. (1.1) using large sets of N-body simulations [25, 26, 27, 28, 29]. This requires performing thousands of N-body simulations in order to obtain sufficient signal-to-noise in the covariance, which makes these estimates extremely costly in terms of computational resources. Estimating the covariance for many different sets of cosmological parameters therefore becomes prohibitive, as does a realistic modeling of baryonic processes.

A complementary approach to simulation estimates is to use perturbation theory. Reference [30] presented a calculation of the trispectrum in the covariance configuration at the 1-loop level, based on the Effective Field Theory (EFT) of large-scale structure (see Ref. [31] for a review). The main limitation of such perturbative approaches is that they are only applicable to sufficiently large scales $k\lesssim k_{\text{NL}}\approx 0.3\>h\,{\rm Mpc}^{-1}$ (at $z=0$), which limits their usefulness in the analysis of data on smaller scales.111We define the nonlinear scale $k_{\text{NL}}$ at a given redshift $z$ through $k_{\text{NL}}^{3}P_{\rm L}(k_{\text{NL}},z)/(2\pi^{2})=1$. There have also been attempts to develop semi-analytical phenomenological models of the covariance matrix on scales $k\ >k_{\text{NL}}$ [32, 33, 10, 1], but these typically involve simplifying assumptions and/or free parameters that need to be tuned to match other covariance estimates, usually simulation-based ones, for any given cosmology considered. Moreover, systematic errors made in these phenomenological estimates are not under rigorous control, and can only be estimated through comparison with simulation-based estimates.

In this paper, our goal is to describe a calculation of the covariance matrix that combines the merits of the simulation- and perturbation theory-based approaches. More concretely, in our approach, we use perturbation theory to identify the mode-coupling terms of the non-Gaussian covariance that can be resummed with simulation-calibrated power spectrum responses (see Ref. [34] for an in-depth discussion and Sec. 2.1 below for an overview). The power spectrum responses measure the fractional change of the local power spectrum in the presence of long-wavelength perturbations, and they can be measured accurately and non-perturbatively with only a few relatively small-volume simulations [35, 36, 37, 38, 39, 40, 41, 42, 43, 44] (with small computational cost compared to that of fully numerical covariance estimates). The types of mode-coupling interactions that are captured by responses are therefore those that describe the coupling between long- and short-wavelength modes. All non-response type terms are calculated using standard perturbation theory (SPT), leaving the whole calculation free from fitting parameters. Moreover, we can use higher-order perturbation theory to estimate the systematic error made in the calculation.

In Ref. [34], we have presented a response-based calculation of the non-Gaussian covariance at tree level in the squeezed regime, i.e., when one of the modes is linear and sufficiently smaller than the other, which can take any other value: $k_{\rm soft}\ll k_{\rm hard}$, $k_{\rm soft}\ll k_{\text{NL}}$, for any $k_{\rm hard}$, where

[TABLE]

This represents an application of the well-known relation between responses and squeezed-limit correlators (in this case, the squeezed trispectrum). In this paper, we go beyond Ref. [34] as we demonstrate how to use responses to resum interaction vertices that involve internal soft-loop momenta, thereby permitting an efficient and accurate evaluation of the covariance for any values of $k_{1}$, $k_{2}$, including (quite crucially) cases in which $k_{1}\approx k_{2}$. This constitutes an example of the use of responses in the calculation of non-squeezed $n$-point functions.

After establishing some notation and summarizing the definitions of power spectrum responses and covariance in Sec. 2, the steps taken in this paper can be outlined as follows:

In Sec. 3, working at tree level in standard perturbation theory, we show how to stitch together the standard perturbation theory and response-based results presented first in Ref. [34] to fully describe the matter covariance in the regime where $k_{\rm soft}\ll k_{\text{NL}}$ and any $k_{\rm hard}$. 2. 2.

Section 4 is devoted to the novel application of responses to calculate loop interactions involving soft loop momenta but fully nonlinear external momenta. Here, we work explicitly at the 1-loop level in perturbation theory, but also describe how to account for higher loops. 3. 3.

We compare our model results to simulation-based estimates of the angle-averaged covariance in Sec. 5. The level of agreement we find across the range of scales probed by the simulations suggests that the calculation presented here (which has no free parameters) captures the majority of the total matter covariance. In Sec. 6, we also look at the prediction for the dependence of $\operatorname{Cov}(\bm{k}_{1},\bm{k}_{2})$ on the angle between the two modes.

The covariance calculation presented here, being based on a well-defined theoretical framework, is particularly useful as it allows us to determine exactly which contributions are being left out at a given point $(k_{1},k_{2})$; in particular, these are higher-loop terms, and certain non-response-type terms. This can be used to estimate the error on the covariance prediction (the error on the error of the matter power spectrum), as well as to guide further developments. These, as well as other concluding remarks are the subject of Sec. 7. In particular, Fig. 8 summarizes which parts of $(k_{1},k_{2})$-space are already completely captured by our calculation, and which parts can benefit from further work.

In Appendix A, we spell out the Feynman rules of cosmological perturbation theory as used in the paper. We collect the expressions to evaluate response functions, as well as the corresponding non-Gaussian covariance terms in Appendices B and C, respectively. The criterion to distinguish between squeezed and non-squeezed configurations is determined in Appendix D. In Appendix E, we demonstrate explicitly the equivalence between 1-loop covariance terms in standard perturbation theory and the response-based description. Finally, in Appendix F, we compare our covariance calculation with the prediction of the model presented in Ref. [1].

In this paper, we assume a flat $\Lambda{\rm CDM}$ cosmology for all numerical results, with the following parameters: $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$, $\sum m_{\nu}=0$. These are the same as those used in Ref. [45] in their estimates of the covariance matrix from simulations, with which we shall compare our results with. Further, in our results below, we use the CAMB code [46] and the Coyote emulator [5] to compute the linear and the nonlinear matter power spectrum, respectively.

2 Definitions and notation

2.1 Power spectrum responses

In this section, we briefly recap the definition and physical content of power spectrum responses, and display the equations that we use in the remainder of the paper. We refer the reader to Ref. [34] for a detailed description of the response formalism. Throughout, we only consider equal-time matter correlators, and will not write the time argument explicitly to ease the notation. The Feynman rules of cosmological perturbation theory (which shall be particularly useful in our considerations below) are summarized in Appendix A. Further, we denote magnitudes of vectors as $k=|\bm{k}|$ and adopt a shorthand notation for the sum of vectors: $\bm{k}_{12\cdots n}=\bm{k}_{1}+\bm{k}_{2}+\cdots+\bm{k}_{n}$. The $n$-th order matter power spectrum response $\mathcal{R}_{n}$ corresponds to the following interaction vertex

[TABLE]

which is interpreted as the response of the nonlinear power spectrum of the small-scale (hard) mode $\bm{k}$ to the presence of $n$ long-wavelength (soft) modes $\bm{p}_{1},...,\bm{p}_{n}$. The dashed blob is thus meant to account for the fully evolved nonlinear matter power spectrum $P_{m}(k)$ and all its possible interactions with the $n$ long wavelength perturbations (including loop interactions — it is thus a resummed vertex). In our notation, $\lim_{\{p_{a}\}\to 0}$ signifies that we only retain the leading contribution in the limit in which all soft momenta are taken to zero. The response $\mathcal{R}_{n}$ depends on the scale $k$, as well as on the angles between the $n$ soft modes and their angles with $\bm{k}$. The response also depends on the ratios of soft wavenumbers, but not on their absolute values.222These responses are to be distinguished from those measured in Refs. [47, 48], which correspond to the derivative of the nonlinear power spectrum with respect to the initial power spectrum (i.e., not to the presence of individual large-scale perturbations).

The diagrammatic representation of $\mathcal{R}_{n}$ helps to understand the connection of the power spectrum response with the squeezed limit of the $(n+2)$-point matter correlation function. Explicitly, attaching power spectrum propagators to the soft momentum lines in Eq. (2.1), we can write

[TABLE]

where the $n!$ factor accounts for the permutations of the $\bm{p}_{a}$. The subscript c denotes connected correlators, while the subscript ${}_{\mathcal{R}_{n}}$ in the $(n+2)$-connected correlator indicates that only certain contributions to the correlation function are captured by $\mathcal{R}_{n}$. There are further response-type contributions which are not included in $\langle\delta(\bm{k})\delta(\bm{k}^{\prime})\delta(\bm{p}_{1})\cdots\delta(\bm{p}_{n})\rangle_{c,\mathcal{R}_{n}}$. These terms are however completely determined by lower order responses, $\mathcal{R}_{m}$, $1\leq m\leq n$, in conjunction with perturbation theory kernels. All other terms that contribute to $\langle\delta(\bm{k})\delta(\bm{k}^{\prime})\delta(\bm{p}_{1})\cdots\delta(\bm{p}_{n})\rangle_{c}$ are small in the squeezed limit.

As described in detail in Ref. [34], the $\mathcal{R}_{n}$ can be expanded in terms of all local gravitational observables associated with the $n$ long-wavelength modes, to a given order in perturbations. These observables, or operators $O$ (which can be constructed using either Lagrangian or Eulerian coordinates) form a basis $\mathcal{K}_{O}$ that unequivocally specifies all of the angular dependence of $\mathcal{R}_{n}$:

[TABLE]

At any given order, there are different equivalent decompositions of the sum in Eq. (2.3), which translate into different expressions for the $\mathcal{K}_{O}$. For instance, Ref. [49] displays an alternative, but mathematically equivalent decomposition at $n=1$ and $n=2$. Here, we will use the Eulerian decomposition described in the main text of Ref. [34]. In this paper, we will only explicitly need the second-order response $\mathcal{R}_{2}\equiv\mathcal{R}_{2}(k,\mu_{1},\mu_{2},\mu_{12},p_{1}/p_{2})$, which is a function of the hard mode $k$ (and time), the cosine angles $\mu_{1}=\bm{p}_{1}\cdot\bm{k}/(p_{1}k)$, $\mu_{2}=\bm{p}_{2}\cdot\bm{k}/(p_{2}k)$, $\mu_{12}=\bm{p}_{1}\cdot\bm{p}_{2}/(p_{1}p_{2})$ and the ratio $p_{1}/p_{2}$. More specifically, for the application to the matter covariance, the relevant kinematic configuration corresponds to $\mu_{1}=\mu$, $\mu_{2}=-\mu$, $\mu_{12}=-1$ and $p_{1}/p_{2}=1$, in which case the expression of $\mathcal{R}_{2}$ can be given as

[TABLE]

where $\mathcal{P}_{\ell}$ is the Legendre polynomial of order $\ell$ and the second equality serves to define the functions $\mathcal{A}(k)$, $\mathcal{B}(k)$ and $\mathcal{C}(k)$, which help to simplify some notation below.

The coefficients $R_{O}(k)$ are called response coefficients and they correspond to the response of the local small-scale power spectrum to specific configurations of the long-wavelength perturbations. At tree level, all $R_{O}(k)$ can be derived by matching the definition of $\mathcal{R}_{2}$ to the squeezed four-point function, in the sense of Eq. (2.2) (see Ref. [34] for the explicit steps of this derivation). In the nonlinear regime of structure formation, the response coefficients must be determined with the aid of N-body simulations. The first three isotropic response coefficients, $R_{1}$, $R_{2}$ and $R_{3}$, where $R_{n}(k)\equiv n!R_{\delta^{n}}(k)$, have already been measured accurately with separate universe simulations [39] (see also Refs. [35, 36, 37, 38]). In these simulations, the presence of an exactly uniform density perturbation in the simulation volume is simulated by using the equivalence to following structure formation in a spatially curved Friedmann-Roberston-Walker spacetime [50]. The remaining coefficients have so far not been measured in N-body simulations due to complications associated with how to model the presence of these anisotropic long-wavelength perturbations [51, 44]. In this paper, we combine the simulation measurements of the isotropic $R_{O}(k)$ with the nonlinear extrapolation of the anisotropic ones put forward in Ref. [34] (see Fig. 1 there for the numerical results). The explicit expressions for all $R_{O}(k)$ used in this paper are given in Appendix B.

2.2 Matter power spectrum covariance

Let $\delta(\bm{x})$ denote the fractional matter density contrast at $\bm{x}$, $\delta(\bm{k})$ its Fourier transform (distinguished by their arguments) and $\hat{P}(\bm{k})$ the estimated power spectrum in a wavenumber bin centered on the Fourier mode $\bm{k}$, in a total survey volume $V$. The matter power spectrum covariance $\operatorname{Cov}(\bm{k}_{1},\bm{k}_{2})$ measures the correlation between the power spectrum of the modes $\bm{k}_{1}$ and $\bm{k}_{2}$ and is defined as (rewriting Eq. (1.1))

[TABLE]

where

[TABLE]

define the matter power spectrum $P_{m}$ and trispectrum $T_{m}$, respectively. The latter contributes to the covariance in the so-called parallelogram configuration, $\bm{k}_{b}=-\bm{k}_{a}$, $\bm{k}_{d}=-\bm{k}_{c}$. Note also that so far we have not restricted ourselves to the covariance of the angle-averaged power spectrum (see Eq. (2.8) below), i.e., we allow for the covariance to depend on the angle between the two wavemodes, $\mu_{12}=\bm{k}_{1}\cdot\bm{k}_{2}/(k_{1}k_{2})$.

As indicated in Eq. (2.5), the two terms in the second line are broadly referred to as the Gaussian and non-Gaussian parts of the covariance, on which we comment further below. Before proceeding however, we note that, for a finite survey, Eq. (2.5) is missing an important additional non-Gaussian contribution. This is the so-called super-sample covariance term [52, 53, 40, 41, 54, 55], which accounts for the coupling of Fourier modes inside the observed surveyed region with density fluctuations whose wavelength is larger than the typical size of the survey. Formally, this term arises from the convolution of the matter trispectrum with the survey window function. The behavior of the super-sample term is well understood and can be described using the first-order power spectrum response $\mathcal{R}_{1}$. Below, we shall compare our covariance results with estimates from standard N-body simulations, which do not include fluctuations on scales larger than the simulation box, and are therefore unable to measure the super-sample term. For this reason, we do not consider the super-sample contribution in our results, but note that its inclusion is straightforward.

For Gaussian initial conditions, the non-Gaussian contribution is only induced by nonlinear structure formation, so that the Gaussian term dominates at early times. This term correlates the power spectra of two modes only if the modes have the same magnitude and are exactly aligned, $\mu_{12}=1$ or anti-aligned $\mu_{12}=-1$. In the literature, the case of angle-averaged power spectra is that which is most commonly considered:

[TABLE]

where the integral is taken over a spherical shell of radius $k_{1}$ and width $\Delta k$, $V_{s}(k_{1})=4\pi k_{1}^{2}\Delta k$, and $V_{f}=(2\pi)^{3}/V$ is the volume of a Fourier cell where $V$ is the total survey volume. In this case, the Gaussian part of the covariance becomes

[TABLE]

where $i,j$ label bins in wavenumber, $N_{k}=V_{s}(k_{i})/V_{f}$ is the number of Fourier modes that are averaged over in a given bin, and the Kronecker delta $\delta_{ij}$ ensures that the Gaussian term contributes only to the diagonal of the angle-averaged covariance matrix, $k_{i}=k_{j}$. Note that the Gaussian covariance depends on the size of the $k$-bins in which the spectra are measured.

The non-Gaussian part of the covariance measures the coupling between Fourier modes that is induced by nonlinear structure formation at late times. This term can also be present in the initial conditions, due to primordial non-Gaussianity, but we do not consider this case here. In the context of standard perturbation theory [56], the non-Gaussian covariance $\operatorname{Cov}^{\rm NG}$ (or the parallelogram matter trispectrum, in the sense of Eq. (2.5)) can be expanded into its tree-, 1-loop and higher-order loop contributions,

[TABLE]

This part of the covariance, which is by far the most challenging to measure and predict, is that which we wish to address specifically in this paper. The main idea behind the calculation that we perform here is that, at a given order in standard perturbation theory, we identify the mode-coupling terms that describe the interactions between hard and soft modes and resum them using power spectrum responses; all other terms can be computed as in standard perturbation theory. Up to 1-loop order, we will see how such a combination of the SPT and response approaches is capable of capturing a substantial part of the total non-Gaussian covariance.

Before proceeding with the more rigorous description of the calculation of $\operatorname{Cov}^{\rm NG}(k_{1},k_{2},\mu_{12})$ in the next sections, we collect here some of the notation that is used throughout. We will use the words standard, response and stitched to refer, respectively, to the SPT-based, response-based and their combined contributions to the total non-Gaussian covariance. More specifically:

Standard tree and standard 1-loop, which we represent as $\operatorname{Cov}^{\rm NG}_{\text{SPT-tree}}$ and $\operatorname{Cov}^{\rm NG}_{\text{SPT-1loop}}$, refers to the standard perturbation theory calculation (at the corresponding tree- or 1-loop levels) that does not employ any response vertices and that loses predictivity whenever any of the external momenta approach $k_{\text{NL}}$. 2. 2.

Response tree and response 1-loop, which we represent as $\operatorname{Cov}^{\rm NG}_{\mathcal{R}\text{-tree}}$ and $\operatorname{Cov}^{\rm NG}_{\mathcal{R}\text{-1loop}}$, refers to the mode-coupling terms between hard and soft modes that exist at tree and 1-loop levels, respectively, and that can be calculated with response vertices. These contributions lose predictivity if the soft modes involved approach $k_{\text{NL}}$, but are otherwise valid for any value of the hard modes, including in the nonlinear regime. 3. 3.

Stitched tree, which we represent as $\operatorname{Cov}^{\rm NG}_{\text{st-tree}}$, refers to a specific combination, or stitching, of the standard and response contributions. The details of this stitching will be clarified in the sections below. Note that while we apply this procedure at tree level here, the stitching can in principle be applied at any order.

It is also useful to organize the angular dependence of the covariance into multipoles as

[TABLE]

The case of the monopole $\ell=0$, $\operatorname{Cov}^{\rm NG,\ell=0}(k_{1},k_{2})$, corresponds to the covariance matrix of angle-averaged spectra, which is the case we shall mostly focus on (with the exception of Sec. 6, where we present the $\ell=2$ and $\ell=4$ predictions).

Finally, throughout we use $k_{\rm soft}$ and $k_{\rm hard}$ to denote the softest and the hardest of the two $k$-modes of the covariance, i.e., $k_{\rm soft}={\rm min}\{k_{1},k_{2}\}$ and $k_{\rm hard}={\rm max}\{k_{1},k_{2}\}$.

3 The stitched non-Gaussian covariance at tree level

In this section, we propose a calculation of the tree-level covariance that combines the standard perturbation theory result, valid only if $k_{1},k_{2}\ll k_{\text{NL}}$, with the response-based description first put forward in Ref. [34], which effectively extends the validity of the calculation to $k_{\rm soft}\ll k_{\text{NL}}$, but any $k_{\rm hard}$, including in the nonlinear regime.

The SPT result for the tree-level non-Gaussian covariance [57] is given by (see e.g. Appendix B of Ref. [34] for explicit expressions for all terms of the general tree-level trispectrum)

[TABLE]

where $F_{2}$ and $F_{3}$ are the symmetrized second- and third-order standard perturbation theory kernels [56]. In the literature, the above equation is sometimes written in a simpler way that anticipates the angle averages that are subsequently taken, but here we opted to remain general. This standard tree level result is only expected to be a good approximation to the full non-Gaussian covariance when both $k_{1}$ and $k_{2}$ are in the linear regime, $\max\{k_{1},k_{2}\}\ll k_{\text{NL}}$. However, as noted already above, if $k_{\rm soft}$ is sufficiently linear and smaller than $k_{\rm hard}$, then it is possible to extend the regime of validity of the tree level calculation to nonlinear values of $k_{\rm hard}$ by making use of the response $\mathcal{R}_{2}$. This follows from noting that, by taking $n=2$, $\bm{p}_{1}=-\bm{p}_{2}=\bm{k}_{\rm soft}$ and $\bm{k}=-\bm{k}^{\prime}=\bm{k}_{\rm hard}$ in Eq. (2.2), one obtains precisely the squeezed limit of the connected $4$-point function (or trispectrum) in the covariance configuration. The following equation provides a schematic picture of this relation:

[TABLE]

where $p_{1},p_{2}$ are understood as much softer than $k,k^{\prime}$. That is, at tree level, the $\mathcal{R}_{2}$ vertex captures the coupling described by one $F_{3}$ kernel in one diagram and two $F_{2}$ kernels in the other. The above equation is shown and used explicitly in Sec. 4.2 of Ref. [34] to derive the shape of $\mathcal{R}_{2}$ at tree level. By replacing the tree-level $\mathcal{R}_{2}$ with its simulation-calibrated shape, then one effectively extends (or resums) the interactions on the right-hand side of Eq. (3.2) to all orders in perturbation theory in the hard mode. Referring the reader to Ref. [34] for more details, here we limit ourselves to showing the final response-based result, which is given by

[TABLE]

where the next-to-leading corrections come from non-response type interactions that are suppressed in the squeezed regime, as well as loop corrections in the soft mode that enter when $k_{\rm soft}$ is no longer much smaller than $k_{\text{NL}}$. We can now put the above two equations together to construct our stitched tree-level covariance as follows333When writing Eqs. (3), (3.3) and (3.4) and all following relations, we have implicitly averaged over the $k$-bins used to estimate $\hat{P}_{m}$. Contrary to the Gaussian case, the non-Gaussian covariance does not depend explicitly on the $k$-bin widths, and hence we omit this averaging to shorten the notation.

[TABLE]

Thus, we use the standard tree-level expression in non-squeezed configurations, but switch to the response tree-level result in squeezed ones. The value of $f_{\text{sq}}$ controls the transition from the non-squeezed to the squeezed regime, and its optimal choice corresponds to a trade-off between two demands. On the one hand, the response prediction is only accurate up to corrections of order $f_{\text{sq}}^{2}$ (Eq. (3.3)), and hence $f_{\text{sq}}$ should be chosen as small as possible. On the other hand, for $k_{\text{hard}}$ that approach or even exceed $k_{\text{NL}}$, the response prediction is more accurate than the SPT-tree prediction, which makes larger values of $f_{\text{sq}}$ beneficial (to maximize the volume in $(k_{1},k_{2})$-space where the response-based result is used). In Appendix D, we describe a procedure to determine the largest value of $f_{\text{sq}}$ that ensures a given accuracy of Eq. (3.4) based on the standard tree-level covariance. From the exercise performed in Appendix D, we take our fiducial choice to be $f_{\text{sq}}=0.5$. From hereon in this paper we shall therefore dub configurations with $k_{\text{soft}}<k_{\text{hard}}/2$ as squeezed.

Our stitched tree-level result of Eq. (3.4) is shown in Fig. 1, for the angle-averaged case ($\ell=0$ in Eq. (2.2)). The matrix is shown as a color plot in the left panel. The two panels on the right show each a slice at fixed $k_{2}$ of the matrix on the left (red), as well as the standard (blue) and response (green) results. The upper right panel represents a slice at $k_{2}<k_{\text{NL}}$, while the lower right shows a slice at $k_{2}>k_{\text{NL}}$. The dashed black lines draw the boundary between the squeezed and non-squeezed limits. When both modes are smaller than $k_{\text{NL}}\approx 0.3\>h\,{\rm Mpc}^{-1}$ (at $z=0$), the color plot displays a smooth transition from the standard- to the response-based results, as it should be by definition. One does note, however, that there are noticeable discontinuities at the junction between the two cases, when both $k$ values are above $k_{\text{NL}}$. This is expected, as the standard tree-level result does not include any nonlinear corrections, while the response prediction does. Furthermore, even if, say, $k_{2}\ll k_{1}$, but $k_{2}$ is of order $k_{\text{NL}}$ or larger, then we do not expect the response tree-level result by itself to be a good description of the covariance. This is because loop contributions become non-negligible in that regime (cf. $\mathcal{O}(k_{\text{soft}}^{2}/k_{\text{NL}}^{2})$ corrections in Eq. (3.3)). We will see below that, in this regime, our tree-level result is negligible compared to the contribution from 1-loop terms, and hence, the unphysical discontinuity in the stitched-tree-level contribution does not affect the much larger total covariance.

4 Non-Gaussian covariance at the 1-loop level

We now extend the calculation of the non-Gaussian part of the covariance by working at 1-loop level (cf. Eq. (2.10)). To do so, we introduce a new concept in large-scale structure perturbation theory which uses responses to describe the coupling between soft internal loop momenta and hard external modes, thereby going beyond the so-far considered application of responses to describe the coupling of soft external modes with hard external modes.

The 1-loop covariance $\operatorname{Cov}^{\rm NG}_{{\rm 1loop}}(k_{1},k_{2},\mu_{12})$ has contributions from nine types of diagrams (see e.g. Fig. 4 of Ref. [30]). In the limit of soft loop momenta $p$, i.e. $p\ll k_{1},k_{2}$, it can be shown that six of these diagrams are linear in $P_{\rm L}(p)$, i.e., they are of the form

[TABLE]

where $k_{i},k_{j},k_{k}$ are of the order of the external modes $k_{1},k_{2}$. Here and below, the dots in square brackets denote perturbation theory kernels $F_{n}$ that we do not write for brevity. On the other hand, the remaining three diagrams contain contributions that involve two powers of $P_{\rm L}(p)$, i.e.,

[TABLE]

These latter three diagrams can be represented exactly as a single diagram that involves two tree-level $\mathcal{R}_{2}$ vertices:

[TABLE]

where the role of the long-wavelength perturbations is played by internal loop momenta. The above equation is written in terms of perturbation theory kernels in Eq. (E) in Appendix E, where we demonstrate explicitly which standard 1-loop covariance terms are captured by the response approach. Based on the Feynman rules augmented with responses (cf. Appendix A), this represents a linking of two diagrammatic representations of the tree-level response $\mathcal{R}_{2}^{\text{tree}}$ (cf. Eq. (3.2)). As it was the case at tree level in the covariance, the generalization of the tree-level $\mathcal{R}_{2}$ to its simulation-calibrated expressions effectively extends the validity of the calculation of these 1-loop terms to nonlinear values of the hard modes $k_{1},k_{2}$. A point that is worth emphasizing is that now both external modes can be of comparable size, i.e., this constitutes an application of the response formalism beyond the commonly used application to describe squeezed-limit correlation functions.

At this point, it may not be clear why capturing the terms in Eq. (4.2) through responses yields a significant advantage. However, as we will comment on in the next subsection, the other 1-loop terms in Eq. (4.1) are suppressed relative to those in Eq. (4.2) if $k_{1},k_{2}$ are sufficiently large. Further, beyond the limit of soft loop momentum, the contributions from the opposite limit, i.e. $p\gg k_{1},k_{2}$, are suppressed due to mass-momentum conservation. By expressing Eq. (4.2) in terms of responses, one can therefore capture a substantial part of the covariance in $(k_{1},k_{2})$-space. Moreover, in the squeezed regime where $k_{\text{soft}}\ll k_{\text{hard}}$, additional response diagrams allow us to capture all 1-loop terms that are leading order in $k_{\rm soft}/k_{\rm hard}$, although we do not calculate these in this paper. We return to this in Sec. 4.1.

By the Feynman rules, the response diagram in Eq. (4) can then be written as

[TABLE]

where $\mu_{1}=\bm{k}_{1}\cdot\bm{p}/(k_{1}p)$, $\mu_{2}=\bm{k}_{2}\cdot\bm{p}/(k_{2}p)$ and the factor of 2 comes from the two possible ways of connecting the loops. In the above equation, we have implicitly fixed the direction of $\bm{k}_{1}$, which means the polar integral of the loop momentum $\bm{p}$ is done w.r.t. $\bm{k}_{1}$.

We note for completeness that Eq. (4), as written, does not correctly describe cases where $\bm{k}_{1}$ and $\bm{k}_{2}$ are very nearly parallel such that $|\bm{p}\pm\bm{k}_{1}\pm\bm{k}_{2}|\approx p$. These terms can be straightforwardly included using responses as well (see Appendix E for more details). However, after angle-averaging, their contribution to low-order multipoles is suppressed by $(p/k_{i})^{2}$ and thus becomes negligible. We therefore do not include them in the main text. After performing the two angle integrals in Eq. (4) one arrives at

[TABLE]

where $\mathcal{A}_{i}\equiv\mathcal{A}(k_{i})\ (i=1,2)$, and similarly for $\mathcal{B}_{i}$ and $\mathcal{C}_{i}$ (cf. Eq. (2.1) for the definition of the $\mathcal{A}(k)$, $\mathcal{B}(k)$ and $\mathcal{C}(k)$ in terms of linear combinations of the response coefficients $R_{O}(k)$).

A key issue to address before evaluating Eq. (4) concerns the value for the maximum loop momentum $p_{\text{max}}$. Our criterion to choose $p_{\text{max}}$ is based on the fact that Eq. (4) is only strictly valid if $p\ll k_{1},k_{2}$, as well as $p<k_{\text{NL}}$, otherwise, the blobs in Eq. (4) would not correspond to response-type interactions. In this paper, we therefore choose the cutoff of the momentum integral to be

[TABLE]

Here, we employ the same fraction $f_{\text{sq}}=0.5$ as used in our stitched tree-level result (cf. Eq. 3.4), which assumes that the departure of the response prediction from the full 1-loop trispectrum away from the soft loop momenta limit scales similarly as in the tree-level case (cf. Fig. 9). This is a reasonable assumption as the relevant interactions at both tree- and 1-loop levels are controlled by $\mathcal{R}_{2}$.

4.1 Additional contributions to the 1-loop covariance

The 1-loop contributions to the covariance that are not captured by Eq. (4) are of two types. One, in the limit of soft loop momentum, is the contribution from the six diagrams that are of the form of Eq. (4.1). The other is the contribution from all 1-loop diagrams with loop momentum $p>p_{\text{max}}$. We comment on both contributions in turn below. We will conclude that the non-response contributions are small compared to the response contributions everywhere except for $k_{1},k_{2}\sim 0.1-0.3\,h\,{\rm Mpc}^{-1}$. These missing terms can nonetheless be included with a stitching procedure analogous to that performed for the tree-level covariance in Sec. 3.

Let us first consider Eq. (4.1). The relative size of these non-response terms compared to the response-type ones (cf. Eq. (4.2)) can be roughly estimated by

[TABLE]

where $k_{\text{soft}}\gtrsim 0.1\>h\,{\rm Mpc}^{-1}$. In the numerator, it makes sense to use the power spectrum evaluated at $k_{\text{soft}}$ because $P_{\rm L}(k_{\text{soft}})>P_{\rm L}(k_{\text{hard}})$ in the regime of interest, so that Eq. (4.7) captures the most relevant terms (we are setting the perturbation theory kernels to unity for this estimate). For our choice of $p_{\text{max}}$, we have for the above ratio

[TABLE]

respectively. This indicates that, at the transition from the linear to the nonlinear regime in soft external momenta, $k_{\rm soft}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$, there are 1-loop terms that are sizeable, but that are not of the type of Eq. (4). When both $k_{1},k_{2}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$ these missing terms can be calculated with standard perturbation theory; if, on the other hand, $k_{\rm soft}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$ but $k_{\rm hard}\gg k_{\rm soft}$, then the missing terms can be evaluated by combining standard perturbation theory and response vertices in the same diagram (see e.g. Eq. (2.12) of Ref. [34]). Note that for values of $k_{\text{soft}}\lesssim 0.1\>h\,{\rm Mpc}^{-1}$, the contribution from the 1-loop term is small, and as a result, it is numerically irrelevant whether the 1-loop contribution is accurate.

We now turn to the second missing 1-loop part, namely the contribution from loop momenta with $p>p_{\text{max}}$. Consider a loop momentum $p\gg k_{1},k_{2}$. This corresponds to mode-coupling interactions in which hard ingoing momenta combine to form outgoing soft momenta. These types of couplings are suppressed by momentum and mass conservation (see Appendix B of Ref. [58] for a more detailed discussion). Specifically, the perturbation theory kernels in this limit scale as $(k_{i}/p)^{2}$ ($i=1,2$), and as a result, the loop integrals in this regime contribute negligibly to the total covariance. This, combined with the fact that the response-type terms dominate for $k_{\text{soft}}\gtrsim k_{\text{NL}}$, restricts the contributions from loop momenta $p>p_{\text{max}}$ to the regime of $k_{1},k_{2}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$, as well.

We shall return to the importance of these missing contributions below, as we analyze the results of our calculation. We stress that the inadequacy of the response-based approach to correctly describe the 1-loop covariance for $k_{1},k_{2}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$ can be circumvented by a stitching to the standard 1-loop calculation (see Ref. [30] for the complete expressions), similar to that implemented in the last section for the tree-level covariance. We leave such a stitching at the 1-loop level (as well as the inclusion of other terms important for $k_{\rm soft}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$ but $k_{\rm hard}\gg k_{\rm soft}$) for future work.

4.2 Estimate of higher-loop contributions

We have argued above that, for sufficiently high $k$, the 1-loop contribution to the covariance is dominated by response-type terms. This does not address, however, the issue of the relevance of higher loops on these scales, which we consider now.

Let us consider the 2-loop contribution to the covariance in the response approach. This corresponds to a single diagram that is the $n=3$ generalization of Eq. (4):

[TABLE]

where the dots in the arguments of $\mathcal{R}_{3}$ represent all the angles involved (omitted for brevity) and the factor of 6 accounts for the permutations of the internal loop momenta. Further, $\int_{\bm{p}}\equiv\int{\rm d}^{3}p/(2\pi)^{3}$. An order-of-magnitude estimate of the relative size of this 2-loop contribution to that of Eq. (4) can be written as

[TABLE]

where

[TABLE]

is the variance of the density field up to the cutoff $p_{\text{max}}$ employed in the loop integrals, and $\langle\rangle_{\hat{\bm{p}}_{i}}$ denotes the angle-average over the responses in Eqs. (4) and (4.9). As a very rough estimate, we now assume that only the isotropic response coefficients $R_{n}/n!$ remain after these angle averages. Note, for instance, that comparing Eqs. (4) and (4) shows that for $\mathcal{R}_{2}$ this is not really correct. Keeping this caveat in mind, we obtain

[TABLE]

By continuing this reasoning to higher loops one obtains

[TABLE]

This order-of-magnitude estimate leaves open the possibility that higher-loop terms in the response approach contribute non-negligibly to the total covariance if $\sigma_{p_{\text{max}}}\gtrsim 1$, which corresponds to $p_{\text{max}}\gtrsim k_{\text{NL}}$. The importance of higher-loop response terms and how their contribution scales with $n$ is dependent also on the details of the shape of the $\mathcal{R}_{n}$, or more accurately, on the specific angle-averages that characterize the corresponding diagrams. These higher-order response functions have however never been fully derived, which prevents us from drawing decisive conclusions here. Interestingly, Ref. [39] found that the Eulerian isotropic response coefficients $R_{n}^{E}(k)$, which measure the response of the power spectrum to evolved isotropic modes, are rapidly suppressed numerically at higher orders $n\geq 2$.

We will return to the potential importance of higher-loop terms as we analyze the results of our covariance calculations below. We stress that, for given $k_{1},k_{2}$, the importance of higher-loop terms should progressively decrease in order to render the response approach to the covariance well-defined and predictive. This highly relevant open issue is left for future investigation.

5 Comparison with simulations: angle-averaged case

We now assess the performance of the matter covariance expressions developed in the previous sections by comparing them to estimates from N-body simulations. In particular, we use the results of Ref. [45], who estimated the covariance matrix of the matter power spectrum by cross-correlating the angle-averaged power spectra from more than 12000 simulation boxes with volume $V=[656.25\ h^{-1}{\rm Mpc}]^{3}$. In this section, we therefore consider only the monopole (angle-averaged) part of our covariance expressions ($\ell=0$ in Eq. (2.2)). In Ref. [45], the authors presented results from two sets of simulations: one called Set A, which consists of $12288$ realizations with $N_{p}=256^{3}$ matter tracer particles; and one called Set B, which is made up of a lower number of realizations, $96$, but at higher resolution $N_{p}=1024^{3}$. Apart from providing an independent estimate of the covariance, the diagonal components of the covariance estimated in Set B are used to derive a correction for the power spectrum as well as the covariance measured from Set A for mass resolution effects (see Ref. [45] for details). In this paper, we show the covariance matrices of Ref. [45] estimated from the spectra of Set B and the spectra of Set A after this correction is applied. Note that the diagonal elements of the covariance of both sets thus agree by definition, and that the correction applied cancels out when considering the correlation coefficient Eq. (5.1).

The cosmological parameters of the simulations of Ref. [45] (cf. end of Sec. 1) are almost the same as those used in Ref. [39] to measure the isotropic response coefficients $R_{1}(k)$ and $R_{2}(k)$, which is why we choose these data to compare our results with. We note that our formalism holds generically for any quintessence-type cosmology, provided the corresponding power spectrum responses are known. Other recent estimates of the covariance matrix using simulations include those of Ref. [41], who account for the super-sample covariance term, as well as those in Ref. [29], which were obtained using over 15000 simulations (see also Refs. [25, 26, 27, 28]).

Before discussing the detailed comparison, we make some cautionary remarks regarding simulation measurements of the power spectrum covariance. As any measurement, they in general have statistical and systematic errors. The statistical errors are due to the finite number of realizations, or total volume, of the simulations. On large scales (small wavenumber) the statistical error on the simulation measurements is dominated by the limited number of modes sampled. For a total simulated volume $V_{t}$, this number is given by $k^{3}V_{t}/(2\pi)^{3}$, and hence it is smallest (largest statistical error) for modes close to the fundamental mode of the individual boxes $k_{\text{fund}}=2\pi/L_{\text{box}}$ (where $L_{\text{box}}$ is the box size). On nonlinear scales and thus higher wavenumber, these sample variance effects become smaller, but the precise error becomes harder to quantify because of mode coupling that effectively correlates the statistical error of the covariance across different wavenumbers. The systematic errors of N-body simulations include the finite resolution due to the number of particles, the subtraction of particle shot noise, and transients from the initial conditions. The first two contributions are expected to be most significant on the smallest scales. Quantifying the systematic error on the estimated power spectrum covariance without numerical convergence tests is very difficult. Further, by definition, simulation-based estimates of the covariance matrix lack the contribution from modes with $k<k_{\text{fund}}=2\pi/L_{\text{box}}$, which can also be seen as systematic error. Note however that, as shown in Ref. [41], these can be included at leading order through the super-sample-variance contribution.

Strictly, for comparison with these simulations, one should also include a minimum value $p_{\text{min}}=k_{\text{fund}}$ in the loop integrals of our calculation. However, for $L_{\text{box}}\sim 650\,h^{-1}\,{\rm Mpc}$, we have found that this makes an entirely negligible numerical difference.

Due to the difficulty of obtaining reliable error estimates on the simulation-based covariance, we will mainly discuss the comparison of our covariance prediction to simulation results in the context of the known deficiencies of the former. As discussed in Secs. 4.1 and 4.2, these are non-response-type terms on quasi-linear scales, and higher-loop response terms on fully nonlinear scales. We deliberately avoid quantifying the exact level of agreement between theory and simulations, since it could be misleading given the above mentioned difficulties in estimating the error on the latter.

5.1 Comparison at $z=0$

The color plots in Figure 2 show the correlation coefficient of the angle-averaged matter power spectrum evaluated at two different wavenumbers $k_{1},k_{2}$, which is defined as

[TABLE]

$r_{\ell}(k_{1},k_{2})$ can take on values between $-1$ and $1$. The upper panels show the contribution from the stitched tree-level non-Gaussian term (cf. Eq. (3.4); upper left) and response 1-loop result (cf. Eq. (4); upper right). The lower left panel shows the total prediction for the angle-averaged covariance, which includes the Gaussian diagonal contribution as well. In the upper panels and in the lower left panel, we use the total covariance in the denominator of $r_{\ell=0}(k_{1},k_{2})$, i.e., the lower left panel is obtained by summing the two upper panels and adding the Gaussian contribution. Figure 3 shows instead a few representative slices at constant $k_{2}$ of the covariance matrix $\operatorname{Cov}^{\ell=0}(k_{1},k_{2})$ (not the correlation coefficient).

The upper panels of Fig. 2 and the panels in Fig. 3 are pedagogical in that they illustrate the kinematical regimes in which the tree level and the 1-loop terms contribute most. In particular, the tree-level result dominates when at least one of the modes is $\lesssim 0.1\>h\,{\rm Mpc}^{-1}$. On the other hand, when both modes are $\gtrsim 0.1\>h\,{\rm Mpc}^{-1}$, then most of the contribution comes from the 1-loop term (recall that we do not include the non-response-type loop contribution). The various panels of Fig. 3 help to visualize the gradual increase in importance of the 1-loop term as $k_{2}$ becomes larger. For instance, in the upper left panel for $k_{2}=0.043\>h\,{\rm Mpc}^{-1}$, the 1-loop contribution is fairly small and almost all of the non-diagonal covariance is captured at tree level (blue line). As $k_{2}$ increases however (left to right, top to bottom), the tree-level result becomes progressively smaller at high $k_{1}$, and is complemented by the growing contribution of the 1-loop term (green line). In light of the relative importance of tree-level and 1-loop contributions, the sharp discontinuities at high $k$ between the two branches of the stitched tree-level result of Eq. (3.4), as well as the extrapolation of the tree-level response to the case of $k_{\text{soft}}\gtrsim k_{\text{NL}}$ are not affecting the total covariance in this regime because the entire stitched tree-level contribution is a small part of the total result.

The lower right panel of Fig. 2 shows $r_{\ell=0}(k_{1},k_{2})$ from the simulation Set A of Ref. [45]. The visual comparison to our prediction does not reveal strong differences in either shape or overall amplitude. A more detailed comparison with simulations is shown in Fig. 4, where we show our total prediction along with the simulation Set A and Set B results of Ref. [45]. Up to the approximation employed in the extrapolation of the anisotropic response coefficients $R_{O}(k)$ (cf. Appendix B), our calculation is guaranteed to capture the total covariance if the soft mode is sufficiently linear, $k_{\text{soft}}\ll k_{\text{NL}}\approx 0.3\>h\,{\rm Mpc}^{-1}$. An interesting application of our calculation in this regime is therefore to test simulation-based estimates of the covariance matrix for systematic errors. Indeed, both simulation sets are in relatively good agreement with our calculation whenever $k_{\text{soft}}\ll k_{\text{NL}}$; this includes roughly the whole $k_{1}$ range in the upper three panels, as well as the low-$k_{1}$ parts of the lower six panels in Fig. 4. The differences between Set A and Set B are likely to be mostly caused by the larger statistical uncertainties in Set B due to the smaller volume covered. Further, the departures seen in the simulations for large-scale modes, $k_{1}\lesssim 0.03\>h\,{\rm Mpc}^{-1}$ (noticeable in the lower three panels of Fig. 4), are likely to be due to insufficient sampling of these modes by both sets of simulations, as we noted already in the beginning of this section.

The lower six panels of Fig. 4 correspond to slices with $k_{2}>0.1\>h\,{\rm Mpc}^{-1}$, in which one notes that our calculation falls short of describing completely the simulation measurements for $k_{1}\gtrsim 0.1\>h\,{\rm Mpc}^{-1}$. As we have discussed in Secs. 4.1 and 4.2, there are two types of terms that are expected to contribute non-negligibly in this regime, presumably accounting for a large fraction of the observed difference between theory and simulations. One corresponds to 1-loop diagrams that cannot be brought into the form of Eq. (4) and that can contribute sizeably when $k_{1},k_{2}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$. The inclusion of these terms via a stitching of standard- and response-based expressions should render the whole calculation fully predictive in this regime. In regimes in which $k_{\rm soft}\sim 0.1-0.3\>h\,{\rm Mpc}^{-1}$ and $k_{\rm hard}>k_{\text{NL}}$ (e.g. $k_{1}=0.2\>h\,{\rm Mpc}^{-1}$ in the lower right panel of Fig. 4), the missing terms can also be added, and in fact, using responses to describe the interactions that involve $k_{\rm hard}$. These terms, however, are not expected to play a major role in cases when $k_{1},k_{2}>0.3\>h\,{\rm Mpc}^{-1}$. In this kinematic regime on the other hand, one expects that 2- and higher-loop terms, which themselves are dominated by response-type terms in this regime (cf. Eq. (4.9)), can account for the missing contribution. An important requirement for the response approach to remain predictive when $k_{1},k_{2}>0.3\>h\,{\rm Mpc}^{-1}$ is, therefore, that higher-loop response-type contributions become progressively smaller. For the time being, we cannot provide a conclusive answer on the exact relative size of higher-loop response terms, and defer that to future work.

It is also instructive to compare predictions for the diagonal of the covariance matrix ($k_{1}=k_{2}$), which is shown in Fig. 5 (solid lines for $z=0$). On large scales, this is dominated by the Gaussian contribution, and the tree-level non-Gaussian contribution is subdominant at all $k$ values. The 1-loop contribution only starts to become important for $k\gtrsim 0.3\>h\,{\rm Mpc}^{-1}$. As a result, for $k\lesssim 0.1\>h\,{\rm Mpc}^{-1}$, there is good agreement between our calculation and the simulations (up to noise), but this is unsurprising because here the result is set by the trivial and well understood Gaussian contribution. For $k\gtrsim 0.3\>h\,{\rm Mpc}^{-1}$, the simulation results at $z=0$ (solid lines) have a higher amplitude than our calculation (better discernible in the less noisy Set A), but, as already mentioned above, this is a regime in which higher-loop terms are expected to contribute non-negligibly, and hence potentially reduce the gap between theory and simulations.

5.2 Comparison at $z=2$

Our prediction for the covariance can be straightforwardly applied to other redshifts as well. To do so, one should use the response coefficients $R_{O}(k)$ at the desired redshift. This includes 1) using the measured isotropic responses from simulations and adjusting appropriately the nonlinear extrapolation of the anisotropic ones (cf. Appendix B); 2) evaluate all spectra at the desired redshift; 3) and using the corresponding value of the nonlinear scale $k_{\text{NL}}$ in setting $p_{\text{max}}$, which increases with redshift.

The color plots in Fig. 6 show the correlation coefficient measured from both simulation sets of Ref. [45] and that of our calculation at $z=2$, as labeled. The correlation coefficients measured from the two simulation sets are noticeably different at this redshift. In Ref. [45], this is attributed to the fact that the total volume of simulations in Set B (96 realizations) is not sufficient to appropriately sample the covariance. Note however that the results from the two simulation sets are in much better agreement at $z=0$, as we have seen above. Interestingly, our prediction agrees markedly better with the result from Set B. This becomes clearer from the lower right panel of Fig. 6, which shows the slice of the covariance matrix at constant $k_{2}=1\>h\,{\rm Mpc}^{-1}$ (the same as the lower right panel of Fig. 3, but for $z=2$). Our calculation underpredicts both simulation set results, but the level of disagreement between our model and Set A is significantly larger than that with Set B. In fact, at $z=2$ the performance of our calculation in reproducing the results from Set B is comparable to the performance of the same in reproducing the results from both Set A and Set B at $z=0$. The lower volume of the higher-resolution simulation Set B unfortunately prevents us from drawing robust conclusions on the significance of its better agreement (compared to Set A) with our theoretical prediction.

This picture becomes different if one focuses only on the diagonal of the covariance matrix at $z=2$, which is shown by the dashed lines in Fig. 5. As mentioned at the beginning of this section, the correction of Set A obtained by matching the diagonal elements to Set B ensures that the two agree on the diagonal, within the noise of the smaller Set B. They both agree well with our prediction. Note that this level of agreement for $k\gtrsim 0.6\>h\,{\rm Mpc}^{-1}$ depends quite crucially on the contribution from the 1-loop term.

6 Angular dependence of the matter power spectrum covariance

We now go beyond the case of the angle-averaged non-Gaussian covariance and analyze its angular dependence, which can be organized into Legendre multipoles according to Eq. (2.2). For the case of the response-based contributions, the well defined analytical dependence on $\mu_{12}$ allows for all multipoles to be evaluated analytically (cf. Appendix C). The angular dependence of the standard perturbation theory contributions (the tree-level one in our case) is more cumbersome and we perform the angle-averages numerically.

Figure 7 displays a few predictions at $z=0$ for the quadrupole ($\ell=2$) and hexadecupole ($\ell=4$), as well as the monopole ($\ell=0$) case studied more extensively in the previous section. The left panel shows the diagonal of these three multipoles. One notes that there is a hierarchy between the diagonal of these terms, with the monopole being the largest and the hexadecupole the smallest. The case depicted by the dashed black line corresponds to the Gaussian result of Eq. (2.9). All multipoles of the Gaussian contribution have this form. On scales $k\lesssim 0.3\>h\,{\rm Mpc}^{-1}$, all multipoles match the Gaussian result because on these scales the fractional size of the non-Gaussian contribution to the diagonal is negligible.

The left two panels of Fig. 7 display slices at constant $k_{2}$ of the multipoles, which show that the hierarchy displayed along the diagonal does not necessarily hold for off-diagonal terms of the covariance. We note also that in the middle panel, the transition between the squeezed and standard tree-level results for $\ell=2$ and $\ell=4$ exhibits a much sharper discontinuity compared to the case for $\ell=0$. This may suggest that our choice of $f_{\text{sq}}=0.5$ may have to be revisited in more careful investigations of the angular dependence of the covariance using our stitched results. We note also that the higher multipoles of the covariance $(\ell>0)$ depend to a much higher degree on the anisotropic response coefficients $R_{O}(k)$, for which we currently only have extrapolations based on the physical reasoning described in Ref. [34].

The higher multipoles of the covariance matrix are much less studied in the literature than the monopole, since most investigations only consider the covariance of the angle-averaged matter power spectrum. One exception is Ref. [27], who measure the angular dependence of the covariance using N-body simulations (see also Ref. [57] for an earlier perturbation theory calculation). Here, we do not attempt to perform detailed comparisons against these simulation results, given the lack of simulation-calibrated anisotropic response coefficients for their cosmology, and defer that to future work. Nevertheless, we point out that the prediction depicted in the left panel of Fig. 7 is in agreement with the hierarchy of the diagonal covariance elements shown in Fig. 10 of Ref. [27].

7 Summary and Discussion

We have described a calculation of the matter power spectrum covariance $\operatorname{Cov}(\bm{k}_{1},\bm{k}_{2})$ based on perturbation theory augmented with specific resummed interaction vertices (the responses), which is applicable in all regimes of structure formation. More specifically, we describe the non-Gaussian part of the matter power spectrum covariance which is equivalent to the parallelogram configuration of the matter trispectrum, $T_{m}(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{2},-\bm{k}_{2})$ (cf. Eq. (2.5)). There are two other important contributions to the total matter covariance, namely the Gaussian diagonal term and the (also non-Gaussian) super-sample contribution, but these are both well understood. Our calculation is built upon the work of Ref. [34], in which the authors have illustrated how the calculation of certain mode-coupling interactions in perturbation theory can be made accurate beyond the perturbative regime with the aid of power spectrum responses.

The $n$-th order power spectrum responses $\mathcal{R}_{n}$ describe the coupling of $n$ long-wavelength modes with the local nonlinear matter power spectrum (cf. Eq. (2.1)), and these responses can be measured accurately with separate universe simulations. The crucial and novel steps of our calculation consist essentially in the identification of the mode-coupling terms in the non-Gaussian covariance that can be described as power spectrum responses (or more technically, that can be resummed to all orders in perturbation theory using responses), thereby enabling efficient and accurate evaluation of these terms in kinematical regimes in which standard perturbation theory breaks down. The well-defined angular structure of the $\mathcal{R}_{n}$ also permits us to straightforwardly determine the angular dependence of the covariance matrix (cf. Sec. 6). Although the formalism presented here still needs as ingredients response measurements from simulations, we stress that the number of simulations that need to be performed for these measurements are orders of magnitude fewer than those that are needed for the direct, fully simulation-based estimation of the power spectrum covariance.

In this paper, we have worked explicitly at tree- and 1-loop-levels in the standard perturbation theory expansion of the non-Gaussian covariance (cf. Eq. (2.10)). At tree level (cf. Sec. 3), we have presented a way to stitch together response-based terms with terms from standard perturbation theory. At the 1-loop level (cf. Sec. 4), we have seen that a response approach is particularly useful because a significant part of the contribution comes from the coupling of soft loop to hard external momenta, $p\ll k_{1},k_{2}$, which are precisely interactions that power spectrum responses are able to capture. We have also pointed out, however, that our response-based 1-loop calculation still leaves important contributions to the covariance uncovered, but which can be added after some additional development (cf. Secs. 4.1 and 4.2).

In order to organize the discussion about which parts of the covariance are already captured by our description, and which require additional work, we can divide the parameter space of $\operatorname{Cov}^{\rm NG,\ell=0}(k_{1},k_{2})$ into 5 kinematic regimes, as illustrated in Figure 8. This division naturally arises when distinguishing three regimes of wavenumber for $k_{1}$ and $k_{2}$: the linear regime ($k_{i}\ll k_{\text{NL}}$), the quasilinear regime $k_{i}\lesssim k_{\text{NL}}$, and the fully nonlinear regime $k\gtrsim k_{\text{NL}}$. The right panels in Fig. 8 show $\operatorname{Cov}^{\rm NG,\ell=0}(k_{1},k_{2})$ as a function of $k_{1}$, while keeping $k_{2}/k_{1}$ constant, i.e., fixed level of squeezing. We now briefly discuss each of these regimes, denoting as throughout $k_{\text{soft}}\equiv\min\{k_{1},k_{2}\}$ and $k_{\text{hard}}\equiv\max\{k_{1},k_{2}\}$.

$\bullet$ I: $k_{1}\ll k_{\text{NL}},\ k_{2}\ll k_{\text{NL}}$. In this regime, with both modes in the linear regime, the covariance is captured completely by the standard tree-level result (cf. Eq. (3)).

$\bullet$ II: $k_{\text{soft}}\ll k_{\text{NL}},\ k_{\text{soft}}\lesssim f_{\text{sq}}k_{\text{hard}}$, for any $k_{\text{hard}}$. In this squeezed regime, with the soft mode being in the linear regime, $\operatorname{Cov}^{\rm NG,\ell=0}(k_{1},k_{2})$ is exactly captured by the second-order power-spectrum response $\mathcal{R}_{2}$ (cf. Eqs. (3.3)), up to corrections that scale as $\mathcal{O}((k_{\text{soft}}/k_{\text{hard}})^{2})$, $\mathcal{O}((k_{\text{soft}}/k_{\text{NL}})^{2})$. Note that the result remains valid for any value of the hard momenta, including $k_{\text{hard}}>k_{\text{NL}}$.

$\bullet$ III: $k_{\text{soft}}\lesssim k_{\text{NL}},\ k_{\text{soft}}\lesssim f_{\text{sq}}k_{\text{hard}}$, for any $k_{\text{hard}}$. This regime is still squeezed, but with quasilinear values of $k_{\text{soft}}$. The response tree-level result that fully determines regime II still contributes, but now loop terms are no longer negligible and become increasingly important with increasing $k_{\text{soft}}$. While in this paper, we have included only the single response-type contribution that is present for generic configurations (cf. Eq. (4)), we expect that all loop contributions in the squeezed regime can be captured by responses, in the sense that non-response-type contributions are suppressed by $(k_{\text{soft}}/k_{\text{hard}})^{2}$.

$\bullet$ IV: $k_{\text{soft}}\sim k_{\text{hard}}\lesssim k_{\text{NL}}$. In this non-squeezed, quasi-linear regime, the standard tree level term still contributes non-negligibly and the 1-loop contributions are important. In this regime, we expect the non-response-type 1-loop terms to be relevant as well. Higher-order loop terms also become increasingly relevant as $k_{1}$ and $k_{2}$ approach $k_{\text{NL}}$.

$\bullet$ V: $k_{\text{soft}}\sim k_{\text{hard}}\gtrsim k_{\text{NL}}$. In this regime, the tree-level and non-response-type 1-loop contributions are negligible and the result is dominated by response-type loop terms. In this paper, we have worked explicitly at 1-loop order, but higher-loop terms (cf. Eq. (4.9) for two-loop) are expected to be significant as well. Crucially, we note that if higher loop contributions are not progressively suppressed, then the approach presented here will not be predictive in this regime.

The discussion points above motivate two immediate steps that can be taken to improve our prediction for $\operatorname{Cov}^{\rm NG}$. One is the inclusion of 1-loop terms that cannot be described with power spectrum responses and which are expected to be important in regime IV. These terms have already been derived and calculated in Ref. [30]. One can therefore include them into our model by following a “stitching” recipe similar to that employed at tree level in this paper. The other improvement is the inclusion of higher-loop response-type contributions, which is expected to result in relevant contributions to regimes III, IV,V. This calculation is also crucial to establish the theoretical consistency of the approach presented here: our prediction on fully nonlinear scales is only robust if the higher-loop contributions can be shown to be progressively suppressed compared to the leading 1-loop contribution derived here. This could happen if the relevant angle-averages of higher-order responses are suppressed.

In regimes I-II, on the other hand, the current calculation captures already the total leading contribution to the covariance. An interesting consequence of this is that comparisons to our calculation in this regime can therefore serve as useful validation checks of simulation-based estimates of the covariance.

In summary, we have paved the way towards the development of a physically motivated framework that enables an efficient calculation of the matter power spectrum covariance without any adjustable free parameters, and which is valid deeply into the nonlinear regime of structure formation where standard perturbative schemes break down. Our calculation can also be generalized to describe matter correlations at different redshift values (by making use of unequal-time power spectrum responses [34]), which is useful for tomographic cosmic shear analyses. Compared to standard ways to estimate the covariance matrix with N-body simulations, an approach combining simulations with analytical results such as the one put forward here has the enormous advantage of requiring far less computational resources. A straightforward consequence of this is that robust and systematic studies of the dependence on cosmology of the covariance can be performed through this approach. The same can be said about the impact of baryonic effects on the covariance matrix. These constitute pieces of information that are very relevant for upcoming observational surveys such as Euclid [22] and LSST [23] whose statistical precision will be at a level sufficient to make these systematic effects on the covariance a pressing concern.

Acknowledgments

We thank Linda Blot for providing the numerical measurements of the power spectrum covariance, and Jean-Michel Alimi, Linda Blot, Pier-Stefano Corasaniti, Joachim Harnois-Déraps, Wayne Hu, Irshad Mohammed, Yann Rasera, Vincent Reverdy, Uro$\check{\rm s}$ Seljak and Zvonimir Vlah 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

In this appendix, we list the Feynman rules we employ to calculate $n$-point functions in cosmological perturbation theory. Our conventions are based on the one used by Ref. [58]. The rules are as follows:

An $n$-point correlation function is represented by a set 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}$, and each such vertex is assigned a factor

[TABLE]

We assign a negative (positive) sign to ingoing (outgoing) momenta. Each ingoing line has to be directly connected to a propagator (linear power spectrum). 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 called loop momenta and are integrated over as

[TABLE]

A tree-level diagram is a diagram without any such loop integrals. 5. 5.

Each diagram is multiplied by the symmetry factor, which accounts for degenerate configurations of the diagram, as well as all nonequivalent labellings of external lines.

To include response-type interactions, we consider one additional rule:

Response-type interaction vertices 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 predictive in the limit where $\sum_{a}p_{a}\ll\min\{k,k_{\text{NL}}\}$, but no restriction is placed on the magnitude of the outgoing momenta, which can be in the nonlinear regime. In our notation, we represent them as dashed blobs. Each such vertex is assigned a factor (cf. Eq. (2.1))

[TABLE]

where the dots in the argument of $\mathcal{R}_{n}$ denote all the relevant cosine angles and soft momenta magnitude ratios that exist at a given order $n$ (described in detail in Sec. 2 of Ref. [34]). The factor $1/2$ serves to cancel the trivial permutation $\bm{k}\leftrightarrow\bm{k}^{\prime}$, which is always present in the response vertex.

Appendix B Response expressions

In this appendix, we show the decomposition of the second-order response $\mathcal{R}_{2}$ as derived in Ref. [34], including the expression of the response coefficients that we use in our covariance calculations.

For general kinematic configurations $\mu_{1}$, $\mu_{2}$, $\mu_{12}$ and $f_{12}=p_{1}/p_{2}$, the second-order matter power spectrum response is given by (in the Eulerian decomposition; see also Ref. [34] for the Lagrangian decomposition and Ref. [49] for another example of the angular decomposition)

[TABLE]

Setting $\mu_{1}=\mu$, $\mu_{2}=-\mu$, $\mu_{12}=-1$ and $f_{12}=1$ yields Eq. (2.1), which is the configuration of $\mathcal{R}_{2}$ which enters the covariance calculations. The nonlinear extrapolation recipe of the $R_{O}(k)$ presented in Ref. [34] results in the following expressions:

[TABLE]

where $G_{n}^{\rm nl}(k)$ correspond to the so called growth-only response functions measured in Ref. [39]. The nonlinear shapes of $R_{1}(k)$ and $R_{2}(k)$ are also those measured directly from the separate universe simulations of Ref. [39]. Reference [39] found that the different isotropic response coefficients ($R_{1}$, $R_{2}$ and $R_{3}$) are similar in shape. In the context of the halo model, this can be explained by the characteristic transition between the 2-halo and 1-halo regimes, which leads to a peak of the response coefficients around $k\sim 0.5\>h\,{\rm Mpc}^{-1}$, with a suppression at higher $k$. The derivation of the above expressions for the $R_{O}(k)$ assumes that these halo model-based considerations extend also to the case of the anisotropic response coefficients.

In the exercise displayed in Fig. 9 below, we use the tree-level expressions of the $R_{O}(k)$. Starting from Eq. (B), these can be obtained with the substitutions (see Refs. [39, 34] for more details) $P_{m}(k)\rightarrow P_{\rm L}(k)$, $G_{1}^{\rm nl}\rightarrow G_{1}^{\rm tree}=26/21$, $G_{2}^{\rm nl}\rightarrow G_{2}^{\rm tree}=3002/1323$ and

[TABLE]

Appendix C Multipoles of the response tree-level and 1-loop covariance

In this appendix, we collect all the Legendre multipoles of the tree-level $\operatorname{Cov}_{\mathcal{R}\text{-tree}}^{\rm NG,\ell}(k_{1},k_{2})$ and 1-loop level $\operatorname{Cov}_{\mathcal{R}\text{-1loop}}^{\rm NG,\ell}(k_{1},k_{2})$ covariance in the response approach, which can be given analytically.

At tree-level, using Eqs. (2.1), (3.3) and (2.2), we can write all non-vanishing multipoles ($\ell=0,2,4$) as

[TABLE]

where, as in the main text, $k_{\text{soft}}={\rm min}\{k_{1},k_{2}\}$ and $k_{\text{hard}}={\rm max}\{k_{1},k_{2}\}$.

Similarly, but for the 1-loop expressions, we can combine Eqs. (2.1), (4) and (2.2) to find

[TABLE]

while all other multipoles vanish exactly (here, $\mathcal{A}_{i}\equiv\mathcal{A}(k_{i})\ (i=1,2)$, and similarly for $\mathcal{B}_{i}$ and $\mathcal{C}_{i}$).

Appendix D Determination of $f_{\text{sq}}$

In this appendix, we describe the procedure used to choose the value of $f_{\text{sq}}$, which controls the transition from the squeezed to the non-squeezed branches in our stitched tree-level result of Eq. (3.4).

The upper triangle of the color plot in Fig. 9 shows the angle-averaged standard tree-level covariance $\operatorname{Cov}_{\text{SPT-tree}}^{\rm NG,\ell=0}(k_{1},k_{2})$ ($\ell=0$ in Eq. (2.2)) using the upper branch of Eq. (3.4). The lower triangle shows the response tree-level result $\operatorname{Cov}_{\mathcal{R}\text{-tree}}^{\rm NG,\ell=0}(k_{1},k_{2})$, but using the tree-level limit of the second-order response $\mathcal{R}_{2}^{\rm tree}$ (cf. Appendix B). By construction, these two results are the same in the limit in which one of the modes is much smaller than the other. A reasonable choice of $f_{\text{sq}}$ should therefore be one that ensures a sufficiently smooth transition between the two branches of Eq. (3.4) at $k_{\text{soft}}=f_{\text{sq}}k_{\text{hard}}$. The upper right panel of Fig. 9 shows the relative difference between the two non-Gaussian tree level expressions, plotted as a function of $k_{\text{soft}}/k_{\text{hard}}$. The several lines shown (grey) correspond to different values of $k_{\text{hard}}$. To guide the eye, the colored curves show the relative difference between the two covariance results along the dashed lines overlaid in the color plot. We consider the two results to be sufficiently close to one another when the relative difference drops below the $25\%$ mark (dashed lines in the upper right panel). The distribution of the values of $f_{\text{sq}}^{25\%}=k_{\text{soft}}/k_{\text{hard}}$ when this happens is shown in the lower right panel of Fig. 9. The mean and median of the distribution are both $\approx 0.5$, and for that reason we choose $f_{\text{sq}}=0.5$. The choice of the $25\%$ figure is somewhat arbitrary, but small variations around it do not affect our tree level results significantly.

Appendix E Explicit relation between the product of two $\mathcal{R}_{2}$ responses and 1-loop covariance terms

In this appendix, we show explicitly which of the 1-loop covariance terms are captured by the response diagram of Eq. (4). In order to compare with the 1-loop expressions in standard perturbation theory, we have to insert the tree-level expression for the response. We then obtain

[TABLE]

where we have used

[TABLE]

which holds in the limit $p\ll k$ as derived in Ref. [34] using the squeezed-limit tree-level trispectrum.444 From the considerations in Sec. 4.2 of Ref. [34], we can write

\displaystyle 2\mathcal{R}_{2}^{\rm tree}(k,\mu_{1},\mu_{2},\mu_{12},p_{1}/p_{2})P_{\rm L}(k^{\prime})P_{\rm L}(p_{1})P_{\rm L}(p_{2})=\lim_{p_{1},p_{2}\to 0}\Big{[}6F_{3}(\bm{k},\bm{p}_{1},\bm{p}_{2})P_{\rm L}(k)P_{\rm L}(p_{1})P_{\rm L}(p_{2})

\displaystyle\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +4F_{2}(-\bm{p}_{1},\bm{k}+\bm{p}_{1})F_{2}(\bm{p}_{2},\bm{k}+\bm{p}_{1})P_{\rm L}(k+p_{1})P_{\rm L}(p_{1})P_{\rm L}(p_{2})+(\bm{k}\leftrightarrow\bm{k}^{\prime})\Big{]},

(E.3) where as everywhere in this paper, the limit means retaining terms at leading order in $p_{1}$ and $p_{2}$, and $\bm{k}^{\prime}=-\bm{k}-\bm{p}_{1}-\bm{p}_{2}$. Specifying to the covariance configuration, $\bm{p}_{1}=\bm{p}$, $\bm{p}_{2}=-\bm{p}$ and approximating $P_{\rm L}(k\pm p)\approx P_{\rm L}(k)$, we arrive at Eq. (E.2).

Our goal is to demonstrate that Eq. (E) is indeed obtained by summing the 1-loop trispectrum terms that are $\mathcal{O}([P(p)]^{2})$ in the covariance configuration. This demonstration can be carried out fairly straightforwardly, despite involving tedious counting of momentum permutations. For completeness, we shall nevertheless lay down a few intermediate steps of the derivation below.

The matter trispectrum that contributes to the covariance at 1-loop (which is directly related to the covariance via a volume factor), $T_{m}^{1-{\rm loop}}(\bm{k}_{1},-\bm{k}_{1},\bm{k}_{2},-\bm{k}_{2})\equiv\bar{T}_{m}^{1-{\rm loop}}(\bm{k}_{1},\bm{k}_{2})$, has contributions from nine types of diagrams (we do not list all of them for brevity, but see e.g. Appendix A of Ref. [30]). Out of these nine diagrams, there are three which contain permutations of momenta that give rise to terms that are quadratic in the power spectrum of the loop momentum $[P(p)]^{2}$, with the remaining six diagrams being linear in $P(p)$ (recall the discussion around Eqs. (4.1) and (4.2)):

[TABLE]

in which the split on the right-hand side separates the diagrams by their order in $P(p)$, as indicated in the subscripts. From hereon in this appendix, we also drop the subscript m to ease the notation. The response-based definition of Eq. (E) describes only the terms that contribute to $T^{1-{\rm loop}}_{[P(p)]^{2}}(\bm{k}_{1},-\bm{k}_{2},\bm{k}_{2},-\bm{k}_{2})$, which are

[TABLE]

where we have adopted the same notation as in Ref. [30]. Next, we write these three diagrams explicitly and demonstrate that their sum yields exactly Eq. (E).

We start with $T_{2222}$, which consists of three terms given by

[TABLE]

Note that out of all possible 24 permutations of four elements, only three correspond to distinct diagrams. In the covariance configuration ($\bm{k}_{2}=-\bm{k}_{1}$, $\bm{k}_{4}=-\bm{k}_{3}$), only two of the permutations give contributions $\propto[P_{\rm L}(p)]^{2}$ for general $\bm{k}_{1},\bm{k}_{3}$.555For instance, in the covariance configuration, the spectra factors of the permutation written explicitly in Eq. (E.6) become $P_{\rm L}(p)[P_{\rm L}(|\bm{p}+\bm{k}_{1}|)]^{2}P_{\rm L}(|\bm{p}+\bm{k}_{1}-\bm{k}_{3}|)$. Noting that we can transform $\bm{p}+\bm{k}_{1}\to\bm{p}$ because we are integrating over $\bm{p}$, we can write $[P_{\rm L}(p)]^{2}P_{\rm L}(|\bm{p}-\bm{k}_{1}|)P_{\rm L}(|\bm{p}-\bm{k}_{3}|)$. The other permutation can also be $\propto[P_{\rm L}(p)]^{2}$, but only if $\bm{k}_{1}=\bm{k}_{3}$, i.e., if the two wavevectors have the same magnitude and are aligned. After angle-averaging, this term contributes negligibly, so we do not consider it further. We note however that this term can be modeled accurately with $\mathcal{R}_{2}$ as well by modifying Eq. (4) to include a term proportional to $\delta_{D}(\bm{k}_{1}-\bm{k}_{3})$. One can then straightforwardly work out that (relabeling $\bm{k}_{3}\rightarrow\bm{k}_{2}$ and with the limit $p\ll k_{1},k_{2}$ understood)

[TABLE]

which matches exactly the contribution in Eq. (E) that is proportional to the product of four $F_{2}$ kernels. The second line is obtained from the first by using the fact that one can transform the loop momentum as $\bm{p}\rightarrow-\bm{p}$ in the integrand. The subscript $\mathcal{R}_{2}$ indicates that we consider only the mode permutations that are $\propto[P_{\rm L}(p)]^{2}$ for general $\bm{k}_{1},\bm{k}_{2}$.

Similarly, the $T_{3311b}$ term can be written as

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

where in total there are now twelve distinct permutations contributing to $T_{3311b}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3},\bm{k}_{4})$. Setting $\bm{k}_{2}=-\bm{k}_{1}$ and $\bm{k}_{4}=-\bm{k}_{3}$, one finds that only the four permutations of $T_{3311b}^{C}$ are $\propto[P_{\rm L}(p)]^{2}$ for any $\bm{k}_{1}$, $\bm{k}_{3}$. Their contribution is obtained straightforwardly as (relabeling $\bm{k}_{3}\rightarrow\bm{k}_{2}$)

[TABLE]

which matches the corresponding term in Eq. (E). Here, similarly to the case of the $T_{2222}$ term above, we have also skipped writing the contribution from terms that are $\propto[P_{\rm L}(p)]^{2}$ when $\bm{k}_{1}=\bm{k}_{2}$ (the four permutations of $T_{3311b}^{A}$), as well as, $\bm{k}_{1}=-\bm{k}_{2}$ (the four permutations of $T_{3311b}^{B}$).

Finally, the $T_{3221c}$ term has also twelve terms, which we can write explicitly as

[TABLE]

with

[TABLE]

[TABLE]

[TABLE]

After setting $\bm{k}_{2}=-\bm{k}_{1}$, $\bm{k}_{4}=-\bm{k}_{3}$, we find four terms that are $\propto[P_{\rm L}(p)]^{2}$ (which come from $T_{3221c}^{B}$ and $T_{3221c}^{C}$) for general $\bm{k}_{1},\bm{k}_{3}$, which we can write as (relabeling $\bm{k}_{3}\rightarrow\bm{k}_{2}$ and interpreting $p\ll k_{1},k_{2}$)

[TABLE]

which is exactly the remaining term in Eq. (E). Once again, we have skipped explicitly writing cases when $\bm{k}_{1}=\pm\bm{k}_{2}$.

As discussed in the main body of the text, the generalization of Eq. (E) to the nonlinear regime (that is, promoting $P_{\rm L}(k_{1})$, $P_{\rm L}(k_{2})$ and $\mathcal{R}_{2}^{\rm tree}$ to their nonlinear versions) then allows one to capture infinitely many higher-loop contributions via the simulation-calibrated nonlinear responses, while always considering only a single soft loop.

Appendix F Comparison to the covariance matrix model of Ref. [1]

In this appendix, we compare our calculation of the covariance with the model put forward in Ref. [1]. In the latter, the total angle-averaged covariance is given by

[TABLE]

The term $\operatorname{Cov}^{\rm G}(k_{1},k_{2})$ is the Gaussian result given by Eq. (2.9). The tree-level contribution $\operatorname{Cov}_{\mbox{\footnotesize{Ref.~{}\cite[cite]{[\@@bibref{Number}{mohammed1}{}{}]}}}}^{\rm NG,tree}(k_{1},k_{2})$ is given by

[TABLE]

i.e., it is obtained by multiplying the standard perturbation theory tree-level result with a correction factor that is the ratio of nonlinear to linear spectra. Concretely, in the results shown in Ref. [1], this correction is realized by comparing the tree-level covariance divided by the linear spectra with the simulation covariance estimates divided by the nonlinear spectra. While there is no obvious physical justification for this rescaling of the tree-level trispectrum, we will see below that it improves the model significantly when compared to simulation measurements.

The 1-loop term is obtained with the aid of a calculation based on functional derivatives of the small-scale power spectrum with respect to the large-scale one. Explicitly, the end result is given by (see Ref. [1] for details)

[TABLE]

with

[TABLE]

In the limit $p/k\rightarrow 0$ we obtain

[TABLE]

In this $p/k\rightarrow 0$ limit, $V(q,k)$ becomes independent of $p$ and Eq. (F) can be directly compared to Eq. (4). One important difference is that in Eq. (4), while performing the angle part of the loop integral, one averages the product of two $\mathcal{R}_{2}$ responses, whereas in Eq. (F), one has the product of two angle-averaged responses. The term $P_{m}(k_{1})P_{m}(k_{2})/P_{\rm L}(k_{1})/P_{\rm L}(k_{2})$ in Eq. (F) is added as a correction term that improves the accuracy on nonlinear scales, in a way similar to the correction employed in Eq. (F.2) for the standard tree level covariance. Finally, the Lorentzian damping term \big{(}1+(p/p_{\text{nl}})^{2}\big{)}^{-2} in Eq. (F), with $p_{\text{nl}}$ being a free parameter, is introduced to help suppress the excess of amplitude of Eq. (F) on nonlinear scales (compared to simulations, as we shall see below). In our response-based description, a similar damping is included without any adjustable parameters by employing the fully nonlinear responses (see Fig. 1 in Ref. [34]). Thus, one can interpret the rescaling by $(P_{m}/P_{\rm L})$ and the phenomenological damping term in Eq. (F) as a rough model of the physical nonlinear response.

Figure 10 compares our predictions with those of Ref. [1], labeled as $Mohammed\ et\ al$. The two panels show the total angle-averaged covariance as a function of $k_{1}$, for the fixed values of $k_{2}$ indicated above each panel. The red curve displays the result from our calculation as described in the main body of this paper. The model of Ref. [1] is shown for two choices of $p_{\text{nl}}$: $p_{\text{nl}}=\infty$ (green) and $p_{\text{nl}}=0.4\>h\,{\rm Mpc}^{-1}$ (magenta). In the left panel, for $k_{2}=0.1\>h\,{\rm Mpc}^{-1}$, the model of Ref. [1] captures the amplitude of the simulation results better, even though the agreement is not perfect. This level of agreement with the simulation results for $k_{1}\gtrsim 0.1\>h\,{\rm Mpc}^{-1}$ depends sensitively on the correction factor $P_{m}(k_{1})P_{m}(k_{2})/P_{\rm L}(k_{1})/P_{\rm L}(k_{2})$ in Eqs. (F.2) and (F). Concretely, dropping the correction term results in a much lower amplitude of the curves whenever $k_{\text{hard}}\gtrsim k_{\text{NL}}$.

On the other hand, the right panel ($k_{2}\approx 0.5\>h\,{\rm Mpc}^{-1}$) shows that the model of Ref. [1], specifically the contribution in Eq. (F), drastically overpredicts the simulation results on small scales (green) if the damping term is not included. This forces one to tune the value of the parameter $p_{\text{nl}}$ to bring the model closer to the simulations (magenta). As we have seen in the main body of the paper, our calculation (red line) still underpredicts the simulation results on nonlinear scales. However, the well-defined physical grounds of the response approach allows us to identify higher loops (which cannot be captured by a term of the type of Eq. (F)) as a likely reason for this discrepancy.

As an exercise, we show as the blue curve the result of our calculation when allowing for a larger loop-momentum cutoff $p_{\text{max}}$. Specifically, the loop momentum is allowed to be as large as $k_{\text{soft}}$ and larger than $k_{\text{NL}}$, at which point the loop integrand is strictly no longer in the response regime. This can be seen as a phenomenological attempt to capture higher-loop contributions that are currently left out. Although it is interesting to observe that the blue curve is capable to describe the simulation results quite well (for these values of $k_{2}$, at least), we stress that such a calculation is not guaranteed to be self-consistent, and that it should not be used to make predictions based on the response approach. Allowing $p_{\text{max}}$ to be a free fitting parameter risks overfitting the simulation measurements, and loses the predictivity of the covariance for other cosmologies. These are issues that generally apply to any model that contains free fitting parameters.