On the asymptotic behaviour of cosmic density-fluctuation power spectra

Sara Konrad∗, Matthias Bartelmann

Institute for Theoretical Physics, Heidelberg University, Germany

*<EMAIL_ADDRESS>
## Abstract

We study the small-scale asymptotic behaviour of the cosmic density-
fluctuation power spectrum in the Zel’dovich approximation. For doing so, we
extend Laplace’s method in arbitrary dimensions and use it to prove that this
power spectrum necessarily develops an asymptotic tail proportional to
$k^{-3}$, irrespective of the cosmological model and the power spectrum of the
initial matter distribution. The exponent $-3$ is set only by the number of
spatial dimensions. We derive the complete asymptotic series of the power
spectrum and compare the leading- and next-to-leading-order terms to derive
characteristic scales for the onset of non-linear structure formation,
independent of the cosmological model and the type of dark matter. Combined
with earlier results on the mean-field approximation for including particle
interactions, this asymptotic behaviour is likely to remain valid beyond the
Zel’dovich approximation. Due to their insensitivity to cosmological
assumptions, our results are generally applicable to particle distributions
with positions and momenta drawn from a Gaussian random field. We discuss an
analytically solvable toy model to further illustrate the formation of the
$k^{-3}$ asymptotic tail.

###### Contents

  1. 1 Introduction
  2. 2 Small-scale asymptotics of the free power spectrum
  3. 3 Implications for cosmic structure formation
  4. 4 Full asymptotic series and analytic toy model
  5. 5 Summary and discussion
  6. A Asymptotics of $d$-dimensional Laplace-Fourier integrals
  7. B Asymptotics of a specific Laplace-Fourier type integral
  8. C Asymptotic expansion of the functions $a_{1,2}(q)$ and $f_{\mu}(q)$
  9. D Analytic toy model – calculations

##  1 Introduction

Numerical simulations reveal that cosmic structures develop universal
properties, for which the radial density profiles of gravitationally-bound
objects are an important example [1, 2, 3, 4, 5, for recent work]. Where does
this universality originate from? This question can be addressed in the
framework of an analytic theory of cosmic structure formation. We use previous
results obtained with the kinetic field theory of cosmic structure formation
(KFT) [6, 7, 8] to study the asymptotic behaviour of the power spectrum of
cosmic density fluctuations in the limit of small scales, or larger wave
numbers, $k\to\infty$.

The central object of KFT is the generating functional for bundles of
classical particle trajectories on the expanding cosmological background
space-time. This generating functional incorporates a probability distribution
for the intial particle distribution in phase space, a Green’s function or
propagator $g_{qp}(t,t^{\prime})$, solving the Hamiltonian equations of
motion, and an interaction term relative to the propagator chosen [9, 10].

Initial particle positions and momenta can be drawn from a velocity potential
$\psi$ because the initial velocity field can be considered curl-free, and
density fluctuations are related to the fluctuations of this potential by the
continuity equation. Assuming that $\psi$ is a Gaussian random field, the
initial particle distribution can be statistically completely characterized by
the power spectrum $P_{\psi}^{\mathrm{(i)}}$ of the velocity potential [6].

One specific and advantageous choice for the particle propagator is given by
the Zel’dovich approximation [11]. It asserts that particles move
approximately along straight lines in a time coordinate set by the linear
growth factor $D_{+}$ of cosmic density fluctuations. Large-scale
gravitational interaction is contained in the Zel’dovich approximation because
initial particle momenta are correlated in such a way that particle streams
are convergent at matter overdensities, and divergent at underdensities. We
have shown in earlier work how particle interactions relative to the
Zel’dovich trajectories can be described in a mean-field approximation [10].

In this paper, we study the asymptotic evolution of the free KFT power
spectrum. By this, we mean the power spectrum of the cosmic density field,
decomposed into classical particles underlying Hamiltonian dynamics,
neglecting particle-particle interactions and modeling particle trajectories
with the Zel’dovich approximation.

We have shown in earlier papers [6, 7] that the free KFT power spectrum is
given by

$\mathcal{P}(k,t)=\mathrm{e}^{-Q}\int\mathrm{d}^{3}q\left(\mathrm{e}^{g_{qp}^{2}(t,0)\,\vec{k}^{\top}\hat{C}_{pp}(\vec{q}\,)\vec{k}}-1\right)\mathrm{e}^{\mathrm{i}\vec{k}\cdot\vec{q}}$
(1)

with the initial auto-correlation matrix

$\displaystyle\hat{C}_{pp}\left(\vec{q}\,\right)=\int_{k^{\prime}}\left(\vec{k}^{\prime}\otimes\vec{k}^{\prime}\,\right)P_{\psi}^{\mathrm{(i)}}(k^{\prime})\,\mathrm{e}^{\mathrm{i}\vec{k}^{\prime}\cdot\vec{q}}$
(2)

of particle momenta. The exponent $Q$ is

$Q=\frac{\sigma_{1}^{2}}{3}\,k^{2}\,g_{qp}^{2}(t,0)\;,$ (3)

where $\sigma_{1}^{2}$ is one of the moments

$\sigma_{n}^{2}=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,k^{2n+2}\,P_{\psi}^{\mathrm{(i)}}(k)$
(4)

of the power spectrum $P_{\psi}^{\mathrm{(i)}}$ of the initial velocity
potential $\psi$. It is important for the present derivation of our results
that $\sigma_{2}$ is finite. If the initial power spectrum does not fall off
steeply enough for large $k$, we introduce an exponential cut-off at a large
wave number $k_{\mathrm{s}}$ to be specified below. Via Poisson’s equation,
this spectrum is related to the initial density-fluctuation power spectrum
$P_{\delta}^{\mathrm{(i)}}$ by

$P_{\psi}^{\mathrm{(i)}}(k)=k^{-4}\,P_{\delta}^{\mathrm{(i)}}(k)\;.$ (5)

Equation (1) has been derived by assuming (i) straight particle trajectories,
(ii) a curl-free initial velocity field with velocity potential $\psi$, (iii)
an initial Gaussian random field, and (iv) taking only the initial momentum-
momentum correlations into account, i.e. neglecting density-density and
density-momentum correlations. The free power spectrum in KFT is equivalent to
the non-linear Zel’dovich power spectrum when we specify $g_{qp}$ to be the
Zel’dovich propagator
$g_{qp}^{\mathrm{(Z)}}(t,0)=D_{+}(t)-D_{+}(t^{\text{(i)}})$, where $D_{+}$ is
the linear growth factor. At the same time, it is the unperturbed zeroth-order
power spectrum in KFT and thus the starting point for calculating the non-
linear power spectrum from KFT perturbation theory or mean-field approaches
[6, 10].

The remainder of this paper has three parts. In Sect. 2, we present our result
for the leading-order, asymptotic behaviour of the free KFT power spectrum
$\mathcal{P}$. Our main result there will be that $\mathcal{P}\sim k^{-3}$ for
$k\to\infty$, where the exponent $-3$ depends neither on the cosmological
model nor on the form of the initial power spectrum. We discuss cosmological
consequences in Sect. 3, deriving a time scale for the onset of non-linear
structure growth on small scales and a maximum amplitude for the free power
spectrum which are both universal in the sense that they depend neither on the
cosmological model nor on the power spectrum of the initial state, and thus on
the dark-matter model.

In Sect. 4, we derive the entire asymptotic series of the free power spectrum.
We show there that $m$-th order terms in this series, $m\in\mathbb{N}_{0}$,
fall off asymptotically like $k^{-3-2m}$, equally independent of cosmology and
the type of dark matter. Also in Sect. 4, we introduce an analytically
solvable toy model which may be useful to clarify the origin of the asymptotic
$k^{-3}$ tail developed by the free power spectrum. In Sect. 5, we summarize
our paper and discuss our conclusions. Detailed calculations are collected in
the Appendix for the main text to be more readable.

##  2 Small-scale asymptotics of the free power spectrum

We derive here the leading-order asymptotics of the free power spectrum for
large wave numbers, $k\to\infty$. Since $k$ appears in the exponent of the
integrand in (1) with a different power than in the oscillating phase factor,
no standard asymptotic method exists that could be applied directly. We thus
develop for this case a modification of Laplace’s method for $d$-dimensional
integrals. We motivate our method in the main text and give a detailed
derivation in Appendix A, involving Morse’s lemma and a resummation of
Laplace’s method.

For evaluating the real exponent in (1), we rotate the coordinate frame such
that the $z$ axis points into $\vec{k}$ direction and write

$g_{qp}^{2}(t,0)\vec{k}^{\top}\hat{C}_{pp}(\vec{q}\,)\vec{k}=-g_{qp}^{2}(t,0)k^{2}f\left(\vec{q}\,\right)\;,$
(6)

defining the function

$f\left(\vec{q}\,\right)=-\int_{k^{\prime}}k_{z}^{\prime
2}P_{\psi}(k^{\prime})\,\mathrm{e}^{\mathrm{i}\vec{k}^{\prime}\cdot\vec{q}}\;.$
(7)

The free power spectrum then reads

$\mathcal{P}(k,t)=\mathrm{e}^{-Q}\left(\int\mathrm{d}^{3}q\,\mathrm{e}^{-g_{qp}^{2}k^{2}f(\vec{q}\,)}\,\mathrm{e}^{\mathrm{i}\vec{k}\cdot\vec{q}}-(2\pi)^{3}\delta_{\mathrm{D}}(\vec{k}\,)\right)\;.$
(8)

We neglect the delta distribution in the following because we focus on $k>0$.

The idea of our modification of Laplace’s method is now the following: For
$k\to\infty$, the integrand sharply peaks where $f(\vec{q}\,)$ has an isolated
global minimum, and it is exponentially suppressed away from it. For the
result of the integral, the most relevant part of the integrand is then the
local environment of the minimum of $f(\vec{q}\,)$. To derive the asymptotics
of $\mathcal{P}(k,t)$ for $k\to\infty$, we thus expand $f(\vec{q}\,)$ up to
quadratic order in $\vec{q}$ and then analytically integrate this new
integral.

Since the integrand in (7) except for the Fourier phase is strictly positive,
the function $f(\vec{q}\,)$ has an isolated global minimum at $\vec{q}=0$,

$\displaystyle\left|f\left(\vec{q}\,\right)\right|$
$\displaystyle\leq\int_{k^{\prime}}k_{z}^{\prime
2}P_{\psi}(k^{\prime})=\frac{1}{(2\pi)^{2}}\int_{0}^{\infty}\mathrm{d}k^{\prime}k^{\prime
4}P_{\psi}(k^{\prime})\int_{-1}^{1}\mathrm{d}\mu^{\prime}\mu^{\prime 2}$
$\displaystyle=\frac{\sigma_{1}^{2}}{3}=-f(0)\;,$ (9)

with a vanishing gradient there, $\vec{\nabla}f(\vec{q}\,)|_{\vec{q}=0}=0$.
The Hessian $A$ of $f(\vec{q}\,)$ at $\vec{q}=0$ has the components

$A_{ij}=\left.\frac{\partial^{2}f\left(\vec{q}\,\right)}{\partial
q_{i}\partial q_{j}}\right|_{\vec{q}=0}=\int_{k^{\prime}}k_{z}^{\prime
2}\,k^{\prime}_{i}k^{\prime}_{j}\,P_{\psi}(k^{\prime})=\frac{\sigma_{2}^{2}}{15}\left(\delta_{ij}+2\delta_{iz}\delta_{jz}\right)$
(10)

in terms of the moment $\sigma_{2}^{2}$. Since $A$ is diagonal and has only
positive entries, it is positive definite and has the determinant

$\det A=3\left(\frac{\sigma_{2}^{2}}{15}\right)^{3}>0\;,$ (11)

and its inverse has the components

$A^{-1}_{ij}=\frac{15}{\sigma_{2}^{2}}\left(\delta_{ij}-\frac{2}{3}\delta_{iz}\delta_{jz}\right)\;.$
(12)

These results show that all requirements listed in Appendix A.3 are met,
except for one. It is left to show that $f(\vec{q}\,)$ is quadratically
integrable in $\mathbb{R}^{3}$. This can be seen as follows: $f(\vec{q}\,)$ is
the inverse Fourier transform of the function
$F(\vec{k})=k_{z}^{2}P_{\psi}(k)$. This function is quadratically integrable
in $\mathbb{R}^{3}$ because it falls of fast enough for large $k$, while
$k^{4}P^{2}_{\psi}(k)\propto k^{2n_{\mathrm{s}}-4}$ with
$n_{\mathrm{s}}\lesssim 1$ for small $k$. Since the Fourier transform maps
quadratically integrable functions into themselves, $f(\vec{q}\,)$ is also
quadratically integrable.

We can thus apply our result (58) of Appendix A with $s=2$ and $d=3$. Using
the expression for $f(0)$ given in (9), the determinant (11) of $A$, the
$zz$-component of its inverse from (12), denoting the linearly evolved $n$-th
moment of the initial velocity potential field by

$\tau^{2}_{n}(t)\coloneqq g_{qp}^{2}(t,0)\sigma_{n}^{2}$ (13)

and introducing the abbreviation

$\Sigma(t)=\frac{5}{2\tau_{2}^{2}(t)}\;,$ (14)

we can state that the leading-order asymptotics of the free power spectrum for
$k\to\infty$ is

$\mathcal{P}(k,t)\sim\frac{3\,(4\pi)^{3/2}}{k^{3}}\Sigma^{3/2}(t)\,\mathrm{e}^{-\Sigma(t)}$
(15)

for $k\to\infty$. This result proves that the free power spectrum inevitably
falls off asymptotically like $k^{-3}$. This result is remarkably independent
of the initial power spectrum. The only condition is that the moment
$\sigma_{2}^{2}$ exists. The apparent exponential damping by $\mathrm{e}^{-Q}$
in (1) is completely canceled by the oscillating integral.

##  3 Implications for cosmic structure formation

In the preceding Section, we proved that the power spectrum for particles on
free trajectories with initially correlated momenta drawn from a Gaussian
random velocity potential always develops a $k^{-3}$ tail for large wave
numbers $k$. This universal power law occurs independently of the cosmological
model and the form of the initial power spectrum at large $k$ values. The
exponent $-3$ is solely determined by the number of spatial dimensions and has
nothing to do with the dark-matter species assumed.

###  3.1 Comparison of analytic to numerical results

To demonstrate the validity of our asymptotic expansion, we proceed to compare
our analytic with numerical results. For the initial density-fluctuation power
spectrum $P^{\mathrm{(i)}}_{\delta}$, we choose the cold dark matter power
spectrum from [12] and multiply it with an exponential regulator with
smoothing scale $k_{\mathrm{s}}$. Figure 1 shows the free power spectrum
(purple lines), together with the leading-order $k^{-3}$ term (15, green
lines) and the $k^{-5}$ term appearing at the next order in the full
asymptotic expansion derived below in (26, blue lines), for two different
redshifts ($z=2$ and $10$) and two different smoothing scales
($k_{\mathrm{s}}=10$ and $1000\,h\,\mathrm{Mpc}^{-1}$). We shall discuss free
power spectra at redshift $z=0$ later. For comparison, we also show the
linearly evolved smoothed power spectra (grey lines),

$P^{\mathrm{(lin)}}_{\delta}(k)=g_{qp}^{2}(t,0)\,P^{\text{(i)}}_{\delta}(k)\;.$
(16)

To relate the propagator to redshift, we choose the Zel’dovich propagator
$g_{qp}^{\mathrm{(Z)}}(t,0)=D_{+}(t)-D_{+}(t^{\text{(i)}})$ and set the
initial time $t^{\mathrm{(i)}}=0$. For calculating the linear growth function
$D_{+}$, we use the cosmological parameters listed in Tab. 1.

Table 1: Cosmological parameters used where needed $\Omega_{\Lambda 0}$ | $\Omega_{\mathrm{m}0}$ | $\Omega_{\mathrm{b}0}$ | $h$ | $\sigma_{8}$ | $n_{s}$
---|---|---|---|---|---
$0.7$ | $0.3$ | $0.04$ | $0.7$ | $0.8$ | $0.96$

As expected, the free power spectrum converges in all four examples to the
$k^{-3}$ asymptotics. A comparison with the linearly evolved power spectrum
cut off at different wave numbers shows that the power on small scales is
enormously increased: all formerly exponentially damped power spectra acquire
a $k^{-3}$ tail. Moreover, the asymptotic behavior is independent of the first
moment of the initial velocity potential $\sigma_{1}^{2}$, i.e. the initial
velocity dispersion.

Figure 1: At large wave numbers, the $k^{-3}$ asymptotics (15) (green lines)
is in perfect agreement with the numerical evaluation of the free power
spectrum (1) (purple lines) for early (top) and late (bottom) times as well as
for large (left) and small (right) smoothing scales. The second order of the
asymptotic expansion (26) (blue lines) predicts the scale of the onset of the
$k^{-3}$ behavior at those wave vectors where the $k^{-5}$ term falls below
the $k^{-3}$ term. The blue dashed line in the top-left panel indicates that
the $k^{-5}$ term is negative. The linearly evolved smoothed power spectra
(16, grey lines) are drawn for comparison. The cosmological parameters used
are listed in Tab. 1.

###  3.2 Time-dependence of the amplitude

The time-dependent amplitude of the leading-order asymptotics in (15),
repeated as $\mathcal{P}^{(0)}(t)$ in (26), is a function of the product
$g_{qp}^{2}(t,0)\sigma_{2}^{2}=\tau_{2}^{2}$. Thus, for straight particle
trajectories, $\tau_{2}$ acts as an effective time coordinate for structure
growth on small scales. For small $\sigma_{2}$, small-scale structures evolve
more slowly, while large values of $\sigma_{2}$ imply faster small-scale
structure growth.

The right-hand panel of Fig. 2 shows the asymptotic amplitude
$\mathcal{P}^{(0)}$ as a function of $\tau_{2}^{2}$, (grey line). Colored
crosses mark the asymptotic amplitude of the power spectra shown in the left
panel. The grey line shows that $\mathcal{P}^{(0)}$ first rises steeply during
the free infall of structures until it reaches a maximum located at

$\tau_{2,\mathrm{max}}=\sqrt{\frac{5}{3}}\;.$ (17)

This maximum indicates stream crossing at small scales because for larger
values of $\tau_{2}$, re-expansion of structures sets in, such that the
asymptotic amplitude $\mathcal{P}^{(0)}$ starts to decrease and eventually
falls off like $\tau_{2}^{-3}$ for late times. Therefore,

$g_{qp,\mathrm{max}}=\sqrt{\frac{5}{3\sigma_{2}^{2}}}\;$ (18)

implicitly sets the time of small-scale stream crossing. It is given in units
of the linear growth factor, if the Zel’dovich propagator is chosen for
$g_{qp}$. The amplitude of the leading-order asymptotic term reaches the
maximum value

$\mathcal{P}_{\mathrm{max}}^{(0)}=3\left(\frac{6\pi}{\mathrm{e}}\right)^{3/2}\approx
54.78\;,$ (19)

which is universal for all cosmological models and all initial power spectra.
The small-scale regime of the power spectrum will never grow beyond this
amplitude when straight particle trajectories and initially correlated momenta
are considered. Just the time for this universal maximum amplitude to be
reached, given in (18) in units of of the linear growth factor, depends on
$\sigma_{2}^{2}$.

Figure 2: Left: The free dimension-less power spectrum (20) for an initial CDM
power spectrum exponentially cut off at
$k_{\mathrm{s}}=1000\,h\,\mathrm{Mpc}^{-1}$ is shown for different redshifts
$z$. The amplitude rises from $z=29$ and reaches the universal maximum (19) at
$z=21$. After the onset of re-expansion, the amplitude decreases while the
$k^{-3}$ behavior persists. Right: The amplitude $\mathcal{P}^{(0)}$ of the
leading-order asymptotic term of $\mathcal{P}$ as a function of the squared
free-streaming time $\tau_{2}$ rises until it reaches the maximum value (19)
at $\tau_{2}=\sqrt{5/3}$, when stream-crossing occurs. After that, re-
expansion leads to a decay of structures, and the amplitude decreases with
$\tau_{2}^{-3}$. Colored points correspond to the asymptotic amplitudes shown
in the left panel. The cosmological parameters used are listed in Tab. 1.

In the left panel of Fig. 2, we show the free dimension-less power spectrum

$\Delta(k)=k^{3}\,\mathcal{P}(k)$ (20)

as a function of $\tau_{2}^{2}(t)$ for
$k_{\mathrm{s}}=1000\,h\,\mathrm{Mpc}^{-1}$ at five different redshifts. The
dimension-less power spectra shown become constant for large values of $k$
because $\Delta(k)\to\mathcal{P}^{(0)}$ (26). Note that the plateaus show that
the $k^{-3}$ asymptotics emerges already before stream crossing ($z=29$ and
$25$ in Fig. 2) and persists after the onset of re-expansion ($z=17$ and $13$
in Fig. 2).

###  3.3 Non-linear Zel’dovich power spectrum at redshift $z=0$

On large scales, the non-linear Zel’dovich power spectrum is known to
correspond to the linear power spectrum, which accurately describes the
structure on the largest scales even today. However, it is also well known
that the Zel’dovich approximation leads to a re-expansion of structures after
stream crossing that goes along with decreasing power at small scales, as
illustrated by Fig. 2.

In order to prevent this re-expansion, techniques like e.g. the adhesion
approximation were applied; [13, for example]. Another method to lift the tail
of the Zel’dovich power spectrum is to use a so-called truncated Zel’dovich
power spectrum, e.g. [14]. This means that the initial power spectrum is
truncated in a way similar to the exponential cut-off. Thus the initial power
is suppressed at wave numbers exceeding $k_{\mathrm{s}}$.

In Fig. 3, we show how this truncation, i.e. the suppression of small-scale
fluctuations at different scales $k_{\mathrm{s}}$ affects the tail of the
Zel’dovich power spectrum at $z=0$. In the top-left panel of Fig. 3, the free
power spectra (colored lines) for four different cut-off scales are shown,
together with the linearly evolved CDM power spectrum (black line). The tail
of all four free power spectra stays below the linearly evolved CDM spectrum.
In the bottom-left panel of Fig. 3, we show the same dimension-less power
spectra that all reach plateau values for $k\gtrsim
100\,h\,\mathrm{Mpc}^{-1}$. In the bottom-right panel of Fig. 3, we show the
amplitude $\mathcal{P}^{(0)}$ as a function of $\tau_{2}^{2}$ (grey line),
where the colored crosses indicate the values for four example power spectra
at $z=0$. The truncation at smaller wave numbers leads to smaller values of
$\sigma_{2}^{2}$ causing the amplitude at scales smaller than $2\pi
k_{\mathrm{s}}^{-1}$ to evolve more slowly. However, since the amplitude is
bound by the maximal value $\approx 54.78$, see (19), there is no truncation
scale such that the free Zel’dovich power spectrum reaches the amplitude of
the linear CDM power spectrum at small scales today.

We conclude that our results agree with the known fact that the truncated
Zel’dovich power spectrum cannot be tuned to yield the amplitude of the linear
CDM power spectrum or even the non-linear power spectrum at $z=0$ [15, for
example]. Additionally, we derived the amplitude that any truncated Zel’dovich
power spectrum can maximally reach at small scales.

Figure 3: The asymptotic amplitude evolves more slowly for smaller smoothing
wave numbers. Top: The free power spectra (1) (colored lines) for different
values of $k_{s}$ and the linearly evolved Bardeen power spectrum (black line)
are shown at today’s redshift $z=0$. The amplitudes at small scales differ by
more than three orders of magnitude. Bottom left: The dimension-less power
spectra as in the upper panel are shown. All free power spectra reach
different plateau values for $k\gtrsim 100$ $h$ Mpc-1. Bottom right: The
asymptotic amplitude $\mathcal{P}^{(0)}$ (grey line) as a function of the free
streaming time $\tau_{2}^{2}$ is shown as in the right panel of Fig. 2.
Colored crosses indicate the plateau values of the dimension-less power
spectra in the bottom left panel. For smaller $k_{s}$, the evolution of
$\mathcal{P}^{(0)}$ proceeds slower, implying a slower growth of small-scale
structures. The cosmological parameters used are listed in Tab. 1.

We finally remark that the exact cancellation of the exponential damping term
$\mathrm{e}^{-Q}$ in (1) visible in the analytic asymptotic expansion requires
particular care in numerical evaluations of the free power spectrum
$\mathcal{P}$. The small-scale behaviour of the exponent in the integrand of
(1) needs to be accurately modelled at least to order $q^{2}$ to avoid
artificial exponential growth or damping.

##  4 Full asymptotic series and analytic toy model

We now proceed to derive an expression of the entire asymptotic series for
$\mathcal{P}$ by applying a resummed version of Laplace’s method to the one-
dimensional integral of the radial coordinate $q$ in (1).

###  4.1 Full asymptotic series

As shown in[16], asymptotic series expansions for $d$-dimensional Laplace
integrals can be obtained without any smoothness conditions on the involved
functions, if only their asymptotic expansions near the critical minimum
exist. The complete expansion series is obtained by first transforming the
integrand to spherical coordinates. Afterwards, the one-dimensional Laplace
method is applied to the integral over the radial coordinate. Finally, the
integrations over the angles occuring in each term in the asymptotic expansion
has to be performed.

We adopt this procedure, but to obtain an asymptotic series for $\mathcal{P}$
as $k\to\infty$, we have to reorder the double series arising by Laplace’s
method. In our final step, we resum the series to arrive at a compact
expression for all coefficients in the asymptotic series.

We return to (1), but this time write the momentum auto-correlation matrix as

$\hat{C}_{pp}\left(\vec{q}\,\right)=-\mathds{1}_{3}a_{1}(q)-\pi_{\parallel}a_{2}(q)\;,$
(21)

where $\pi_{\parallel}$ is the projector

$\pi_{\parallel}=\frac{\vec{q}\otimes\vec{q}}{q^{2}}$ (22)

parallel to $\vec{q}$ and the momentum auto-correlation functions

$a_{1}(q)=\frac{\xi_{\psi}^{\prime}(q)}{q}\;,\quad
a_{2}(q)=\xi_{\psi}^{\prime\prime}(q)-a_{1}(q)$ (23)

are determined by derivatives of the auto-correlation function

$\xi_{\psi}(q)=\int_{k}\,P_{\psi}^{\mathrm{(i)}}(k)\,\mathrm{e}^{\mathrm{i}\vec{k}\cdot\vec{q}}$
(24)

of the initial velocity potential $\psi$.

We align the wave vector $\vec{k}\neq 0$ in (1) with the $z$-axis and
transform to spherical coordinates $(q,\varphi,\mu)$, where $\mu=q_{z}/q$ is
the cosine of the angle enclosed by $\vec{q}$ and $\vec{k}$. Since the
integration kernel is independent of the azimuthal angle $\varphi$, we can
write

$\mathcal{P}(k,t)\sim
2\pi\mathrm{e}^{-Q}\int_{-1}^{1}\mathrm{d}\mu\int_{0}^{q_{\mathrm{max}}}\mathrm{d}q\,q^{2}\mathrm{e}^{-g_{qp}^{2}k^{2}\left[a_{1}(q)+\mu^{2}a_{2}(q)\right]}\mathrm{e}^{\mathrm{i}k\mu
q}\quad\mbox{as}\quad k\to\infty\;.$ (25)

The $q$ integration can stop at the arbitrary upper bound $q_{\mathrm{max}}>0$
because the asymptotic behaviour of the integrand is determined by the
exponent in the immediate neighbourhood of its minimum. As explained above, we
can neglect the delta distribution resulting from integrating the Fourier
phase over the unity subtracted from the integrand in (1).

The resulting asymptotic series is derived in Appendix B. The result is

$\mathcal{P}(k,t)\sim\sum_{m=0}^{\infty}\frac{\mathcal{P}^{(m)}(t)}{k^{3+2m}}\quad\mbox{as}\quad
k\to\infty\;,$ (26)

with the time-dependent coefficients

$\mathcal{P}^{(m)}(t)\coloneqq\frac{\pi^{3/2}}{(2m)!}\lim_{x\to
0}\,\frac{\mathrm{d}^{2m}}{\mathrm{d}x^{2m}}\left\\{\left(\frac{x^{2}/g_{qp}^{2}(t,0)}{a_{1}(x)+\frac{\sigma_{1}^{2}}{3}}\right)^{m+1}\left.\frac{\mathrm{d}^{m}}{\mathrm{d}y^{m}}\frac{y^{m}\sqrt{R(x,y)}}{\mathrm{e}^{R(x,y)}}\right|_{y=y_{0}(x)}\right\\}\;.$
(27)

containing the definitions

$R(x,y)=\frac{x^{2}/g_{qp}^{2}(t,0)}{a_{2}(x)}\frac{y}{1+y}\quad\mbox{and}\quad
y_{0}(x)=\frac{a_{2}(x)}{a_{1}(x)+\sigma_{1}^{2}/3}\;.$ (28)

The coefficients $\mathcal{P}^{(m)}(t)$ are therefore given by compact
expressions. Using the asymptotic expansions of the functions $a_{1}(q)$ and
$a_{2}(q)$ for $q\to 0$ derived in Appendix C, the first two coefficients turn
out to be

$\displaystyle\mathcal{P}^{(0)}(t)$
$\displaystyle=3\,(4\pi)^{3/2}\,\Sigma^{3/2}(t)\,\mathrm{e}^{-\Sigma(t)}\;,$
$\displaystyle\mathcal{P}^{(1)}(t)$
$\displaystyle=\frac{(4\pi)^{3/2}}{28}\,\frac{\sigma_{3}^{2}}{\sigma_{2}^{2}}\,\Sigma^{5/2}(t)\,\mathrm{e}^{-\Sigma(t)}\left[123-132\Sigma(t)+20\Sigma^{2}(t)\right]\;.$
(29)

In general, since the asymptotic series of $a_{1}(q)$ and $a_{2}(q)$ in the
limit $q\to 0$ are even in $q$ with the $n$-th coefficients being proportional
to $\sigma_{n+1}^{2}$, (29) implies that $\mathcal{P}^{(m)}(t)$ depends only
on moments of the initial velocity potential that are of order $m+2$ and
lower. Furthermore,

$\mathcal{P}^{(m+1)}(t)\propto\sigma_{m+3}^{2}\quad\mbox{for all}\quad
m=0,1,2,\ldots\;,$ (30)

showing that higher orders of the asymptotics series for the free power
spectrum $\mathcal{P}$ carry higher orders of the initial velocity potential.

###  4.2 Analytic toy model

The free power spectrum can be calculated analytically for the initial density
power spectrum

$P_{\delta}^{\mathrm{(i)}}(k)=Akj_{3}(k)\;,$ (31)

with amplitude $A$ and the dimension-less wave number $k$ referred to some
arbitrary scale $k_{0}$ set to unity. This spectrum is shown in the left panel
of Fig. 4. Clearly, this example is unphysical for cosmic structure formation
because real fields, like the cosmic density fluctuations, cannot have
negative-valued power spectra. Nevertheless, we find this example very
instructive from a mathematical point of view, and worth discussing.

Figure 4: Left: The toy-model power spectrum (31) increases $\propto k^{4}$
for small $k\ll k_{0}$. For large $k\gg k_{0}$, it oscillates like a sine
wave, phase shifted by $-3\pi/2$. Dashed lines indicate negative values.
Right: On large scales (small wave numbers), the three terms $T_{1,2,3}$
cancel to several orders. At intermediate scales, $\mathcal{P}$ approximately
equals $T_{3}$, while $T_{1}$ dominates on small scales.

We show in App. D that the spectrum (31) has the moment

$\frac{\sigma_{1}^{2}}{3}=\frac{A}{8\pi}\;,$ (32)

and the correlation functions $a_{1}(q)$ and $a_{2}(q)$ are

$a_{1}(q)=\frac{\sigma_{1}^{2}}{3}\left(q^{2}-1\right)\Theta\left(1-q\right)\;,\quad
a_{2}(q)=\frac{2\sigma_{1}^{2}}{3}\,q^{2}\,\Theta\left(1-q\right)$ (33)

with the Heaviside step function $\Theta$. The radial coordinate $q$ is
dimension-less as well, expressed in units of $k_{0}^{-1}$. With these
expressions, aligning the wave vector once more with the $z$ axis, we compute
the free power spectrum

$\displaystyle\mathcal{P}(k,t)$
$\displaystyle=2\pi\int_{-1}^{1}\mathrm{d}\mu\int_{0}^{k_{0}^{-1}}\mathrm{d}q\,q^{2}\left[\mathrm{e}^{-Q(1+2\mu^{2})q^{2}}-\mathrm{e}^{-Q}\right]\mathrm{e}^{\mathrm{i}k\mu
q}$ $\displaystyle=T_{1}(k)+T_{2}(k)+T_{3}(k)$ (34)

with $Q$ from (3). With the definition

$\Sigma=\left(4\tau_{1}^{2}\right)^{-1}\;,$ (35)

the three terms $T_{1,2,3}(k)$ resulting from the integration in (51) are

$\displaystyle T_{1}(k)$
$\displaystyle=\frac{3(4\pi\Sigma)^{3/2}}{k^{3}}\mathrm{e}^{-\Sigma}\operatorname{Re}\mathop{\mathrm{erf}}\left(\tau_{1}k+\frac{\mathrm{i}}{2\tau_{1}}\right)\;,$
$\displaystyle T_{2}(k)$
$\displaystyle=-\frac{2(6\pi\Sigma)^{3/2}}{k^{3}}\mathrm{e}^{-3\Sigma/2-Q}\operatorname{Re}\mathop{\mathrm{erf}}\left[\frac{1}{\sqrt{6}}\left(2\tau_{1}k+\frac{3\mathrm{i}}{2\tau_{1}}\right)\right]\;,\quad\mbox{and}$
$\displaystyle T_{3}(k)$
$\displaystyle=-\frac{4\pi}{k}\,\mathrm{e}^{-Q}\,j_{1}(k)\;.$ (36)

The term $T_{3}$ falls off exponentially for $k\to\infty$. For $T_{1}$ and
$T_{2}$, we recall $\mathop{\mathrm{erf}}z^{*}=\mathop{\mathrm{erf}}^{*}z$ for
$z\in\mathbb{C}$ to write
$2\operatorname{Re}\mathop{\mathrm{erf}}z=\mathop{\mathrm{erf}}z+\mathop{\mathrm{erf}}z^{*}$
and use the asymptotic expansion

$\mathop{\mathrm{erf}}(ak+\mathrm{i}b)+\mathop{\mathrm{erf}}(ak-\mathrm{i}b)=2+\mathcal{O}\left(\frac{\mathrm{e}^{-a^{2}k^{2}}}{k}\cos(2abk)\right)\;,$
(37)

from which we can conclude that $T_{2}$ is also exponentially suppressed. For
$T_{1}$, we get

$T_{1}(k)\sim\frac{3(4\pi\Sigma)^{3/2}}{k^{3}}\,\mathrm{e}^{-\Sigma}\quad\mbox{as}\quad
k\to\infty\;,$ (38)

which is in perfect agreement with the asymptotics that we derived above.

Figure 4 illustrates what happens to this power spectrum: Free streaming first
suppresses power on small scales and leads to an exponential decrease of the
power spectrum. The collective motion of the particles with correlated momenta
wipes out the oscillatory pattern in the initial power spectrum and builds up
the power on small scales.

##  5 Summary and discussion

In Sect. 2 of this paper, we derived the leading-order asymptotic behaviour on
small scales, $k\to\infty$, of the free density-fluctuation power spectrum
$\mathcal{P}$ in KFT, given by (1). By applying rigorous results for
$d$-dimensional integrals of the Fourier-Laplace type, derived in Appendix A,
we proved that the free power spectrum always develops a $k^{-3}$ tail,
irrespective of the cosmological model and the shape of the power spectrum
characterizing the initial particle distribution. This small-scale behavior
has been suggested before by [15] for the non-linear Zel’dovich power
spectrum. Independently of our work, the $k^{-3}$ asymptotics has recently
been derived in the framework of Lagrangian perturbation theory by [17],
finding the same amplitude that we derive here.

After discussing cosmological consequences, we then derived the complete
asymptotic series of $\mathcal{P}$. We proved that terms of order $m$ fall off
asymptotically proprotional to $k^{-3-2m}$, derived a closed equation for the
coefficients $\mathcal{P}^{(m)}$ of this series, and gave specific expressions
for these coefficients in leading and next-to-leading order.

An important next step will be to study whether the $k^{-3}$ behavior at
leading order, which we proved here for free trajectories, also persists when
more gravitational interactions are included, for example via perturbation or
mean-field theory. Additionally, we only took initial momentum-momentum
correlations into account, neglecting initial density-density and initial
density-momentum correlations. At large scales, it has been shown that this is
a valid approach. At small scales, however, this assumption also has to be
studied more thoroughly.

It is an important implication of our work that the exponential damping factor
in front of the integral (1) is exactly cancelled by the asymptotic expansion.
Freely streaming, classical particles initially drawn from a Gaussian random
density and velocity field will thus not lead to an exponential damping of
small-scale structures, but to structures with a dimension-less power spectrum
$k^{3}\mathcal{P}$ developing a constant level for $k\to\infty$. Apart from
its physical significance, this cancellation of the exponential damping term
also has implications for the implementation of codes for evaluating the free
power spectrum. It is crucial for the evaluation at large wave numbers that
the small-scale behavior of the initial momentum-correlation function is
implemented asymptotically correctly up to order $q^{2}$. The results will
otherwise diverge artificially and exponentially.

We further showed that the time evolution of the amplitude of the leading-
order, $k^{-3}$ asymptotic tail depends only on the product
$g_{qp}^{2}(t,0)\sigma_{2}^{2}$, showing that the evolution of the small-scale
structure proceeds more slowly when the initial power spectrum is normalized
to have a smaller moment $\sigma_{2}^{2}$. At early times, the amplitude rises
steeply, indicating the convergence of particle streams due to their initially
correlated momenta. The leading-order asymptotic term then reaches a universal
maximum value $\approx 54.78$ when $g_{qp}^{2}(t,0)\sigma_{2}^{2}=5/3$. This
maximum amplitude marks small-scale stream crossing, where particle
trajectories cross and the particles subsequently move apart. This leads to
the re-expansion of the small-scale structures that is known to occur in the
Zel’dovich approximation and that goes along with the decreasing amplitude at
small scales.

With this, we lay the foundation for future work to explore the possible
imprint of these scales in observables and to derive the precise relations for
specific dark matter candidates.

## Acknowledgements

We thank Manfred Salmhofer for numerous helpful discussions and for reviewing
a draft of this paper. Constructive comments by Robert Lilow, Ricardo Waibel,
Leif Seute and Hannes Riechert contributed to clarifying the line of reasoning
and improving numerical algorithms. This work is funded by the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy EXC 2181/1 - 390900948 (the Heidelberg STRUCTURES
Excellence Cluster).

##  Appendix A Asymptotics of $d$-dimensional Laplace-Fourier integrals

We consider $d$-dimensional Laplace-Fourier integrals

$P(k)=\int_{\Omega}\mathrm{e}^{-|k|^{s}f(x)}\mathrm{e}^{\mathrm{i}k\cdot
x}\mathrm{d}^{d}x$ (39)

where the vector $k\in\mathbb{R}^{d}$ appears with different powers in the
real and the imaginary exponent of the integrand. For this type of integrals,
we wish to study the asymptotics in the limit $|k|\to\infty$. To this end, we
need to prepare some concepts and notation. Here and below, $D$ denotes the
derivative of a function with respect to its arguments.

###  A.1 Preliminary remarks

Let $\Omega\subset\mathbb{R}^{d}$ be an open subset of $\mathbb{R}^{d}$ and
$x_{0}\in\Omega$ a non-degenerate critical point of the smooth function
$f:\Omega\to\mathbb{R}$ on $\Omega$ with $f\in C^{\infty}$. Then, by Morse’s
lemma, neighbourhoods $U$, $V$ of the points $y=0$ and $x=x_{0}$ and a
diffeomorphism $h:U\to V$ exist such that

$\left(f\circ h\right)(y)=f(x_{0})+\frac{1}{2}y^{\top}Qy\;,\quad
Q=\mathrm{diag}\left(\mu_{1},\ldots,\mu_{d}\right)\;,$ (40)

where the Jacobian $H:=Dh$ of $h$ has unit determinant at $y=0$, $\det
H_{0}=1$ with $H_{0}:=H(0)$. Then, by the chain rule and $(Df)(x_{0})=0$,
$H_{0}$ diagonalises the Hessian matrix

$A:=\left(D^{2}f\right)\Bigr{|}_{x=x_{0}}$ (41)

of $f$ in $x_{0}$ such that

$Q=H^{\top}_{0}AH_{0}\;.$ (42)

With the same chart $h:U\to V$, the Jacobian $H=Dh$ and a smooth function
$g:\Omega\to\mathbb{R}$ with $g\in C^{\infty}$, we define the function
$G:U\to\mathbb{R}$ by $G(y)=(g\circ h)(y)\det H(y)$. We further introduce
multi-indices $\alpha=(\alpha_{1},\alpha_{2},\ldots,\alpha_{d})$ and agree on
the notation

$\displaystyle|\alpha|$
$\displaystyle:=\alpha_{1}+\alpha_{2}+\ldots+\alpha_{d}\;,$
$\displaystyle\alpha!$
$\displaystyle:=\alpha_{1}!\alpha_{2}!\cdots\alpha_{d}!\;,$
$\displaystyle\Gamma(\alpha)$
$\displaystyle:=\Gamma(\alpha_{1})\Gamma(\alpha_{2})\cdots\Gamma(\alpha_{d})\;,$
$\displaystyle\mu^{\alpha}$
$\displaystyle:=\mu_{1}^{\alpha_{1}}\mu_{2}^{\alpha_{2}}\cdots\mu_{d}^{\alpha_{d}}\;,$
(43)

where $\mu=(\mu_{1},\mu_{2},\ldots,\mu_{d})$ is a $d$ dimensional real-valued
vector. We finally define the symbol

$\delta(\alpha)=\begin{cases}1&\mbox{all $\alpha_{j}$ even}\\\
0&\mbox{else}\end{cases}$ (44)

and the derivative operator

$D^{\alpha}G(0)=\frac{\partial^{|\alpha|}}{\partial^{\alpha_{1}}y_{1}\ldots\partial^{\alpha_{d}}y_{d}}G(y)\Biggr{|}_{y=0}\;.$
(45)

With these preparations, we can now continue with the following theorem for
the asymptotic behaviour of a class of integrals on $\Omega$ [18].

###### Theorem 1

Let $J(\lambda)$ be the integral

$J(\lambda)=\int_{\Omega}\mathrm{e}^{-\lambda f(x)}g(x)\,\mathrm{d}x$ (46)

with $x\in\mathbb{R}^{d}$ and $\Omega\subset\mathbb{R}^{d}$, whose asymptotic
behaviour we wish to determine for $\lambda\to\infty$. We assume that

  1. 1.

$f,g\in C^{\infty}$ on $\Omega$;

  2. 2.

$J(\lambda)$ converges absolutely for sufficiently large $\lambda$;

  3. 3.

$f$ has a minimum at, and only at, $x_{0}\in\Omega$ such that

$\rho(\varepsilon):=\inf_{\Omega\setminus
B_{\varepsilon}(x_{0})}f(x)-f(x_{0})>0$ (47)

for all $\varepsilon>0$, where $B_{\varepsilon}(x_{0})$ is the open ball with
radius $\varepsilon$ around $x_{0}$; and that

  4. 4.

the Hessian matrix $A$ of $f$ in $x_{0}$ is positive definite.

Then, the integral $J(\lambda)$ has the asymptotic expansion

$J(\lambda)\sim\mathrm{e}^{-\lambda
f(x_{0})}\sum_{n=0}^{\infty}\frac{c_{n}}{\lambda^{d/2+n}}$ (48)

for $\lambda\to\infty$, with

$c_{n}=\sum_{|\alpha|=2n}\delta(\alpha)\left(\frac{2}{\mu}\right)^{(\alpha+1)/2}\Gamma\left(\frac{\alpha+1}{2}\right)\frac{D^{\alpha}G(0)}{\alpha!}\;.$
(49)

The $\mu_{i}$, $1\leq i\leq d$, are the (real and positive) eigenvalues of
$A$.

###  A.2 Laplace-Fourier integrals in $d$ dimensions

We now apply the result (48) with the coefficients (49) to Laplace-Fourier-
type integrals of the form

$J(\lambda,k)=\int_{\Omega}\mathrm{e}^{-\lambda
f(x)}\mathrm{e}^{\mathrm{i}kx}\mathrm{d}^{d}x\;,$ (50)

whose asymptotic behaviour we wish to determine for $\lambda\to\infty$, where
$|k|$ is assumed to be a large parameter at the same time. Without loss of
generality, we set $x_{0}=0$ here. For integrals of this type, the function
$G$ defined above is

$G(y)=\mathrm{e}^{\mathrm{i}kh(y)}\det H(y)\;.$ (51)

Since only the derivatives of $G$ acting on the exponential factor increase
the powers of $k$, and since additionally $G(0)=1$, we can state

$D^{\alpha}G(0)\sim\prod_{j=1}^{d}\left[\mathrm{i}\left(H^{\top}_{0}k\right)_{j}\right]^{\alpha_{j}}=\mathrm{i}^{|\alpha|}\prod_{j=1}^{d}\left(H^{\top}_{0}k\right)_{j}^{\alpha_{j}}\quad\mbox{for}\quad|k|\to\infty$
(52)

and move on to write

$c_{n}(k)\sim\sum_{|\alpha|=2n}\delta\left(\alpha\right)\left(-1\right)^{n}\prod_{j=1}^{d}\left(\frac{2}{\mu_{j}}\right)^{(\alpha_{j}+1)/2}\Gamma\left(\frac{\alpha_{j}+1}{2}\right)\frac{\left(H^{\top}_{0}k\right)_{j}^{\alpha_{j}}}{\alpha_{j}!}\quad\mbox{for}\quad|k|\to\infty\;.$
(53)

Since the factor $\delta(\alpha)$ ensures that only even $\alpha_{j}$ can
occur in the sum, it is convenient to define a multi-index $n$ with elements
$n_{j}=\alpha_{j}/2$. Evaluating the Gamma functions for odd arguments in (53)
and noting that $\prod\mu_{j}=\det A$, we arrive at

$c_{n}(k)\sim\sqrt{\frac{(2\pi)^{d}}{\det
A}}\left(-\frac{1}{2}\right)^{n}\sum_{|n|=n}\prod_{j=1}^{d}\frac{1}{n_{j}!}\frac{\left(H^{\top}_{0}k\right)_{j}^{2j}}{\mu_{j}^{n_{j}}}\quad\mbox{for}\quad|k|\to\infty\;.$
(54)

Since (42) implies that $H_{0}Q^{-1}H_{0}^{\top}=A^{-1}$, the remaining sum
equals

$\sum_{|n|=n}\prod_{j=1}^{d}\frac{1}{n_{j}!}\frac{\left(H^{\top}_{0}k\right)_{j}^{2j}}{\mu_{j}^{n_{j}}}=\frac{1}{n!}\left[\left(H^{\top}_{0}k\right)^{\top}Q^{-1}\left(H^{\top}_{0}k\right)\right]^{n}=\frac{\left(k^{\top}A^{-1}k\right)^{n}}{n!}$
(55)

allowing us to write (48) in the form

$J(\lambda,k)\sim\mathrm{e}^{-\lambda
f(0)}\sqrt{\frac{(2\pi)^{d}}{\lambda^{d}\det
A}}\exp\left(-\frac{k^{\top}A^{-1}k}{2\lambda}\right)\quad\mbox{for}\quad\lambda\to\infty\;.$
(56)

For integrals of the type

$P(k)=\int_{\Omega}\mathrm{e}^{-|k|^{s}f(x)}\mathrm{e}^{\mathrm{i}kx}\mathrm{d}^{d}x$
(57)

with $s\geq 2$, we substitute again $\lambda\to|k|^{s}$ in (51). We then
obtain immediately from (56)

$P(k)\sim\mathrm{e}^{-|k|^{s}f(0)}\sqrt{\frac{(2\pi)^{d}}{|k|^{sd}\det
A}}\exp\left(-\frac{k^{\top}A^{-1}k}{2|k|^{s}}\right)\quad\mbox{for}\quad|k|\to\infty\;,$
(58)

which proves the result (3).

###  A.3 Extension to $\mathbb{R}^{d}$

We finally wish to extend the integration domain $\Omega$ to $\mathbb{R}^{d}$.
We begin with the conditions imposed in Theorem 1, but relax conditions (1)
and (2) to

  1. $1^{\prime}$.

the Hessian $A=D^{2}f$ exists in $x_{0}=0$ and

  2. $2^{\prime}$.

$f\in L_{2}(\mathbb{R}^{d})$, i.e. $f$ is quadratically integrable in
$\mathbb{R}^{d}$

and add as a further condition:

  1. 3b.

The infimum

$\sigma(\varepsilon):=\inf_{\mathbb{R}^{d}\setminus
B_{\varepsilon}(0)}|f(0)|-|f(x)|>0$ (59)

for all $\varepsilon>0$.

Note that conditions (3) and (3b) together imply that $f(0)<0$. Then, the
integral

$P(k)=\int_{\mathbb{R}^{d}}\mathrm{e}^{-|k|^{s}f(x)}\mathrm{e}^{\mathrm{i}kx}\mathrm{d}^{d}x$
(60)

has the asymptotic expansion (58).

To prove this statement, we consider the integral

$\displaystyle\bar{P}(k):$
$\displaystyle=\int_{\mathbb{R}^{d}}\left[\mathrm{e}^{-|k|^{s}f(x)}-1+|k|^{s}f(x)\right]\mathrm{e}^{\mathrm{i}kx}\mathrm{d}^{d}x$
$\displaystyle=P(k)-(2\pi)^{N}\delta_{\mathrm{D}}(k)+|k|^{s}\tilde{f}(k)$
$\displaystyle=P(k)+|k|^{s}\tilde{f}(k)$ (61)

for $k>0$, where $\tilde{f}$ denotes the Fourier transform of $f$. We further
split the integration domain as

$\bar{P}(k)=\left(\int_{\Omega}+\int_{\mathbb{R}^{d}\setminus\Omega}\right)\left[\mathrm{e}^{-|k|^{s}f(x)}-1+|k|^{s}f(x)\right]\mathrm{e}^{\mathrm{i}kx}\mathrm{d}^{d}x=:\bar{P}_{\Omega}(k)+\bar{P}_{\bar{\Omega}}(k)\;.$
(62)

For the asymptotic expansion of $\bar{P}_{\Omega}(k)$, we estimate

$\bar{P}_{\Omega}(k)\leq\int_{\Omega}\mathrm{e}^{-|k|^{s}f(x)}\mathrm{e}^{\mathrm{i}kx}\mathrm{d}^{d}x+\left[1+|k|^{s}|f(0)|\right]V_{\Omega}\;,$
(63)

where $V_{\Omega}$ is the volume of the domain $\Omega$. Since the second term
on the right-hand side of (63) is exponentially suppressed compared to the
asymptotic expansion (58) of the first term (recall that $f(0)<0$), we can
conclude that

$\bar{P}_{\Omega}(k)\sim\mathrm{e}^{-|k|^{s}f(0)}\sqrt{\frac{(2\pi)^{d}}{|k|^{sd}\det
A}}\exp\left(-\frac{k^{\top}A^{-1}k}{2|k|^{s}}\right)\quad\mbox{for}\quad|k|\to\infty\;.$
(64)

For $\bar{P}_{\bar{\Omega}}(k)$, we estimate

$\displaystyle\bar{P}_{\bar{\Omega}}(k)$
$\displaystyle\leq\sum_{n=2}^{\infty}\frac{|k|^{ns}}{n!}\int_{\mathbb{R}^{d}\setminus\Omega}\left|f(x)\right|^{n}\mathrm{d}^{d}x$
$\displaystyle\leq\sum_{n=2}^{\infty}\frac{|k|^{ns}}{n!}\left[|f(0)|-\sigma(\varepsilon)\right]^{n-2}\underbrace{\int_{\mathbb{R}^{d}}|f(x)|^{2}\mathrm{d}^{d}x}_{=:\,C\,<\,\infty}$
$\displaystyle=\frac{C}{\left[|f(0)|-\sigma(\varepsilon)\right]^{2}}\left\\{\mathrm{e}^{-|k|^{s}\left[f(0)+\sigma(\varepsilon)\right]}-1+|k|^{s}\left[f(0)+\sigma(\varepsilon)\right]\right\\}$
$\displaystyle\sim\frac{C}{\left[|f(0)|-\sigma(\varepsilon)\right]^{2}}\mathrm{e}^{-|k|^{s}\left[f(0)+\sigma(\varepsilon)\right]}\quad\mbox{for}\quad|k|\to\infty\;.$
(65)

Again, this result is exponentially suppressed by the factor
$\exp(-|k|^{s}\sigma(\varepsilon))$ compared to the asymptotic expansion (58).
Since $\tilde{f}\in L_{2}(\mathbb{R}^{d})$, we can conclude

$\displaystyle P(k)$
$\displaystyle=\bar{P}(k)-|k|^{s}\tilde{f}(k)\sim\bar{P}(k)\sim\bar{P}_{\Omega}(k)$
$\displaystyle\sim\mathrm{e}^{-|k|^{s}f(0)}\sqrt{\frac{(2\pi)^{d}}{|k|^{sd}\det
A}}\exp\left(-\frac{k^{\top}A^{-1}k}{2|k|^{s}}\right)\quad\mbox{for}\quad|k|\to\infty\;,$
(66)

as claimed.

##  Appendix B Asymptotics of a specific Laplace-Fourier type integral

We now wish to study the asymptotic behavior of the specific Laplace-Fourier
type integral

$P(k)=\int_{-1}^{1}\mathrm{d}\mu\int_{0}^{\infty}\mathrm{d}q\;q^{2}\mathrm{e}^{-k^{2}g_{qp}^{2}\left[a_{1}(q)+\mu^{2}a_{2}(q)\right]}\mathrm{e}^{\mathrm{i}k\mu
q}$ (67)

for $k\to\infty$. To this end, we begin with an asymptotic expansion of the
$q$ integral, for which we modify Erdélyi’s theorem [19, 20, 21] for Laplace
integrals, which states:

######  Theorem 2 (Erdélyi)

Let $I(\lambda)$ be an integral of the form

$I(\lambda)=\int_{a}^{b}\mathrm{e}^{-\lambda f(x)}g(x)\,\mathrm{d}x\;,$ (68)

where $f(x)$ is a real function of the real variable $x$, while $g(x)$ may be
real or complex. Then, if

  1. 1.

$f(x)>f(a)$ for $x\in(a,b)$ and

$\inf_{[a+\delta,b)}f(x)-f(a)>0$ (69)

for $\delta>0$;

  2. 2.

$f^{\prime}(x)$ and $g(x)$ are continuous in a neighbourhood of $a$, except
possibly at $a$;

  3. 3.

$f$ and $g$ admit asymptotic expansions

$\displaystyle f(x)$ $\displaystyle\sim
f(a)+\sum_{k=0}^{\infty}a_{k}(x-a)^{k+\alpha}\;,$ $\displaystyle g(x)$
$\displaystyle\sim\sum_{k=0}^{\infty}b_{k}(x-a)^{k+\beta-1}\;,$ (70)

$f$ can be term-wise differentiated,

$f^{\prime}(x)\sim\sum_{k=0}^{\infty}a_{k}(k+\alpha)(x-a)^{k+\alpha-1}$ (71)

for $x\to a^{+}$, where $\alpha>0$ and $\mathrm{Re}\,\beta>0$; and

  4. 4.

$I(\lambda)$ converges absolutely for sufficiently large $\lambda$; then

the integral $I(\lambda)$ has the asymptotic expansion

$I(\lambda)\sim\mathrm{e}^{-\lambda
f(a)}\sum_{n=0}^{\infty}\frac{\Gamma(\nu)\,c_{n}}{\lambda^{\nu}}$ (72)

for $\lambda\to\infty$, where $\nu:=(n+\beta)/\alpha$. The coefficients
$c_{n}$ can be expressed by $a_{n}$ and $b_{n}$ as

$c_{n}=\frac{1}{\alpha
a_{0}^{\nu}}\sum_{m=0}^{n}\frac{b_{n-m}}{m!}\,d_{m,n}\quad\mbox{with}\quad
d_{m,n}=\lim_{x\to
0}\frac{\mathrm{d}^{m}}{\mathrm{d}x^{m}}\left(1+\sum_{j=1}^{\infty}\frac{a_{j}}{a_{0}}\,x^{j}\right)^{-\nu}\;.$
(73)

For applying Erdelyi’s theorem, we set $\lambda:=k^{2}g_{qp}^{2}$ and define

$f_{\mu}(q):=a_{1}(q)+\mu^{2}a_{2}(q)\quad\text{and}\quad
g_{\mu,k}(q):=q^{2}\mathrm{e}^{\mathrm{i}k\mu q}\;.$ (74)

As we shall see in the next section, this function has the asymptotic
expansion

$f_{\mu}(q)\sim-\frac{\sigma_{1}^{2}}{3}+\sum_{m=0}^{\infty}a_{2m}(\mu)\,q^{2m+2}$
(75)

with

$a_{2m}(\mu)=\frac{(-1)^{m+2}\sigma_{m+2}^{2}}{(5+2m)(3+2m)!}\left(1+2(m+1)\mu^{2}\right)\;,$
(76)

and the asymptotic expansion of the function $g_{\mu,k}(q)$ is

$g_{\mu,k}(q)\sim\sum_{m=0}^{\infty}b_{m}(\mu,k)\,q^{m+3-1}\quad\mbox{with}\quad
b_{m}(\mu,k)=\frac{(\mathrm{i}k\mu)^{m}}{m!}\;.$ (77)

We identify $\alpha=2$ and $\beta=3$ and find from (68) and (72)

$\displaystyle I_{\mu,k}(\lambda)$
$\displaystyle:=\int_{0}^{\infty}\mathrm{d}q\,q^{2}\mathrm{e}^{-\lambda
f_{\mu}(q)}\mathrm{e}^{\mathrm{i}k\mu q}$
$\displaystyle\sim\int_{0}^{b}\mathrm{d}q\;\mathrm{e}^{-\lambda
f_{\mu}(q)}g_{\mu,k}(q)\sim\mathrm{e}^{-\lambda
f_{\mu}(0)}\sum_{n=0}^{\infty}\frac{\Gamma\left(\frac{n+3}{2}\right)c_{n}(\mu,k)}{\lambda^{(n+3)/2}}\;.$
(78)

The coefficients $c_{n}$ are

$c_{n}(\mu,k)=\frac{a_{0}^{-(n+3)/2}(\mu)}{2}\sum_{m=0}^{n}\frac{(\mathrm{i}k\mu)^{n-m}}{m!(n-m)!}d_{m,n}(\mu)$
(79)

with

$\displaystyle d_{m,n}(\mu)$ $\displaystyle=\lim_{x\to
0}\frac{\mathrm{d}^{m}}{\mathrm{d}x^{m}}\left(1+\sum_{j=0}\frac{a_{2j}(\mu)}{a_{0}(\mu)}x^{2j}\right)^{-\frac{n+3}{2}}$
$\displaystyle=a_{0}^{(n+3)/2}(\mu)\lim_{x\to
0}\frac{\mathrm{d}^{m}}{\mathrm{d}x^{m}}\left[\frac{f_{\mu}(x)-f_{\mu}(0)}{x^{2}}\right]^{-\frac{n+3}{2}}\;.$
(80)

Since $f_{\mu}(q)$ is even in $q$ (see also Sect. C below), all odd
derivatives vanish, implying $d_{2m+1,n}(\mu)=0$. Thus in the coefficients
$c_{n}(\mu,k)$, only even values of $m$ contribute to the sum. We further add
that, due to the integration over the parameter $\mu$, it is sufficient to
consider only the real part, which is even in $\mu$, i.e. only terms for even
numbers of $n$. This brings us to

$\displaystyle\operatorname{Re}I_{\mu,k}(\lambda)$
$\displaystyle\sim\frac{\mathrm{e}^{\lambda\frac{\sigma_{1}^{2}}{3}}}{2}\,\sum_{n=0}^{\infty}\frac{\Gamma\left(n+\frac{3}{2}\right)}{\lambda^{n+3/2}}\,\sum_{m=0}^{n}\frac{(\mathrm{i}k\mu)^{2n-2m}}{(2m)!(2n-2m)!}\lim_{x\to
0}\frac{\mathrm{d}^{2m}}{\mathrm{d}x^{2m}}\left[\frac{f_{\mu}(x)-f_{\mu}(0)}{x^{2}}\right]^{-n-\frac{3}{2}}\;.$
(81)

This is an asymptotic expansion of the integral for $\lambda\to\infty$. Note
that this is not an asymptotic series for $k\to\infty$, which can be seen as
follows. Each term in the sum over $n$ is a finite sum over $m$. Considering
for each $n$ the $m=0$ term in the second sum, we observe that terms of order
$k^{3}$ are contained in each $n$ term. If, however, the asymptotic series
above was an asymptotic series for $k\to\infty$, then the terms of the first
sum would have to contain terms of descending order in $k$. To get the first
term of the asymptotic series in $k$, we therefore have to to collect all
terms in the total expression that are of order $k^{-3}$. For the second term,
we have to collect all $k^{-5}$-terms and so on. The result of this ordering,
together with substituting $\lambda\to g^{2}_{qp}k^{2}$, is the asymptotic
series

$\operatorname{Re}I_{\mu,k}\left(g_{qp}^{2}k^{2}\right)\sim\frac{\mathrm{e}^{g_{qp}^{2}k^{2}\frac{\sigma_{1}^{2}}{3}}}{2k^{3}}\sum_{m=0}^{\infty}\frac{(-k^{2}\mu^{2})^{-m}}{(2m)!}\tilde{c}_{m}(\mu^{2})$
(82)

for $k\to\infty$ with the coefficients

$\tilde{c}_{m}(\mu^{2}):=\lim_{x\to
0}\frac{\mathrm{d}^{2m}}{\mathrm{d}x^{2m}}\left[\frac{x^{2}/g_{qp}^{2}}{f_{\mu}(x)-f_{\mu}(0)}\right]^{\frac{3}{2}}\sum_{n=0}^{\infty}\frac{\Gamma\left(n+m+\frac{3}{2}\right)}{(2n)!}\left[\frac{-\mu^{2}x^{2}/g_{qp}^{2}}{f_{\mu}(x)-f_{\mu}(0)}\right]^{n+m}\;.$
(83)

We now proceed by performing the $\mu$-integral for each individual term,

$\int_{-1}^{1}\mathrm{d}\mu\,\left(-\mu^{2}\right)^{-m}\tilde{c}_{m}(\mu^{2})=\sqrt{\pi}\lim_{x\to
0}\frac{\mathrm{d}^{2m}}{\mathrm{d}x^{2m}}\left(\frac{x^{2}/g_{qp}^{2}}{a_{1}(x)+\frac{\sigma_{1}^{2}}{3}}\right)^{m+1}\frac{\mathrm{d}^{m}}{\mathrm{d}y^{m}}y^{m}\frac{\sqrt{R(x,y)}}{\mathrm{e}^{R(x,y)}}\;,$
(84)

to find the asymptotic expansion for the power spectrum

$\displaystyle P(k)$
$\displaystyle=\int_{-1}^{1}\mathrm{d}\mu\int_{0}^{\infty}\mathrm{d}q\,q^{2}\mathrm{e}^{-k^{2}g_{qp}^{2}f_{\mu}(q)}\mathrm{e}^{\mathrm{i}k\mu
q}$
$\displaystyle\sim\sqrt{\pi}\frac{\mathrm{e}^{k^{2}\frac{\tau_{1}^{2}}{3}}}{2k^{3}}\sum_{m=0}^{\infty}\frac{k^{-2m}}{(2m)!}\lim_{x\to
0}\frac{\mathrm{d}^{2m}}{\mathrm{d}x^{2m}}\left(\frac{x^{2}/g_{qp}^{2}}{a_{1}(x)+\frac{\sigma_{1}^{2}}{3}}\right)^{m+1}\frac{\mathrm{d}^{m}}{\mathrm{d}y^{m}}y^{m}\frac{\sqrt{R(x,y)}}{\mathrm{e}^{R(x,y)}}\;.$
(85)

The function $R(x,y)$ is defined in (28).

##  Appendix C Asymptotic expansion of the functions $a_{1,2}(q)$ and
$f_{\mu}(q)$

The functions $a_{1,2}(q)$ defined in (4) in terms of the correlation function
$\xi_{\psi}(q)$ (5) of the velocity potential can be written as

$a_{1}(q)=-\frac{1}{2\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,P_{\delta}^{\mathrm{(i)}}(k)\,\frac{j_{1}(kq)}{kq}\;,\quad
a_{2}(a)=qa_{1}^{\prime}(q)\;.$ (86)

Inserting the Taylor series for the first-order spherical Bessel function
$j_{1}(kq)$ leads to the asymptotic series

$\displaystyle a_{1}(q)$
$\displaystyle\sim\sum_{n=0}^{\infty}\frac{(-1)^{n+1}\sigma_{n+1}^{2}}{(3+2n)(1+2n)!}\,q^{2n}\quad\mbox{and}$
$\displaystyle a_{2}(q)$
$\displaystyle\sim\sum_{n=1}^{\infty}\frac{(-1)^{n+1}\,2n\,\sigma_{n+1}^{2}}{(3+2n)(1+2n)!}\,q^{2n}\quad\mbox{as}\quad
q\to 0\;.$ (87)

The function $f_{\mu}(q)$ from (74) thus has the asymptotic series

$f_{\mu}(q)=-\frac{\sigma_{1}^{2}}{3}+\sum_{m=0}^{\infty}\frac{(-1)^{m+2}\sigma_{m+2}^{2}}{(5+2m)(3+2m)!}\left(1+2(m+1)\mu^{2}\right)\,q^{2m+2}\;.$
(88)

Comparing to (70), we can read off $\alpha=2$ and

$a_{0}(\mu)=\frac{\sigma_{2}^{2}}{30}\left(1+2\mu^{2}\right)\;,\quad
a_{1}(\mu)=0\;,\quad
a_{2}(\mu)=-\frac{\sigma_{3}^{2}}{840}\left(1+4\mu^{2}\right)\;.$ (89)

##  Appendix D Analytic toy model – calculations

We use the identity

$\int_{0}^{\infty}\mathrm{d}x\,x^{\nu-\mu+1/2}\sqrt{xy}\,J_{\mu}(ax)\,J_{\nu}(xy)=\frac{2^{\nu-\mu+1}y^{\nu+1/2}}{\Gamma(\mu-\nu)a^{\mu}}\left(a^{2}-y^{2}\right)^{\mu-\nu-1}\,\Theta(a-y)$
(90)

valid for $-1<\operatorname{Re}\nu<\operatorname{Re}\mu$ and $a>0$ to
calculate the correlation functions $a_{1,2}(q)$ for the power spectrum (31).
Setting $\nu=3/2$ or $\nu=5/2$, $\mu=7/2$, $a=k_{0}^{-1}$ and $y=q$ in (90)
results in

$a_{1}(q)=-\frac{1}{2\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,P_{\delta}^{(i)}(k)\frac{j_{1}(kq)}{kq}=\frac{A}{8\pi}\,\left(q^{2}-1\right)\,\Theta(1-q)$
(91)

and

$a_{2}(q)=\frac{1}{2\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,P_{\delta}^{(i)}(k)j_{2}(kq)=\frac{A}{4\pi}\,q^{2}\,\Theta(1-q)\;.$
(92)

Using

$\lim_{q\to 0}\,\frac{j_{1}(kq)}{kq}=\frac{1}{3}$ (93)

in (91), we immediately conclude

$\frac{\sigma_{1}^{2}}{3}=-a_{1}(q=0)=\frac{A}{8\pi}\quad\Rightarrow\quad
A=\frac{8\pi\sigma_{1}^{2}}{3}\;.$ (94)

## References

  * [1] V. Springel, S. D. White, A. Jenkins, C. S. Frenk, N. Yoshida, L. Gao, J. Navarro, R. Thacker, D. Croton, J. Helly _et al._ ,  _Simulations of the formation, evolution and clustering of galaxies and quasars_ ,  nature 435(7042), 629 (2005).
  * [2] D. Xu, V. Springel, D. Sluse, P. Schneider, A. Sonnenfeld, D. Nelson, M. Vogelsberger and L. Hernquist,  _The inner structure of early-type galaxies in the illustris simulation_ ,  MNRAS469(2), 1824 (2017),  10.1093/mnras/stx899,  1610.07605.
  * [3] V. Springel, R. Pakmor, A. Pillepich, R. Weinberger, D. Nelson, L. Hernquist, M. Vogelsberger, S. Genel, P. Torrey, F. Marinacci _et al._ ,  _First results from the illustristng simulations: matter and galaxy clustering_ ,  Monthly Notices of the Royal Astronomical Society 475(1), 676 (2018).
  * [4] K. T. E. Chua, A. Pillepich, M. Vogelsberger and L. Hernquist,  _Shape of dark matter haloes in the illustris simulation: effects of baryons_ ,  MNRAS484(1), 476 (2019),  10.1093/mnras/sty3531,  1809.07255.
  * [5] M. Vogelsberger, F. Marinacci, P. Torrey and E. Puchwein,  _Cosmological simulations of galaxy formation_ ,  Nature Reviews Physics 2(1), 42 (2020).
  * [6] M. Bartelmann, F. Fabis, D. Berg, E. Kozlikin, R. Lilow and C. Viermann,  _A microscopic, non-equilibrium, statistical field theory for cosmic structure formation_ ,  New Journal of Physics 18(4), 043020 (2016),  10.1088/1367-2630/18/4/043020.
  * [7] M. Bartelmann, F. Fabis, E. Kozlikin, R. Lilow, J. Dombrowski and J. Mildenberger,  _Kinetic field theory: effects of momentum correlations on the cosmic density-fluctuation power spectrum_ ,  New Journal of Physics 19(8), 083001 (2017),  10.1088/1367-2630/aa7e6f.
  * [8] M. Bartelmann, E. Kozlikin, R. Lilow, C. Littek, F. Fabis, I. Kostyuk, C. Viermann, L. Heisenberg, S. Konrad and D. Geiss,  _Cosmic structure formation with kinetic field theory_ ,  Annalen der Physik 531(11), 1800446 (2019).
  * [9] M. Bartelmann,  _Trajectories of point particles in cosmology and the zel’dovich approximation_ ,  Physical Review D 91(8), 083524 (2015).
  * [10] M. Bartelmann, J. Dombrowski, S. Konrad, E. Kozlikin, R. Lilow, C. Littek, C. Pixius and F. Fabis,  _Kinetic field theory: Non-linear cosmic power spectra in the mean-field approximation_ ,  SciPost Physics 10(6), 153 (2021),  10.21468/SciPostPhys.10.6.153,  2011.04979.
  * [11] Y. B. Zel’dovich,  _Gravitational instability: an approximate theory for large density perturbations_ ,  A&A500, 13 (1970).
  * [12] J. M. Bardeen, J. Bond, N. Kaiser and A. Szalay,  _The statistics of peaks of gaussian random fields_ ,  The Astrophysical Journal 304, 151 (1986).
  * [13] D. Weinberg and J. Gunn,  _Largescale structure and the adhesion approximation_ ,  Monthly Notices of the Royal Astronomical Society 247, 260 (1990).
  * [14] P. Coles, A. L. Melott and S. F. Shandarin,  _Testing approximations for non-linear gravitational clustering_ ,  Monthly Notices of the Royal Astronomical Society 260(4), 765 (1993),  10.1093/mnras/260.4.765.
  * [15] P. Schneider and M. Bartelmann,  _The power spectrum of density fluctuations in the zel’dovich approximation_ ,  Monthly Notices of the Royal Astronomical Society 273(2), 475 (1995).
  * [16] W. Fulks and J. Sather,  _Asymptotics. ii. laplace’s method for multiple integrals_ ,  Pacific Journal of Mathematics 11(1), 185 (1961).
  * [17] S.-F. Chen and M. Pietroni,  _Asymptotic expansions for the large scale structure_ ,  arXiv preprint arXiv:2002.11357 (2020).
  * [18] R. Wong,  _Asymptotic approximations of integrals_ ,  No. ARRAY(0x560f39c829c0) in Classics in applied mathematics. Society for Industrial and Applied Mathematics, Philadelphia [u.a.], unabridged, corrected republ edn.,  ISBN 0-89871-497-4 and 978-0-89871-497-5,  This SIAM ed. is an unabridged, corrected republication of the work first publ. by Acad. Press, Boston [u.a.], 1989 (2001).
  * [19] G. Doetsch, ed.,  _Anwendungen der Laplace-Transformation : 1. Abt_ ,  No. ARRAY(0x56379abe9680) in Lehrbücher und Monographien aus dem Gebiete der exakten Wissenschaften. Birkhäuser, Basel ; Stuttgart (1955).
  * [20] A. Erdélyi,  _Asymptotic expansions_ ,  Dover Publ., New York, NY,  ISBN 0-486-60318-0 and 978-0-486-60318-6 (1956).
  * [21] A. Erdélyi,  _General asymptotic expansions of laplace integrals_ ,  Archive for Rational Mechanics and Analysis 7(1), 1 (1961),  10.1007/BF00250746.