Separate Universe Simulations with IllustrisTNG: baryonic effects on power spectrum responses and higher-order statistics

TL;DR

This study uses separate universe simulations with IllustrisTNG to measure baryonic effects on power spectrum responses and higher-order lensing statistics, finding that baryonic physics minimally impacts the response functions at certain scales and redshifts.

Contribution

It introduces a method to quantify baryonic effects on power spectrum responses using separate universe simulations with the IllustrisTNG model, highlighting the minimal impact of baryons on response functions at specific scales.

Findings

Baryonic effects do not significantly alter the growth-only response function G_1(k,z) at z<3.

Higher-order lensing statistics show varying sensitivity to baryonic effects, with the squeezed bispectrum being least affected.

Ignoring baryonic effects on lensing covariances slightly overestimates parameter errors, with negligible bias impact.

Abstract

We measure power spectrum response functions in the presence of baryonic physical processes using separate universe simulations with the IllustrisTNG galaxy formation model. The response functions describe how the small-scale power spectrum reacts to long-wavelength perturbations and they can be efficiently measured with the separate universe technique by absorbing the effects of the long modes into a modified cosmology. Specifically, we focus on the total first-order matter power spectrum response to an isotropic density fluctuation $R_1(k,z)$, which is fully determined by the logarithmic derivative of the nonlinear matter power spectrum ${\rm dln}P_m(k,z)/{\rm dln}k$ and the growth-only response function $G_1(k,z)$. We find that $G_1(k,z)$ is not affected by the baryonic physical processes in the simulations at redshifts $z < 3$ and on all scales probed ($k \lesssim 15h/{\rm Mpc}$, i.e. length scales $\gtrsim 0.4 {\rm Mpc}/h$). In practice, this implies that the power spectrum fully specifies the baryonic dependence of its response function. Assuming an idealized lensing survey setup, we evaluate numerically the baryonic impact on the squeezed-lensing bispectrum and the lensing super-sample power spectrum covariance, which are given in terms of responses. Our results show that these higher-order lensing statistics can display varying levels of sensitivity to baryonic effects compared to the power spectrum, with the squeezed-bispectrum being the least sensitive. We also show that ignoring baryonic effects on lensing covariances slightly overestimates the error budget (and is therefore conservative from the point of view of parameter error bars) and likely has negligible impact on parameter biases in inference analyses.

Figures

Tables

Table 1. Table 1: Parameters of the cosmologies simulated in this paper. The SepHigh  and SepLow cosmologies mimic, respectively, the effect of δ L ​ 0 = 0.05 subscript 𝛿 𝐿 0 0.05 \delta_{L0}=0.05 and δ L ​ 0 = − 0.05 subscript 𝛿 𝐿 0 0.05 \delta_{L0}=-0.05 long-wavelength density perturbations in the Fiducial cosmology. The comoving box size used is quoted in units of Mpc / h Mpc ℎ {\rm Mpc}/h in the corresponding cosmology. The value of the spectral index n s = 0.967 subscript 𝑛 𝑠 0.967 n_{s}=0.967 is the same for the three cosmologies, which share also the same amplitude of the initial power spectrum: that which corresponds to σ 8 ​ ( z = 0 ) = 0.816 subscript 𝜎 8 𝑧 0 0.816 \sigma_{8}(z=0)=0.816 in the Fiducial cosmology. We have performed simulations with two particle/mass element numbers N p = 625 3 subscript 𝑁 𝑝 superscript 625 3 N_{p}=625^{3} (called TNG300-3), N p = 1250 3 subscript 𝑁 𝑝 superscript 1250 3 N_{p}=1250^{3} (called TNG300-2) and with (dubbed Hydro) and without (dubbed Gravity) hydrodynamical processes (for the Hydro runs, the particle number is actually twice the quoted values: N p subscript 𝑁 𝑝 N_{p} gas cells and N p subscript 𝑁 𝑝 N_{p} dark matter mass elements). The same random seed was used to generate the initial conditions of all the simulations. The original IllustrisTNG simulations comprise also a higher-resolution set called TNG300-1 ( N p = 2500 3 subscript 𝑁 𝑝 superscript 2500 3 N_{p}=2500^{3} ), whose spectra (Springel et al., 2018 ) can be used to draw considerations on the convergence of our TNG300-2 results.

Name	${\Omega }_{m0}$	${\Omega }_{b0}$	${\Omega }_{\Lambda 0}$	$h$	$L_{\mathrm{box}}\left[\frac{\mathrm{Mpc}}{h}\right]$
Fiducial	$0.3089$	$0.0486$	$0.6911$	$0.6774$	$205$
SepHigh	$0.3194$	$0.0502$	$0.7146$	$0.6662$	$201.608$
SepLow	$0.2991$	$0.0471$	$0.6691$	$0.6884$	$208.337$

Equations

$$P_m(\mathbf{k},t|\mathcal{O}(\mathbf{x},t))=P_m(k,t)\left[1+\sum _{\mathcal{O}}{R}_{\mathcal{O}}(k,t)\mathcal{O}(\mathbf{x},t)\right],$$

$$P_m(\mathbf{k},t|{\delta }_L(\mathbf{x},t))$$

$$R_1(k,t)=\frac{1}{{\delta }_L(\mathbf{x},t)}\left[\frac{P_m(\mathbf{k},t|{\delta }_L(\mathbf{x},t))}{P_m(k,t)}-1\right],$$

$$R_1(k,t)=1-\frac{1}{3}\frac{\mathrm{dln}{P}_m(k,t)}{\mathrm{dln}k}+G_1(k,t),$$

$$R_K(k,t)=G_K(k,t)-\frac{\mathrm{dln}{P}_m(k,t)}{\mathrm{dln}k},$$

$${\tilde{\rho }}_m(t)={\rho }_m(t)\left[1+{\delta }_L(t)\right],$$

$$\frac{{\tilde{\Omega }}_{m0}{\tilde{h}}^2}{\tilde{a}(t{)}^3}=\frac{{\Omega }_{m0}{h}^2}{a(t{)}^3}\left[1+{\delta }_L(t)\right],$$

$${\tilde{\Omega }}_{m0}{\tilde{h}}^2={\Omega }_{m0}{h}^2.$$

$$1+{\delta }_L(t)={\left[1+{\delta }_a(t)\right]}^{-3},$$

$$H^2(t)={\left(\frac{\dot{a}(t)}{a(t)}\right)}^2=\frac{8\pi G}{3}\left[{\rho }_m(t)+{\rho }_{\Lambda }\right],$$

$${\tilde{H}}^2(t)={\left(\frac{\dot{\tilde{a}}(t)}{\tilde{a}(t)}\right)}^2=\frac{8\pi G}{3}\left[{\tilde{\rho }}_m(t)+{\rho }_{\Lambda }\right]-\tilde{K}\tilde{a}(t{)}^{-2},$$

$$\tilde{K}a(t{)}^{-2}$$

$$\frac{\tilde{K}}{H_0^2}=\frac{5}{3}{\Omega }_{m0}\frac{{\delta }_{L0}}{D(t_0)},$$

$$\ddot{D}(t)+2H(t)\dot{D}(t)-4\pi G{\rho }_m(t)D(t)=0.$$

$${\tilde{H}}_0$$

$${\tilde{\Omega }}_{m0}$$

$${\tilde{\Omega }}_{\Lambda 0}$$

$$1$$

$$t(a_{\mathrm{out}})\equiv {\int }_0^{a_{\mathrm{out}}}\frac{\mathrm{d}a}{aH(a)}={\int }_0^{a_{\mathrm{out}}^{\ast }}\frac{\mathrm{d}\tilde{a}}{\tilde{a}\tilde{H}(\tilde{a})}.$$

$$G_1(k,z)$$

$$G_1^{\text{SepHigh}}(k,z)$$

$$G_1^{\text{SepLow}}(k,z)$$

$$\kappa (\mathbf{\theta })={\int }_0^{{\chi }_S}\mathrm{d}\chi g(\chi ){\delta }_m(\mathbf{x}=\chi \mathbf{\theta },z),$$

$$C_{\kappa }(\ell )={\int }_0^{{\chi }_S}\mathrm{d}\chi \frac{g(\chi {)}^2}{{\chi }^2}{P}_m(k_{\ell },z),$$

$$B_m({\mathbf{k}}_h,{\mathbf{k}}_h,{\mathbf{k}}_s;z)$$

$$B_{\kappa }({\mathbf{\ell }}_1,{\mathbf{\ell }}_2,{\mathbf{\ell }}_3)={\int }_0^{{\chi }_S}\mathrm{d}\chi \frac{g(\chi {)}^3}{{\chi }^4}{B}_m({\mathbf{k}}_{{\ell }_1},{\mathbf{k}}_{{\ell }_2},{\mathbf{k}}_{{\ell }_3};z),$$

$$B_{\kappa }({\ell }_h,{\ell }_h,{\ell }_s)$$

$${\mathrm{Cov}}_{\kappa }({\ell }_1,{\ell }_2)={\mathrm{Cov}}_{\kappa }^G({\ell }_1,{\ell }_2)+{\mathrm{Cov}}_{\kappa }^{SSC}({\ell }_1,{\ell }_2)+{\mathrm{Cov}}_{\kappa }^{cNG}({\ell }_1,{\ell }_2).$$

$${\mathrm{Cov}}_{\kappa }^G({\ell }_1,{\ell }_2)=\frac{4\pi }{{\Omega }_W{\ell }_1\Delta {\ell }_1}{\left[C({\ell }_1)+C_{\mathrm{shape}}\right]}^2{\delta }_{ij},$$

$${\mathrm{Cov}}_{\kappa }^{SSC}({\ell }_1,{\ell }_2)$$

$$P_m(k,z)=P_m^{Emu.}(k,z){\mathcal{I}}^{\mathrm{baryon}}(k,z),$$

Full text

††pagerange: Separate Universe Simulations with IllustrisTNG: baryonic effects on power spectrum responses and higher-order statistics–Separate Universe Simulations with IllustrisTNG: baryonic effects on power spectrum responses and higher-order statistics

1 Introduction

The large-scale distribution of matter in the Universe is one of the best tools to address some of the deepest open problems in fundamental physics. These include studies of the nature of dark energy and dark matter, tests of the law of gravity on large scales, tests of early-Universe inflationary models, as well as constraints on the absolute value of neutrino masses. Within large-scale structure analyses, most attention has been devoted to studies of the clustering pattern of galaxies, which trace the underlying clustering of matter in a biased manner, and to studies of background galaxy shapes that are coherently distorted (cosmic shear) by intervening gravitational potentials via weak gravitational lensing. Observational analyses of galaxy clustering (Alam et al., 2017; Beutler et al., 2017; Sánchez et al., 2017; Gil-Marín et al., 2017) and weak lensing data (Hildebrandt et al., 2017; Abbott et al., 2018; Joudaki et al., 2018; Hikage et al., 2018) have provided tighter and tighter constraints on various cosmological parameters over the years, and are expected to continue to do so as next-generation surveys like DESI (Levi et al., 2013), LSST (The LSST DESC et al., 2018) and Euclid (Laureijs et al., 2011) come online.

The most popular way to characterize the statistical information encoded in large-scale structure is via the two-point correlation function, or its Fourier transform, the power spectrum. In galaxy clustering studies, the main theoretical modeling steps involve descriptions of galaxy bias (Desjacques et al., 2018) and redshift-space distortions. The most frequently used methods are based on perturbation theory (Bernardeau et al., 2002), whose limited accuracy in the nonlinear regime of structure formation has been restricting most analyses to sufficiently large-scales $k\lesssim 0.3\ h/{\rm Mpc}$ ($\gtrsim 30\ {\rm Mpc}$). Weak-lensing analyses, on the other hand, are sensitive to the total clustering of matter, which is better understood in the nonlinear regime and permits extending the analysis to smaller length scales. Here, an important factor in setting the smallest length scale used in observational analyses is associated with the uncertain impact that baryonic physical processes can have on the small-scale power spectrum. A number of different hydrodynamical simulations (van Daalen et al., 2011; Hellwing et al., 2016; Springel et al., 2018; Chisari et al., 2018) have consistently reported a suppression of the total clustering power on scales $k\gtrsim 1h/{\rm Mpc}$, but the precise shape and magnitude of the effect depend on the exact implementation of the baryonic physical processes, most notably feedback effects by active galactic nuclei (AGN) and the corresponding impact on halo dark matter and gas profiles (Schneider et al., 2018). Current weak lensing analyses circumvent this problem by either simply discarding the data in regimes that may be affected by baryons, or modeling the impact of baryons with a number of nuisance parameters that are subsequently marginalized over. Given that a significant part of the constraining power in weak-lensing surveys comes from small scales, the development of strategies to incorporate or mitigate the impact of baryonic effects is a pressing issue to address. This has consequently become the subject of recent active work on lensing power spectrum analyses (Semboloni et al., 2011a; Schneider & Teyssier, 2015; Eifler et al., 2015; Mead et al., 2015; Harnois-Déraps et al., 2015; Chisari et al., 2018; Huang et al., 2018; Schneider et al., 2018).

The fact that the late-time distribution of matter in the Universe is appreciably non-Gaussian distributed immediately leads, however, to the observation that the power spectrum (despite being the most popular summary statistic) cannot describe all of the available information. This naturally motivates incorporating higher-order $N$-point correlation functions in observational analyses as a way to gain access to the information that is missed by the power spectrum, and consequently, obtain improved constraints on cosmological parameters. For instance, the gains that can be achieved by incorporating the bispectrum (the Fourier transform of the three-point function) in cosmological analyses are by now well established (Cooray & Hu, 2001; Dodelson & Zhang, 2005; Sefusatti et al., 2006; Sato & Nishimichi, 2013; Kayo & Takada, 2013; Gil-Marín et al., 2017; Chan & Blot, 2017; Coulton et al., 2018; Rizzato et al., 2018; Yankelevich & Porciani, 2018). The body of work on higher-order $N$-point functions in the nonlinear regime is however much smaller compared to that on the power spectrum. The reason for this can be traced back to the fact that the estimators of the $N$-point functions become rapidly numerically intensive with increasing $N$ (Sabiu et al., 2019), as well as the fact that these $N$-point functions live in higher-dimensional parameter spaces that require large-volume, high-resolution cosmological simulations to obtain sufficiently precise measurements. Note also that knowledge of higher-order statistics is important even for power spectrum analysis since its covariance is given by a specific configuration of the four-point correlation function (Scoccimarro et al., 1999).

The larger numerical resources needed to study higher-order correlation functions become even more critical for the case of hydrodynamical simulations, which are appreciably more demanding to run than equivalent gravity-only counterparts. This helps explain the relatively small number of results on the impact of baryonic physics on higher-order lensing correlation functions, which remains a largely unexplored subject in the literature, despite a few interesting existing studies. For instance, Semboloni et al. (2013) used results from the OWLS suite of simulations (Schaye et al., 2010) to predict the impact of baryonic physics on the three-point lensing correlation function. Further, Osato et al. (2015) investigated the impact of baryonic effects on weak-lensing peaks and Minkowski functionals, which are statistics sensitive to higher-order $N$-point functions (see also Yang et al. (2013), and Castro et al. (2018) for a study of baryonic effects on the probability distribution function of weak-lensing maps). Studies such as these can provide interesting insights on the physics of nonlinear structure formation and are therefore crucial to fully exploit the constraining power of existing and future datasets.

In this paper, we take a number of steps forward in understanding the impact of baryonic effects on higher-order $N$-point correlation functions by using the response approach to perturbation theory developed by Barreira & Schmidt (2017b, a) as an extension of previous work done by Chiang et al. (2014); Wagner et al. (2015b); Chiang et al. (2015). Power spectrum responses are functions of scale and time that encode how the small-scale power spectrum responds (or reacts) to the presence of long-wavelength density and tidal field perturbations. These response functions can be measured efficiently in the nonlinear regime using so-called separate universe simulations (Li et al., 2014a, b; Wagner et al., 2015a; Dai et al., 2015; Baldauf et al., 2016; Wagner et al., 2015b; Schmidt et al., 2018), which are simulations that absorb the impact of the long-wavelength modes into a redefinition of the cosmology that is simulated. In a perturbation theory sense, these response functions describe the coupling between small- and large-scale Fourier modes, quite importantly, for nonlinear values of the small-scale modes. In practice, this permits one to straightforwardly evaluate a number of higher-order correlation functions in the nonlinear regime (Barreira & Schmidt, 2017b). Importantly, accurate and precise measurements of the response functions can be obtained with the same particle resolution and simulation volume as for the matter power spectrum, thereby permitting systematic studies of higher-order correlation functions without exacerbated demand for numerical resources.

Concretely, in this paper we apply the separate universe simulation technique using the IllustrisTNG galaxy formation model (Weinberger et al., 2017; Pillepich et al., 2018b) to measure for the first time the matter power spectrum responses in the presence of baryonic physical processes. All past matter power spectrum response measurements (Li et al., 2014a; Wagner et al., 2015b; Baldauf et al., 2016; Schmidt et al., 2018) have been carried out using gravity-only separate universe simulations. Once the dependence of the power spectrum and its response functions on baryonic effects is determined, the machinery of the response approach to perturbation theory (Barreira & Schmidt, 2017b) can then be used to predict the corresponding impact of baryons on a number of higher-order correlation functions. We shall do so explicitly for the case of the squeezed weak-lensing convergence bispectrum, as well as for the super-sample covariance (SSC; Takada & Hu (2013); Barreira et al. (2018b)) part of the covariance matrix of the weak lensing convergence power spectrum (which is a four-point function).

The rest of this paper is organized as follows. In Sec. 2, we outline the main concepts behind matter power spectrum responses and the separate universe simulations we performed. In Sec. 3, we present and discuss our numerical measurements of the power spectrum response functions with the IllustrisTNG model. In Sec. 4, we use the response approach to perturbation theory to calculate the impact of baryonic effects on the squeezed lensing bispectrum and lensing power spectrum covariances. We summarize and conclude in Sec. 5.

2 Power spectrum responses and Separate Universe simulations

In this section, we outline the main concepts behind power spectrum response functions and separate universe simulations. These have been discussed at length in past works and we thus refer the interested reader to the cited literature for more details.

2.1 Power spectrum responses

Power spectrum responses are functions that describe how the local, nonlinear matter power spectrum responds to the presence of long-wavelength perturbations. Formally, we can write (Barreira & Schmidt, 2017b)

[TABLE]

where $P_{m}(\bm{k},t|\mathcal{O}(\bm{x},t))$ is the nonlinear matter power spectrum at physical time $t$ measured in some local volume around position $\bm{x}$ that is embedded in a long-wavelength perturbation $\mathcal{O}$, $P_{m}(k,t)$ is the global matter power spectrum measured at cosmic mean density and $R_{\mathcal{O}}(k,t)$ are response coefficient functions that specify the response of the power spectrum to each specific long-wavelength perturbation $\mathcal{O}$ (see also Bertolini & Solon (2016)). The meaning of the local power spectrum is that it is measured in a region of size $L$ that is sufficiently smaller than the wavelength of the perturbation $\mathcal{O}$, i.e., $1/k<L\ll 1/p$, where $p$ is the Fourier mode of the perturbation $\mathcal{O}$.

The total matter power spectrum response is specified by listing all possible long-wavelength perturbations $\mathcal{O}$ the local power spectrum can depend on. The leading order ones (see Barreira & Schmidt (2017b) for a discussion) are those that involve second spatial derivatives of the gravitational potential $\Phi$, such as, isotropic density perturbations $\mathcal{O}(\bm{x},t)=\delta(\bm{x},t)\sim\nabla^{2}\Phi(\bm{x},t)$ or tidal field perturbations $\mathcal{O}(\bm{x},t)=\hat{k}^{i}\hat{k}^{j}K_{ij}(\bm{x},t)=\hat{k}^{i}\hat{k}^{j}\left[\nabla_{i}\nabla_{j}/\nabla^{2}-\delta_{ij}/3\right]\delta(\bm{x},t)$ (note how for tidal fields there is an explicit dependence on the direction of the small-scale mode $\bm{k}$, which is why the local power spectrum on the left-hand side of Eq. (1) can be anisotropic in general (Akitsu et al., 2017; Akitsu & Takada, 2018; Li et al., 2018; Schmidt et al., 2018)). In this paper, we focus on isotropic density fluctuations, for which Eq. (1) specializes to

[TABLE]

i.e., the response functions $R_{n}(k,t)$ can be interpreted as the coefficients of a Taylor expansion of the local power spectrum in the long-wavelength isotropic overdensity $\delta_{L}(\bm{x},t)$ (the subscript L emphasizes the long-wavelength nature of the perturbation). For small enough $\delta_{L}$ values, we can linearize Eq. (2) and solve for the first-order isotropic power spectrum response

[TABLE]

which is the response function we wish to measure below with separate universe simulations of the IllustrisTNG model (see Wagner et al. (2015b) for measurements of higher-order isotropic responses, $R_{n}$ $(n>1)$, in gravity-only simulations).

The presence of the long-wavelength perturbation $\delta_{L}(x,t)$ induces three physically distinct effects on the small-scale power spectrum (Li et al., 2014a; Wagner et al., 2015b):

Reference density effect. This accounts for the fact that measurements of the local power spectrum are performed with respect to a modified mean density $\rho_{m}(\bm{x},t)=\rho_{m}(t)\left[1+\delta_{L}(\bm{x},t)\right]$, where $\rho_{m}(t)$ is the cosmic mean density w.r.t. which the global power spectrum is defined 111We distinguish background-averaged quantities by omitting their dependence on the spatial coordinate $\bm{x}$.. 2. 2.

Dilation effect. This describes the rescaling of the local scale factor $a(\bm{x},t)$, which evolves differently than the global one $a(t)$ because of the long-wavelength mode. In other words, the expansion of spacetime inside the perturbation $\delta_{L}(\bm{x},t)$ takes place at a different rate than the global expansion of the Universe. In practice, this induces a rescaling of the physical scales inside the perturbation. 3. 3.

Growth-only effect. This describes the actual mode-coupling between the linear, long-wavelength mode and nonlinear small-scale power spectrum modes.

Accounting for these effects permits one to write the first-order response $R_{1}(k,t)$ as

[TABLE]

where $G_{1}(k,t)$ is called the growth-only response function that describes effect (iii) above in isolation. Strictly speaking, this is the only term for which separate universe simulations are needed, with the remainder of the response function $R_{1}$ being completely specified by the logarithmic derivative of the nonlinear matter power spectrum ${\rm dln}P_{m}(k,t)/{\rm dln}k$.

In some of our results below in Sec. 4, we will also use the first-order response to a long-wavelength tidal field, which admits the following decomposition

[TABLE]

where $G_{K}(k,t)$ is the growth-only response to a tidal field (Schmidt et al., 2018).

2.2 Separate Universe formalism

The following observation forms the basis of the implementation of the separate universe formalism to measure power spectrum responses (Li et al., 2014a, b; Wagner et al., 2015a; Dai et al., 2015; Baldauf et al., 2016; Wagner et al., 2015b; Schmidt et al., 2018): local structure formation in a given cosmology in the presence of a long-wavelength perturbation is equivalent to global structure formation in an appropriately modified cosmology. For the particular case of isotropic density perturbations $\delta_{L}$ that we focus on in this paper, this corresponds to a modified cosmology with a modified mean background density

[TABLE]

where, here and throughout, quantities with and without a tilde are defined in the modified and fiducial cosmologies, respectively. In Eq. (6), we have dropped the coordinate argument $\bm{x}$ as we assume that the long-wavelength perturbation $\delta_{L}$ is effectively constant inside the local volume where one wishes to measure the power spectrum; this is why the effect can be mimicked by rescaling the mean matter density everywhere in the modified cosmology. Equation (6) can be written as

[TABLE]

where $\Omega_{m0}$ is the fractional matter density parameter and $h$ is the dimensionless Hubble constant $H_{0}=100h\ {\rm km/s/Mpc}$. Owing to the different mean background densities, the expansion of spacetime takes place at different rates in the two cosmologies, i.e., $a(t)\neq\tilde{a}(t)$. We adopt the usual convention that $a(t_{0})=1$ (where $t_{0}$ is the present-day physical time in the fiducial cosmology), and at some sufficiently early epoch $t_{i}$, when the amplitude of the long-mode can be neglected, we take $a(t_{i})=\tilde{a}(t_{i})$. In this early time limit, Eq. (7) becomes

[TABLE]

This equation can be plugged back into Eq. (7) to derive

[TABLE]

where $\tilde{a}(t)=a(t)\left[1+\delta_{a}(t)\right]$.

In this paper, our fiducial cosmology is given by a spatially-flat $\Lambda{\rm CDM}$ model, for which the Friedmann equation is given by

[TABLE]

where an overdot denotes a derivative w.r.t. physical time (we ignore the impact of radiation and massive neutrinos on the expansion rate). For the modified cosmology, the Friedmann equation follows as

[TABLE]

where we have allowed for non-zero spatial curvature and the physical energy density of the cosmological constant $\rho_{\Lambda}$ is not affected by the long-wavelength density mode. The value of $\tilde{K}$ can be derived by taking the difference of Eqs. (10) and (11)

[TABLE]

Given that $\tilde{K}$ is a constant, we can derive its value at any time using Eq. (12). At early times, the expansion of the Universe is matter-dominated, $H^{2}(t)\approx\Omega_{m0}H_{0}^{2}a(t)^{-3}$, and we can linearize Eq. (9) to get $\delta_{a}(t)\approx-\delta_{L}(t)/3$; further noting that during matter domination $\delta_{L}(t)\propto a(t)$, we have $\dot{\delta}_{a}(t)=H(t)\delta_{a}(t)$. We can thus write

[TABLE]

where $\delta_{L0}$ is the value of the long-wavelength overdensity today, and $D(t)$ is the linear growth factor of the fiducial cosmology that obeys

[TABLE]

For the definition of cosmological parameters such as $H_{0}$ and $\Omega_{m0}$, $N$-body/hydro codes define their numerical values at the time the scale factor is equal to unity, but in the modified and fiducial cosmologies this occurs at different physical times. The parameters that are used in the code should thus be evaluated at a time $\tilde{t}_{0}$ defined as $\tilde{a}(\tilde{t}_{0})=1$:

[TABLE]

where the quantity $\delta_{H}$ is defined by Eq. (15) and is used subsequently in Eqs. (16) and (17). The fractional curvature density is given by $\tilde{\Omega}_{K0}=-\tilde{K}/\tilde{H}_{0}^{2}$. What is left to specify is the value of $\delta_{H}$, which can be obtained from the present-day value of the Friedmann equation in the modified cosmology:

[TABLE]

where we have used that the fiducial cosmology is spatially flat $\Omega_{m0}+\Omega_{\Lambda 0}=1$. This finalizes the derivation of the cosmological parameters of the modified cosmology given the parameters of the fiducial cosmology, as well as the present-day value of the long-wavelength mode $\delta_{L0}$.

2.3 Separate Universe simulations

The simulation results we present in this paper were obtained using the moving-mesh code AREPO (Springel, 2010; Pakmor et al., 2016) with the IllustrisTNG hydrodynamical galaxy formation model (Weinberger et al., 2017; Pillepich et al., 2018b), which is an upgraded version of the original Illustris model (Genel et al., 2014; Vogelsberger et al., 2014). In this paper, we focus on the impact of baryons on matter power spectrum responses, and hence, our results are complementary to those of Springel et al. (2018) who investigated already the corresponding impact on the power spectrum itself, for the same implementation of the TNG model. Here, for brevity, we skip describing the details of the IllustrisTNG model and refer the reader to the cited literature for more details, including Marinacci et al. (2018); Pillepich et al. (2018a); Naiman et al. (2018); Springel et al. (2018); Nelson et al. (2018a) for the first results from the IllustrisTNG simulations and Nelson et al. (2018b) for an overview of the publicly available simulation data.

The cosmological parameters of the fiducial and separate universe simulations we consider are summarized in Table 1. The Fiducial cosmology parameters are the same as in the main IllustrisTNG runs; the two separate universe runs SepHigh and SepLow correspond to $\delta_{L0}=0.05$ and $\delta_{L0}=-0.05$, respectively. This choice for the amplitude of the long-wavelength modes is motivated by having a strong enough perturbation capable of inducing sizeable changes to the power spectrum, but keeping higher order corrections $\mathcal{O}(\delta_{L}^{2})$ in Eq. (2) negligible. Within this linear regime ($\delta_{L0}\ll 1$), any two values of $\delta_{L0}$ are sufficient to make one measurement of the first-order response by finite differencing the corresponding power spectra. Here, we opt to utilize three values values ($\delta_{L0}=0,\pm 0.05$) to be able to make a measurement using an overdensity and an underdensity (note however that the result should be independent of the sign of $\delta_{L0}$); cf. Eqs. (20)-(22) below). The simulation box size of the Fiducial cosmology is $L_{\rm box}=205{\rm Mpc}/h\approx 300{\rm Mpc}$ and we consider two particle/mass element resolutions: $N_{p}=625^{3}$ and $N_{p}=1250^{3}$, to which we refer below as TNG300-3 and TNG300-2, respectively.222This is the same terminology as that used in the main IllustrisTNG runs. In the latter, for this box size, there is a higher-resolution set called TNG300-1 with $N_{p}=2500^{3}$. In this paper, we are interested in the clustering of matter on scales $k\lesssim 10\ h/{\rm Mpc}$, on which Springel et al. (2018) has shown that the TNG300-1 and TNG300-2 runs display $\%$-level agreement. To isolate the impact of baryonic physics on the power spectrum responses we have run standard IllustrisTNG hydrodynamical simulations (dubbed Hydro below), as well as gravity-only simulations without baryonic processes (dubbed Gravity below).

The initial conditions were generated with the N-GenIC code (Springel, 2015) using the Zel’dovich approximation at the same initial redshift $z_{i}=127$ for all the cosmologies (this ignores the fact that $z(t)\neq\tilde{z}(t)$, but which is a small effect at early times).333The use of the Zel’dovich approximation, although not as accurate as second-order perturbation theory methods (2LPT), is not expected to have a visible impact on our results since the simulations start at sufficiently high redshift and our main results consist of power spectra ratios, and not their absolute values. The shape of the linear matter power spectrum is the same in the three cosmologies because the physical matter densities are the same (cf. Eq. (8)). We generate the initial power spectrum file using the CAMB code (Lewis et al., 2000; Lewis & Challinor, 2011) at $z=0$ and rescale it to $z_{i}$ using the growth factor $D(t)$ of the fiducial cosmology. An extra step that is needed for the modified cosmologies is to change the units of the wavenumbers and power spectrum from $h/{\rm Mpc}$ and ${\rm Mpc}^{3}/h^{3}$ to $\tilde{h}/{\rm Mpc}$ and ${\rm Mpc}^{3}/\tilde{h}^{3}$; effectively, this is done by multiplying the $k$ values by $\left(1+\delta_{H}\right)^{-1}$ and the power spectra values by $\left(1+\delta_{H}\right)^{3}$.

Another noteworthy consideration concerns the choice of the box size for the SepHigh and SepLow simulations. Given that the measurement of the responses is effectively the ratio of the spectra in the modified and fiducial cosmologies (cf. Eq. (3)), it is in one’s interest to use the same random seed in the generation of the initial conditions and then ensure that the phases of the density fluctuations coincide in between the two cosmologies as this significantly reduces cosmic variance in the measurement. Specifically, if $\tilde{L}_{\rm box}$ is the comoving size of the box of the modified cosmology in units ${\rm Mpc}/\tilde{h}$ and $L_{\rm box}$ is the comoving size of the fiducial simulation box in units ${\rm Mpc}/h$, then equating the two yields $\tilde{L}_{\rm box}=L_{\rm box}\tilde{h}/h$ and ensures that the density fluctuations match on comoving scales and at all times. In other words, the comoving fundamental modes of the simulations of the modified and fiducial cosmologies are the same at all times $\tilde{k}_{F}=(2\pi/\tilde{L}_{\rm box})\tilde{h}/{\rm Mpc}=(2\pi/{L}_{\rm box}){h}/{\rm Mpc}=k_{F}$. The box size values we use are listed in Table 1. It is important to note that with this choice, the ratio of the spectra in the simulations of the modified and fiducial cosmologies at the same numerical value of $k$ yields actually a direct measurement of the growth-only response $G_{1}(k,t)$ (see e.g. Sec. III. B. of Li et al. (2014a) for a detailed explanation). The full isotropic response can, however, always be obtained with Eq. (4) using the derivative of the matter power spectrum.444An alternative approach is to choose to match the physical scales, $\tilde{a}(t)\tilde{L}_{\rm box}{\rm Mpc}/\tilde{h}=a(t)L_{\rm box}{\rm Mpc}/h$, which directly measures the total response $R_{1}(k,t)$, but with the price that now the density fluctuations have matching scales only at a single physical time. That is, the measurement of the response at each redshift requires a set of separate universe simulations with different box sizes to be performed.

We wish to compare the fiducial and modified cosmologies at the same physical time, but this happens at different numerical values of the scale factor, which is the coordinate that AREPO and most $N$-body/hydro codes use to specify the output times of the particle/mass element data. Hence, if $a_{\rm out}$ represents the scale factor of the Fiducial cosmology we are interested in, then we output the particle/mass element data in the modified cosmologies at a scale factor $a_{\rm out}^{*}$, defined as

[TABLE]

We will show results for $a_{\rm out}=0.25,0.33,0.5,0.66,1.0$, which correspond to $z_{\rm out}=3,2,1,0.5,0$. In the SepHigh and SepLow simulations, these correspond, respectively, to $a_{\rm out}^{*}\approx 0.249,0.331,0.495,0.658,0.983$ and $a_{\rm out}^{*}\approx 0.251,0.336,0.505,0.675,1.017$.

Finally and important for the interpretation of our results is the fact that the hydrodynamical physics model is kept the same in between the fiducial and separate universe runs. Our results should thus be interpreted as solving nonlinear structure formation in the presence of a long-wavelength density perturbation, in a Universe with the hydrodynamical physical processes as specified by the IllustrisTNG model.555We note for completeness that the simulations make use of tabulated values of UV background radiation as a function of redshift, which should be converted to physical time for the separate universe runs. We skipped adjusting for these slight shifts in timing in the reionization history, which is slightly inconsistent, but which is not expected to have a visible impact on our total matter power spectrum measurements.

3 Response measurements

We measure the power spectrum in the simulations with the publicly-available nbodykit code (Hand et al., 2018). We assign the mass distribution of the particle/mass elements in the simulation onto a regular grid with $1024$ cells on a side, which is then Fourier-transformed to measure the power spectrum (for the case of a gas cell, we approximate all of its mass to lie at its geometrical center, before interpolating it to the regular grid.) In the Hydro runs, the total mass is distributed across different mass element types: gas, dark matter, stars and black holes. Unless otherwise specified, our power spectrum measurements in the Hydro runs correspond to the total mass in all mass components. We also always consider power spectra with the mass-weighted shot-noise power subtracted $P_{\rm shot}=L_{\rm box}^{3}/N_{\rm eff}$, where $N_{\rm eff}=(\sum_{i}m_{i})^{2}/\sum_{i}m_{i}^{2}$. Further, when we quote redshift values and units with $h$ retained, we always assume the Fiducial cosmology.

Given the measured power spectrum, the first-order growth-only isotropic response function is calculated as

[TABLE]

with

[TABLE]

where $\delta_{L}(z)=\delta_{L0}D(z)/D(z=0)$ and the numerical values of the comoving modes $k$ are the same in units of ${\rm Mpc}^{-1}$ in the Fiducial, SepHigh and SepLow simulations (the superscripts indicate which cosmology each measurement corresponds to). In theory, $G_{1}^{\text{SepHigh}}(k,z)=G_{1}^{\text{SepLow}}(k,z)$, and hence, any discrepancy between the two results can be used to signal inadequacies in our numerical measurements; the difference between $G_{1}^{\text{SepHigh}}$ and $G_{1}^{\text{SepLow}}$ can also be used as a rough guide for the precision of our response measurements.

3.1 Response of the total matter field

Figure 1 shows $G_{1}(k,z)$ measured from the Hydro (red) and Gravity (blue) simulations of the TNG300-2 box. The solid lines show the result as given by Eq. (20), with the shaded area indicating the region bracketed by $G_{1}^{\text{SepHigh}}$ and $G_{1}^{\text{SepLow}}$. On large scales, both the Hydro and Gravity runs recover the expected redshift- and scale-independent linear theory result $G_{1}^{\rm linear}(k,z)=26/21$. Furthermore and interestingly, the Hydro and Gravity response measurements continue to agree with one another to very good precision for all of the scales and epochs shown. This includes scales $k\gtrsim 5\ h/{\rm Mpc}$ deep in the nonlinear regime of structure formation where baryonic processes display already a marked impact on the power spectrum (cf. Fig. 4 below). More specifically, the relative difference between the Hydro and Gravity results is below $4\%$ and it does not exhibit any trend with redshift or scale. There are two exceptions at $z=0.5,k\in\left[2,4\right]\ h/{\rm Mpc}$ and $z=1,k\in\left[5,10\right]\ h/{\rm Mpc}$, in which the difference becomes of order $6\%$ and $9\%$, respectively (not clear from the scale in Fig. 1). We attribute these exceptions to statistical fluctuations associated with having only one set of phases of the initial conditions. We note also that in these exceptions, the Hydro and Gravity results remain within the shaded area of one another. We thus report that, within the precision of our measurements with the IllustrisTNG model, there is no evidence of any baryonic physics imprint on the first-order growth-only power spectrum response function. In other words, the physical processes through which baryons alter the clustering of matter are largely insensitive to structure formation taking place at cosmic mean or in slightly overdense/underdense regions.

For reference, Fig. 1 shows also the result obtained by Wagner et al. (2015b) using gravity-only simulations with the Gadget code (Springel, 2005). The small differences between their result and ours can be attributed to different cosmological parameter values. Their curve is also smoother as it is an average over 16 realizations of the initial conditions.

Figure 2 compares the $G_{1}(k,z)$ measurements from the TNG300-2 (solid curves) and TNG300-3 (dashed curves and shaded areas) boxes. The results are shown for $z=0$ and $z=1$ and they illustrate that (i) the Gravity runs of both resolutions are in good agreement for all of the scales shown, but (ii) the Hydro results are somewhat discrepant on scales $k\gtrsim 2h/{\rm Mpc}$ ($z=1$) and $k\gtrsim 0.3h/{\rm Mpc}$ ($z=0$). The reduced accuracy of the Hydro TNG300-3 simulations in measuring the response function is not entirely surprising given its poorer performance in resolving as realistically the intricate baryonic physical processes that take place on small length scales (Pillepich et al., 2018a; Springel et al., 2018; Nelson et al., 2018a). The width of the red shaded area in Fig. 2 is also telling of this poorer resolution as it indicates that the power spectrum is responding in a qualitatively different way to an overdense and underdense perturbation, when this response should be the same.

The lack of convergence between the Hydro runs of the TNG300-2 and TNG300-3 boxes could raise the question of whether or not the Hydro TNG300-2 results are actually converged. We can list, however, a number of reasons that convince us about the robustness of our $G_{1}(k,z)$ measurements at TNG300-2 resolution. First, Fig. 6 of Springel et al. (2018) shows that the clustering of matter with the TNG300-2 resolution agrees to better than $2\%$ with the higher resolution simulation TNG300-1 ($N_{p}=2500^{3}$), on the scales $k\lesssim 15h/{\rm Mpc}$ we probe in this paper. Second, compared to the TNG300-3 case, the appreciably smaller difference between the $G_{1}^{\text{SepHigh}}$ and $G_{1}^{\text{SepLow}}$ curves (depicted by the shaded area) is also indicative of a faithful capture of the response of the power spectrum to the long-wavelength perturbation by the TNG300-2 resolution. Additionally, the nearly perfect overlap between the Gravity and Hydro TNG300-2 results shown in Fig. 1 (which is absent at TNG300-3 resolution) can also be invoked to strengthen further the conclusion that the result is converged. Alternatively, this would mean that any errors due to insufficient resolution were canceling almost exactly some actual physical baryonic impact on $G_{1}(k,z)$ on all scales and all redshifts, which we deem to be highly unlikely. Finally, the degree of numerical convergence in hydrodynamical simulations is in general redshift dependent (being better at higher redshift), but our results show no evidence of a redshift dependent difference between the $G_{1}(k,z)$ measured with the TNG300-2 Gravity and Hydro simulations; this is not the case at TNG300-3 resolution, for which the Gravity and Hydro results are closer at $z=1$ than at $z=0$.

3.2 Response of individual mass components

The results discussed thus far for the Hydro simulations correspond to the response of the total matter power spectrum computed using the total mass distribution in gas, dark matter, stars and black holes. It is also interesting, however, to measure the response of the power spectrum of the individual mass components. The result is shown in Fig. 3 for the gas (green) and dark matter (black) components, as well as for the mass distribution comprising gas cells, stars and black holes (magenta). The results correspond to the TNG300-2 box and the total power spectrum (red) is repeated from Fig. 1 for comparison. The figure shows that the dark matter power spectrum responds effectively in the same way as the total matter power spectrum; this is as expected as the total clustering of matter is dominated by the (more abundant) dark matter component.

The response of the power spectrum of the gas component displays however a slightly different behavior. This is better noticed at redshifts $z>1$, where the additional hydrodynamical processes in which the gas participates seem to make its power spectrum respond less strongly to the presence of long-wavelength modes. These hydrodynamical phenomena include additional pressure forces felt by the gas cells, as well as the fact that the mass in gas cells is not conserved, as stars form and black holes grow in the simulations. The scales at which the gas response deviates from that of the dark matter is also redshift dependent: $k\gtrsim 3h/{\rm Mpc}\ (z=1)$, $k\gtrsim 0.8h/{\rm Mpc}\ (z=2)$, and effectively all scales shown at $z=3$, including large linear scales where the gas response plateaus at a slightly different value. Interestingly, the response of the gas+stars+black holes power spectrum is closer to that of the total power spectrum response; this is particularly noticeable at $z=3$. The mass in the gas+stars+black holes is conserved in the simulations, which suggests that the non-conservation of mass in the gas component plays an important role in the corresponding response function. A reasonable physical picture is that the enhanced star formation and black hole growth caused by the long-wavelength density perturbation leaves less mass in the gas component, effectively making its clustering less responsive to the perturbation. Once mass conservation is restored by including the mass in stars and black holes, then the corresponding response function is brought closer to that of the dark matter component (which is conserved in the TNG model). This physical picture is also in line with the fact that star formation rates peak between $z=2$ and $z=3$, which coincides with the epochs when the gas response differs the most from the others in Fig. 3.

Regarding the responses of the gas and gas+stars+black holes components, we note that the difference between $G_{1}^{\rm SepHigh}$ and $G_{1}^{\rm SepLow}$ (shaded areas) can be non-negligible on small length scales, $k\gtrsim 4h/{\rm Mpc}$ at $z=0$, $z=0.5$ and $z=3$, which may signal some lack of convergence of the TNG300-2 box for these components there. On larger scales, however, the width of the shaded areas is smaller (i.e., $G_{1}^{\rm SepHigh}\approx G_{1}^{\rm SepLow}$), which suggests that our results are faithful representations of the expected behavior of these response functions. A stronger case for the convergence of our gas response results at TNG300-2 resolution on scales $k\lesssim 4\ h/{\rm Mpc}$ can be made by noting that the clustering of the gas component itself is reasonably converged at this resolution. This can be seen in Fig. 5 of Springel et al. (2018), which shows that the 2-point correlation function of the gas component of the TNG300-2 resolution agrees with that of the higher resolution TNG300-1 simulation to better than $5\%$ on length scales $r\gtrsim 0.5{\rm Mpc}/h$. Figure A2 of Pillepich et al. (2018b) also shows that the gas fraction in haloes with mass $2\times 10^{12}-2\times 10^{13}\ M_{\odot}/h$ found at TNG300-2 resolution (labeled TNG25-128 there) agrees to better than a few percent with that found in higher resolution runs. We note, however, that the clustering of stars and black holes at TNG300-2 resolution is not as converged as the clustering of gas and dark matter (cf. Fig. 5 of Springel et al. (2018)). Consequently, a more quantitative study of the physical picture outlined above to explain the differences between the responses of the gas and gas+stars+black holes components would certainly benefit from higher resolution response measurements.

We leave as future work a more in-depth investigation of the response of the gas component. Such a study can provide interesting theoretical insights to Lyman-alpha forest flux power spectrum analyses (McDonald et al., 2006; Slosar et al., 2013; Chabanier et al., 2018), which probe the line-of-sight distribution of neutral hydrogen clouds via the Lyman-alpha absorption lines they imprint on background quasar spectra. Similarly, line intensity mapping studies (Kovetz et al., 2017) (including 21cm line emission (Pritchard & Loeb, 2012; Furlanetto et al., 2006; Villaescusa-Navarro et al., 2018)) can also benefit from knowing how the distribution of emitting gas changes as a function of long-wavelength matter density perturbations. In fact, a number of works in these directions have taken place already with Chiang et al. (2017) studying the response of the Lyman-alpha power spectrum to fluctuations in the quasar distribution and with Giri et al. (2019) investigating the response of the 21cm power spectrum signal.

4 Baryonic effects on higher-order lensing statistics

The response approach to perturbation theory developed by Barreira & Schmidt (2017b) is a formalism that uses response functions to rigorously evaluate a number of higher-order $N$-point functions in the nonlinear regime of structure formation. In this section, we use the response approach and the measurements of the baryonic impact (or lack thereof) on the power spectrum response functions discussed in Sec. 3 to predict the impact of baryonic effects on higher-order weak-lensing statistics. We will focus in particular on the case of the squeezed lensing bispectrum and the lensing super-sample power spectrum covariance.

In the considerations below, we will limit ourselves to writing down the relevant equations without derivations; the interested reader is referred to the cited literature for details.

4.1 Lensing equations

4.1.1 The lensing convergence power spectrum

The lensing convergence field $\kappa$ is defined as (see e.g. Hoekstra & Jain (2008); Kilbinger (2015) for reviews on weak lensing)

[TABLE]

where $\bm{\theta}$ is a two-dimensional position on the sky, $\chi$ is the comoving distance (we always leave implicit that $z\equiv z(\chi)$), $\delta_{m}$ is the three-dimensional total matter density contrast, $g(\chi)=(3/2c^{2})H_{0}^{2}\Omega_{m}(1+z)(\chi_{S}-\chi)\chi/\chi_{S}$, $c$ is the speed of light and $\chi_{S}$ is the comoving distance to a single source galaxy redshift $z_{S}$. Under the Limber approximation (which is valid for the multipoles $\ell>20$ we consider below), the angular power spectrum of the lensing convergence is given in terms of the three-dimensional total matter power spectrum as

[TABLE]

where $\ell$ is the multipole moment associated with $\bm{\theta}$ and $k_{\ell}=\ell/\chi$. The above equation makes apparent how the baryonic effects on the three-dimensional matter power spectrum propagate to the two-dimensional lensing one.

4.1.2 The lensing convergence squeezed bispectrum

The three-dimensional total matter bispectrum is defined as $B_{m}(\bm{k}_{1},\bm{k}_{2},\bm{k}_{3})=(2\pi)^{3}\delta_{D}(\bm{k}_{1}+\bm{k}_{2}+\bm{k}_{3})\langle\tilde{\delta}(\bm{k}_{1})\tilde{\delta}(\bm{k}_{2})\tilde{\delta}(\bm{k}_{3})\rangle_{c}$, where $\tilde{\delta}$ is the Fourier transform of the matter density contrast (not to be confused with the meaning of tilded quantities in Sec. 2), $\delta_{D}$ is the Dirac-delta function and the subscript c indicates the correlation function is a connected one; note that the bispectrum depends in general on the relative orientation of the Fourier modes. In the so-called squeezed-limit configuration, the matter bispectrum can be written in terms of responses as

[TABLE]

where $\bm{k}_{h}$ is a small-scale (or hard) Fourier mode, $\bm{k}_{s}$ is a large-scale (or soft) Fourier mode and $\mu_{\bm{k}_{h},\bm{k}_{s}}$ is the cosine angle between the two wavevectors; we have also specialized to isosceles configurations for simplicity, $\bm{k}_{1}=\bm{k}_{2}=\bm{k}_{h}$. Equation (25) is strictly only valid for linear values of $k_{s}$ and $k_{s}\ll k_{h}$666More precisely, there are corrections to Eq. (25) which are of order $(k_{s}/k_{h})^{2}$. Note also that $P_{m}(k_{s},z)$ is effectively the linear matter power spectrum as $k_{s}$ is kept in the linear regime., but the hard mode $k_{h}$ can be deep in the nonlinear regime. Again under the Limber approximation, we can write the bispectrum of the lensing convergence in terms of the three-dimensional matter one as (Cooray & Hu, 2001; Dodelson & Zhang, 2005)

[TABLE]

where $\bm{k}_{\ell_{i}}=\bm{\ell}_{i}/\chi$ ($i=1,2,3$). Specializing to squeezed (isosceles) configurations, $\bm{\ell}_{1}=\bm{\ell}_{2}=\bm{\ell}_{h}$, $\bm{\ell}_{3}=\bm{\ell}_{s}$, as well as averaging over the directions of the modes $\bm{\ell}_{h}$ and $\bm{\ell}_{s}$ permits one to write

[TABLE]

where $R_{\perp}=R_{1}+R_{K}/6$ and the equation is valid up to corrections that scale as $\left(\ell_{s}/\ell_{h}\right)^{2}$. This equation shows that the impact of baryonic effects on the squeezed-lensing bispectrum is controlled by the impact of baryonic effects on the small-scale total matter power spectrum $P_{m}$, as well as on the responses $R_{1}$ and $R_{K}$. In Sec. 3, we have seen that the growth-only response $G_{1}(k,z)$ is effectively independent of baryonic effects, which implies that the baryonic effects on $R_{1}$ are specified by the modifications to $P_{m}$ via the logarithmic derivative term ${\rm dln}P_{m}(k,z)/{\rm dln}k$ in Eq. (4). Regarding the tidal response $R_{K}$ of Eq. (5), we have not explicitly checked whether the corresponding growth-only response $G_{K}$ depends on baryonic effects, but the lack of any impact on $G_{1}$ suggests any such modification on $G_{K}$ should be negligible as well.777The tidal response $R_{K}$ can also be evaluated with separate universe methods, but this requires modifying the N-body/hydro codes to allow them to follow the anisotropic expansion of the local spacetime induced by the long-wavelength tidal field; see Schmidt et al. (2018) for the details of such an implementation in gravity-only simulations. We proceed with this assumption on $G_{K}$, in which case the way $R_{K}$ depends on baryonic effects is specified by the power spectrum according to Eq. (5). Consequently, the whole dependence of the squeezed-lensing bispectrum on baryonic effects can be predicted using the same dependence on the power spectrum.

4.1.3 The lensing covariance matrix

Another interesting and powerful application of the response approach is in the calculation of the covariance matrix of the lensing convergence power spectrum (Barreira et al., 2018b, a). For a given power spectrum estimator $\hat{C}(\ell)$, the covariance is defined as $\operatorname{Cov}_{\kappa}(\ell_{1},\ell_{2})=\langle\hat{C}(\ell_{1})\hat{C}(\ell_{2})\rangle-\langle\hat{C}(\ell_{1})\rangle\langle\hat{C}(\ell_{2})\rangle$, and it can be split into three physically distinct contributions known as the Gaussian term (G), the super-sample covariance term (SSC) and the connected non-Gaussian term (cNG):

[TABLE]

The covariance matrix is a crucial ingredient in parameter inference analyses using weak lensing data, which makes it observationally relevant to ask how baryonic effects on the covariance may affect parameter constraints.

Ignoring a number of real-life complications, such as masking or non-uniform survey depth, the Gaussian term can be written as

[TABLE]

where $\Delta\ell_{1}$ is the size of the multipole bin in which the power spectrum is measured, $\delta_{ij}$ is a Kronecker delta that ensures the result is only non-zero if the two modes $\ell_{1},\ell_{2}$ are in the same bin, $\Omega_{W}=4\pi f_{\rm sky}$ with $f_{\rm sky}$ the surveyed sky fraction and $C_{\rm shape}$ is a shot noise component that arises due to the unknown (i.e. unlensed) shapes of the source galaxies. A key point to note is that the impact of baryonic effects on the Gaussian term is uniquely determined by the way $C(\ell_{1})$ depends on baryonic effects, which we can readily predict using Eq. (24). When we evaluate this term below, we do so for 50 $\ell$ bins equally distributed in log-scale between $\ell_{\rm min}=20$ and $\ell_{\rm max}=5000$, $f_{\rm sky}=0.363$ and $C_{\rm shape}=\sigma_{e}^{2}/2/\bar{n}_{\rm gal}$, with a RMS source galaxy ellipticity of $\sigma_{e}=0.37$ and $\bar{n}_{\rm gal}=30/{\rm arcmin}^{2}$ (these values are inspired by the expected specifications for surveys like Euclid).

Under the Limber approximation, the lensing SSC term can be written as (Takada & Hu, 2013; Barreira et al., 2018b)

[TABLE]

where $\tilde{W}(\bm{\ell})$ is the Fourier transform of the survey footprint on the sky.888Barreira et al. (2018b) has actually shown that the Limber approximation can underestimate the amplitude of the small-scale SSC term by $\approx 10\%$ for $f_{\rm sky}\approx 0.3$. Here, we opt to display and evaluate the equations in the Limber approximation because they are slightly simpler and because they depend on the exact same physical ingredients. Our conclusions on the corresponding impact of baryonic effects thus hold for both derivations. This term is part of the connected four-point function that defines the covariance matrix of a two-point function and its physical meaning is that it describes the coupling between the nonlinear sub-survey modes, $\ell_{1},\ell_{2}$, with linear super-survey modes $\bm{\ell}$ that are integrated over. Similarly to the squeezed lensing bispectrum case, knowing how the power spectrum and its first-order responses $R_{1}$ and $R_{k}$ depend on baryonic effects permits one to straightforwardly evaluate how the SSC term is modified by the same effects. In the numerical results below we consider a disk-shaped survey window function for which $|\tilde{W}(\bm{\ell})|^{2}=\Omega_{W}^{2}\left[2J_{1}(\ell\theta_{W})/(\ell\theta_{W})\right]^{2}$, where $\theta_{W}=\sqrt{\Omega_{W}/\pi}$ and $J_{1}$ is the first-order spherical Bessel function.

Finally, the cNG term corresponds to the rest of the connected four-point function that does not belong to the SSC term; effectively, it describes the coupling between the sub-survey modes $\ell_{1},\ell_{2}$ that arises as the density field becomes less Gaussian distributed. Barreira et al. (2018a) showed recently that for the survey specifications of future (as well as current) surveys, this term is likely to play a negligible role in cosmological analysis using weak-lensing data, and for that reason, we skip evaluating it numerically in this paper. We note for completeness that the majority of this term can be evaluated with second-order power spectrum response functions (Barreira & Schmidt, 2017a).

4.2 Numerical results

4.2.1 Baryonic effects on higher-order lensing statistics

As mentioned above, the dependence of the squeezed-lensing bispectrum and of the lensing SSC term on baryonic effects is encapsulated by the corresponding dependence of the power spectrum. This happens via the power spectrum $P_{m}(k,z)$ itself that enters the above equations, as well as via the logarithmic derivative ${\rm dln}P_{m}(k,z)/{\rm dln}k$ that enters $R_{1}(k,z)$ and $R_{K}(k,z)$ (cf. Eqs. (4) and (5)). Evaluating the latter derivative numerically using the power spectrum measured from a simulation would however yield too noisy a result. To circumvent this, in our numerical results, we evaluate power spectra as

[TABLE]

where $P_{m}^{\rm Emu.}(k,z)$ is the gravity-only power spectrum evaluated using the CosmicEmu emulator (Heitmann et al., 2014, 2016) and $\mathcal{I}^{\rm baryon}(k,z)$ is a baryon boost factor defined as the ratio of the power spectrum from the IllustrisTNG Hydro and Gravity runs

[TABLE]

This ratio is appreciably smoother than the power spectra measurements themselves, which together with the also smooth emulator curves yields sufficiently noise-free derivatives ${\rm dln}P_{m}(k,z)/{\rm dln}k$. We use the TNG300-2 results to evaluate $\mathcal{I}^{\rm baryon}(k,z)$ at the available redshifts, from which we interpolate to carry out the line-of-sight integrals that characterize the calculation of the lensing quantities. We have checked that the CosmicEmu and TNG300-2 gravity-only spectra agree to better than $4\%$ at $z=0,1$ on scales $k<5h/{\rm Mpc}$. We note, however, that the accuracy of the gravity-only spectra is not critical for the interpretation of our results, which focus on the relative impact of baryonic effects and not on absolute spectra values.

The scale-dependence of $\mathcal{I}^{\rm baryon}(k,z)$ is shown in the upper panel of Fig. 4 for $z=0,0.5,1,2$, as labeled. The result shown displays the well known suppression effect due to AGN feedback that becomes important on scales typically inside dark matter haloes $k\gtrsim 1h/{\rm Mpc}$. For instance, on scales $k=3-5h/{\rm Mpc}$, the baryonic effects in the TNG300-2 simulations induce a suppression of power that is of order $10-15\%$ at $z=0.5,1$. Note also the non-trivial interplay between the scale- and redshift-dependence of $\mathcal{I}^{\rm baryon}(k,z)$; for instance, at $k\approx 10h/{\rm Mpc}$ there is a stronger effect at $z=0$ than $z=2$, but the situation is reversed at $k\approx 2-3h/{\rm Mpc}$. These quantitative results are in line with those reported in Springel et al. (2018). The lower panel of Fig. 4 shows the corresponding impact of the baryonic effects on ${\rm dln}P_{m}(k,z)/{\rm dln}k$ ($z=0.5$), which is the quantity that enters explicitly the first-order responses $R_{1}(k,z)$ and $R_{K}(k,z)$ in Eqs. (4) and (5), respectively.

Figure 5 shows the imprint that baryonic effects cast on the lensing power spectrum $C_{\kappa}(\ell_{h})$ (blue), the squeezed-lensing bispectrum $B_{\kappa}(\ell_{h},\ell_{h},\ell_{s})$ (red and green for two values of $\ell_{s}$) and the diagonal of the SSC term $\operatorname{Cov}_{\kappa}^{SSC}(\ell_{h},\ell_{h})$, for a single source redshift $z_{S}=1$, as labeled. The figure shows that the suppression effects on the small-scale three-dimensional power spectrum shown in Fig. 4 manifest themselves also in the line-of-sight-projected lensing quantities shown in Fig. 5. There are some differences in the exact size and scale-dependence of the effect, whose origin can be traced back to two main factors. One is the fact that the line-of-sight integrands of $C_{\kappa}$, $B_{\kappa}$ and $\operatorname{Cov}^{SSC}_{\kappa}$ peak at different redshifts, and hence, are affected differently by the time-dependent baryonic effects on $P_{m}(k,z)$; for example, the integrands of the SSC term tend to peak at lower redshifts (not shown). The other is that the calculation of $B_{\kappa}$ and $\operatorname{Cov}^{SSC}_{\kappa}$ depends also on the responses $R_{1}(k,z)$ and $R_{K}(k,z)$, which are actually enhanced (not suppressed) by the baryonic effects. This is because baryons make the logarithmic derivative ${\rm dln}P_{m}(k,z)/{\rm dln}k$ more negative (cf. Fig. 4), but this derivative enters $R_{1}(k,z)$ and $R_{K}(k,z)$ with a negative sign (cf. Eqs. (4) and (5)). Quantitatively, for our idealized lensing setup, Fig. 5 shows that the SSC term is that which is most strongly affected by the baryonic effects, whose size exceeds $1\%$ ($5\%$) on multipoles $\gtrsim 900$ ($\gtrsim 2000$). The squeezed-bispectrum is that which is the least affected with the same effects exceeding $1\%$ ($5\%$) on multipoles $\gtrsim 2000$ ($\gtrsim 4000$). At a fixed $\ell=3000$, the suppression induced by baryons on the lensing power spectrum, squeezed-bispectrum and SSC term is at the $6\%$, $2-4\%$ and $12\%$ levels, respectively.

We have explicitly checked (although not shown) that these quantitative results for $z_{S}=1$ do not change critically if we choose other source redshift values. For higher $z_{S}$, the impact of baryons gets reduced because the lensing integrands peak at higher redshift, where the impact of baryons on the three-dimensional power spectrum is generically weaker (cf. Fig. 4). This effect is however not very large because of the broad lensing integrands. For instance, the scales above which the impact on the SSC term exceeds $1\%$ are $\ell\approx 500,900,1100$, for $z_{S}=0.5,1,2$, respectively; for the lensing squeezed-bispectrum, the same values are $\ell\approx 1100,2000,2500$. Further, the fact that the lensing SSC term is more strongly affected than the power spectrum, and that the latter is in turn more strongly affected than the lensing squeezed-bispectrum holds equally for $z_{S}=0.5,1,2$. The behavior of the baryonic effects within $z_{S}\in\left[0.5,2\right]$ should be indicative of what is expected from going beyond the idealized single source redshift scenario we adopted here to the redshift distributions of current and future surveys (which span roughly this redshift range).

It is interesting to compare our squeezed lensing bispectrum results to those of Semboloni et al. (2013), who investigated the impact of baryonic effects on lensing three-point statistics as well. They work with the third-order $\langle M^{3}_{\rm ap}(\theta)\rangle$ aperture mass statistic (Pen et al., 2003; Jarvis et al., 2004; Schneider et al., 2005), which is sensitive to general configurations of the lensing bispectrum (not just squeezed, as the case here using the response approach).999See Semboloni et al. (2011b); Simon et al. (2013); Fu et al. (2014) for a number of observational analyses of $\langle M^{3}_{\rm ap}(\theta)\rangle$, none of which explicitly take baryonic effects into account however. They evaluate the three-dimensional bispectrum using an approximate model based on perturbation theory and explore a number of different implementations of the baryonic physical processes present in the OWLS simulations. They also consider more realistic galaxy source distributions, compared to our idealized single source redshift scenario. Semboloni et al. (2013) reported that the baryonic effects can have an impact on the three-point statistic that may be up to $\approx 2$ times larger than that on the two-point statistic (see e.g. Fig. 1 in Semboloni et al. (2013)), and that this distinct sensitivity could be useful in the design of baryon mitigation strategies in real data analyses. Our results in Fig. 5 display, on the other hand, a weaker baryonic impact on $B_{\kappa}(\ell_{h},\ell_{h},\ell_{s})$ than on $C_{\kappa}(\ell_{h})$. We note however that the differences in methodology mentioned above between the two works make it hard to establish a robust comparison.

4.2.2 Baryonic effects on parameter inference using weak-lensing data

The result shown in Fig. 6 is useful to understand the impact that baryonic effects can have on parameter inference analysis using weak-lensing data. The left panel shows the G+SSC lensing covariance matrix (cf. Eq. (4.1.3)) for the idealized lensing setup described in Sec. 4.1; specifically, it shows the ratio of the matrix calculated with baryonic effects to that without baryonic effects taken into account. The upper right panel shows the impact that baryonic effects on the covariance calculation leave on the cumulative signal-to-noise ratio

[TABLE]

for fixed $C_{\kappa}$ calculation with baryonic effects included. Figure 6 shows that taking baryonic effects into account lowers the amplitude of the covariance matrix entries on small-scales, which effectively increases the signal-to-noise of the analysis. For the case of our idealized single source redshift lensing setup, the increase is modest, $\approx 5\%$ for $\ell_{\rm max}=5000$, but an interesting way to interpret this result is that ignoring baryonic effects on the weak-lensing covariance calculation represents an overestimation of the total statistical error budget, and it is thus conservative from the point of view of inferred parameter error bars.

In addition to controlling the size of the error bars on parameters, a miscalculation of the covariance matrix can in general also bias the central value of the parameter constraints. To estimate this effect we have drawn $1000$ realizations of a convergence power spectrum data vector $C^{{\rm data}}_{\kappa}(\ell_{i})$ (where $i$ runs over the $50$ multipole bins) from a multivariate Gaussian distribution with a mean $C_{\kappa}(\ell)$ and $\operatorname{Cov}_{\kappa}(\ell_{1},\ell_{2})$ (G+SSC) evaluated with baryonic effects taken into account. For each sampled data vector we have evaluated the $\chi^{2}$ quantity

[TABLE]

for three distinct cases:

Baryons in signal and covariance: both $C_{\kappa}$ and $\operatorname{Cov}_{\kappa}$ are calculated with baryonic effects taken into account; 2. 2.

No baryons in the covariance: the calculation of $C_{\kappa}$ takes baryonic effects into account, but the calculation of $\operatorname{Cov}_{\kappa}$ does not; 3. 3.

No baryons in the signal: the calculation of $C_{\kappa}$ does not take baryonic effects into account, but the calculation of $\operatorname{Cov}_{\kappa}$ does.

The lower right panel of Fig. 6 shows the histogram of the difference between the $\chi^{2}$ of cases (ii) and (iii) above to that of case (i). Ignoring the impact that baryonic effects can have on the signal (red) induces an appreciable degradation of the goodness-of-fit $\Delta\chi^{2}/{dof}\approx 22$. In real data analysis, this would trigger shifts in the cosmological parameter values away from their true values to attempt minimizing the $\chi^{2}$ value. This well-known result is what justifies the many efforts currently being undertaken to incorporate baryonic physics in parameter inference analysis using weak-lensing spectra (Semboloni et al., 2011a; Schneider & Teyssier, 2015; Eifler et al., 2015; Mead et al., 2015; Harnois-Déraps et al., 2015; Chisari et al., 2018; Huang et al., 2018; Schneider et al., 2018). On the other hand, neglecting baryonic effects on the covariance matrix results in an appreciably smaller change in the goodness-of-fit $\Delta\chi^{2}/{dof}\approx-0.02$. This is depicted by the distribution shown in black in the lower right panel of Fig. 6, which is plotted with a different x-axis scale to permit visualizing the shape of the distribution.

We note that these considerations, although indicative of a small practical impact of baryonic effects in weak-lensing covariance matrices, refer to our idealized lensing setup. More robust, quantitative and survey-specific conclusions should be drawn from more realistic simulated likelihood analyses including lensing tomography, realistic source galaxy distributions and systematic effects (intrinsic alignments, photo-$z$, etc.), as well as measuring the impact of baryonic effects at the level of parameter constraints. The fact that baryonic effects, in particular AGN feedback, are expected to lower the small-scale entries of lensing covariance matrices is however robust to the exact choice of the source redshift distribution.

5 Summary and Conclusions

In this paper, we have used the separate universe simulation technique applied to the IllustrisTNG hydrodynamical galaxy formation model to measure the impact of baryonic physical processes on matter power spectrum response functions. The responses specify the scale- and redshift-dependence of the response of the small-scale matter power spectrum to the presence of long-wavelength density and tidal field fluctuations. The separate universe formalism is a technique that incorporates these long-wavelength modes via a redefinition of the cosmology that is actually simulated (cf. Sec. 2.2).

We focused on measurements of the first-order isotropic growth-only power spectrum response function $G_{1}(k,z)$, which together with the matter power spectrum $P_{m}(k,z)$, fully specifies the first-order isotropic response $R_{1}(k,z)$ (cf. Eq. (4)). Armed with the knowledge of the baryonic effects on $P_{m}(k,z)$, $R_{1}(k,z)$, as well as on the first-order tidal response $R_{K}(k,z)$, one can then use the machinery of the response approach to perturbation theory to straightforwardly evaluate how baryonic effects impact a number of higher-order matter/lensing correlation functions in the nonlinear regime of structure formation (cf. Sec. 4.1).

We have carried out separate universe simulations on the TNG300 box ($L_{\rm box}=205h/{\rm Mpc}\approx 300{\rm Mpc}$) at two particle/mass element resolutions: $N_{p}=1250^{3}$ (TNG300-2) and $N_{p}=625^{3}$ (TNG300-3). For each, we have run hydrodynamical simulations with the full TNG physics model (which we dubbed Hydro runs), as well as gravity-only counterparts (which we dubbed Gravity runs). While the cosmological parameters of the simulation are modified to mimic the presence of long-wavelength density perturbations, the remaining parameters of the TNG galaxy formation model are held fixed. Our measurements describe therefore the response of the matter power spectrum at fixed galaxy formation physics model.

Our main results can be summarized as follows:

• •

The growth-only total matter power spectrum response $G_{1}(k,z)$ is independent of baryonic effects on redshifts $z<3$ and scales $k\lesssim 15h/{\rm Mpc}$ (cf. Fig. 1). This means that the dependence of the full isotropic response $R_{1}(k,z)$ on baryonic effects is totally encoded in the term $\propto{\rm dln}P_{m}(k,z)/{\rm dln}k$ (cf. Fig. 4 and Eq. (4)).

• •

The $G_{1}(k,z)$ response of the power spectrum of the dark matter component alone is nearly indistinguishable of the response of the total matter fluid (cf. Fig. 3). Our results suggest, however, that the hydrodynamical physical processes felt by the gas cells (e.g. loss of mass to form stars or fuel black hole growth) lower the amplitude of its $G_{1}(k,z)$ response w.r.t. those of the total and dark matter components. These differences are also more pronounced at higher redshift (cf. Fig. 3).

• •

The exact size and scale-dependence of the impact of baryonic physics on the lensing squeezed-bispectrum and super-sample covariance term can differ from that of the power spectrum (cf. Fig. 5). For instance, for a single source redshift of $z_{S}=1$ and at $\ell=3000$, baryonic effects suppress the amplitude of the lensing power spectrum by $\approx 6\%$, the lensing squeezed-bispectrum by $\approx 2-4\%$ and the diagonal of the SSC term by $\approx 12\%$.

• •

Skipping the incorporation of baryonic effects in the calculation of weak-lensing covariance matrices overestimates its amplitude on small scales, and is thus conservative from the point of view of inferred parameter error bars. In other words, gravity-only covariances yield artificially larger parameter errors, even if only by small amounts (cf. Fig. 6). Further, ignoring baryonic effects on the covariance has a negligible impact on the $\chi^{2}$ goodness-of-fit, compared to the impact from ignoring the same effects on the modeling of the signal.

The steps that we have taken in this paper to obtain our results can prove useful in the design of strategies to incorporate or mitigate the impact of baryonic physics in real data analyses. This is the case, for instance, at the level of the modeling of the signal in cosmological inferences using the lensing bispectrum, as well as at the level of current observational data analyses using the lensing power spectrum in the evaluation of its covariance matrix. A noteworthy advantage of using the response approach to perturbation theory is that it permits one to study higher-order statistics in the nonlinear regime with nearly the same numerical costs as two-point function studies.

Finally, beyond the study of first-order matter power spectrum responses, there are a number of additional interesting investigations that can be carried out with separate universe simulations of the IllustrisTNG model. An example is the exploration of larger values of the amplitude $\delta_{L0}$ of the long-wavelength mode to measure higher-order power spectrum response functions; these find applications in the calculation of the cNG term of the matter power spectrum covariance (Barreira & Schmidt, 2017a) and matter bispectrum covariance (Barreira, 2019). Another example is the study of galaxy bias as the response of the galaxy number density to the presence of long-wavelength modes (Lazeyras et al., 2016; Li et al., 2016; Baldauf et al., 2016); this would be complementary to the galaxy bias measurements performed already by Springel et al. (2018) using the large-scale limit of galaxy clustering. Similarly, the study of galaxy power spectrum responses with realistic galaxy formation processes becomes also possible and it can provide interesting insights on galaxy power spectrum covariances, as well as on the squeezed galaxy bispectrum that is a relevant observable in studies of primordial non-Gaussianity.

Acknowledgements

We would like to thank Anz̆e Slosar for useful comments and discussions. The simulations used in this work were run on the Cobra supercomputer at the Max Planck Computing and Data Facility (MPCDF) in Garching near Munich. AB and FS acknowledge support from the Starting Grant (ERC-2015-STG 678652) "GrInflaGal" of the European Research Council.