Cosmological Model Parameter Dependence of the Matter Power Spectrum Covariance from the DEUS-PUR $Cosmo$ Simulations

Future galaxy surveys will provide accurate measurements of the matter power spectrum across an unprecedented range of scales and redshifts. The analysis of these data will require to accurately model the imprint of non-linearities on the matter density field, which induces a non-Gaussian contribution to the data covariance. As the imprint of non-linearities is cosmology dependent, a further complication arises from accounting for the cosmological dependence of the non-Gaussian part of the covariance. Here, we study this using a dedicated suite of N-body simulations, the Dark Energy Universe Simulation - Parallel Universe Runs (DEUS-PUR) $Cosmo$. These consist of 512 realizations for 10 different cosmologies where we vary the matter density $\Omega_m$, the amplitude of density fluctuations $\sigma_8$, the reduced Hubble parameter $h$ and a constant dark energy equation of state $w$ by approximately $10\%$. We use these data to evaluate the first and second derivatives of the power spectrum covariance with respect to a fiducial $\Lambda$CDM cosmology. We find that the variations can be as large as $150\%$ depending on the scale, redshift and model parameter considered. Using a Fisher matrix approach, we evaluate the impact of using a covariance estimated at a fiducial model rather than the true underlying cosmology. We find that the estimated $1\sigma$ errors are affected at approximately $5\%$, $20\%$, $50\%$ and $120\%$ level when assuming non-fiducial values of $h$, $w$, $\Omega_m$ and $\sigma_8$ respectively. These results suggest that the use of cosmology-dependent covariances is key for precision cosmology.


INTRODUCTION
The upcoming generation of galaxy surveys will provide accurate measurements of the clustering of matter across an unprecedented range of scales and redshifts (e.g. LSST Science Collaboration et al. 2009;Laureijs et al. 2011;DESI Collaboration et al. 2016;Akeson et al. 2019). Precise estimates of the matter power spectrum from measurements of the spatial distribution of galaxies and the weak gravitational lensing shear will enable to test models beyond the standard ΛCDM scenario and investigate the nature of dark energy. These datasets will be sensitive to the imprints of the non-linear regime of gravitational collapse of matter, as E-mail: lblot@mpa-garching.mpg.de such they need to be accurately modelled if one aims to infer unbiased cosmological parameter constraints.
In the past years, this has sparked a major theoretical and numerical effort to provide accurate predictions of galaxy clustering observables and the associated data covariances on quasi-linear and non-linear scales. On large scales the matter density field is Gaussian. Consequently, the matter power spectrum covariance has a diagonal structure and is simply proportional to the square of the power spectrum itself. Finite survey volume effects induce a non-Gaussian contribution also known as super-sample covariance (see e.g. Hamilton et al. 2006;Takada & Hu 2013), while the nonlinearities of the matter density field that develop at small scales induce mode correlations that further contribute to the non-Gaussian structure of the covariance, making the matrix non-diagonal (Meiksin & White 1999;Scoccimarro et al. 1999). Analytical approaches to estimate these effects have been developed in a vast literature (see e.g. Mohammed & Seljak 2014;Bertolini et al. 2016; Barreira & Schmidt 2017). Nevertheless, N-body simulations remain the primary tool to investigate the imprint of non-linearities, while providing the necessary benchmark to test the validity of analytical model predictions (see e.g. Takahashi et al. 2009;Kiessling et al. 2011;Blot et al. 2015;Villaescusa-Navarro et al. 2019). Estimating the covariance with the level of accuracy that is required to correctly analyse future galaxy survey data demands sampling the matter power spectrum from a very large suite of N-body simulations. As an example, in Blot et al. (2015) we have estimated the covariance using ∼ 10 4 independent N-body simulations and shown that non-linearities induce significant deviations from the Gaussian prediction on modes k 0.25 h Mpc −1 and at redshift z < 0.5. Furthermore, by taking advantage of the large simulation suite, Blot et al. (2016) have shown that more than > 5000 realisation are necessary to reduce the impact of sample covariance errors on the estimated cosmological parameter uncertainties to sub-percent level.
In these studies the power spectrum covariance has been evaluated for a fixed fiducial cosmological model. However, the imprint of non-linearities on the matter power spectrum is cosmology dependent (see e.g. Ma et al. 1999;Casarini et al. 2009;Alimi et al. 2010). Hence, it is natural to expect that such dependence extends to the non-Gaussian part of the matter power spectrum covariance. Neglecting the variation of the covariance with the cosmological model parameters can introduce spurious systematic errors in the parameter inference analysis. This has been investigated in the past in the context of weak lensing shear power spectrum measurements especially in relation to the super-sample covariance (Eifler et al. 2009;Labatie et al. 2012;Harnois-Déraps et al. 2019) and several methodologies have been developed to extrapolate the cosmological dependence from a finite set of simulations (see e.g. Morrison & Schneider 2013;White & Padmanabhan 2015;Reischke et al. 2017).
Here, we aim to specifically investigate the cosmological dependence of the matter power spectrum covariance due to small scale non-linearities. We will make use of a large ensemble of N-body simulations for several cosmological parameter configurations to compute the first-and second-order derivatives of the power spectrum covariance. Our intent is to determine the amplitude of such derivatives and perform a preliminary evaluation of their impact on cosmological parameter inference through a forecast analysis. Moreover, to facilitate further progress in the analytical modelling of the cosmological dependence of the power spectrum covariance, we have made publicly available the numerical simulation data used in the study presented here.
The paper is organized as follows. In Section 2 we describe the simulations set and the covariance estimator. In Section 3 we present our results on the cosmological dependence of the covariance and its impact on cosmological parameter inference analyses. In Section 4 we present our conclusions.

Numerical Codes & Simulation Pipeline
Building upon the automated pipeline developed for the Dark Energy Universe Simulations -Parallel Universe Runs (DEUS-PUR) project (Blot et al. 2015), we have realized a large suite of N-body simulations for different sets of cosmological parameters. We refer to this simulation suite as DEUS-PUR Cosmo.
In the following, we will briefly describe the simulation pipeline and we refer the interested readers to Blot et al. (2015) for a more detailed description. The cosmological parameters for the different runs are provided by the user through namelist files, while the sequential call to the various codes, from the computation of the input tables containing the cosmological functions to the post-processing of the simulations, is entirely automatised. For a given cosmological model the first step consists in computing the linear matter power spectrum using the code CAMB (Lewis et al. 2000) and solving the Friedmann equations using a dedicated code called NewDarkCosmos. The respective output tables are input to the code generating the initial conditions in the former case and the N-body solver in the latter case. Then Gaussian initial conditions are generated with an optimized version of the code MPGRAFIC (Prunet et al. 2008). The simulations are run using an improved version of the Adaptive Mesh Refinement (AMR) N-body code RAM-SES (Teyssier 2002), which uses a multigrid Poisson solver (Guillet & Teyssier 2011). Finally, halos are detected with the halo finder code pFoF (Roy et al. 2014), which is based on the friends-of-friends algorithm, while power spectra are computed using an optimized version of the code POWER-GRID (Prunet et al. 2008), which uses a Fast Fourier Transform (FFT) algorithm. The matter density field is estimated on a cartesian grid (twice thinner than the coarse AMR grid) with a Cloud-in-Cell mass assignment scheme. To minimize the effect of aliasing, we exclude all modes beyond half the Nyquist frequency of the density grid.

Cosmological Models & Simulation Characteristics
For each of the 10 cosmological models listed in Table 1 we have run a set of 512 cosmological N-body simulations that share the same initial phases across different cosmologies such as to reduce the numerical noise in the computation of the derivatives of the matter power spectrum and its covariance. Each simulation consists of a cubic box of length size L box = 328.125 h −1 Mpc with 512 3 particles (corresponding to a particle mass mp = 1.88 × 10 10 M /h for the fiducial cosmology). The cosmological parameters have been chosen such that we vary one parameter of interest at a time in a symmetric way with respect to the fiducial value. This allows us to obtain a very accurate estimate of the first and second derivative of any quantity in the vicinity of the reference cosmological model. Our fiducial cosmology corresponds to model 2 in Table 1, which is a flat ΛCDM model calibrated on WMAP-7 year data (Komatsu et al. 2011). We have set the baryon density Ω b h 2 = 0.02258 and the scalar spectral index ns = 0.963 consistently with the values of the WMAP-7 cosmological analysis, while the amount of radiation and  Figure 1. Variance of the matter power spectrum as a function of wave number for different cosmologies normalized to the linear theory (top panels) and the Gaussian expectation (bottom panels). The continuous line indicates the fiducial cosmology case while the shaded area represents the variation when the parameter indicated in the title of each panel is varied. Colours from dark brown to yellow corresponds to decreasing redshifts: z = 2, z = 1, z = 0.5 and z = 0. As we can see non-linearities play an important role for k > 0.1 − 0.2 h Mpc −1 and are cosmology-dependent. year data, while all other models differ for a variation of one of the parameter values (in bold). Models 5-6 are characterized by a ±13% variation of σ 8 with respect to the fiducial value, models 7-8 by a ±20% variation of Ωm and models 9-10 by a ±7% variation of h. Ωm, while models 9 and 10 are associated with variations of the reduced Hubble constant h. We have also run as test case a set of simulations for model 1, which is an Einstein-de Sitter model (Eds) without dark energy.
It is worth remarking that the parameter variations considered here are of the order of the parameter uncertainties from state-of-the-art cosmological experiments (Planck Collaboration et al. 2018;Abbott et al. 2018). Moreover, there is mounting evidence of tensions between the values of a number of cosmological parameters inferred from different probes within the ΛCDM model. This is indeed the case of the Hubble parameter H0 as measured from low redshift probes which is found to be in ∼ 5σ tension with the one inferred from early universe probes such as the Cosmic Microwave Background (CMB) anisotropy power spectra from Planck (see e.g. Verde et al. 2019, and references therein), resulting in a variation of the best-fit value of the order of 10%. Secondly, the parameter combination S8 = σ8(Ωm/0.3) 0.5 measured from weak lensing probes is consistently lower than the value inferred from the Planck data (Hildebrandt et al. 2020;Abbott et al. 2018;Hikage et al. 2019). Thus, the interest of exploring the cosmological dependence of the matter power spectrum covariance on a similar range of parameter values.

Covariance Estimator and Parameter Derivatives
We compute the matter power spectrum P (k) of each realization in band powers of size ∆k = 2π/L box . We evaluate the covariance between two different modes using the unbiased sample covariance estimator: where Ns is the number of realization, Pi(k) is the matter power spectrum of the i−th realization andP (k) = Ns i=1 Pi(k)/Ns is the sample mean.
We estimate the first and second derivatives of the power spectrum covariance with respect to the cosmological parameters for each pair of modes using the finite difference approximation: whereθ is the fiducial cosmological parameter value and ∆θ the finite variation of its value.

Variance of the Matter Power Spectrum
We evaluate the variance of the matter power spectrum (i.e. the diagonal part of the covariance matrix). This is shown in Figure 1 for the different cosmological models (panel left to right) and redshifts (lines from yellow to dark brown). The top panels show the variance normalized to the linear prediction, where N modes is the number of modes in the bin ∆k and P lin (k) is the linear matter power spectrum; the lower panels show the variance normalized to the Gaussian prediction whereP is the average non-linear matter power spectrum from the 512 independent N-body realizations. We may notice that at small wavenumbers (k < 0.1 h Mpc −1 ) the estimated variance is consistent with the linear Gaussian prediction, which validates the results of our simulations in this regime. On the other hand at larger wavenumbers, we observe a strong departure that increases as function of the wavenumber and for decreasing redshifts. The amplitude of this deviation reaches up to a factor of 10 4 greater than the linear prediction (top panels) at k ≈ 3 h Mpc −1 at z = 0, and a factor 25 with respect to the Gaussian case (bottom panels). In the first case this is partly due to the fact that the linear power spectrum significantly underestimates P (k) at these scales and redshifts, while in the latter case this is due to the well-known nongaussian contribution from the non-linear regime of matter clustering (see e.g. Blot et al. 2015, and references therein). In all cases, we can see that the amplitude and slope of the departures from the linear and Gaussian expectations depend on the cosmological parameters in a non-trivial way. This motivates a detailed study of the cosmological dependence of the covariance which we discuss next.

Matter Power Spectrum Covariance Derivatives
We evaluate the first and second derivatives of the power spectrum covariance with respect to the cosmological parameters which we plot in the top and bottom panels of Fig. 2 respectively. Panels from left to right show the redshift evolution at z = 2, 1, 0.5 and 0. The intensity mapping is set such that positive (negative) derivatives are shown in red (blue), while vanishing matrix elements are shown in white. Panels from top to bottom correspond to derivatives with respect to Ωm, σ8, h and w respectively. First, we may notice that the sign of the derivatives follows from that of the matter power spectrum. As an example, the first derivative of the covariance with respect to Ωm is negative. This is because at a fixed value of σ8, a positive variation of Ωm decreases the overall amplitude of the matter power spectrum. Hence the covariance decreases, which results in a negative derivative. Conversely, a positive variation of σ8 at a fixed Ωm value increases the overall amplitude of matter power spectrum, thus resulting in a positive derivative. We can also see that the first-order derivative of the covariance increases in absolute value from high to low redshift. Moreover, at a given redshift the largest elements are those corresponding to pairs of modes consisting of a large mode coupled to a small one, which is indicative of the onset of non-linearities that grow at lower redshifts while propagating to larger scales. This leads to a characteristic off-diagonal structure of the first-order derivative of the covariance that similar to that of the covariance itself (see e.g. Fig. 3 in Blot et al. 2015). It is worth noticing that the first-order derivative of the covariance is larger for Ωm and σ8 and smallest for w, which follows from the dependence of the matter power spectrum on these parameters. We observe a similar trend in the case of the second-order derivatives, shown in Fig. 2, which all have positive values except for the case of h.
We can use these derivatives to infer an understanding of the dependence of the matter power spectrum covariance on the cosmological parameters. In particular, we can consider a Taylor expansion of the covariance up to second-order around the fiducial cosmology: the validity of this approximation depends on the expansion coefficients, i.e. the covariance derivatives normalised to the fiducial covariance, to be O(1). We show these ratios in Fig. 3 for the first-order (top panel) and second-order (bottom panel) terms respectively. In the first-order case we can see that the largest matrix elements are less than unity for Ωm, h and w, though still large enough to cause a slow convergence of the Taylor expansion along these parameter directions. Instead, in the case of σ8 the largest matrix elements exceed unity even on quasi-linear scales corresponding to modes k 1 h Mpc −1 . This suggests that the dependence of the covariance on σ8 may be highly non-linear. Using the same colour coding, we can see that the second-order terms are much smaller than the first-order ones. Moreover, most of the second-order contributions are smaller than unity, meaning that an important part of the information about the cosmological dependence of the covariance is encoded in the first two derivatives. In any case, it is striking that a variation of order 10% of the cosmological parameters induces a change of the covariance matrices between 10% and 150% depending on cosmology, redshift and scale.

Cosmological Parameter Inference Forecast
In order to assess the impact of the cosmological parameter dependence of the covariance on the cosmological parameter inference we perform a simple Fisher matrix analysis.
In principle, given the non-Gaussian structure that such dependence of the covariance may induce on the cosmological parameter posterior, a more rigorous approach would be to perform a Markov Chain Monte Carlo analysis of a set of synthetic matter power spectra for our fiducial cosmology and let the covariance vary along the sampling of the cosmological parameters. However, given the limited number of parameter configurations for which we have evaluated the covariance and the potentially non-linear nature of its parameter dependence, we are unable to perform such an analysis. Moreover, even though this would be the correct Bayesian way to infer model parameter constraints from datasets with parameter dependent covariances (see e.g. Ma et al. 2016), this is not how current matter power spectrum analyses are performed, since the data covariance is fixed at a fiducial cosmology. If the fiducial model is far from the true cosmology, then the constraints can be biased and the errors misestimated. Therefore, to evaluate the extent to  which this affects the cosmological parameter inference we compute the Fisher matrix: where we have considered the parameter vector θ = {Ωm, σ8, w, h}. We compute the derivative of the matter power spectrum for the fiducial cosmology specified by the vector of valuesθ = {0.2573, 0.801, −1, 0.72} with the finite difference approximation using the spectra from the simulations, while we use the covariance matrix of the simulated cosmologies specified by the vector of values θ * from Table 1. We evaluate Eq. (5) assuming 15 uncorrelated redshift bins in the range 0.15 < z < 1.75 corresponding to the redshift of the snapshots of our simulation suite. In Figure 4 we show two examples of the 1 and 2σ contours inferred from the evaluation of the Fisher matrix assuming kmax = 1 h Mpc −1 and obtained using the covariance evaluated for different non-fiducial values of σ8 (left panel) and w (right panel). The contours inferred with the covariance evaluated at the fiducial cosmology are shown as solid line. We can see that with respect to this case, setting the covariance to a different cosmological model leads to a modification of the area within the confidence regions as well as the angle of the degeneracies between different pairs of parameters.
In Figure 5 we show the variation of the area (left panel) and the angle (right panel) of the 1σ contours for different combination of parameters as function of kmax when using a covariance computed for the non-fiducial cosmological models from Table 1 with different values of Ωm (blue lines), h (orange lines), σ8 (green lines) and w (pink lines). The thick solid lines mark the largest deviations from the case with fiducial covariance. We can see that the largest deviation of the contour area occurs in the case of the covariance being computed for non-fiducial values of σ8 and Ωm, while difference are smaller for h and w. Quite importantly, deviations are of the order of 50% on quasi-linear scales corresponding to kmax ∼ 0.1 − 0.2 h Mpc −1 , which are already probed by current galaxy surveys (see e.g. Beutler et al. 2017). The effect on the angle of the parameter degeneracies is smaller with maximal deviations not exceeding the 5% level up to kmax ∼ 2 h Mpc −1 . For higher kmax the largest impact is associated with σ8, while the effect remains smaller for the other parameters.
In Figure 6 we show the variation of the 1σ parameter errors as a function of kmax when using a covariance in a cosmology with one of the parameters set to a non-fiducial value (indicated in the title of each panel). As already noted above, the most dramatic variations occur for non-fiducial values of σ8 and Ωm. Quite remarkably, all the parameter errors considered here are already affected at low kmax values. In a more realistic data analysis setting one would sample the likelihood over a large number of points in the cosmological parameter space, while re-computing the whole data covariance at each evaluation. However, this may be unfeasible. To ease this problem one possibility is to model the cosmological dependence of the variance (i.e. the diagonal elements of the covariance) so that it can be varied during the sampling, while keeping the off-diagonal structure of the covariance fixed to a given cosmology. We explore this idea by repeating the Fisher forecast with an approximate covariance given by: where r k 1 ,k 2 is the correlation coefficient: In this case only the off-diagonal structure ofC is computed in a non-fiducial cosmology θ * , while the diagonal uses the correct parameter valuesθ. This mimics the scenario in which the variance is computed at the cosmology of the sampling points while the correlation coefficient is kept fixed at a given cosmology. The results are shown in Figures 7-8, where we can see that the impact on cosmological parameter errors is now significantly reduced, especially for low kmax values. We would like to stress that in a realistic galaxy clustering analysis, the impact of the cosmological dependence of the covariance on the cosmological parameter inference might have a smaller effect than what we have found here. This is because such an analysis will propagate uncertainties on the galaxy bias as well as the effect of shot noise. Whether such nuisance parameters can reduce or absorb the impact of the cosmological dependence of covariance goes beyond the scope of this work and we leave to a future study.

CONCLUSION
In this work we have investigated the cosmological dependence of the matter power spectrum covariance. To this purpose we have used the DEUS-PUR Cosmo set of simulations consisting of a large set of independent N-body realizations for different cosmological models characterized by different values of the cosmic matter density Ωm, the amplitude of matter density fluctuations σ8, the reduced Hubble parameter h and the dark energy equation of state w. This dataset has enabled us to estimate the covariance matrix for different cosmological parameter values and evaluate its first and second derivatives around a fiducial cosmological model. We found that the non-gaussian part of the covariance from the non-linear clustering of matter exhibits a varying degree of dependence on the different parameters at different redshifts. In particular σ8 and w have the largest impact at high redshift, while Ωm and h at low redshift. The analysis of the covariance derivatives indicates that the convergence of a second order Taylor expansion around the fiducial cosmology to approximate the cosmological dependence of the covariance is rather slow since the first order coefficients of the expansion are of order of unity. In the case of σ8, the first order coefficients is larger than unity at high redshift, potentially indicating a non-linear dependence of the covariance on this parameter. The different cosmological model parameters considered here span ∼ 10% variation around the fiducial cosmology, and yet it can lead to important differences in the power spectrum covariance up 100% level on some scales and redshifts. We have evaluated the impact of the cosmological parameter dependence of the covariance on cosmological parameter inference through a Fisher matrix approach. In particular, we have estimated the parameter uncertainties assuming a non-fiducial model covariance as function of the maximum mode kmax probed by a galaxy survey. We found significant differences with respect to the case with the covariance set to the fiducial cosmology. The largest effect occurs for non-fiducial values of σ8 and Ωm with deviations larger than 50% level already at modest kmax ∼ 0.1 − 0.2 h Mpc −1 . On the other hand, the impact on the degeneracy between pair of parameters is less significant, exceeding the 10% level only at kmax > 1 h Mpc −1 .
These results suggest that the cosmological parameter dependence of the non-Gaussian part of the covariance may impact the cosmological analyses from future surveys of the large scale structures. It is worth emphasizing that the quasilinear and non-linear scales over which the power spectrum covariance exhibits such a large dependence on the cosmo-logical parameters are also probed by cosmic shear measurements. Hence, it is reasonable to expect that the effects we have found in our analysis may also impact the parameter inference from weak lensing observations.
The current approach of keeping the covariance fixed to the fiducial cosmology when sampling the likelihood is likely to alter the shape of the posterior and consequently introduce systematic uncertainties on the cosmological parameter inference. Since it is not possible to run thousands of simulations to evaluate the covariance for each point in the parameter space that is explored by the likelihood sampling, the cosmological dependence of the covariance need to be modelled. Here, we have explored the possibility of modelling such dependence by fixing the off-diagonal part of the covariance matrix to the fiducial cosmology, while letting only the diagonal part vary with cosmology. This significantly reduces the misestimation of the parameter errors from the Fisher analysis, particularly at scales probed by galaxy clustering measurements. The dataset from DEUS-PUR Cosmo simulations provides an ideal benchmark to test models of the cosmological dependence of the covariance. To this purpose we have made the power spectra used in this work publicly available.
Further investigation is indeed necessary for a more robust assessment of the potential bias induced on the parameter estimation beyond the Fisher matrix approach. Our evaluation of the first-and second-order derivatives of the covariance can provide the foundation for a study that accounts for the non-Gaussian structure of the likelihood, for example using the so called Derivative Approximation for Likelihoods (DALI, Sellentin et al. 2014) method. We leave this investigation to a future work.

DATA AVAILABILITY
The power spectra from the DEUS-PUR Cosmo simulations can be downloaded at https://cosmo.obspm.fr/public-d atasets/.