TY - JOUR
AU1 - Mohammed,, Irshad
AU2 - Seljak,, Uroš
AU3 - Vlah,, Zvonimir
AB - Abstract We evaluate the covariance matrix of the matter power spectrum using perturbation theory up to dominant terms at 1-loop order and compare it to numerical simulations. We decompose the covariance matrix into the disconnected (Gaussian) part, trispectrum from the modes outside the survey (supersample variance) and trispectrum from the modes inside the survey, and show how the different components contribute to the overall covariance matrix. We find the agreement with the simulations is at a 10 per cent level up to k ∼ 1 h Mpc−1. We show that all the connected components are dominated by the large-scale modes (k < 0.1 h Mpc−1), regardless of the value of the wave vectors k, k΄ of the covariance matrix, suggesting that one must be careful in applying the jackknife or bootstrap methods to the covariance matrix. We perform an eigenmode decomposition of the connected part of the covariance matrix, showing that at higher k, it is dominated by a single eigenmode. The full covariance matrix can be approximated as the disconnected part only, with the connected part being treated as an external nuisance parameter with a known scale dependence, and a known prior on its variance for a given survey volume. Finally, we provide a prescription for how to evaluate the covariance matrix from small box simulations without the need to simulate large volumes. cosmology: theory, large-scale structure of Universe 1 INTRODUCTION The distribution of matter in the Universe contains a wealth of information about the energy content of the Universe, its properties and evolution. Initial distribution is thought to be a Gaussian random field, but as a result of the gravitational instability, the tiny fluctuations in the initial matter distribution evolve non-linearly and produce non-Gaussian correlations. The simplest statistic for data analysis is the two-point correlation function, or its Fourier transform, the power spectrum. This contains significant cosmological information since it is sensitive to many parameters, and much of the information comes from deeply in the non-linear regime. The two-point analysis is the prime focus of many cosmological observables like weak lensing (WL) or galaxy clustering, which are the key observable probes in the current and future generation surveys like Dark Energy Survey (DES),1 Dark Energy Spectroscopic Instrument (DESI),2 Large Synoptic Survey Telescope (LSST),3Euclid4 etc. Our ability to extract useful constraints on the cosmological parameters from these surveys depend on our ability to model the statistical properties of the distribution of matter in the Universe, particularly the matter power spectrum, in the non-linear regime (k > 0.2 h−1 Mpc). In recent years, there have been many efforts in providing accurate estimates of the matter power spectrum for various cosmological models and redshifts using perturbation theory (Bernardeau et al. 2002), the halo model (Cooray & Sheth 2002), simulations (Smith et al. 2003; Heitmann et al. 2009, 2010; Lawrence et al. 2010; Takahashi et al. 2012) and semi-analytic models (Mohammed & Seljak 2014; Seljak & Vlah 2015). The current precision from simulations is at 1 per cent up to k ∼ 1 h Mpc−1, and can potentially be improved further with simulations, if necessary. There are additional complications such as the baryonic effects (Semboloni et al. 2011; van Daalen et al. 2011; Semboloni, Hoekstra & Schaye 2013; Mohammed et al. 2014), which redistribute the matter within the halo centres and change the power spectrum at high k, and it is likely that WL will not be able to extract much reliable information at high k (or ℓ) without including these effects. The same also applies to massive neutrinos, which suppress structure formation as they free stream. We will not model these processes here and focus on the dark matter part only. While the matter power spectrum predictions, in the absence of baryonic effects, are under control, for a complete analysis, one also needs its covariance matrix. The covariance matrix of the matter power spectrum is important in order to perform any statistical inference analysis on the cosmological data. Therefore, future surveys would require an accurate estimation of the covariance matrix of the matter power spectrum in order to perform the cosmological parameters estimation. An accurate quantification of the covariance matrix is crucial in order to derive constraints on the cosmological parameters using observables modelled on the matter power spectrum. Failing to do so will mislead the interpretation of the data. In the linear regime, computing the covariance matrix is reasonably straightforward, even if it can be computationally expensive for surveys with complicated masks. Because of the non-linear evolution of the matter density field, the Fourier modes become correlated, and this generates a non-Gaussian component to the covariance matrix (Scoccimarro, Zaldarriaga & Hui 1999; Cooray & Hu 2001; Hu & White 2001), which we also call the connected part. It affects both diagonal and non-diagonal elements of the covariance matrix. Its characterization has been developed in terms of the physically motivated halo model (Takada & Hu 2013; Mohammed & Seljak 2014), or numerical simulations (Takahashi et al. 2009; Sato et al. 2011; Harnois-Déraps & Pen 2012; Dodelson & Schneider 2013; Li, Hu & Takada 2014a; Blot et al. 2016, 2015). Another approach is using perturbation theory (PT), which is the focus of the work here. Our approach is to identify the perturbative terms that dominate the covariance matrix, and PT can be used to understand better the structure of the covariance matrix. This paper is organized as follows. In Section 2, we review the approach and describe the perturbation theory up to partial 1-loop calculations for the full covariance matrix of the matter power spectrum. In Section 3, we outline a comparison of our model with two sets of cosmological simulations. In Section 4, we decompose the covariance matrix based on simulations and the analytic model into a vector that fully describes the full covariance matrix using a principal component analysis (PCA). In Section 4.1, we study the degeneracies between covariance matrix model and other cosmological parameters and perform a Fisher information matrix analysis. We also address the question of how to evaluate covariance matrix from small box simulations. Finally, we conclude with a discussion in Section 5. During the completion of our work, an independent analysis including all PT terms up to 1-loop has been presented by Bertolini et al. (2016a), with results comparable to ours. 2 COVARIANCE MATRIX IN PERTURBATION THEORY The simplest of the statistics is the two-point correlation function of the matter density perturbations δ, or its Fourier transform matter power spectrum P(k), \begin{equation} \langle \tilde{\delta }(\boldsymbol {k}) \tilde{\delta }^{\star }(\boldsymbol {k}^{\prime }) \rangle = (2\pi )^3 \delta _{\rm D}(\boldsymbol {k}-\boldsymbol {k}^{\prime }) P(k), \end{equation} (1) where δD is the Dirac delta function, k is the magnitude of the wave vector |$\boldsymbol {k}$| and |$\tilde{\delta }(\boldsymbol {k})$| is the Fourier transform of the matter density perturbations δ, defined as \begin{equation} \tilde{\delta }(\boldsymbol {k}) = \int \delta (\boldsymbol {x}) {\rm e}^{{\rm i} \boldsymbol {k\cdot x}} {\rm d}^3{\boldsymbol x}. \end{equation} (2) For a given survey volume V, one can compute |$\tilde{\delta }(\boldsymbol {k})$| using equation (2). Then the power spectrum can be estimated by dividing the survey volume into shells centred at ki as \begin{equation} P(k_i) = V_{\rm f}\int _{k_i} \frac{{\rm d}^3\boldsymbol {k}}{V_{\rm s}(k_i)} \tilde{\delta }(\boldsymbol {k})\tilde{\delta }(-\boldsymbol {k}), \end{equation} (3) where |$V_{\rm s}(k_i) = 4\pi k_i^2 \Delta k$| is the volume of the ith shell, Δk is the width of the k shell and Vf = (2π)3/V is the volume of the fundamental cell in Fourier space. The matter power spectrum is a two-point function, its statistical properties defining its uncertainties and statistical properties, precision to which it can be measured, and correlations at various scales can be quantified in terms of a covariance matrix Cov(ki, kj), which is a four-point function and defined as \begin{equation} {\rm Cov}(k_i, k_j) \equiv \langle P({k_i}) P({k_j}) \rangle - \langle P({k_i}) \rangle \langle P({k_j}) \rangle , \end{equation} (4) where the angle brackets represent an ensemble average. We interchangeably use the notation Cov(ki, kj) and Covij. Similarly, we also interchangeably use P(ki) or Pi for the matter power spectrum. A full covariance matrix can be decomposed into \begin{equation} {\rm Cov}_{ij}^{\rm Full} = {\rm Cov}_{ij}^{\rm G} + {\rm Cov}_{ij}^{\rm NG}. \end{equation} (5) The first contribution is the disconnected contribution, also known as the Gaussian contribution, and is always present due to the random phases of the modes of the density field. As the matter power spectrum is computed by averaging over those Fourier modes, a smaller number of modes give larger statistical uncertainty over the mean power spectrum. Therefore, this contribution is inversely proportional to the number of modes. As we have only one survey volume to observe, the largest scale modes are very few, and hence the covariance matrix is dominated by this part for low k. When noise can be ignored, it is known as the sampling variance. It dominates along the diagonal elements of the covariance matrix: for a finite volume of the survey the modes are mixed, but if we bin the modes with binning Δk > 2π/R, where R3 = V is the survey volume, then the window is diagonal. The Gaussian contribution can be estimated as \begin{equation} {\rm Cov}_{ij}^{\rm G} = \frac{2}{N_k} \delta _{ij} P(k_i) P(k_j), \end{equation} (6) where Nk = Vs/Vf is the total number of k modes in the corresponding shell, and δij is the Kronecker delta function which is unity for i = j, and zero otherwise. It scales inversely with the fourth power of the growth factor for its redshift evolution. The second contribution is the connected part, or the non-Gaussian part, and can be expressed in terms of the trispectrum as \begin{equation} {\rm Cov}_{ij}^{\rm NG} = \frac{\bar{T}({k_1,-k_1,k_2,-k_2})}{V}, \end{equation} (7) where |$\bar{T}$| is the bin averaged trispectrum, the fourth-order connected moment of the density perturbation, given by \begin{eqnarray} \bar{T}({k_1,-k_1,k_2,-k_2}) = \int _{V_{k_i}} \int _{V_{k_j}} \frac{{\rm d}^3 {\boldsymbol k}_1}{V_{k_i}} \frac{{\rm d}^3 {\boldsymbol k}_2}{V_{k_j}} T({\boldsymbol k}_1,-{\boldsymbol k}_1,\boldsymbol {k}_2,-\boldsymbol {k}_2), \nonumber\\ \end{eqnarray} (8) where the two integral on the right are across the ith and jth shell. The Gaussian part of the covariance depends on the binning scheme in the k shells, whereas the non-Gaussian part is fairly independent of any binning scheme, except for a small dependence for low k (see Section 2.2 for details). However, both terms scale inversely with the survey volume. We split the non-Gaussian contribution into the parts that come from modes outside the survey (SSC) and inside the survey, and we further split the latter into tree-level and 1-loop terms, \begin{equation} {\rm Cov}_{ij}^{\rm NG} = {\rm Cov}_{ij}^{\rm SSC}+ {\rm Cov}_{ij}^{{\rm tree-level}} + {\rm Cov}_{ij}^{{\rm 1-loop}} \end{equation} (9) where the terms on the right are the supersample covariance (SSC), tree-level and 1-loop contributions, respectively. The first two terms are both tree-level contributions, and their split is motivated by the different calculational approach to the two terms. To this, we add a specific 1-loop term because it significantly improves the results. In principle, we should include all 1-loop terms as well as add counter-terms to 1-loop, but we will not explore these extensions here (see Bertolini et al. 2016a,b). 2.1 Supersample covariance The first NG contribution is known as the supersample covariance (SSC). It is caused by the coupling between all k modes that are larger than the survey, and the modes inside the survey. It is caused by the survey window effects mixing the modes and does not show up in periodic box simulations. Its leading effect can be derived as a tree level trispectrum in the squeezed limit. This contribution was first pointed out by Hamilton, Rimes & Scoccimarro (2006), and further studied by Sefusatti et al. (2006), Sato et al. (2009), Takahashi et al. (2009), Kayo, Takada & Jain (2013), Takada & Bridle (2007), Takada & Jain (2009) and de Putter et al. (2012). Because this contribution comes from the mode that is constant across the survey, it can be viewed as a local curvature term within the survey, which can further be mimicked by a change in the background density δb. Therefore, this term can be modelled completely by the response of the matter power spectrum to the change in background density (Takada & Hu 2013; Li et al. 2014a; Li, Hu & Takada 2014b). Further, this response can be accurately calculated using the separate universe simulations, where the same effect can be mimicked by a change in cosmological parameters (Sirko 2005; Baldauf et al. 2011, 2016b; Gnedin, Kravtsov & Rudd 2011; Wagner et al. 2015). This term is the dominant contribution in the quasi-linear regime, and can be modelled as \begin{equation} {\rm Cov}_{ij}^{\rm SSC} = \sigma _{\rm b}^2 \frac{\mathrm{\partial} P(k_i)}{\mathrm{\partial} \delta _{\rm b}} \frac{\mathrm{\partial} P(k_j)}{\mathrm{\partial} \delta _{\rm b}}, \end{equation} (10) where ∂P(ki)/∂δb is the response of the matter power spectrum to the change in the background density δb, and |$\sigma ^2_{\rm b}$| is the variance of the mean density field in the survey window, defined as \begin{equation} \sigma _{\rm b}^2 \equiv \langle \delta _{\rm b}^2\rangle = \frac{1}{V^2} \int \frac{{\rm d}^3\boldsymbol {q}}{(2\pi )^3} \tilde{W}(\boldsymbol {q})^2 P(q), \end{equation} (11) where |$\tilde{W}(\boldsymbol {q})$| is the survey window function (with units of volume). The responses depend on whether we divide by the local density (SSC-local), in which case the low k limit is \begin{equation} \frac{{\rm d}\ln P(k)}{{\rm d}\delta _{\rm b}}_{{\rm SSC-local}} ={5 \over 21}-{1 \over 3}{{\rm d} \ln P(k) \over {\rm d} \ln k}, \end{equation} (12) or we do not divide by the local density (SSC-global), which gives low k limit of \begin{equation} \frac{{\rm d}\ln P(k)}{{\rm d}\delta _{\rm b}}_{{\rm SSC-global}} ={47 \over 21}-{1 \over 3}{{\rm d} \ln P(k) \over {\rm d} \ln k}. \end{equation} (13) We will use numerical derivative obtained from separate universe simulations (Li et al. 2014a). In general, SSC term is well understood and can easily be computed by running two separate universe simulations. 2.2 Tree-level covariance In the tree-level perturbation theory, and ignoring the window function effects that lead to SSC term above, the full trispectrum can be written as \begin{eqnarray} &&{T(\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3,\boldsymbol {k}_4)} \nonumber\\ &&{\quad = 4\left[ F_2(\boldsymbol {k}_{12},-\boldsymbol {k}_1) F_2(\boldsymbol {k}_{34},\boldsymbol {k}_4) P(k_1)P(k_{34})P(k_3) + \rm {cyc.}\right] }\nonumber\\ &&{ \qquad +\, 6\left[F_3(\boldsymbol {k}_1,\boldsymbol {k}_2,\boldsymbol {k}_3) P(k_1)P(k_2)P(k_3) + \rm {cyc.}\right],} \end{eqnarray} (14) where |$\boldsymbol {k}_{ij} = \boldsymbol {k}_i+\boldsymbol {k}_j$|⁠, and cyc. denotes the cyclic permutations of the arguments. The kernels F2 and F3 are the solutions to the equations of motion of gravitational instability to second and third order, respectively. Although this is the full expression for the trispectrum, only parallelogram configurations contribute to the non-Gaussian covariance, as in equation (8). Therefore, the parallelogram configuration of the trispectrum that contributes to the non-Gaussian covariance is given by (Scoccimarro et al. 1999) \begin{eqnarray} \bar{T}(k_i, k_j) = \int _{k_i} \frac{{\rm d}^3 k_1}{V_{\rm s}(k_i)}\int _{k_j} \frac{{\rm d}^3 k_2}{V_{\rm s}(k_j)} \nonumber\\ {\times} (12 F_3^{({\rm s})}(\boldsymbol k_1, -\boldsymbol k_1, \boldsymbol k_2) P_{\rm L}(\boldsymbol k_1)^2 P_{\rm L}(\boldsymbol k_2) \nonumber \\ &&+\,8 F^{({\rm s})}_{2}(\boldsymbol k_1 - \boldsymbol k_2, \boldsymbol k_2)^2 P_{\rm L}(\boldsymbol k_2)^2 P_{\rm L}(\boldsymbol k_2 - \boldsymbol k_1) \nonumber \\ &&+\,8 F^{({\rm s})}_{2}({\boldsymbol k}_1 - \boldsymbol k_2,{\boldsymbol k}_2)F^{({\rm s})}_{2}({\boldsymbol k}_2 - {\boldsymbol k}_1,{\boldsymbol k}_1) P_{\rm L}({\boldsymbol k}_2) P_{\rm L}({\boldsymbol k}_1) P_{\rm L} \nonumber \\ {\times}\,({\boldsymbol k}_2 - {\boldsymbol k}_1) + \lbrace {\boldsymbol k}_1\leftrightarrow {\boldsymbol k}_2 \rbrace ), \end{eqnarray} (15) where |$F^{({\rm s})}_2$| and |$F^{({\rm s})}_3$| are the standard Eulerian PT kernels (see e.g. Bernardeau et al. 2002). The diagonal part can now be explicitly written as \begin{eqnarray} \frac{\bar{T}(k_i, k_i)}{ \left[P_{\rm L}(k_i)\right]^3} = -\frac{44}{63} + \frac{1}{49}\int _{-1}^1 {\rm d}\mu \left( 3+10\mu \right)^2 \frac{P_{\rm L}\left( k \sqrt{2(1-\mu )} \right)}{P_{\rm L}(k_i)}. \nonumber\\ \end{eqnarray} (16) Using the approximation suggested in Scoccimarro et al. (1999), one gets |$\bar{T}(k_i, k_i) / \left[P_{\rm L}(k_i)\right]^3\simeq 232/441$|⁠, although somewhat better approximation (at low k) can be obtained |$\bar{T}(k_i, k_i) / \left[P_{\rm L}(k_i)\right]^3\simeq 454/441 - 62/343 \,n_{\rm s}\simeq 2620/3087$|⁠. Fig. 1 shows the trispectrum contribution to the covariance for few different binning schemes. Binning makes essentially no effect, as expected. This would appear to be in conflict with Takahashi et al. (2008), where they studied binning dependence of the ratio between non-linear and linear mode, finding strong dependence. However, this is caused by division with the linear mode, which induces Gaussian fluctuations in the ratio. Figure 1. Open in new tabDownload slide The trispectrum contribution to the covariance for different binning schemes in different colours. Blue, green, red and black correspond to a linear bin width (Δk) of 0.01, 0.025, 0.05 and 0.001, respectively. Figure 1. Open in new tabDownload slide The trispectrum contribution to the covariance for different binning schemes in different colours. Blue, green, red and black correspond to a linear bin width (Δk) of 0.01, 0.025, 0.05 and 0.001, respectively. 2.3 1-Loop covariance In this section, we describe the calculations of the 1-loop non-Gaussian covariance contribution. A full calculation has been recently presented in Bertolini et al. (2016a,b). Here, we will instead do a simplified calculation, where we first derive a functional derivative of high k power to the low q power. In the next step, we compute the fluctuations of low q power due to the finite volume effects, i.e. the sampling variance fluctuations. Finally, we obtain the covariance by combining the two. The 1-loop PT on the power spectrum has been intensively studied in previous works such as Juszkiewicz (1981), Vishniac (1983), Juszkiewicz, Sonoda & Barrow (1984), Coles (1990), Suto & Sasaki (1991), Makino, Sasaki & Suto (1992), Jain & Bertschinger (1994), Baugh & Efstathiou (1994), Lokas et al. (1996), Scoccimarro & Frieman (1996) and Bernardeau et al. (2002). Here, we briefly review these results. In perturbation theory framework, one can describe the full matter power spectrum as \begin{equation} P(k,z) = P^{(0)}(k,z) + P^{(1)}(k,z) +\cdots + P^{(n)}(k,z), \end{equation} (17) where the superscript n denotes the n-loop contribution. P(0)(k, z), or the 0-loop contribution, is just the linear power spectrum, such that P(0)(k, z) = PL(k)D2(z). The 1-loop contribution, P(1)(k, z), consists of two terms and can be written as \begin{eqnarray} P^{(1)}(k,z) = P_{13}(k,z) + P_{22}(k,z) \nonumber \\ = 6\int F_3^{({\rm s})}(\boldsymbol {k}, \tilde{\boldsymbol q},-\tilde{\boldsymbol q})P_{\rm L}(k,z)P_{\rm L}(\tilde{q},z) {\rm d}^3\tilde{q} \nonumber\\ &&+\, 2\int \left( F_2^{({\rm s})}(\boldsymbol {k}-\tilde{\boldsymbol q},\tilde{\boldsymbol q})\right)^2 P_{\rm L}(|\boldsymbol {k}-\tilde{\boldsymbol q}|,z)P_{\rm L}(\tilde{q},z) {\rm d}^3 \tilde{q}. \nonumber\\ \end{eqnarray} (18) Here Pij denotes the amplitude given by a connected diagram representing the contribution from 〈δiδj〉c to the power spectrum. We have assumed Gaussian initial conditions, for which Pij vanishes if i + j is odd. The two contributions, P13 and P22, have different structure. P13 is, in general, negative and describes the decorrelation of the propagator, while P22 is positive definite and describes the effects of mode coupling between modes with wave vectors |$\boldsymbol {k}-\boldsymbol {q}$| and |$\boldsymbol {q}$|⁠. The functional derivatives of these function with respect to the linear power spectrum yields \begin{eqnarray} N \frac{\delta P_{13}( \boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} = 3 P_{\rm L}(\boldsymbol k) \int {\rm d}^3 \tilde{q}\, F_3^{\rm (s)}(\boldsymbol {k},\tilde{ \boldsymbol q },-\tilde{ \boldsymbol q }) \delta ^{\rm D}(\boldsymbol q-\tilde{ \boldsymbol q})\nonumber\\ &&+\, 3 \delta ^{\rm D} (\boldsymbol k - \boldsymbol q) \int \frac{{\rm d}^3 \tilde{q}}{(2\pi )^3} \, F_3^{\rm (s)}(\boldsymbol {k},\tilde{ \boldsymbol q },-\tilde{ \boldsymbol q }) P_{\rm L}(\tilde{ \boldsymbol q}) \nonumber \\ = 3 F_3^{\rm (s)}(\boldsymbol {k},\boldsymbol {q},-\boldsymbol {q}) P_{\rm L}(\boldsymbol k) \!+\! 3 P_{13}(\boldsymbol k)/P_{\rm L}(\boldsymbol k) \delta ^{\rm D} (\boldsymbol k - \boldsymbol q), \nonumber \\ N \frac{\delta P_{22}(\boldsymbol k)}{\delta P_{\rm L}(\boldsymbol q)} = 4 \int {\rm d}^3 \tilde{q} \left[ F_2^{\rm (s)}(\boldsymbol {k}-\tilde{ \boldsymbol q },\tilde{ \boldsymbol q })\right]^2 P_{\rm L}(\boldsymbol {k}-{\boldsymbol q}) \delta _{\rm D} (\boldsymbol q-\tilde{\boldsymbol q}) \nonumber \\ = 4 \left( F_2^{\rm (s)}(\boldsymbol {k}\!-\!\boldsymbol {{q}},\boldsymbol {{q}})\right)^2 P_{\rm L}(\boldsymbol {k}-\boldsymbol {q}), \end{eqnarray} (19) where N is the normalization pre-factor to be chosen below. Dropping the δ function part of the result we have for the total 1-loop result: \begin{eqnarray} N \frac{\delta P_{{\rm 1-loop}}(\boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} = N \frac{\delta P_{22}(\boldsymbol k)}{\delta P_{\rm L}(\boldsymbol q)} + 2 N \frac{\delta P_{13}( \boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} \nonumber \\ = 4 \left( F_2^{\rm (s)}(\boldsymbol {k}-\boldsymbol {{q}},\boldsymbol {{q}})\right)^2 P_{\rm L}(|\boldsymbol {k}-\boldsymbol {q}|) \nonumber\\ &&+\, 6 F_3^{\rm (s)}(\boldsymbol {k},\boldsymbol {q},-\boldsymbol {q})P_{\rm L}(k). \end{eqnarray} (20) Performing the angle averaging |$\langle X(\boldsymbol q) \rangle _\Omega = \int \frac{{\rm d} \Omega _{\boldsymbol q}}{4\pi } X(\boldsymbol q)$|⁠, we get \begin{eqnarray} 2 N \left\langle \frac{\delta P_{{\rm 1-loop}}(\boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} \right\rangle _\Omega = 2 \int {\rm d}\mu \,\left( F_2^{\rm (s)}(\boldsymbol {k}-\boldsymbol {{q}},\boldsymbol {{q}})\right)^2 P_{\rm L}(|\boldsymbol {k}-\boldsymbol {q}|)\nonumber\\ &&+\, 3 \int {\rm d} \mu \, F_3^{\rm (s)}(\boldsymbol {k},\boldsymbol {q},-\boldsymbol {q})P_{\rm L}(k). \end{eqnarray} (21) It is instructive to look at the limiting cases. From above, it readily follows (Sherwin & Zaldarriaga 2012): \begin{eqnarray} 2 N \left\langle \frac{\delta P_{{\rm 1-loop}}(\boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} \right\rangle _\Omega &\sim & \frac{2519}{2205} P_{\rm L}( \boldsymbol k) - \frac{47}{105} k P_{\rm L}^{\prime }( \boldsymbol k) \nonumber\\ &&+ \frac{1}{10} k^2 P_{\rm L}^{\prime \prime }( \boldsymbol k),\,{\rm as}\,\,k \gg q, \end{eqnarray} (22) and similarly \begin{eqnarray} 2 N \left\langle \frac{\delta P_{{\rm 1-loop}}(\boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} \right\rangle _\Omega &\sim & \frac{9 k^4}{49 q^4} P_{\rm L}( \boldsymbol q) - \left( \frac{61 k^2}{315 q^2} - \frac{4 k^4}{105 q^4} \right)\nonumber\\ {\times} \,P_{\rm L}( \boldsymbol k),\,\,{\rm as}\,\,k \ll q. \end{eqnarray} (23) The result of the fill one-loop functional derivative given in equation (21) and the corresponding k ≫ q expansion is shown in Fig. 2. Figure 2. Open in new tabDownload slide The 1-loop functional derivatives |$2N \langle \delta P_{{\rm 1-loop}} (k)/\delta P_{\rm L}(q) \rangle _\Omega$| (black dashed lines) given in equation (21) at redshift z = 0.0. Values of q are from left to right are { 0.01, 0.018, 0.032, 0.056, 0.10, 0.18, 0.32, 0.56, 1.0, 1.8, 3.2}[h−1 Mpc]. Red lines correspond to contributions in the limiting cases k ≫ q in equation (22). Figure 2. Open in new tabDownload slide The 1-loop functional derivatives |$2N \langle \delta P_{{\rm 1-loop}} (k)/\delta P_{\rm L}(q) \rangle _\Omega$| (black dashed lines) given in equation (21) at redshift z = 0.0. Values of q are from left to right are { 0.01, 0.018, 0.032, 0.056, 0.10, 0.18, 0.32, 0.56, 1.0, 1.8, 3.2}[h−1 Mpc]. Red lines correspond to contributions in the limiting cases k ≫ q in equation (22). We can introduce the normalized functional derivative, |${\boldsymbol {V}(\boldsymbol {q},\boldsymbol {k})}$|⁠, by dividing the expression above by the power per mode Δ2(q) = 4πq3P(q)/(2π)3 as in Nishimichi, Bernardeau & Taruya (2014), \begin{equation} \boldsymbol {V}(q,k) = {P_{\rm L}(q) \over \Delta ^2(q)}\left\langle \frac{\delta P_{{\rm 1-loop}}(\boldsymbol k)}{\delta P_{\rm L}( \boldsymbol q)} \right\rangle _\Omega. \end{equation} (24) In the low q limit, i.e. q ≪ k, the functional derivative becomes independent of q: \begin{equation} \lim _{q/k\rightarrow 0} \boldsymbol {V}(q,k) = W(k) P(k), \end{equation} (25) where \begin{equation} \boldsymbol {W}(k) = \frac{2519}{2205} - \frac{47}{105} \frac{{\rm d}\ln P(k)}{{\rm d}\ln k} + \frac{1}{10} \frac{{\rm d}^2\ln P(k)}{{\rm d}\ln k^2}. \end{equation} (26) The expressions above give the mean response at a given k to the power at a given q. At low q, the power at a given q will fluctuate due to the finite volume. The Gaussian part of the covariance scales as 2P2/Nk and number of modes Nk scales as Nk = Vk3Δln k/(2π)3. The full expression for 1-loop covariance contribution can now be written as \begin{equation} {\rm Cov}_{ij}^{{\rm 1-loop}} = \left(\frac{1}{V\pi ^2} \int P^2_{\rm L}(q) q^2 \boldsymbol {V}(q,k_i) \boldsymbol {V}(q,k_j) {\rm d}q\right). \end{equation} (27) Equation (27) is an integral of the term |$P^2_{\rm L}(q)q^2$| and the functional derivatives. We can define \begin{equation} S=\left( \frac{1}{V\pi ^2} \int P^2_{\rm L}(q) q^2{\rm d}q \right). \end{equation} (28) Fig. 3 shows the integral of S to qmax, showing that the scales at which this contribution is important are mostly linear (q < qnl ∼ 0.2 h Mpc−1).This suggests the non-linear corrections are likely to be small, especially at higher redshifts. Figure 3. Open in new tabDownload slide The 1-loop variance S as a function of qmax (defined in equation 28) for a volume of 1(h−1 Gpc)3. The integral converges to a few per cent at qmax ∼ 0.3 h Mpc−1. We also show the effect of damping for z = 0 (qnl = 0.3 h Mpc−1) and for z = 1 (qnl = 0.7 h Mpc−1). Figure 3. Open in new tabDownload slide The 1-loop variance S as a function of qmax (defined in equation 28) for a volume of 1(h−1 Gpc)3. The integral converges to a few per cent at qmax ∼ 0.3 h Mpc−1. We also show the effect of damping for z = 0 (qnl = 0.3 h Mpc−1) and for z = 1 (qnl = 0.7 h Mpc−1). In the low q limit, the 1-loop contribute separates into \begin{equation} {\rm Cov}_{ij}^{{\rm 1-loop}} = S\boldsymbol {W}(k_i)P(k_i) \boldsymbol {W}(k_j)P(k_j). \end{equation} (29) We will also explore simple extensions beyond the 1-loop. Since the integral in equation (28) extends into the non-linear regime, its high q contribution may not be reliable. We expect the response to high q modes to be suppressed in functional derivatives, so their effect is reduced (Nishimichi et al. 2014). The transition is governed by the non-linear scale qnl where |$\Delta ^2(q)=\int q^2{\rm d}qP(q)W_R^2(q)/2\pi ^2 = 1.35^2$|⁠, where WR is a Gaussian window WR = exp [−(0.46q/qnl)2/2], so the expectation is that qnl ∼ 0.3 h Mpc−1 at z = 0. Following Nishimichi et al. (2014). we will assume a Lorentzian damping form and explore the 1-loop expression: \begin{eqnarray} {\rm Cov}_{ij}^{\rm NL} = \left( \frac{1}{V\pi ^2} \int {P^2_{\rm L}(q)\over (1+ (q/q_{{\rm nl}})^2)^2} q^3 \boldsymbol {V}(q,k_i) \boldsymbol {V}(q,k_j) {\rm d}\ln q \right). \nonumber\\ \end{eqnarray} (30) Fig. 3 also shows the integral of S along with the damping for z = 0. We see that damping slightly reduces the value of S at z = 0. The damping effect is even smaller for higher redshifts. In the subsequent figures, we refer the non-linear model to the sum of all contributions to covariance where 1-loop term is given by equation (30). 3 COMPARISON WITH SIMULATIONS We used two sets of cosmological simulations, from Blot et al. (2015, hereafter B15) and Li et al. (2014a, hereafter L14), to compare with our model. For both simulations, the cosmology is the flat Λ cold dark matter (ΛCDM) with particular cosmological parameters values listed in Table 1. The simulation suite from L14 also contains an explicit calculation of the SSC term measured using a corresponding separate universe simulation. B15 have a much larger volume, with a total effective volume greater than 3000 h−3 Gpc3. L14 gives 3584 matter power spectra at redshift zero, which we use to generate a full covariance matrix. B15 provides their covariance matrix at four different redshifts: 0.0, 0.5, 1.0, 2.0. Table 1. Specifications for the two sets of the cosmological simulations used for comparison. Parameter . Number of simulations . Box size (h−1 Mpc) . Ωm . σ8 . h . ns . Ωb . w . L14 3584 500 0.286 0.82 0.7 0.96 0.047 −1.0 B15 12 288 650 0.257 0.801 0.72 0.963 0.044 −1.0 Parameter . Number of simulations . Box size (h−1 Mpc) . Ωm . σ8 . h . ns . Ωb . w . L14 3584 500 0.286 0.82 0.7 0.96 0.047 −1.0 B15 12 288 650 0.257 0.801 0.72 0.963 0.044 −1.0 Open in new tab Table 1. Specifications for the two sets of the cosmological simulations used for comparison. Parameter . Number of simulations . Box size (h−1 Mpc) . Ωm . σ8 . h . ns . Ωb . w . L14 3584 500 0.286 0.82 0.7 0.96 0.047 −1.0 B15 12 288 650 0.257 0.801 0.72 0.963 0.044 −1.0 Parameter . Number of simulations . Box size (h−1 Mpc) . Ωm . σ8 . h . ns . Ωb . w . L14 3584 500 0.286 0.82 0.7 0.96 0.047 −1.0 B15 12 288 650 0.257 0.801 0.72 0.963 0.044 −1.0 Open in new tab Figs 4 –7 show the comparison of our model for the covariance matrix to B15 covariance matrix at redshifts 0.0, 0.5, 1.0, 2.0, respectively. Each panel shows various contributions (1-loop, tree-level, Gaussian) to the covariance, a full model and the corresponding simulation output. Figure 4. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 0.0. The total curve (red) is the sum of all contributions with 1-loop given by equation (27), whereas non-linear curve (dashed green) is the sum of all contributions with 1-loop given by equation (30). Similar convention applies to the subsequent figures. Figure 4. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 0.0. The total curve (red) is the sum of all contributions with 1-loop given by equation (27), whereas non-linear curve (dashed green) is the sum of all contributions with 1-loop given by equation (30). Similar convention applies to the subsequent figures. Figure 5. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 0.5. For non-linear curve (dashed green), we used qmax = 0.6. Figure 5. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 0.5. For non-linear curve (dashed green), we used qmax = 0.6. Figure 6. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 1.0. Figure 6. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 1.0. Figure 7. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 2.0. Figure 7. Open in new tabDownload slide Comparing analytic model with B15 covariance using 1-loop exact functional derivatives at redshift 2.0. At redshift zero, there is a good agreement in quasi-linear regime up to k ∼ 0.2 h−1 Mpc, but as we go to the non-linear regime, the 1-loop term becomes large and overestimates the covariance. We have explored the non-linear model of equation (30) with one free parameter and found that with qnl = 0.2 we find a good agreement even for higher k. A similar trend can also be seen for a comparison with L14 covariance at redshift zero in Fig. 8. Figure 8. Open in new tabDownload slide Comparing analytic model with L14 covariance using 1-loop exact functional derivatives at redshift 0.0. Figure 8. Open in new tabDownload slide Comparing analytic model with L14 covariance using 1-loop exact functional derivatives at redshift 0.0. As we go higher in redshift, the 1-loop term scales inversely with the eighth power of the growth factor whereas both tree level terms scale inversely with the sixth power of the growth factor. Therefore, the 1-loop contribution to the covariance decreases more rapidly than other contributions. Fig. 5 also shows the non-linear model of equation (30) with qnl = 0.6. In addition, we expect the correction to S from non-linear scales to become less important, because qnl increases. As a result, we find a much better agreement between our model and the simulations up to much higher k. In Figs 9 and 10, we show the comparison of the covariance of L14 data set, with local and global SSC term, respectively, added to the calculation. For large scales, nearly quasi-linear scales, global SSC term dominates, and the 1-loop contribution is subdominant. Therefore, the agreement with the simulations is better than without SSC contribution even without additional non-linear correction. For higher redshifts, the 1-loop term is suppressed relative to SSC, and we expect an even better agreement with simulations. Figure 9. Open in new tabDownload slide Comparing analytic model with L14 covariance using 1-loop exact functional derivatives at redshift 0.0 with local SSC term. The SSC term contributes to simulations (thin black line), total (solid red line) and non-linear model (dashed green line). Figure 9. Open in new tabDownload slide Comparing analytic model with L14 covariance using 1-loop exact functional derivatives at redshift 0.0 with local SSC term. The SSC term contributes to simulations (thin black line), total (solid red line) and non-linear model (dashed green line). Figure 10. Open in new tabDownload slide Comparing analytic model with L14 covariance using 1-loop exact functional derivatives at redshift 0.0 with global SSC term. The SSC term contributes to simulations (thin black line), total (solid red line) and non-linear model (dashed green line). Figure 10. Open in new tabDownload slide Comparing analytic model with L14 covariance using 1-loop exact functional derivatives at redshift 0.0 with global SSC term. The SSC term contributes to simulations (thin black line), total (solid red line) and non-linear model (dashed green line). 4 EIGENMODE DECOMPOSITION OF THE COVARIANCE MATRIX The SSC part of the covariance is a rank 1 tensor (equation 10), which means that it has a single non-zero eigenvalue. The 1-loop term expression in equation (27) also suggests a similar form. While the tree level trispectrum does not decompose like that, it is mostly important at low k, where the Gaussian contribution dominates. It is thus interesting to ask about the eigenvalue structure of the entire covariance matrix (Harnois-Déraps & Pen 2012). Here, we use a PCA on the non-Gaussian part of the covariance. We start by calculating the covariance matrix of global SSC version of L14 (total 3584 power spectra) data set using equation (4). We perform a PCA of the matrix Cij, such that \begin{equation} C_{ij} = \frac{{\rm Cov}_{ij}^{\rm full} - {\rm Cov}_{ij}^{\rm Gauss}}{P(k_i)P(k_j)}. \end{equation} (31) The eigenvalues of the matrix Cij are shown in the left-hand panel of Fig. 11. Clearly, one eigenvalue is much larger than all of the others. Therefore, the full matrix Cij can be well approximated using a single principal component |$d_1 = \sqrt{\lambda _1}v_1$|⁠, where λ1 is the largest eigenvalue and v1 is the corresponding eigenvector, such that |$C_{ij}=d_id_j^{\dagger}$|⁠. The right-hand panel of Fig. 11 shows the first two principal components (d1 and |$d_2 = \sqrt{\lambda _2}v_2$|⁠). Figure 11. Open in new tabDownload slide Showing eigenvalues and eigenvectors (only two largest) for the covariance matrix (with SSC) and the reduced covariance matrix of L14 data set. Dotted horizontal lines are the reference to the flatness. Figure 11. Open in new tabDownload slide Showing eigenvalues and eigenvectors (only two largest) for the covariance matrix (with SSC) and the reduced covariance matrix of L14 data set. Dotted horizontal lines are the reference to the flatness. In order to show the dominance of the first component explicitly, we perform another exercise by removing the first eigenvector from the underlying data set. To do so, we subtract a contribution from each of the power spectra, such that the new power spectra (⁠|$P^{\prime }_i(k)$|⁠) are \begin{equation} P^{\prime }(k_i) = P(k_i) - \alpha d_1(k_i) \langle P(k_i)\rangle. \end{equation} (32) Here α is the best-fitting coefficient when fitting Pi(k) to (1 + αd1(k))〈P(k)〉. Fig. 12 shows the distribution of α, which is very Gaussian. This is not surprising, since it is dominated by the SSC term variance on the scale of the survey box, which is Gaussian distributed. We computed the covariance matrix of the new reduced data set P΄(ki). The corresponding eigenvalues can be seen in the left-hand panel of Fig. 11. It can be noticed that the largest eigenvalue of the reduced covariance matrix is very close to the second largest eigenvalue of the original covariance. Similar trend can also be seen in the eigenvectors in the right-hand panel of Fig. 11. Figure 12. Open in new tabDownload slide Distribution of α as in equation (32). Figure 12. Open in new tabDownload slide Distribution of α as in equation (32). We also performed the diagonal decomposition of the analytic matrix (using equation 30), i.e. using the covariance from our analytic model. Fig. 13 shows the comparison of the analytic and numerical covariance matrix decomposition, particularly for the first eigenvector and the diagonal part of the covariance. The left-hand panel shows the comparison without the SSC term whereas the right-hand panel shows the comparison with global SSC term. As expected, the agreement improves for the full covariance matrix with SSC term. Fig. 14 shows a comparison of the single eigenvector model of the full covariance matrix with the L14 simulations, as well as the diagonal. We see that analytic model makes about 20 per cent error at k ∼ 1 h Mpc−1 without SSC, and about 10 per cent error with SSC. The non-linear model where high q modes are damped reduces this error to a few per cent only in the SSC case, in agreement with previous figures. Figure 13. Open in new tabDownload slide The first eigenvector d1 (in dashed lines) and the diagonal component of the covariance matrix for analytic (1-loop and damped 1-loop, which we call non-linear) versus numerical approach, at z = 0. The left-hand panel shows the covariance without SSC term, the middle panel with global SSC term. We see that analytic model makes about 20 per cent error at k ∼ 1 h Mpc−1 without SSC, and about 10 per cent error with SSC. The non-linear model reduces this error to a few per cent only in the SSC case, in agreement with previous figures. The right-hand panel shows the relative deviations between two different realizations of the power spectrum (different pair in different colours). This shows that the shape of the eigenvector can be recovered with only two simulations at higher k. Figure 13. Open in new tabDownload slide The first eigenvector d1 (in dashed lines) and the diagonal component of the covariance matrix for analytic (1-loop and damped 1-loop, which we call non-linear) versus numerical approach, at z = 0. The left-hand panel shows the covariance without SSC term, the middle panel with global SSC term. We see that analytic model makes about 20 per cent error at k ∼ 1 h Mpc−1 without SSC, and about 10 per cent error with SSC. The non-linear model reduces this error to a few per cent only in the SSC case, in agreement with previous figures. The right-hand panel shows the relative deviations between two different realizations of the power spectrum (different pair in different colours). This shows that the shape of the eigenvector can be recovered with only two simulations at higher k. Figure 14. Open in new tabDownload slide Comparing analytic model with L14 covariance (global SSC case). We compare the single eigenvector model (dashed green) and total 1-loop predictions (solid thick red) against simulations (thin black). The residual covariance after subtracting the single eigenvector component is shown with dotted line. Figure 14. Open in new tabDownload slide Comparing analytic model with L14 covariance (global SSC case). We compare the single eigenvector model (dashed green) and total 1-loop predictions (solid thick red) against simulations (thin black). The residual covariance after subtracting the single eigenvector component is shown with dotted line. Having the decomposition of the analytic form of the covariance matrix gives us a handle to perform the redshift evolution of the principle component without explicitly computing the covariance at each redshift. Fig. 15 shows the first eigenvector at five different redshifts (with global SSC). Figure 15. Open in new tabDownload slide Redshift evolution of the first eigenvector of the full covariance matrix (analytic) using L14 cosmological model including SSC covariance. Figure 15. Open in new tabDownload slide Redshift evolution of the first eigenvector of the full covariance matrix (analytic) using L14 cosmological model including SSC covariance. 4.1 Degeneracies with cosmological parameters The analysis above suggests that instead of using a full covariance matrix, one can use the Gaussian covariance, and replace the non-Gaussian part with a fictitious external parameter α, whose response is given by the first eigenvector d1(k)P(k), as in equation (32). The parameter α has zero expected mean, with 〈α2〉 determined by the covariance matrix calculations above. However, one can also try to determine it from the data. In this section, we explore the degeneracies between the first eigenvector of the full covariance matrix with SSC contribution to the response of various cosmological parameters to the matter power spectrum. A related analysis with just SSC has been performed in Li et al. (2014b). We use the eigenvector of the full covariance matrix of the L14 data set. The left-hand panel of Fig. 16 shows the comparison of this response to the cosmological parameters responses. At low k, there is a degeneracy between the d1(k) and the amplitude of the fluctuations (σ8) to the matter power spectrum. In this regime, the covariance matrix calculation of the first eigenvalue gives a prior, shown as the horizontal dashed line on the right-hand side. However, this degeneracy is broken in the highly non-linear regime (k > 0.3 h Mpc−1), suggesting that the data may determine the covariance matrix amplitude if one uses high k information. Figure 16. Open in new tabDownload slide Left: the response of the cosmological parameters and the covariance parameter to the matter power spectrum. Right: showing the forecast errors on each parameter for fixed α (in dashed lines) and marginalizing over α (in solid lines). Figure 16. Open in new tabDownload slide Left: the response of the cosmological parameters and the covariance parameter to the matter power spectrum. Right: showing the forecast errors on each parameter for fixed α (in dashed lines) and marginalizing over α (in solid lines). To explore the effects of this degeneracy on cosmological parameter estimation, we performed a Fisher matrix analysis for the two cases: (i) when α is fixed, and hence, only the response of the cosmological parameters plays a role, (ii) while marginalizing over α, incorporating the response of the α. The full fisher matrix can be written as \begin{equation} F_{\mu \nu } = \frac{\mathrm{\partial} P(k_i)}{\mathrm{\partial} p_{\mu }} C^{\rm Gauss}_{ij} \frac{\mathrm{\partial} P(k_j)}{\mathrm{\partial} p_{\nu }}, \end{equation} (33) where p is a set containing cosmological parameters and α. We used three cosmological parameters for this analysis – amplitude of the fluctuations σ8, spectral index ns and Hubble constant in the units 100 km s−1 Mpc−1 h. Therefore, \begin{equation} p \equiv \lbrace \sigma _8, n_{\rm s}, h, \alpha \rbrace. \end{equation} (34) For this particular exercise, we used only the Gaussian part of the covariance matrix. We carried out this exercise for different kmax (from 0.1 to 1.2) and computed the expected errors on each parameter. The right-hand panel of Fig. 16 shows the expected errors on the parameters. The solid lines represent the expected errors on the parameters when marginalized over α and dashed lines when α kept fixed. Because of the strong degeneracy between α and σ8, the marginalization errors are large for small k. However, as we increase kmax this degeneracy is broken, giving improved constraints on all cosmological parameters including α. We note that this analysis is still somewhat idealized since in a typical WL survey a given angular scale receives contributions from a projection over a wide range of scales, and their relative contribution changes with angle. Our analysis is most directly applicable to a tomographic analysis of WL data, where one isolates the lensing effects from a localized volume so that projection effects are minimized. 4.2 Convergence of the covariance matrix and the required volume of simulations As we argued the non-Gaussian covariance is dominated by a single eigenmode, whose k dependence is fixed, and only its amplitude α varies from a realization to a realization. By the central limit theorem, we expect the distribution of αs to be close to a Gaussian. This is shown in Fig. 12, normalized to its variance. We see that the distribution is indeed very close to a Gaussian. In the context of this model, we can ask how many simulations are needed to fully determine the covariance matrix. To determine the k dependence of the eigenvector d1(k), we just need two simulations: the difference in P(ki) between the two simulations will be proportional to d1(ki)P(ki). The proportionality factor is not relevant if we use the method described in the previous subsection, where k dependence is used to determine the amplitude and remove the non-Gaussian part of the covariance matrix. The right-hand panel of Fig. 13 shows the relative deviations of the two power spectra from two simulations for few different pair of simulations. All curves look similar except for their normalization and agree with the shape of the eigenvector itself, but the realization noise at lower k is significant. If we wish to determine the full covariance matrix including its normalization, we also need to determine 〈α2〉. Since α has a Gaussian distribution, we can obtain a fractional error of (2/N)1/2 on its variance with N simulations. While Gaussian distribution is the best case scenario, it is still a relatively slow convergence: if one wishes to determine covariance matrix to 10 per cent one needs 200 simulations, while 1 per cent requires 20 000 simulations. As all our results are shown as Covij/(PiPj), any resolution issue for different simulations should cancel out. However, we can also exploit the scalings of individual terms of covariance matrix with volume. The Gaussian term is analytic and the SSC term can be well determined by a response to a constant curvature (separate universe methods), and the value of |$\sigma _{\rm b}^2$| on the scale of the survey given by equation (11). Note that this term does not scale simply as 1/V, it also depends on the power spectrum P(q). We are thus left with the modes inside the volume of the survey. Equation (7) states that the covariance is given by the trispectrum divided by the volume V. One can compute the trispectrum using a small volume simulations and scale by the survey volume. But one must include all of the trispectrum contributions. The mode contribution to the simplified version of the 1-loop integral is shown in Fig. 3. Modes with q < 0.02 h Mpc−1 contribute about 1 per cent of the total integrand. If we are willing to tolerate 1 per cent error, then the required box size to get all the contributing modes is V ∼ (300 h−1 Mpc)3. In general, we, therefore, expect that the minimal volume per simulation is roughly of this size. If the desired error of the covariance is 10 per cent, then we can run 200 simulations of V ∼ (300 h−1 Mpc)3, i.e. a total volume of V ∼ 5 (h−1 Gpc)3. If the desired error is 1 per cent, then the required volume is V ∼ 500 (h−1 Gpc)3. To reduce the disconnected part of the covariance matrix, which is effectively noise for the connected part, one can simply take broader bins of k. As shown in Fig. 1, the covariance is expected to be very insensitive to the binning scheme. Note that the effective volume of future surveys is of the order of tens to hundreds of (h−1 Gpc)3, so the required volumes are of the order of a single survey volume. But the simulations still need to be done on periodic boxes, to eliminate the SSC term, which would otherwise dominate the covariance of sub-boxes (although see below for an alternative method where this is not needed). SSC term is added separately to the total covariance matrix, using the correct value of |$\sigma _{\rm b}^2$| from equation (11). To test this, we compare the two covariance matrices at our disposal, which have been simulated at different volumes. We rescale them by the volume ratio, which is (650 h−1 Mpc)3 for B15 and (500 h−1 Mpc)3 for L14, a factor of 2.2. We see in Fig. 17 that this gives a very good agreement between the two. There are residual differences at the level of 10 per cent. At low k these are probably due to noise fluctuations. It is unclear if the remaining differences are due to the differences in the cosmological model, or due to some differences in the simulations. Figure 17. Open in new tabDownload slide Comparing L14 covariance (solid black) to that of B15 (red dashed), where B15 covariance has been renormalized to the same volume as L14, i.e. by a factor (650/500)3. Unrescaled B15 is shown as red solid. Figure 17. Open in new tabDownload slide Comparing L14 covariance (solid black) to that of B15 (red dashed), where B15 covariance has been renormalized to the same volume as L14, i.e. by a factor (650/500)3. Unrescaled B15 is shown as red solid. Finally, we mention that one may be able to get the covariance matrix simply from the survey data itself. As we argued above, the volumes of typical future surveys are tens of (h−1 Gpc)3, which is not too different from the volume needed to reliably simulate covariance matrix. One can divide a survey into many subvolumes of a typical size of (300 h−1 Mpc)3 and measure the power spectrum in each subvolume. To reduce the fluctuations for the disconnected part one should use relatively broad k bins, and remove the disconnected part (including any possible window function effects that induce mixing between the bins). The main remaining issue is that dividing an observed volume into subvolumes also induces SSC term on the subvolume (same is true for jackknife or bootstrap methods that use a similar subdivision). This requires thus to estimate and remove the SSC term from the subvolumes. This can be done by correlating the spatially dependent subvolume power spectrum with the mean density on the subvolume, which is a form of a squeezed limit bispectrum (Chiang et al. 2014). Once the mean SSC response is determined, one can multiply it with the subvolume mean density and subtract it out of the subvolume power spectrum, thereby removing SSC term from the subvolumes. From the remaining power spectra, one can now compute the covariance matrix whose connected term will be dominated by the modes inside the subvolume, and which is expected to scale as 1/Vsubvolume. Once the connected part of the covariance matrix is determined, one can scale it as Vsubvolume/V to get the corresponding covariance of the whole survey from modes inside the subvolume. To this, one needs to add SSC term by computing |$\sigma _{\rm b}^2$| from equation (11) using the window function defined by the whole survey, multiplied by the response as derived from the squeezed limit bispectrum. Finally, one adds the disconnected part of the covariance matrix, including the survey window function effects. It is a straightforward technique in order to estimate the covariance directly from data, however, yet to be verified. We leave this analysis to future work. 4.3 Connecting perturbation theory to the halo model The halo model has been a very successful model for non-linear clustering of the dark matter power spectrum (see e.g. Cooray & Sheth 2002 for a review). It decomposes the non-linear power spectrum into the correlations between haloes, the so-called two-halo term, and correlations within the haloes, the so-called one-halo term. The two-halo term has been originally modelled as a linear power spectrum smoothed on the virial radius scale of the haloes. This does not capture effects like the smoothing of baryonic acoustic peaks. A simple generalization has been presented in Seljak & Vlah (2015), where the two-halo term is instead replaced by a Zel'dovich approximation (ZA) power spectrum. The difference between the full power spectrum and ZA is defined to be the one-halo term. This difference has to vanish for k → 0, meaning that the one-halo term needs to be compensated and is not a constant, which happens at very low k. The one-halo term on large scales is given by Poisson fluctuations at zero lag and can be calculated as the second moment of mass integrated over the halo mass function. An alternative description pursued here is to use perturbation theory. PT has been applied to the power spectrum and at 1-loop order, it is expected to do well only at very low k, while at higher k it requires modest EFT corrections, which make 1-loop sufficiently accurate up to k ∼ 0.2 h Mpc−1 (Baldauf, Schaan & Zaldarriaga 2016a; Foreman, Perrier & Senatore 2016; Vlah et al. 2016). In Seljak & Vlah (2015), it was argued that we can define a regime where both halo model and 1-loop PT (with EFT) are equally valid so that one can calibrate the one-halo term on PT, and then use the halo model to extend it to higher k. This allows a smooth continuation from PT to a high k regime where the halo model is valid. When compared to the halo model predictions in terms of the second moment integrated over the halo mass function, this requires a halo mass definition that is somewhat higher than the standard virial mass definition M200. A similar discussion applies to the covariance matrix. In the halo model, the covariance at high k arises from Poisson fluctuations. A particularly simple form has been developed in Mohammed & Seljak (2014), where the connected term has been derived as \begin{equation} {\rm Cov}_{ij}^{\rm connected}=P(k_i)P(k_j)V^{-1}\delta _{A_0}^2. \end{equation} (35) In the halo model, |$\delta _{A_0}^2$| can be related to the fourth moment of mass integrated over the halo mass function. This calculation somewhat overpredicts the true value when compared to simulations (Mohammed & Seljak 2014), but it is reasonably close. One can see that the form of the halo model prediction is identical to PT calculation in the limit where there is a single eigenvector that does not depend on k, which is approximately valid, as seen in Fig. 13. In PT, the variance is given by the value of S in equation (28). This is replaced by |$\delta _{A_0}^2$| in the halo model. If we insist that the two descriptions agree with each other in the regime of overlap (k ∼ 0.2 h Mpc−1), then this determines the value of |$\delta _{A_0}^2$| from PT. Hence, the fourth moment of the halo mass function is determined by equation (28), just as the second moment of halo mass function is determined by 1-loop PT (Seljak & Vlah 2015). This is a consistency that restricts the form of the halo mass function. Note that the agreements are not perfect: while for the power spectrum, the standard one-halo term calculation underpredicts PT (unless the halo mass is increased), it slightly overpredicts it for the covariance, but these differences could also be caused by uncertainties in the halo mass function and halo mass definitions. This picture allows a convenient connection of PT to the halo model on small scales. For covariance, we expect the small scale clustering to be dominated by one-halo term from lower mass haloes, which are less rare, and hence we expect |${\rm Cov}_{ij}^{\rm connected}/P(k_i)P(k_j)$| to decrease as we go to higher k. Ultimately, this agreement is just a justification for the halo model being applicable to both the power spectrum and its covariance. While the one halo term has to be the dominant term at high k in the power spectrum, essentially by definition of the halo mass function, for higher order moments the halo model needs to be justified case by case. We have seen that it works well for covariance. For variance of covariance, the halo model prediction is given by the eighth moment of the mass integrated over the halo mass function, and predicts very large relative fluctuations, of order 30 per cent for a volume of 1 (Gpc h−1)3 (Mohammed & Seljak 2014). The corresponding 1-loop PT calculation is given by 112π3∫P(q)4d3q/[∫P(q)2d3q]2/V and it gives two orders of magnitude lower variance of covariance, suggesting Poisson fluctuations are not very important, as also explicitly verified by simulations in Fig. 12 (note that this figure includes SSC term). This suggests that the halo model cannot be reliably applied to the calculations of the eighth moment. 5 CONCLUSIONS In this paper, we developed a perturbative model for the covariance matrix of the power spectrum, going up to 1-loop in PT, but without including all the 1-loop terms (see Bertolini et al. 2016a for the full 1-loop calculation). In addition, we go beyond 1-loop by including the non-linear damping of response functions (Nishimichi et al. 2014), which improves our results at higher k. Overall, our approach predicts the results from simulations to about 10 per cent accuracy in the quasi-linear regime. The largest contribution to the covariance is SSC, which arises due to the coupling of the modes to the ones which are larger than the survey scale and which has been extensively studied previously (Li et al. 2014a,b). This effect comes from the modes outside the survey and cannot be captured by the jackknife or bootstrap methods, which subdivide the full volume into many smaller volumes. In fact, these methods overestimate the covariance terms because they generate SSC from the smaller regions. We propose an alternative approach which removes this problem. On large scales, the second largest contribution comes from tree-level terms from modes inside the survey, while on smaller scales 1-loop terms dominate. The dominant 1-loop term, and the only one we include, comes from the sampling variance fluctuations of large-scale modes within the survey volume, which induce a coherent response by small-scale modes which we model using 1-loop power spectrum. When we allow for a damped non-linear response to modes that are highly non-linear (Nishimichi et al. 2014), we obtain some additional improvements, especially at lower redshifts. Damping does not affect results at higher redshifts because all of the modes that have significant sampling fluctuations are linear. We explore the structure of the covariance matrix and find that its non-Gaussian part is strongly dominated by a single eigenmode. This suggests that the non-Gaussian response has always the same shape, only its amplitude varies. We analyse the probability distribution of this amplitude and find that it is close to a Gaussian. Thus, the convergence rate of covariance matrix simulations scales as a Gaussian, with (2/N)1/2 giving the relative error after N simulations. One possible alternative approach is to ignore the non-Gaussian covariance in the analysis and instead include the eigenvector response as a fictitious external parameter in the analysis. This parameter can, in principle, be determined from the data itself, but it is quite degenerate with other cosmological parameters. While the analysis given in this paper provides several important insights into the nature of the covariance matrix of the two-point correlators, our results cannot be directly applied to the data. For WL observations, which are measuring the dark matter correlations, one needs to perform the projection along the line of sight. This leads to a decorrelation of the non-Gaussian covariance because different l correspond to different effective redshift, and hence different effective volume. This differs from our analysis where all modes are maximally correlated because they are coupled to the same long wavelength modes within the given volume. So far, this has only been addressed in the context of the halo or SSC model (Schaan, Takada & Spergel 2014; Krause & Eifler 2016). The 3D analysis performed here is more likely to be of immediate application to the galaxy clustering, but this would require adding biasing and redshift space distortions to the model, and it is unclear whether successful analytic models can be built. Nevertheless, the decomposition of the covariance matrix into the three components, disconnected, connected but generated from modes outside the survey and connected and generated from modes inside the survey, allows one, in principle, to build the full covariance matrix from relatively small simulation volumes, and possibly even from the data itself, by subdividing into subvolumes and determining each of the three components separately. It is worth pursuing this further to determine the optimal approach that delivers the highest accuracy with this technique. We thank Matias Zaldarriaga for useful discussions. US acknowledges support from NASA grant NNX15AL17G. ZV is supported in part by the U.S. Department of Energy contract to SLAC no. DE-AC02-76SF00515. IM is supported by Fermi Research Alliance, LLC under Contract No. De-AC02-07CH11359 with the United States Department of Energy. 1 " http://www.darkenergysurvey.org 2 " http://desi.lbl.gov 3 " https://www.lsst.org 4 " http://www.euclid-ec.org REFERENCES Baldauf T. , Seljak U., Senatore L., Zaldarriaga M., 2011 , J. Cosmol. Astropart. Phys. , 10 , 031 Crossref Search ADS Baldauf T. , Schaan E., Zaldarriaga M., 2016a , J. Cosmol. Astropart. Phys. , 3 , 007 Baldauf T. , Seljak U., Senatore L., Zaldarriaga M., 2016b , J. Cosmol. Astropart. Phys. , 09 , 007 Baugh C. M. , Efstathiou G., 1994 , MNRAS , 270 , 183 Crossref Search ADS Bernardeau F. , Colombi S., Gaztañaga E., Scoccimarro R., 2002 , Phys. Rep. , 367 , 1 Crossref Search ADS Bertolini D. , Schutz K., Solon M. P., Walsh J. R., Zurek K. M., 2016a , Phys. Rev. D , 93 , 123505 Crossref Search ADS Bertolini D. , Schutz K., Solon M. P., Zurek K. M., 2016b , J. Cosmol. Astropart. Phys. , 06 , 052 Crossref Search ADS Blot L. , Corasaniti P. S., Alimi J.-M., Reverdy V., Rasera Y., 2015 , MNRAS , 446 , 1756 (B15 ) Crossref Search ADS Blot L. , Corasaniti P. S., Amendola L., Kitching T. D., 2016 , MNRAS , 458 , 4462 Crossref Search ADS Chiang C.-T. , Wagner C., Schmidt F., Komatsu E., 2014 , J. Cosmol. Astropart. Phys. , 5 , 048 Crossref Search ADS Coles P. , 1990 , MNRAS , 243 , 171 Crossref Search ADS Cooray A. , Hu W., 2001 , ApJ , 554 , 56 Crossref Search ADS Cooray A. , Sheth R., 2002 , Phys. Rep. , 372 , 1 Crossref Search ADS de Putter R. , Wagner C., Mena O., Verde L., Percival W. J., 2012 , J. Cosmol. Astropart. Phys. , 4 , 019 Crossref Search ADS Dodelson S. , Schneider M. D., 2013 , Phys. Rev. D , 88 , 063537 Crossref Search ADS Foreman S. , Perrier H., Senatore L., 2016 , J. Cosmol. Astropart. Phys. , 5 , 027 Crossref Search ADS Gnedin N. Y. , Kravtsov A. V., Rudd D. H., 2011 , ApJS , 194 , 46 Crossref Search ADS Hamilton A. J. S. , Rimes C. D., Scoccimarro R., 2006 , MNRAS , 371 , 1188 Crossref Search ADS Harnois-Déraps J. , Pen U.-L., 2012 , MNRAS , 423 , 2288 Crossref Search ADS Heitmann K. , Higdon D., White M., Habib S., Williams B. J., Lawrence E., Wagner C., 2009 , ApJ , 705 , 156 Crossref Search ADS Heitmann K. , White M., Wagner C., Habib S., Higdon D., 2010 , ApJ , 715 , 104 Crossref Search ADS Hu W. , White M., 2001 , ApJ , 554 , 67 Crossref Search ADS Jain B. , Bertschinger E., 1994 , ApJ , 431 , 495 Crossref Search ADS Juszkiewicz R. , 1981 , MNRAS , 197 , 931 Crossref Search ADS Juszkiewicz R. , Sonoda D. H., Barrow J. D., 1984 , MNRAS , 209 , 139 Crossref Search ADS Kayo I. , Takada M., Jain B., 2013 , MNRAS , 429 , 344 Crossref Search ADS Krause E. , Eifler T., 2016 , preprint (arXiv:1601.05779) Lawrence E. , Heitmann K., White M., Higdon D., Wagner C., Habib S., Williams B., 2010 , ApJ , 713 , 1322 Crossref Search ADS Li Y. , Hu W., Takada M., 2014a , Phys. Rev. D , 89 , 083519 (L14 ) Crossref Search ADS Li Y. , Hu W., Takada M., 2014b , Phys. Rev. D , 90 , 103530 Crossref Search ADS Lokas E. L. , Juszkiewicz R., Bouchet F. R., Hivon E., 1996 , ApJ , 467 , 1 Crossref Search ADS Makino N. , Sasaki M., Suto Y., 1992 , Phys. Rev. D , 46 , 585 Crossref Search ADS Mohammed I. , Seljak U., 2014 , MNRAS , 445 , 3382 Crossref Search ADS Mohammed I. , Martizzi D., Teyssier R., Amara A., 2014 , MNRAS , preprint (arXiv:1410.6826) Nishimichi T. , Bernardeau F., Taruya A., 2014 , preprint (arXiv:1411.2970) Sato M. , Hamana T., Takahashi R., Takada M., Yoshida N., Matsubara T., Sugiyama N., 2009 , ApJ , 701 , 945 Crossref Search ADS Sato M. , Takada M., Hamana T., Matsubara T., 2011 , ApJ , 734 , 76 Crossref Search ADS Schaan E. , Takada M., Spergel D. N., 2014 , Phys. Rev. D , 90 , 123523 Crossref Search ADS Scoccimarro R. , Frieman J. A., 1996 , ApJ , 473 , 620 Crossref Search ADS Scoccimarro R. , Zaldarriaga M., Hui L., 1999 , ApJ , 527 , 1 Crossref Search ADS Sefusatti E. , Crocce M., Pueblas S., Scoccimarro R., 2006 , Phys. Rev. D , 74 , 023522 Crossref Search ADS Seljak U. , Vlah Z., 2015 , Phys. Rev. D , 91 , 123516 Crossref Search ADS Semboloni E. , Hoekstra H., Schaye J., van Daalen M. P., McCarthy I. G., 2011 , MNRAS , 417 , 2020 Crossref Search ADS Semboloni E. , Hoekstra H., Schaye J., 2013 , MNRAS , 434 , 148 Crossref Search ADS Sherwin B. D. , Zaldarriaga M., 2012 , Phys. Rev. D , 85 , 103523 Crossref Search ADS Sirko E. , 2005 , ApJ , 634 , 728 Crossref Search ADS Smith R. E. et al. , 2003 , MNRAS , 341 , 1311 Crossref Search ADS Suto Y. , Sasaki M., 1991 , Phys. Rev. Lett. , 66 , 264 Crossref Search ADS PubMed Takada M. , Bridle S., 2007 , New J. Phys. , 9 , 446 Crossref Search ADS Takada M. , Hu W., 2013 , Phys. Rev. D , 87 , 123504 Crossref Search ADS Takada M. , Jain B., 2009 , MNRAS , 395 , 2065 Crossref Search ADS Takahashi R. et al. , 2008 , MNRAS , 389 , 1675 Crossref Search ADS Takahashi R. et al. , 2009 , ApJ , 700 , 479 Crossref Search ADS Takahashi R. , Sato M., Nishimichi T., Taruya A., Oguri M., 2012 , ApJ , 761 , 152 Crossref Search ADS van Daalen M. P. , Schaye J., Booth C. M., Dalla Vecchia C., 2011 , MNRAS , 415 , 3649 Crossref Search ADS Vishniac E. T. , 1983 , MNRAS , 203 , 345 Crossref Search ADS Vlah Z. , Seljak U., Yat Chu M., Feng Y., 2016 , J. Cosmol. Astropart. Phys. , 3 , 057 Crossref Search ADS Wagner C. , Schmidt F., Chiang C.-T., Komatsu E., 2015 , MNRAS , 448 , L11 Crossref Search ADS © 2016 The Authors Published by Oxford University Press on behalf of the Royal Astronomical Society
TI - Perturbative approach to covariance matrix of the matter power spectrum
JF - Monthly Notices of the Royal Astronomical Society
DO - 10.1093/mnras/stw3196
DA - 2017-04-01
UR - https://www.deepdyve.com/lp/oxford-university-press/perturbative-approach-to-covariance-matrix-of-the-matter-power-aUW8FZniBb
SP - 780
VL - 466
IS - 1
DP - DeepDyve
ER -