Response function of the large-scale structure of the universe to the small scale inhomogeneities

In order to infer the impact of the small-scale physics to the large-scale properties of the universe, we use a series of cosmological $N$-body simulations of self-gravitating matter inhomogeneities to measure, for the first time, the response function of such a system defined as a functional derivative of the nonlinear power spectrum with respect to its linear counterpart. Its measured shape and amplitude are found to be in good agreement with perturbation theory predictions except for the coupling from small to large-scale perturbations. The latter is found to be significantly damped, following a Lorentzian form. These results shed light on validity regime of perturbation theory calculations giving a useful guideline for regularization of small scale effects in analytical modeling. Most importantly our result indicates that the statistical properties of the large-scale structure of the universe are remarkably insensitive to the details of the small-scale physics, astrophysical or gravitational, paving the way for the derivation of robust estimates of theoretical uncertainties on the determination of cosmological parameters from large-scale survey observations.

Wide field galaxy surveys are widely considered for unveiling the detailed geometrical properties or energy content of the universe [1]. Large-scale projects, such as the EUCLID mission [26], are planned in the coming decade, aiming at the determination of these properties with an unprecedented accuracy. Such measurements rely to a large extent on the use of the statistical properties of the large-scale cosmic structures up to scales entering the weakly non-linear regime, where the sole linear theory cannot be used. Such a scientific program could then only be achieved if the properties of the large-scale cosmological structure can be safely predicted either from numerical simulations or from analytical investigations for any given cosmological model. In particular it is important that such observables are shielded from the details of small scale astrophysics and gas physics at galactic or sub-galactic scales.
One way to reformulate this question is to quantify how small-scale structures can impact the growth on large scales as soon as modes are entering the nonlinear regime. Perturbation theory (PT) of the structure formation is a powerful framework to precisely predict the nonlinear gravitational dynamics of the cosmic fluid from the first principle, at least when gravity only is at play (see [2] for a review). The importance of such methods has been heightened after the detection of the baryon acoustic oscillations (BAOs) in the clustering of galaxies at late times (e.g., [3]), making precise predictions of the nonlinear matter power spectrum crucially important.
PT calculations show precisely that mode couplings between different scales are unavoidable. This makes these calculations in general difficult to develop in a controlled manner. We propose here to quantify such couplings with the use of a two-variable kernel function [27], defined as the linear response at wave mode k with respect to an initial perturbation of the linear power spectrum at wave mode q. In the context of PT calculations, Ref. [4] showed progressive broadening of the kernel function as increasing the PT order, and speculated that a regularization scheme in the UV domain is required to give a realistic estimate of the high-order PT contributions. The recent paper by [5] also pointed out the unsuccessful convergence of the PT series at late times and proposed a simple ansatz based on the Padé approximation to suppress the strong UV sensitivity seen in the standard PT (SPT).
If the broadness of the kernel at late times suggested from PT is true, physics at very small scale can influence significantly the matter distribution on large scales, where the acoustic feature is prominent. It also poses a question to the reliability of simulations, with which we can follow the evolution of Fourier modes only in a finite dynamic range. We here present direct measurement of the kernel structure from cosmological N -body simulations. We show that this allows a direct test of regularization schemes employed in analytical models.
Definition and methodology.-Here we wish to introduce a well-defined kernel function and investigate it at fully nonlinear level. We consider the nonlinear power spectrum as a functional of the linear counterpart, i.e., P nl = P nl [P lin ], and define the kernel function as its functional derivative: We omit the explicit dependence on z from the arguments in what follows. The normalization for K is chosen such that it contributes to the change in the nonlinear power spectrum with uniform weights per decade. We can construct a simple estimator of the kernel function from numerical simulations based on this definition. We prepare two initial conditions with small modulations in the linear power spectrum over a finite interval of wave mode q, evolve them to a late time, and take the difference of the power spectra measured from the two. That iŝ where the two perturbed linear spectra are given by In the above, the index j runs over the initial wave-mode bins, and we choose log-equal binning, ln q j+1 − ln q j = ∆ ln q. The other index i is used for the final wave-mode bins, which we set identically to the initial one. It is straightforward to show that the estimatorK approaches to the kernel function K defined in Eq. (1), when ∆ ln q and ∆ ln P lin are small. The definition (1) is advantageous in that it allows the measurement in this way at the fully nonlinear level [28]. Note that the kernel function in a similar matrix form was first discussed in Ref. [6]. The usefulness of such an object was shown in understanding the mode transfer under various local transformations of the cosmic density field. Numerical analysis: We adopt a flat-ΛCDM cosmology consistent with the five-year observation by the WMAP satellite [7] with parameters (Ω m , Ω b /Ω m , h, A s , n s ) = (0.279, 0.165, 0.701, 2.49 × 10 −9 , 0.96), which are the current matter density parameter, baryon fraction, the Hubble constant in units of 100km/s/Mpc, the scalar amplitude normalized at k 0 = 0.002Mpc −1 and its power index, respectively. The matter transfer function is computed with these parameters using the CAMB code [8].
We run three sets of simulations with different volume and number of particles as listed in Table I. They are intended to test the convergence of the measured kernel function. The initial conditions are created using a parallel code developed in [9,10] based on the second-order Lagrangian PT (e.g., [11,12]) at the optimal redshifts following Ref. [13]. We evolve the matter distribution using Gadget2 [14] with the Tree-PM calculation. We finally measure the power spectrum by fast Fourier transform of the Cloud-in-Cell (CIC) density estimates on 1024 3 mesh with the CIC kernel deconvolved in Fourier space.
For each set of simulations, we prepare multiple initial conditions with linear power spectra perturbed by ±1% over q j ≤ q < q j+1 . The q-bins start at q 1 = 0.006h Mpc −1 (0.012h Mpc −1 ) for L10-N9 (L9-N9 and L9-N8) with the bin width ∆ ln q = ln( √ 2). We consider 15 or 13 bins depending on the simulation set as listed in"bins" column of Table I. We run four random realizations for L9-N9 and L9-N8 to estimate the statistical scatter, and the initial conditions with perturbed spectra at different bins are created with exactly the same random phases for each realization. The total numbers of runs used in this analysis are also shown in Table I. Shape of the kernel and comparison with PT.-We are now in position to present the kernel function measured from N -body simulations. The combination K(k, q)P lin (q) is plotted at a fixed k shown by the vertical arrow as a function of q in Fig. 1. The heavy overlap among the three sets of simulations ensures that the results are converged against the resolution and volume. We hereafter discuss the results of L9-N9, which has the best spatial resolution. At high redshifts, we can see a strong peak at k = q arising from the trivial linear calculation. Nonlinear coupling then gradually grows with time and the peak feature gets less significant. One of the key features that can be observed here is that a large contribution comes from small wave modes (q < k) suggesting that the growth of structure is dominated by mode flows from large to small scales. Not surprisingly, the formation of structure is more effectively amplified when it is part of a larger structure than when it contains small scale features. showing the perturbative order of the calculation. We show a negative sign in the legend when the contribution is negative. Note that we ignore terms proportional to the Dirac delta function at k = q, which is meaningful only when binning is considered.
Such findings are fully in line with expectations from PT calculations. We show the analytical calculation in Fig. 2 up to the two-loop level (ignoring at this stage binning effects). We present the contribution from P ij (k) ∝ δ (i) δ (j) , where δ (i) is the ith-order term in the PT expansion. The terms in the same loop order cancel at the IR domain (q < k) due to the extended galilean invariance of the motion equations as shown and analyzed in e.g., [15][16][17][18][19]. On the other hand, the UV domain is entirely dominated by P 13 (k) and P 15 (k) at one and two loops, respectively. Such terms can be alternatively described as the correction to the density propagator. They have been shown to dominate the behavior of the UV domain at any fixed order in SPT.
We then rescale the kernel at various redshifts as T (k, q) = [K(k, q) − K lin (k, q)]/[qP lin (k)], where K lin is the trivial linear contribution, and plot them in Fig. 3. They are compared with the one-loop PT calculation (solid), which is now time-independent. The simulation data indeed shows little time evolution at q k in striking agreement with the PT predictions, reproducing the expected q dependence and amplitude [29] obtained from the one-loop calculation and the change of sign one expects between the IR and the UV domain. The small but non-negligible z-dependence at k ∼ q is further accounted for by the two-loop calculation (see the figure legend for line types). Note that at the final wave mode plotted here (i.e., k = 0.161 h Mpc −1 ), the two-loop SPT prediction for the nonlinear power spectrum agrees with simulations within 1% at z 1 and the agreement gets worse at lower redshift reaching to ∼ 5% at z = 0 (see e.g., [5]).
At q 0.3 h Mpc −1 however, the measured kernel function is observed to be damped compared to perturbation theory predictions at one or two-loop order. As can be seen on Fig. 3, the one-loop SPT (solid line) predicts the kernel function to reach a constant [30]; at the two-loop order, it is expected to grow in amplitude with time. The numerical measurements show on the other hand that the scaled kernel function is strongly damped with decreasing redshift. It is such that the couplings between scales take place effectively between modes of similar wavelengths. This effect is particularly important at late time. At redshift zero, the departure between two-loop predictions and numerical results is striking. Furthermore analysis of the kernel structure at three and higher loop order (see e.g., [4]) suggests that SPT calculations, taken at any finite order, predict an even larger amplitude of the kernel function in the high q region. It strongly suggests that this anomaly is genuinely non-perturbative.
We propose an effective description of this observed behavior. As illustrated on Fig. 4 we found that this can be accounted for with a Lorentzian shape: T eff. (k, q) = T 1−loop (k, q) + T 2−loop (k, q) 1 1 + (q/q 0 ) 2 (4) characterized by a time-dependent critical wave mode, q 0 (z) = 0.3D −2 + (z)h/Mpc, where D + is the linear growth rate of the fluctuations. Note that, as it can be checked on Fig. 4, q 0 shows no dependence on k preserving the k dependence of the UV part of the kernel function. This dependence is in full agreement with PT predictions. Discussion-. What is the origin of this behavior? We have not yet identified its origin. In particular it is not clear whether it can be explained in the context of the physics of single stream dynamics or if it is associated with genuine shell crossing effects. The latter seems the more natural explanation. It could then be accounted for with the help of Effective Field Theory (EFT) approaches advocated initially in [20,21] (see also [22]). It is possible however that such damping effect originates from simpler mechanisms in single-stream physics. It has been shown in particular that the nonlinear density propagator, which expresses how a given wave mode is evolving with time, is exponentially damped by the large-scale displacements. This is the standard result on which the Renormalized Perturbation Theory approach is based [23,24]. As explicitly shown in [25] equal time power spectra are however insensitive to displacements of the global system, that is originating from wave modes smaller than k. Displacements at intermediate scales are however expected to induce some effective damping for equal time spectra. The physical idea behind such a mechanism is that the force driving the collapse of a large-scale perturbation (e.g., a cluster of galaxies) is affected by the small scale inhomogeneities within the structure (say galaxies), but that this dependence might be damped when such small scale inhomogeneities are actually moving within the structure. This is however beyond the scope of this presentation to evaluate the importance of this effect.
Summary-. We have presented the direct measurement of the kernel function that governs the dependence of the nonlinear power spectrum on the initial spectrum in the course of cosmic structure formation. This measurement was done using a large ensemble of N -body simulations that differ slightly in their initial conditions. It was found to be robust with respect to the simulation resolution so we believe it reveals genuine physical processes that take place along the development of cosmic gravitational instabilities.
Our results are to a large extent in agreement with SPT predictions at one or two loop order. However we find an anomalous strong damping of the function at the UV domain which cannot be explained from such calculations. Although we do not have identified the origin of such an effect, we point out several explanations that could be investigated in the future, either from shell crossings or relative displacement effects.
Irrespectively of its actual origin, our finding supports the idea that gravity is localized in scale, i.e., Fig. 3 shows that the amplitude of the power spectrum of wave mode k depends only on the initial power spectrum at wave modes q close to k. It supports the idea that the nonlinear density power spectrum is an accurate tracer of the linear density power spectrum at similar scales, with the small scale physics -which is also driven by complex non-gravitational processes and therefore cannot be controlled as effectively -being screened away through this damping mechanism. It eventually helps us justify the use of power spectra as probes of fundamental cosmological parameters and the use of either perturbation theory calculations or large-scale N-body simulations to calibrate such probes.
We thank Patrick Valageas for fruitful discussions on analytical calculations of the kernel function. This works is supported in part by grant ANR-12-BS05-0002 of the French Agence Nationale de la Recherche. TN is supported by JSPS. AT is supported by a Grant-in-Aid for Scientific Research from the Japan Society for the Promotion of Science (JSPS) (No. 24540257). Numerical calculations for the present work have been carried out on Cray XC30 at Center for Computational Astrophysics, CfCA, of National Astronomical Observatory of Japan.