Fourier imaging of non-linear structure formation

TL;DR

This paper introduces a Fourier space decomposition technique to analyze non-linear structure formation in cosmology, revealing mode coupling, power transfer, and vorticity effects from N-body simulations, aiding model calibration.

Contribution

It presents a novel Fourier decomposition method to study non-linear mode coupling and power transfer in cosmological structure formation from N-body simulations.

Findings

Identifies power transfer from large to small scales.

Detects the stable clustering regime.

Highlights the importance of vorticity at certain scales.

Abstract

We perform a Fourier space decomposition of the dynamics of non-linear cosmological structure formation in LCDM models. From N-body simulations involving only cold dark matter we calculate 3-dimensional non-linear density, velocity divergence and vorticity Fourier realizations, and use these to calculate the fully non-linear mode coupling integrals in the corresponding fluid equations. Our approach allows for a reconstruction of the amount of mode coupling between any two wavenumbers as a function of redshift. With our Fourier decomposition method we identify the transfer of power from larger to smaller scales, the stable clustering regime, the scale where vorticity becomes important, and the suppression of the non-linear divergence power spectrum as compared to linear theory. Our results can be used to improve and calibrate semi-analytical structure formation models.

Tables

Table 1. Table 1 : The table shows parameters and initialization methods for the N 𝑁 N -body simulations used in this work. R box subscript 𝑅 box R_{\rm box} denotes the box size in units of Mpc / h Mpc ℎ {\rm Mpc}/h , N part subscript 𝑁 part N_{\rm part} denotes the number of N 𝑁 N -body particles, z i subscript 𝑧 𝑖 z_{i} is the initial N 𝑁 N -body starting redshift, where the initial conditions are generated with either the Zel’dovich Approximation (ZA) or with an added term from second order Lagrangian perturbation theory (2LPT).

Sim	$R_{\mathrm{box}}$	$N_{\mathrm{part}}$	$z_i$	ZA/2LPT
A	64	${512}^3$	199	ZA
B	64	${512}^3$	49	ZA
C	64	${512}^3$	49	2LPT
D	256	${512}^3$	49	ZA
E	1024	${512}^3$	49	ZA
F	4096	${512}^3$	49	ZA
G	256	${1024}^3$	49	ZA

Equations

$$d{s}^2=g_{\mu \nu }d{x}^{\mu }d{x}^{\nu }=a^2(\tau )[-(1+2\psi )d{\tau }^2+(1-2\varphi ){\delta }_{ij}d{x}^{i}d{x}^j].$$

$$T_{\nu ;\mu }^{\mu }={\partial }_{\mu }{T}_{\nu }^{\mu }+{\Gamma }_{\lambda \mu }^{\lambda }{T}_{\nu }^{\mu }-{\Gamma }_{\mu \nu }^{\lambda }{T}_{\lambda }^{\mu }=0.$$

$$T_0^0$$

$$T_0^i$$

$$T_j^i$$

$$\dot{\delta }=$$

$$-\theta \delta -u^i{\partial }_i\delta$$

$$+3\left(1+\frac{\delta P}{\delta \rho }\right)\dot{\varphi }\delta -\frac{\delta P}{\delta \rho }\theta \delta -u^i{\partial }_i\left(\frac{\delta P}{\delta \rho }\delta \right)$$

$$-({\partial }_i\psi -3{\partial }_i\varphi )\left(1+\delta +w+\frac{\delta P}{\delta \rho }\delta \right)u^i.$$

$$\dot{\delta }=-(1+\delta )\theta -u^i{\partial }_i\delta ,$$

$${\dot{u}}^i=$$

$$-\frac{1}{\bar{\rho }\left(1+\delta +w+\frac{\delta P}{\delta \rho }\delta \right)}\left[{\delta }^{ik}({\partial }_j+{\partial }_j\psi -3{\partial }_j\varphi )+{\delta }_j^k{\delta }^{il}{\partial }_l\varphi \right]T_k^j.$$

$${\tilde{T}}_j^i=$$

$$(\rho +P){\sigma }_j^i=T_j^i-{\tilde{T}}_j^i.$$

$${\dot{u}}^i=-\frac{\dot{a}}{a}{u}^i-u^j{\partial }_j{u}^i-{\delta }^{ij}{\partial }_j\psi -\frac{{\delta }^{ik}}{\rho }{\partial }_j\left(\rho {\sigma }_k^j\right).$$

$$q_{\theta }\equiv$$

$$q_w^i\equiv$$

$$\dot{\theta }=-\frac{\dot{a}}{a}\theta -{\partial }_i{u}^j{\partial }_j{u}^i-u^i{\partial }_i\theta -{\nabla }^2\psi -q_{\theta }.$$

$${\dot{w}}^i=-\frac{\dot{a}}{a}{w}^i-{\delta }_{km}{\epsilon }^{ilm}{\partial }_l(u^j{\partial }_j{u}^k)-q_w^i.$$

$$\delta (\mathbf{k})$$

$$\delta (\mathbf{x})$$

$$\alpha (\mathbf{x})\beta (\mathbf{x})={\mathbf{F}}^{-1}\left\{\mathbf{F}[\alpha (\mathbf{x})]\ast \mathbf{F}[\beta (\mathbf{x})]\right\},$$

$$\alpha (\mathbf{k})\ast \beta (\mathbf{k})=\int d^3{k}^{\prime }\alpha ({\mathbf{k}}^{\prime })\beta (\mathbf{k}-{\mathbf{k}}^{\prime })=\int d^3{k}^{\prime }\alpha (\mathbf{k}-{\mathbf{k}}^{\prime })\beta ({\mathbf{k}}^{\prime }),$$

$$\mathbf{F}[I_{\delta }(\mathbf{x})]=$$

$$=$$

$$\dot{\delta }(\mathbf{k})=-\theta (\mathbf{k})-I_{\delta }(\mathbf{k}),$$

$$I_{\theta \delta }(\mathbf{k})=\int d^3{k}^{\prime }\frac{\mathbf{k}\cdot {\mathbf{k}}^{\prime }}{k^{\prime 2}}\theta ({\mathbf{k}}^{\prime })\delta ({\mathbf{k}}^{\prime \prime }),$$

$$I_{w\delta }(\mathbf{k})=-\int d^3{k}^{\prime }\frac{{\mathbf{k}}^{\prime \prime }\times {\mathbf{k}}^{\prime }}{k^{\prime 2}}\cdot \mathbf{w}({\mathbf{k}}^{\prime })\delta ({\mathbf{k}}^{\prime \prime }).$$

$$\dot{\theta }(\mathbf{k})=-\frac{\dot{a}}{a}\theta (\mathbf{k})+k^2\psi (\mathbf{k})-q_{\theta }(\mathbf{k})-I_{\theta }(\mathbf{k}),$$

$$I_{\theta \theta }(\mathbf{k})=\frac{1}{2}{k}^2\int d^3{k}^{\prime }\frac{{\mathbf{k}}^{\prime }\cdot {\mathbf{k}}^{\prime \prime }}{k^{\prime 2}{k}^{\prime \prime 2}}\theta ({\mathbf{k}}^{\prime })\theta ({\mathbf{k}}^{\prime \prime }).$$

$$I_{\theta \mathbf{w}}(\mathbf{k})=-\int d^3{k}^{\prime }\frac{(\mathbf{k}+{\mathbf{k}}^{\prime \prime })\cdot {\mathbf{k}}^{\prime }}{k^{\prime 2}{k}^{\prime \prime 2}}\theta ({\mathbf{k}}^{\prime }){\mathbf{k}}^{\prime }\times {\mathbf{k}}^{\prime \prime }\cdot \mathbf{w}({\mathbf{k}}^{\prime \prime }),$$

$$I_{\mathbf{ww}}(\mathbf{k})=\int d^3{k}^{\prime }\frac{{\mathbf{k}}^{\prime \prime }\times {\mathbf{k}}^{\prime }\cdot \mathbf{w}({\mathbf{k}}^{\prime })}{k^{\prime 2}}\frac{{\mathbf{k}}^{\prime }\times {\mathbf{k}}^{\prime \prime }\cdot \mathbf{w}({\mathbf{k}}^{\prime \prime })}{k^{\prime \prime 2}}.$$

$$i_{\alpha }^{\beta }(k,k^{\prime })=\frac{\int d\Omega s^{\beta }(\mathbf{k})\cdot \int d{\Omega }^{\prime }{\tilde{i}}_{\alpha }(\mathbf{k},{\mathbf{k}}^{\prime })}{2\int d\Omega },$$

$$s^{\beta}(\mathbf{k})=\left(\begin{array}[]{c}{\rm sgn}(\beta_{re}(\mathbf{k}))\\ {\rm sgn}(\beta_{im}(\mathbf{k}))\end{array}\right),$$

Full text

1 Introduction

Cosmological structure formation can in principle be tracked exactly by solving the coupled Einstein-Boltzmann equations for all species. For species which are non-relativistic (such as cold dark matter and baryons) it is advantageous to solve the hierarchy of momentum integrated Boltzmann equations, i.e. the continuity equation, the Euler equation, and corresponding equations for the higher velocity moments, and using the fact that for non-relativistic species the higher velocity moments are strongly suppressed. For cold dark matter it is a very good approximation to truncate the hierarchy equations and simply solve the continuity and Euler equations.

In linear perturbation theory these equations are easily solved, but on smaller scales (d\mbox{\;\raisebox{1.29167pt}{<}!!!!!\raisebox{-3.87495pt}{\sim}\;}100-200 Mpc) structures have evolved into the non-linear regime and perturbation theory breaks down. This means that terms which are quadratic in the fluid variables must be taken into account. Methods based on diagram resummation [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] or higher order perturbation theory [12, 13, 14, 15, 16, 17, 18] have been shown to be fairly accurate beyond the scales where linear theory breaks down. However, they are only applicable in the semi-linear regime and fail once the truly non-linear regime is entered. Alternatively, semi-analytic methods such as HALOFIT [19] use knowledge of the qualitative behaviour of quantities such as the power spectrum in the linear and non-linear regimes as well as a number of fitting parameters to interpolate between the two regimes. These parameters in turn must be calibrated using $N$-body simulations, and using the best current methods an accuracy of a few percent can be achieved for $\Lambda$CDM-like models (see e.g. [20]).

Even in the simplest possible case of a pure Cold Dark Matter (CDM) density field, the problem of non-linear structure formation is highly challenging. Formally, the problem is to solve the Boltzmann equation for a collisionless fluid in a metric provided by the Einstein equations. On scales much smaller than the horizon this extremely hard problem can be well approximated by Newtonian theory in an expanding FRW background. Since CDM is assumed to be cold, it can be approximated by a non-viscous fluid with $c_{s}=0$ and the Boltzmann equation can then be reduced to the continuity and Euler equations, i.e. the first two moment-equations of the Boltzmann hierarchy.

A fully non-linear calculation directly in Fourier space requires one to perform a 3-dimensional integral for each wavenumber. The time required to perform structure formation simulations with this method is several orders of magnitude larger than what can currently be considered feasible. The most successful approach to solving these equations is instead to use $N$-body simulations (i.e. a coarse grained particle representation of phase-space) and simply solve the Newtonian equation of motion for all particles. This approach can be mainly thought of as a real space solution of the fluid equations, although fast Fourier transforms are used for calculating long-range gravitational forces.

The next generation of large-scale structure surveys such as EUCLID will achieve a measurement of observables like the shear power spectrum to almost cosmic variance precision over a wide range of scales. In order to extract information on the underlying cosmological parameters this in turn requires quantities such as the matter power spectrum to be calculated to 1% level precision (see e.g. [21] for a recent discussion). While this is achievable for a few different models it is currently impossible to perform scans over multi-dimensional parameter spaces using $N$-body simulations. The theoretical understanding and numerical implementation of non-linear cosmological structure formation models is therefore significantly challenged by the progress in observational cosmology. To confront this challenge, one needs a deep understanding of the dynamics of cosmological mode coupling and thereby how non-linear structures form in the Universe.

The purpose of the present study is to gain insights into the structure of the continuity and Euler equations in Fourier space. Using $N$-body simulations we calculate exactly the non-linear source terms in these equations and identify the various regimes, from purely linear structure formation on large scales to the stable clustering regime governing the smallest scales. We propose that these insights can be used to construct new semi-analytic methods for calculating quantities such as the matter power spectrum.

Since we neglect higher moments in the Euler equation, our calculations do not contain the full picture of structure formation at scales close to virialization, $k\sim 1h/{\rm Mpc}$ at $z\sim 0$. At these scales, anisotropic stress must be taken into account.

The high dimensionality of the problem requires $N$-body simulations to follow the time-evolution, but it is possible at specified redshifts to perform a full non-linear calculation in Fourier space, and thereby probe the dynamics of cosmological mode coupling.

We note that [22] and [23] also investigate the effect of mode coupling using $N$-body simulations, but in a very different way from the one presented in this paper.

The paper is structured as follows. In section 2 we derive the non-linear continuity and Euler equations in Fourier space from relativistic energy-momentum conservation. In section 3 we then explain how we have implemented the calculation of the mode coupling integrals in these equations, and in section 4 we present our results. Finally, section 5 contains our conclusions.

2 The non-linear fluid equations

The following theoretical outline is performed in the conformal Newtonian gauge [24], where deviations from flat space-time are parametrized in terms of the metric potentials $\phi$ and $\psi$, $a$ is the scale factor, and $\tau$ is the proper time. The line element is given by

[TABLE]

The relativistic evolution of the energy-momentum tensor, $T^{\mu}_{\nu}$, is determined by requiring its covariant derivative to vanish

[TABLE]

The 10 independent components in the symmetric energy-momentum tensor can be used to define 10 fluid variables, $\delta$, $u^{i}$, $\delta P/\delta\rho$ and $\Sigma^{i}_{j}$, as follows

[TABLE]

where $w=\bar{P}/\bar{\rho}$, $\bar{\rho}$ and $\bar{P}$ is the average density and pressure, respectively, $\delta^{i}_{j}$ is the Kronecker delta and $\Sigma_{i}^{i}=0$.

2.1 Real space equations

Below we linearize the equations with respect to the metric perturbations, $\psi$ and $\phi$, since they are very small for all scales relevant to cosmology. The 10 fluid variables are not assumed to be small, and the equations are therefore non-linear in these quantities. After deriving the full general relativistic equations we identify the terms relevant for studying mode coupling for a cold fluid in a Newtonian setting.

The continuity equation in real space and its Newtonian limit

Setting $\nu=0$ in Eq. (2.2) and defining the velocity divergence $\theta=\partial_{i}u^{i}$ we get the continuity equation

[TABLE]

The terms in the first line are linear terms, the ones in the second line are Newtonian non-linear terms, while the third line includes pressure and non-linear horizon terms. The terms in the last line can be neglected since $(\partial_{i}\psi-3\partial_{i}\phi)u^{i}=\partial_{i}(\psi u^{i}-3\phi u^{i})-(\psi-3\phi)\theta\ll\theta$. Here we have assumed that the metric potentials are much less than unity, and that, e.g., $\partial_{i}(\psi u^{i})\ll\partial_{i}u^{i}$, since the general relativistic version of the Poisson equation constrains the relation between spatial derivatives of the metric potentials and the fluid variables.

The non-linear Newtonian limit of the continuity equation is given by

[TABLE]

where pressure has been neglected, as well as time-derivatives of the metric perturbations.

The Euler equation in real space and its Newtonian limit

Setting $\nu=i$ in the energy-momentum conservation equations, Eq. (2.2), one arrives at the Euler equation

[TABLE]

As was the case for the continuity equation, we assume that the metric potentials do not vary extremely rapidly, and therefore that the flow velocities are small compared to the speed of light. With this assumption, one can neglect the $\partial_{j}\phi$ terms in the second line, but one must keep the $\partial_{j}\psi T^{j}_{k}$ term for species with isotropic pressure.

The energy-momentum tensor for a perfect fluid is given by $T^{\mu}_{\nu}=P\delta^{\mu}_{\nu}+(\rho+P)U^{\mu}U_{\nu}$. Since $U^{\mu}=dx^{\mu}/\sqrt{-ds^{2}}$ and $u^{i}=dx^{i}/d\tau$, $U^{i}=u^{i}/a$ to lowest order in the metric potentials.

Throughout this work we will assume the Pressureless Perfect Fluid approximation and set $P=0$ (i.e. our results will be applicable to any cold dark matter simulation but not to cases where e.g. neutrinos are important). The spatial part of $T^{\mu}_{\nu}$ can then be written as

[TABLE]

and from it we define the intrinsic stress tensor $\sigma^{i}_{j}$ as follows

[TABLE]

Note that this definition is different from the stress tensor $\Sigma^{i}_{j}$.

The non-linear Newtonian limit of the Euler equation can then be found by neglecting time-derivatives of the metric potentials and terms involving $P$. Using Eq. (2.5), Eq. (2.6) then reduces to

[TABLE]

From the intrinsic stress tensor we can then define divergence and curl quantities as follows

[TABLE]

where $\varepsilon^{ijk}$ is the Levi-Civita symbol. Throughout this paper we will use the terms curl and vorticity interchangeably. Applying $\partial_{i}$ to Eq. (2.9) one finds the evolution equation for the velocity divergence, namely the Euler equation

[TABLE]

We now define the vorticity of the velocity field, $w^{i}=(\nabla\times\mathbf{u})^{i}$, and construct its equation of motion by applying the curl operator, $\varepsilon^{ijk}\partial_{j}$, to Eq. (2.9)

[TABLE]

Since the curl of a gradient is zero the equation of motion is not directly sourced by the gravitational field. Vorticity can therefore only be generated by non-linear evolution through derivatives of the term $u^{j}\partial_{j}u^{i}$.

2.2 Fourier space equations

Define the forward, $\mathbf{F}$, and inverse, $\mathbf{F}^{-1}$, Fourier transforms as

[TABLE]

Rewriting a product of variables in real space in terms of their Fourier components can be achieved by using the convolution theorem

[TABLE]

with

[TABLE]

i.e., the convolution operator is commutative. Note that the sum of the arguments equals $\mathbf{k}$, i.e. when multiplying two Fourier modes the resulting wave has a wavelength equal to the sum of the arguments of the two Fourier modes.

Using a Helmholtz decomposition the velocity field, $\mathbf{u}$, can be split into a curl-free part, $\nabla\times(\nabla\Upsilon)=0$, and a divergence-free part, $\nabla\cdot(\nabla\times\mathbf{A})=0$, as follows $\mathbf{u}=-\nabla\Upsilon+\nabla\times\mathbf{A}$. Using the definitions of $\theta$ and $\mathbf{w}$ this can also be written as $\mathbf{u}=\nabla^{-2}\nabla\theta+(\nabla\times)^{-2}\nabla\times\mathbf{w}$, where $\nabla^{2}=\delta^{ij}\partial_{i}\partial_{j}$ is the comoving Laplace operator. The inverse Laplacian, $\nabla^{-2}$, applied to the exponential, $e^{i\mathbf{k}\mathbf{x}}$, has the property $\nabla^{-2}e^{i\mathbf{k}\mathbf{x}}=-k^{-2}e^{i\mathbf{k}\mathbf{x}}$, since $\nabla^{2}e^{i\mathbf{k}\mathbf{x}}=-k^{2}e^{i\mathbf{k}\mathbf{x}}$, where $k^{2}=\delta^{ij}k_{i}k_{j}$. Since $\mathbf{w}$ is divergence free, $\nabla\times\nabla\times\mathbf{w}=-\nabla^{2}\mathbf{w}$, so that $(\nabla\times)^{-2}(e^{i\mathbf{k}\mathbf{x}}\mathbf{w})=k^{-2}(e^{i\mathbf{k}\mathbf{x}}\mathbf{w})$.

The continuity equation in Fourier space

Defining $I_{\delta}=\theta\delta+\mathbf{u}\cdot\nabla\delta$, and using that the Fourier coefficient of $\mathbf{u(\mathbf{x})}$ is given by $\frac{i}{k^{2}}[-\theta(\mathbf{k})\mathbf{k}+\mathbf{k}\times\mathbf{w}(\mathbf{k})]$, we get

[TABLE]

where $\mathbf{k}^{\prime\prime}\equiv\mathbf{k-k^{\prime}}$.

The Newtonian Fourier space version of the continuity equation, Eq. (2.5), is then given by

[TABLE]

where we define $I_{\delta}(\mathbf{k})\equiv I_{\theta\delta}(\mathbf{k})+I_{w\delta}(\mathbf{k})$, with the subscripts referring to the fluid variables in the integral. The individual components are given by

[TABLE]

and

[TABLE]

Since the convolution operator is commutative we could just as well interchange $\mathbf{k}^{\prime}$ and $\mathbf{k}^{\prime\prime}$ and thereby get a different appearance for the integrand of $I_{\delta}(\mathbf{k})$. The two different ways of writing the integrand gives the same value of the integral, but specific features in the integrand will appear at different values of $\mathbf{k}^{\prime}$ in the two different representations. But, importantly, features at a specific value of $\mathbf{k}^{\prime}$ measured relative to the argument of $\delta$ is an invariant. We will use this fact below to attach physical meaning to various parts of the integrand.

The Euler equation in Fourier space

The Euler equation in Fourier space is found by Fourier transforming Eq. (2.12)

[TABLE]

with $I_{\theta}(\mathbf{k})=I_{\theta\theta}(\mathbf{k})+I_{\theta\mathbf{w}}(\mathbf{k})+I_{\mathbf{w}\mathbf{w}}(\mathbf{k})$.

$I_{\theta\theta}$ is found from the following considerations. The contribution from the real space Euler equation to $I_{\theta\theta}$ is symmetric under exchange of the two $\theta(\mathbf{x})$ variables. The convolution theorem breaks this symmetry, but it can be recovered by applying the convolution theorem twice, with the roles of the two $\theta$s exchanged, and then taking the average of the result. Finally, one arrives at a Fourier space version that is symmetric under exchange of the arguments of the two $\theta$s, namely $\mathbf{k^{\prime}}$ and $\mathbf{k^{\prime\prime}}$

[TABLE]

In the two cases mentioned above, $\frac{1}{2}k^{2}$ should be replaced with $k^{2}-\mathbf{k}\cdot\mathbf{k}^{\prime}$ and $\mathbf{k}\cdot\mathbf{k}^{\prime}$. The version with $\frac{1}{2}k^{2}$ is the one quoted in [25, 26, 27, 22, 28]. Note that the different ways of writing the integrand introduces an ambiguity in attaching physical meaning to the non-linear contribution from a given value of $\mathbf{k}^{\prime}$ or $\mathbf{k^{\prime\prime}}$, but that the version in Eq. (2.22) is the only one that preserves the symmetry of the real space Euler equation.

The mixed term, with a kernel symmetry broken by the explicit dependence on $\mathbf{k}$, is given by

[TABLE]

and finally the symmetric curl-curl term is

[TABLE]

3 Numerical implementation details

3.1 Cosmology and simulations

In all cases we perform our calculations for a flat cosmology with $\Omega_{m}=0.3$, $h=0.7$, $n_{s}=1$, $A_{s}=2.3\cdot 10^{-9}$, and $\Omega_{b}=0.05$ (our results do not depend on this specific choice and could equally well have been performed with e.g. Planck 2015 best-fit parameters [29]). The evolution of density and velocity fluctuations in linear theory is followed using CAMB [30], and the $N$-body particle initial conditions are found from a weighted sum of the CDM and baryon transfer functions. We have used the initial conditions generator first described in [31], and our pure CDM $N$-body simulations evolved with gadget-2 [32] are listed in Table 1. The initial conditions are created with either the Zel’dovich Approximation [33] or with a correction added from second order Lagrangian perturbation theory (2LPT) [34, 35, 36].

3.2 Assignment schemes

When calculating the power spectra of real space fields, it is necessary to assign $N$-body particle masses and velocities to a grid. There are various ways of performing this assignment. A deconvolved Clouds-in-cell (CIC) algorithm works well for the density field, but this method has severe drawbacks for the velocity field. Assume that there are no $N$-body particles in the 8 cubes surrounding a given grid point. For the density field, this will result in $\delta(\mathbf{x})=-1$, which can be considered a fair estimate of the “real” density. But the velocity field will be assigned the value 0, which is a very poor estimate for, e.g., the large-scale bulk flows in voids. These missing values will significantly affect the velocity power spectrum statistics, leading to less power on large scales, and more power on small scales. With the CIC scheme, this error can only be reduced by using a very coarse assignment grid.

Other assignment schemes can also be used, e.g., the Delaunay tesselation was used in [22, 37]. In this paper we use the adaptive smoothing length interpolation scheme of [38] to calculate velocity field statistics. In this method each $N$-body particle is assigned a smoothing length calculated as the distance to its 33rd nearest neighbour ($\pm 2$). This smoothing length is then used to interpolate the particle velocities to all grid points within its smoothing length. As a result, the velocity field is well sampled, also in very low density void regions.

We have investigated the effect of varying the number of neighbours used to find the smoothing length, from 33 to 22 and 66, and found no noticeable effect. Since the smoothing length interpolation kernel varies in space, the interpolation kernel is not deconvolved. The velocity field is found on a $1024^{3}$ grid, but in our subsequent mode coupling calculations we only use the largest scale modes from a $128^{3}$ sub-grid. The effect of any deconvolution should be very small on this sub-grid.

For definiteness, we use the deconvolved CIC method to calculate the density field Fourier modes, and the smoothing length method to calculate velocity divergence and curl Fourier modes. In all cases, the modes are calculated on a $1024^{3}$ regular grid. The divergence and curl quantities are calculated in $k$-space from the velocity Fourier modes.

3.3 A correlation measure

The effect of the mode coupling integrals must be measured relative to the linear sourcing terms in the fluid equations. $I_{\delta}(\mathbf{k})$ should therefore be compared to the sign of $\theta(\mathbf{k})$, and $I_{\theta}(\mathbf{k})$ should be contrasted with the sign of $-\psi(\mathbf{k})\propto\delta(\mathbf{k})$. Defining $i_{\alpha}(\mathbf{k,k^{\prime}})$ from $I_{\alpha}(\mathbf{k})=\int d^{3}k^{\prime}i_{\alpha}(\mathbf{k,k^{\prime}})$, we construct the following statistics as a measure of the amount of mode coupling between different scales

[TABLE]

where

[TABLE]

and

[TABLE]

where $\beta$ can be either of $\delta$, $\theta$ and $w$.

It is the quantity in Eq. 3.1 which is presented in the figures in this paper. The correlation measure quantifies the average correlation between two different wavenumbers in the mode coupling integrals found in the continuity and Euler equations.

3.4 Practical calculation of the mode coupling integrals

In Eq. (3.1) one must perform a sum over $\Omega^{\prime}$ and an average over $\Omega$. Since the $\delta(\mathbf{x})$ and $\theta(\mathbf{x})$ fields are real, the corresponding complex Fourier coefficients are packed in the volume with $k_{z}\geq 0$. For $k_{z}=0$ and $k_{z}=k_{N}$, where $k_{N}$ is the Nyquist frequency of the grid, redundant coefficients are stored.

Now, the deconvolution theorem requires one to perform an integration over all of $k$-space, and not just the part with $k_{z}\geq 0$. This entails that $\int d\Omega^{\prime}$ must be performed over all angles. The $\int d\Omega$ integral then estimates an average quantity for a given value of $k$. In this case it is required to perform an isotropic sampling, but this does not require that $\int d\Omega$ should be performed over all angles. Only one half of the angles must be sampled, and we choose the ones with $k_{z}\geq 0$ for simplicity. The sampling with $k_{z}<0$ gives information that is related to the information for $k_{z}\geq 0$ through complex conjugation.

In practice the integration over $\Omega$ is shifted into one over $\Omega^{\prime\prime}$ since $\mathbf{k}^{\prime\prime}$ is an argument of the perturbation variables in the mode coupling integrals in the continuity and Euler equations.

To create a binned version in $k$ and $\tilde{k}$, $\tilde{k}\in(k^{\prime},~{}k^{\prime\prime})$, performed by the operator $\mathcal{B}$, see below, we need to perform a 6-dimensional integration. We use the definition that ${}^{\prime}()^{\prime}$ around integration limits means that redundant symmetry elements must be removed when performing the integral. To avoid cumbersome notation, the integration limits are denoted with $N$, the integer corresponding to the Nyquist frequency, $k_{N}$.

The practical implementation of the integrator uses the following rewrite of the 6-dimensional integral

[TABLE]

In the above expression $\mathcal{B}$ assigns the first term to $\mathbf{k=k^{\prime\prime}+k^{\prime}}$ and the second term to $\mathbf{k=k^{\prime\prime}-k^{\prime}}$, and can bin in either $\tilde{k}=k^{\prime}$ or $\tilde{k}=k^{\prime\prime}$. For the continuity equation $\alpha\in\{\theta\delta,w\delta\}$ and $\beta=\theta$, whereas $\alpha\in\{\theta\theta,\theta w,ww\}$ and $\beta=\delta$ when calculating the amount of correlation in the Euler equation. Note however that in this paper we do not calculate the $\alpha\in\{\theta w,ww\}$ terms, since they are subdominant on the scales shown.

3.5 Convergence of results

We have compared the two simulations with $512^{3}$ and $1024^{3}$ $N$-body particles in a $256{\rm Mpc}/h$ box. The difference in $i_{\theta\delta}^{\theta}(k^{\prime\prime};k)$ is smaller than the thickness of the lines in the figures presented in the following section.

For a given value of the one dimensional grid size, $N$, the Fourier modes are sampled isotropically out to $\tilde{k}_{I}=\frac{2\pi}{R_{\rm box}}\frac{N}{2}$ for $\tilde{k}\in\{k^{\prime},k^{\prime\prime}\}$, i.e. the arguments of the Fourier variables. For $N=128$ and $R_{\rm box}=64{\rm Mpc}/h$ this gives $\tilde{k}_{I}=2\pi$. This means that, e.g., $i_{\alpha}^{\beta}(k^{\prime\prime};k)$ is only exactly sampled for $k^{\prime\prime}<2\pi-k$. This limit is very conservative and the accuracy only decreases slowly with $k^{\prime\prime}$.

To assess this issue, we have compared the results for the $64{\rm Mpc}/h$ box for grid sizes of $64^{6}$ and $128^{6}$ at $z=0$. We have found that the $64^{6}$ and $128^{6}$ cases give results that deviate by $\sim$ a line thickness in the figures in the next section.

4 Results

4.1 The physical picture and its power spectra statistics

Fig. 1 shows a slice of the $N$-body simulation volume. The figure displays $\delta$, $\mathbf{v}$, $\theta$ and $\mathbf{w}$. Several things are directly visible. First, the velocity field changes direction perpendicular to the largest filaments, which indicates that $N$-body particles are streaming into the filaments. Second, the velocity field does not display any significant small-scale structure, which mathematically can be understood from the linear theory relation $v\propto\delta/k$. Physically, this means that the gravitational field of the smaller clusters and filaments is too small to significantly alter the velocity of nearby structures. These structures are instead moved by fluctuations on even larger scales.

The velocity divergence field traces the density field exactly in linear theory. But in non-linear theory, the divergence falls behind the density field. Physically this can be understood as follows: The larger scales, $k\sim 0.1h/{\rm Mpc}$, are dominated by voids, with mass flowing out of these voids today. As the voids become more empty, less mass can flow out of them, a physical fact which is not captured in linear theory, and explains why the linear theory divergence power spectrum is higher than its non-linear counterpart. At smaller scales, $k\gtrsim 0.5h/{\rm Mpc}$, this effect happened at higher redshift. So even though these scales are dominated by clusters at low redshift, the extra suppression of the divergence power spectrum in non-linear theory is a consequence of mass flowing out of under-density regions at higher redshift.

The vorticity shows virtually no large scale perturbations, but is solely concentrated inside the largest clusters and filaments. Mathematically this can be explained by the fact that the vorticity evolution equation, Eq. (2.13), is sourced by non-linear perturbations. The vorticity changes sign (from blue to red), over a short distance scale, which shows that the clusters are rotating. Using perturbation theory [39] finds the same shape as we do for the vorticity power spectrum on larger scales, while [22] also predict the turn-over in the vorticity power spectrum.

The physical picture from Fig. 1 is captured by the power spectra statistics shown in Fig. 2. This figure shows the linear theory power spectra for $\delta$ ($\propto\theta$ in linear theory), as well as the non-linear $P_{\delta}$, $P_{\theta}$ and $P_{w}$ power spectra created by combining spectra from 4 different simulations (B, D, E and F). The divergence power spectra have been normalized so that they match the density power spectrum in each image at the largest scale simulated. The divergence and vorticity power spectra are directly comparable. $P_{\delta}$ is found using the deconvolved CIC mass assignment scheme, whereas $P_{\theta}$ and $P_{w}$ are found by using the adaptive smoothing length (SL) method of [38]. Since the shape of $P_{\delta}$ and $P_{\theta}$ are virtually identical at $z=49$, this shows that the SL method, and thereby the velocity power spectra are accurate on the scales shown.

At lower redshift, the non-linear divergence power spectra falls below the linear one. This effect begins at small scales, but as time progresses, it can be seen on much larger scales as well. At $z=0$ the linear and non-linear divergence power spectra only become similar for $k\lesssim 0.02h/{\rm Mpc}$.

The vorticity power spectrum is completely negligible at large scales and high redshift, and only becomes comparable to $P_{\theta}$ for $k\gtrsim 1h/{\rm Mpc}$ at $z\lesssim 0.5$. The vorticity power spectra peaks at $k\sim 1.5h/{\rm Mpc}$ today, at $k\sim 2.5h/{\rm Mpc}$ for $z=1$ and $k\sim 5h/{\rm Mpc}$ at $z=2.33$. It is noticeable that the divergence power spectrum has a kink in its slope at these same scales and redshifts.

Note that the curl power spectrum is non-zero at the $N$-body starting redshift, $z=49$. This initial curl field is generated by the Zel’dovich Approximation.

It is worth noting that the $\delta$ power spectrum changes slope at $k=1.5h/{\rm Mpc}$ at $z=0$. This signifies the stable clustering regime. Finally, at $z=0$ the non-linear density power spectrum lies below linear theory for $k\sim 0.02-0.1h/{\rm Mpc}$. This suppression of density power is driven by the non-linear suppression in the velocity divergence.

4.2 Initial mode coupling

Fig. 3 shows $i_{\theta\delta}^{\theta}(k^{\prime\prime};k)$, Eq. (3.4), i.e. a function of $k^{\prime\prime}$ for specific values of $k$. The amplitude in each figure is normalised in such a way that the maximum positive amplitude is unity. Results from 3 different simulations are shown. One (simulation A) initialized with ZA at $z=199$ (black and red lines) and then evolved until $z=49$ (green line), another (simulation B) initialized with ZA at $z=49$ (blue line) and finally a simulation (C) initialized with 2LPT at $z=49$ (magenta line). Focusing on the $z=49$ results, it can in all cases be seen that power is moved from larger scales to smaller scales relative to a particular mode ($k$). There is in fact a very sharp transition at $k$ (the value of $k$ is shown with the vertical dashed line). This interpretation is based on the physical understanding, that the non-linear terms in the continuity equation are driven by collapsing density perturbations, and therefore that a physical understanding of the various terms in the mode coupling integral should be based on a binning in the argument of $\delta$, namely $k^{\prime\prime}$.

In the upper panel this effect is not significantly present at $z=199$ due to noise, whereas the power transfer is present in the lower panel. But in both cases the power transfer from larger to smaller scales is a pattern that evolves during the $N$-body simulation. For this to occur, a certain fraction of modes must have an initial pattern that moves power from large to small scales, otherwise the coherent buildup of this pattern over all modes will not evolve over time.

The transfer of power from large to small scales can also be generated by the ZA alone, by applying it to a linear theory transfer function and uncorrelated random numbers. This means that the Zel’dovich Approximation by itself generates mode coupling. Adding 2LPT to ZA leads to additional mode coupling, and it can be seen that this additional mode coupling is reproduced by the simulation initialized at $z=199$ and evolved until $z=49$.

4.3 The redshift evolution of mode coupling

Figs. 4 and 5 show the redshift evolution of $i_{\theta\delta}^{\theta}(k^{\prime\prime};k)$ divided by the square of the linear growth factor.111In practice we approximate the square of the growth factor by $\sqrt{P_{\delta}(3k_{f})P_{\theta}(3k_{f})}$ for the continuity equation and $P_{\theta}(3k_{f})$ for the Euler equation, where $k_{f}$ is the fundamental frequency of the simulation volume.

Focusing on Fig. 4, it can be seen that there is a particular pattern that builds up, and that a given mode mostly couples to density perturbations which have roughly the same magnitude of the wavenumber. It can also clearly be seen that statistically, power is constantly moved from larger to smaller scales.

From Fig. 4 it can be seen that the absolute amplitude of mode coupling, divided by the linear growth factor squared, grows as a function of redshift. For $k=0.74h/{\rm Mpc}$ the mode coupling levels off to a constant shape and amplitude for $z\lesssim 1$. For this redshift range, the perturbations in this wavenumber can therefore be calculated with the linear theory growth rate multiplied with a factor that is time-$independent$.

For $k>1h/{\rm Mpc}$, see Fig. 5, the direction of power transfer becomes redshift dependent. For $k\sim 1.5h/{\rm Mpc}$, the amplitude begins to decrease for $z\lesssim 1$, while for $k$ larger than $\sim 2h/{\rm Mpc}$ the pattern has changed sign at low redshift. At $k\sim 3h/{\rm Mpc}$ the pattern first changes sign, and then decreases its amplitude once again. Notice that this sign change is local, in the sense that it does not happen for scales where $k^{\prime\prime}$ and $k$ are disparate.

The scale and redshift where this pattern shows up, matches the scale of virialisation. The change in sign can be explained by the fact that during infall smaller halos merge and transfer power from smaller to larger scales. Furthermore, fluctuations can get stretched by tidal forces as they fall into larger structures.

It is worth stressing that this pattern in Fourier space is transient and that the phenomena is followed by the stable clustering regime.

Instead of binning in $k^{\prime\prime}$ one could bin in $k^{\prime}$. But since $k^{\prime}$ is not the argument of the variable driving the physics in the continuity equation, namely $\delta$, this binning method does not give significant physical insight.

If the binning is instead performed in the angle between $\mathbf{k}$ and $\mathbf{k^{\prime\prime}}$, one finds that the largest contributions to the mode coupling integral come from small angles. From a variance point of view, this can easily be understood from the fact that, in this case, $k^{\prime}$ is reduced, thereby increasing the expected absolute value of the terms $\theta(\mathbf{k^{\prime}})$ and $\mathbf{k}\cdot\mathbf{k^{\prime}}/k^{\prime 2}$.

Figs. 6 and 7 show the effect of curl at $z=0$ and $z=1$, respectively, calculated from the simulation with a $64{\rm Mpc}/h$ box. It can be seen that the effect of curl is limited, and matters the most for scales while they are in the process of virialisation.

From the figures it can be seen that before virialisation, the effect of curl is to reduce non-linear structure formation, since the effect of curl is negative in the upper panels at both redshifts. During virialisation, middle panels, the pattern is more complex. But after virialisation, lower panels, the effect of curl is opposite the effect of the divergence, i.e., the curl works in the opposite direction of the total flow of power.

The effect of curl appears to be small since we show $i_{w\delta}^{\theta}$, i.e. the curl contribution to the mode coupling integral when correlating with the divergence. One should not expect curl to correlate significantly with the divergence, but to facilitate a comparison as a function of $k^{\prime\prime}$ we must display $i_{\theta\delta}^{\theta}$ and $i_{w\delta}^{\theta}$, i.e. both quantities must be correlated with the same fluid variable.

Fig. 8 shows the evolution of the dominant mode coupling term in the Euler equation, namely $i_{\theta\theta}^{\delta}(k^{\prime\prime};k)$. The top panel shows the time-evolution for $k=0.023h/{\rm Mpc}$, which is around the scale where the linear and non-linear divergence power spectra begin to diverge today. Despite the presence of noise, it can clearly be seen that the mode coupling leads to divergence being ’moved’ towards smaller scales, i.e. with the interpretation that the divergence decreases.

The panel with $k=0.49h/{\rm Mpc}$ shows a pattern that is reminiscent of the one found from the $i_{\theta\delta}^{\theta}$ integrand in the continuity equation, except for the fact that the negative amplitude is larger in absolute terms than the positive amplitude. Again, this accounts for the fact that the non-linear divergence is smaller than the linear one. The difference between the mode coupling integrals in the continuity and Euler equations at semi-linear scales is the presence of the extra factor of $-1$ in the kernel of $i_{\theta\theta}^{\delta}$ as compared to $i_{\theta\delta}^{\theta}$ (the $k^{\prime\prime-2}$ factor only leads to an even larger difference, compare Eqs. (2.19) and (2.22) with the substitution $\mathbf{k^{\prime\prime}=k-k^{\prime}}$). It is this extra factor that accounts for the suppression of the velocity divergence in non-linear theory as compared to linear theory.

The two lower panels in Fig. 8 show the same qualitative behaviour as was seen from the continuity equation. Interestingly, as virialisation occurs, the sum over $i_{\theta\theta}^{\delta}$ tends towards positive values, so that the non-linear divergence power spectrum grows faster than its linear counterpart. The physical driver is the $k^{2}\psi(\mathbf{k})$ term at very non-linear scales, and the effect can be seen as a change in slope for $P_{\theta}$ at low redshift and $k\gtrsim 1$ (see Fig 2).

We have not added the $\theta$-curl and curl-curl terms to the $\theta$-$\theta$ term, since the latter term is the dominant term, except on very small scales, $k\gtrsim 2-5h/{\rm Mpc}$. We have likewise not added the velocity dispersion term, $q_{\theta}$, but as indicated in [22], $q_{\theta}$ might become relevant for $k\gtrsim 1h/{\rm Mpc}$ at $z=0$. But we stress that adding these sub-dominant terms, do not affect the qualitative conclusions presented in this paper.

To sum up, at $z=0$ the scale $k\sim 1.5h/{\rm Mpc}$ (length scale $\sim 6~{}{\rm Mpc}$) stands out as the scale where the $\delta$ and $\theta$ power spectra change their slope and where the vorticity power spectrum attains its maximum. It is also the scale where the mode coupling pattern in Fig. 5, upper panel, and Fig. 8, third panel, becomes roughly featureless. It is likewise the scale where the effect of vorticity is largest, see Fig. 6, middle panel. The above considerations are also valid at $z=1$, if one translates to a length scale of $\sim 4~{}{\rm Mpc}$ corresponding to $k\sim 2.5h/{\rm Mpc}$.

5 Conclusions

We have presented a novel way of studying non-linear structure formation, namely by imaging the fully non-linear mode coupling integrals appearing in the Fourier space continuity and Euler equations using $N$-body simulations. This method allows for a direct reconstruction of the coupling between any two wavenumbers and can be easily visualised.

In accordance with expectations we find that in the linear regime there is no mode coupling. In the mildly non-linear regime the mode coupling integrals are dominated by transfer of power from large to small scales, primarily coupling modes which are close in $k$-space. This behaviour corresponds to the infall regime where structure formation is dominated by the initial collapse of structures, prior to the onset of virialisation.

On smaller scales (corresponding to $k\sim 2-3\,h$/Mpc at $z=0$) where virialisation sets in and the stable clustering regime is entered, this simple behaviour becomes much more complex. Around this scale vorticity starts to build up and structure formation slows down. In the mode coupling integrals this can be seen as a temporary reversal of the power transfer (lack of small scale resolution in our simulations prevents us from studying this effect in the extremely small scale limit). This change in behaviour sets in exactly where the curl contribution to the velocity power spectrum becomes comparable to the divergence part, and qualitatively the effect can be thought of as the build-up of pressure once virialisation sets in and the local velocity tensor is isotropised.

The results presented here are a step towards understanding the behaviour of the mode coupling integrals in the fully non-linear regime and our results can potentially be used to calibrate and improve semi-analytic models of structure formation, either by having these models calculate the statistics presented in this paper, or by extending the tools presented in this paper to calculate quantities which naturally arise in various perturbation theory models. In this work we have not investigated how the mode coupling function shown in Figs. 4 and 5 depends on the cosmological model used, but this will be the focus of a future study. If the functional form turns out to be close to universal it will be extremely useful for devising new, semi-analytic models of non-linear structure formation based on the Fourier space fluid equations. Since we have not included anisotropic stress in our calculations, calibration of semi-analytic models to our results cannot be performed reliably below the virialization scale, $k\gtrsim 1h/{\rm Mpc}$ at $z\sim 0$.

Acknowledgements

We acknowledge computing resources from the Danish Center for Scientific Computing (DCSC), and thank Y. Y. Y. Wong and T. Tram for comments on the manuscript. This work was supported by the Villum Foundation.