min}$ instead of $s>0$. This is the case of eq. (6) with $\alpha=0$ (thus
$h=\mathcal{R}$), which reads

$\displaystyle\frac{s_{\rm
min}}{\sigma_{2}}=\frac{9}{8}\frac{(a_{m}H_{m})^{2}}{\sigma_{2}}\delta_{c}\,,$
(84)

Eq. (81) becomes

$\displaystyle n_{\rm max}(\mathcal{R},s_{\rm
min})=\frac{8}{\pi^{2}\sigma_{1}^{2}\sigma_{2}^{2}\sqrt{\sigma_{2}^{2}\sigma_{0}^{2}-\sigma_{1}^{4}}}\int_{s_{\rm
min}}^{\infty}ds\int_{0}^{s}dr\,r(s^{2}-r^{2})\exp\left[-\frac{(\mathcal{R}^{2}\sigma_{2}^{2}+4s^{2}\sigma_{0}^{2}-4s\mathcal{R}\sigma_{1}^{2})}{2(\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4})}-\frac{4r^{2}}{\sigma_{2}^{2}}\right]\,,$
(85)

where, according to the previous argument, we set to $s_{\rm min}$ the lower
limit of integration over $s$ and define $n_{\rm max}(\mathcal{R},s_{\rm
min})$ such that $n_{\rm max}(\mathcal{R},s_{\rm min}=0)=n_{\rm
max}(\mathcal{R})$. The integration over $r$ gives

$\displaystyle n_{\rm max}$ $\displaystyle(\mathcal{R},s_{\rm min})=$ (86)
$\displaystyle\underbrace{\frac{\sigma_{2}^{2}}{4\pi^{2}\sigma_{1}^{2}\sqrt{\sigma_{2}^{2}\sigma_{0}^{2}-\sigma_{1}^{4}}}\int_{s_{\rm
min}}^{\infty}ds\left(\frac{4s^{2}}{\sigma_{2}^{2}}+e^{-\frac{4s^{2}}{\sigma_{2}^{2}}}-1\right)\exp\left\\{-\frac{\sigma_{0}^{2}\sigma_{2}^{2}}{2(\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4})}\left[\frac{4s^{2}}{\sigma_{2}^{2}}-\frac{4s}{\sigma_{2}}\left(\frac{\sigma_{1}^{2}}{\sigma_{2}\sigma_{0}}\right)\frac{\mathcal{R}}{\sigma_{0}}+\frac{\mathcal{R}^{2}}{\sigma_{0}^{2}}\right]\right\\}}_{\equiv\int_{s_{\rm
min}}^{\infty}ds\,\bar{n}_{\rm max}(\mathcal{R},s)}\,,$

where

$\displaystyle\bar{n}_{\rm
max}(\mathcal{R},s)\equiv\frac{\sigma_{2}^{2}}{4\pi^{2}\sigma_{1}^{2}\sqrt{\sigma_{2}^{2}\sigma_{0}^{2}-\sigma_{1}^{4}}}\left(\frac{4s^{2}}{\sigma_{2}^{2}}+e^{-\frac{4s^{2}}{\sigma_{2}^{2}}}-1\right)\exp\left\\{-\frac{\sigma_{0}^{2}\sigma_{2}^{2}}{2(\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4})}\left[\frac{4s^{2}}{\sigma_{2}^{2}}-\frac{4s}{\sigma_{2}}\left(\frac{\sigma_{1}^{2}}{\sigma_{2}\sigma_{0}}\right)\frac{\mathcal{R}}{\sigma_{0}}+\frac{\mathcal{R}^{2}}{\sigma_{0}^{2}}\right]\right\\}\,,$
(87)

can be interpreted as the number density of maxima with field value in the
range $\mathcal{R}$ to $\mathcal{R}+d\mathcal{R}$ and curvature
$-\triangle\mathcal{R}$ in the range $2s$ to $2(s+ds)$. We note that the
argument of the exponential function in $\bar{n}_{\rm max}(\mathcal{R},s)$ is
invariant under the exchange $\mathcal{R}/\sigma_{0}\leftrightarrow
2s/\sigma_{2}$.

From the structure of the argument of the exponential function in
$\bar{n}_{\rm max}(\mathcal{R},s)$, it is natural to introduce the
dimensionless parameter $\gamma\equiv\sigma_{1}^{2}/\sigma_{2}\sigma_{0}$
which is completely determined, as we shall explain in a moment, by the
properties of the power spectrum. Furthermore, it is easy to see that we have
$0<\gamma<1$. In turn, this condition implies
$\sigma_{2}^{2}\sigma_{0}^{2}-\sigma_{1}^{4}>0$ so that the square root in eq.
(87) is always real valued. Let us verify the non-trivial condition
$\sigma_{2}^{2}\sigma_{0}^{2}-\sigma_{1}^{4}>0$; from the definition in eq.
(73), we have

$\displaystyle\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4}=\left[\int\frac{dk}{k}\mathcal{P}_{\mathcal{R}}(k)\right]\left[\int\frac{dk^{\prime}}{k^{\prime}}\mathcal{P}_{\mathcal{R}}(k^{\prime})k^{\prime\,4}\right]-\left[\int\frac{dk}{k}\mathcal{P}_{\mathcal{R}}(k)k^{2}\right]^{2}=\int\frac{dk}{k}\frac{dk^{\prime}}{k^{\prime}}\mathcal{P}_{\mathcal{R}}(k)\mathcal{P}_{\mathcal{R}}(k^{\prime})\left(k^{\prime\,4}-k^{\prime\,2}k^{2}\right)\,.$
(88)

To conclude that $\sigma_{0}^{2}\sigma_{2}^{2}-\sigma_{1}^{4}>0$, we only need
to show that $(k^{\prime\,4}-k^{\prime\,2}k^{2})>0$ since the power spectrum
is positive definite and the integrals over $k$ cover the positive real axis.
In the first double-integral, we can change $k\leftrightarrow k^{\prime}$
without altering the result, and this means that we can also substitute
$k^{\prime\,4}\to(k^{4}+k^{\prime\,4})/2$; if we do this transformation, the
factor $(k^{\prime\,4}-k^{\prime\,2}k^{2})$ becomes
$(k^{\prime\,2}-k^{2})^{2}/2$ which is always positive. This concludes the
proof.

From eq. (87) we see that the value of $\gamma$ controls the amount of
correlation between $\mathcal{R}$ and $s$. If we take $\gamma=0$, the two
variables are completely uncorrelated (because $\sigma_{1}^{2}=0$). On the
contrary, $\gamma=1$ corresponds to the case in which they are maximally
correlated.

To proceed further, we assume the following analytical expression for the
power spectrum

$\displaystyle\mathcal{P}_{\mathcal{R}}(k)=\frac{A_{g}}{\sqrt{2\pi}v}\exp\left[-\frac{\log^{2}(k/k_{\star})}{2v^{2}}\right]\,.$
(89)

As far as the spectral moments in eq. (73) are concerned, we find the
analytical result

$\displaystyle\sigma_{j}^{2}=A_{g}k_{\star}^{2j}e^{2j^{2}v^{2}}\,,$ (90)

which implies $\sigma_{0}^{2}=A_{g}$ and $\gamma=e^{-2v^{2}}$, where we see
that $0<\gamma<1$ as expected. The three parameters $\\{A_{g},v,k_{\star}\\}$
fully specify our problem. The physical picture is the following.

  * $\circ$

The scale $k_{\star}$ represents the comoving wavenumber at which the power
spectrum peaks. This quantity is related to the mass of the black holes
produced after the collapse of the regions where the overdensity field is
above threshold. We take $k_{\star}=O(10^{14})$ Mpc-1 corresponding to $M_{\rm
PBH}=O(10^{18})$ g.

  * $\circ$

The parameter $v$ controls the broadness of the power spectrum and, in turn,
the broadness of the mass distribution of the black holes. The case $v=0.1$,
for instance, corresponds to a very narrow power spectrum. This, in turn, will
generate a very narrow mass distribution of black holes.

  * $\circ$

The amplitude of the power spectrum $A_{g}$ is related to the abundance of
dark matter in the present-day Universe in the form of black holes, and
typical values are of order $A_{g}/\sqrt{2\pi}v=O(10^{-2})$. To fix ideas, to
get $\mathcal{P}_{\mathcal{R}}(k_{\star})=10^{-2}$ one needs $A_{g}=2.5\times
10^{-3}$ for $v=0.1$.

Let us now pause for a moment to clarify the rationale of the computation that
we are doing. We are interested in regions of space where the field
$\mathcal{R}$ has large curvature $-\triangle\mathcal{R}$ since in this case
the overdensity field (which is proportional to $-\triangle\mathcal{R}$) takes
large values. In eq. (87) we can select regions with large curvature by
implementing (as done in eq. (86)) a lower limit of integration over $s$.
However, eq. (87) always associates, by construction, regions with large
curvature (consequently, peaks of the overdensity field) with local maxima of
$\mathcal{R}$. This association can be analytically justified as follows.
Consider the Taylor expansion in eq. (66) which we rewrite in two spatial
dimensions taking a local maximum as stationary point of $\mathcal{R}$

$\displaystyle\mathcal{R}(\vec{x})=\mathcal{R}_{\rm
M}+\frac{1}{2}\sum_{i,j=1}^{2}\mathcal{R}_{ij}(\vec{x}_{\rm
M})(\vec{x}-\vec{x}_{\rm M})_{i}(\vec{x}-\vec{x}_{\rm M})_{j}\,.$ (91)

This quadratic equation can be written in a canonical form (that is without
cross terms) if we do a coordinate transformation
$\vec{x}\to\vec{x}^{\,\prime}\equiv(x^{\prime},y^{\prime})$ such that the new
axes are aligned along the eigenvectors of the matrix
$\mathcal{R}_{ij}(\vec{x}_{\rm M})$. We are free to do this rotation because
of isotropy. In such case, eq. (91) takes the form

$\displaystyle\mathcal{R}(\vec{x}^{\,\prime})=\mathcal{R}_{\rm
M}-\frac{1}{2}(\lambda_{1}x^{\prime\,2}+\lambda_{2}y^{\prime\,2})\,,$ (92)

where $\lambda_{i=1,2}>0$ are the (minus) eigenvalues of
$\mathcal{R}_{ij}(\vec{x}_{\rm M})$. Notice that in eq. (92) we used
homogeneity to shift the position of the maximum $\vec{x}_{\rm M}$ to the
origin of the new coordinate system. From eq. (92) we see that an iso-density
surface with constant
$\mathcal{R}_{\vec{x}^{\prime}}\equiv\mathcal{R}(\vec{x}^{\,\prime})$ is an
ellipse with canonical equation

$\displaystyle\frac{\lambda_{1}}{2(\mathcal{R}_{\rm
M}-\mathcal{R}_{\vec{x}^{\prime}})}x^{\prime\,2}+\frac{\lambda_{2}}{2(\mathcal{R}_{\rm
M}-\mathcal{R}_{\vec{x}^{\prime}})}y^{\prime\,2}=1\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}a_{i}\equiv\bigg{[}\frac{2(\mathcal{R}_{\rm
M}-\mathcal{R}_{\vec{x}^{\prime}})}{\lambda_{i}}\bigg{]}^{1/2}\,,$ (93)

and semi-axes $a_{i=1,2}$. The key point is that the actual magnitude of the
eigenvalues $\lambda_{i}$ depends on the steepness of the field $\mathcal{R}$
around the position of its maximum. This follows from eq. (92) if we apply the
Laplacian with respect to the coordinates $\vec{x}^{\,\prime}$ since we find

$\displaystyle\lambda_{1}+\lambda_{2}=-\triangle\mathcal{R}(\vec{x}^{\,\prime})\,.$
(94)

Suppose now that the point $\vec{x}^{\,\prime}$ coincides with a peak of the
overdensity field, $\vec{y}_{\rm pk}$. Using eq. (6), we find

$\displaystyle\lambda_{1}+\lambda_{2}=-\triangle\mathcal{R}(\vec{y}_{\rm
pk})\gtrsim\frac{9}{4}(aH)^{2}\delta_{c}=2\sigma_{2}\underbrace{\frac{9}{8}\frac{(aH)^{2}}{\sigma_{2}}\delta_{c}}_{\gg
1}~{}~{}~{}~{}~{}~{}\Longrightarrow~{}~{}~{}~{}~{}~{}\frac{\lambda_{i}}{\sigma_{2}}\gg
1\,.$ (95)

In eq. (95) we used the fact that we typically expect $(aH)^{2}/\sigma_{2}\gg
1$; we will comment in more detail about this estimate at the end of this
section. Notice that in eq. (95) we assumed that the two eigenvalues
$\lambda_{i}$ have the same magnitude. This is because we have separately
$\lambda_{1}=-\mathcal{R}_{xx}(\vec{y}_{\rm pk})$ and
$\lambda_{2}=-\mathcal{R}_{yy}(\vec{y}_{\rm pk})$, and the second derivatives
$\mathcal{R}_{xx}$ ad $\mathcal{R}_{yy}$ have the same covariance (see eq.
(72)). We can now do the same expansion in eq. (91) but with respect to the
first derivatives of $\mathcal{R}$ at $\vec{x}_{\rm M}$. The stationary
condition $\mathcal{R}_{i}(\vec{x}_{\rm M})=0$ reads

$\displaystyle\mathcal{R}_{i}(\vec{x})-\sum_{j=1}^{2}\mathcal{R}_{ij}(\vec{x}_{\rm
M})(\vec{x}-\vec{x}_{\rm
M})_{j}=0~{}~{}~{}~{}~{}~{}\Longrightarrow~{}~{}~{}~{}~{}~{}\left\\{\begin{array}[]{c}\mathcal{R}_{x}(\vec{x}^{\,\prime})+\lambda_{1}x^{\prime}=0\\\
\\\
\mathcal{R}_{y}(\vec{x}^{\,\prime})+\lambda_{2}y^{\prime}=0\end{array}\right.$
(99)

where in the last line we introduced, as done before, the eigenvalues
$\lambda_{i=1,2}$. We identify again the point $\vec{x}^{\,\prime}$ with a
peak of the overdensity field so that we can use the estimate in eq. (95).
Furthermore, $\vec{x}^{\,\prime}$ represents now, by construction, the
distance between the local maximum of $\mathcal{R}$ and the peak of the
overdensity field. We can estimate this distance by means of eq. (95) and eq.
(99). In eq. (99), $\mathcal{R}_{x}(\vec{x}^{\,\prime})$ and
$\mathcal{R}_{y}(\vec{x}^{\,\prime})$ are not equal to zero (since we moved
away from the local maximum of $\mathcal{R}$) and their magnitude can be
estimated (in a probabilistic sense) by means of the covariance matrix. We
have
$\mathcal{R}_{x}(\vec{x}^{\,\prime})\approx\mathcal{R}_{y}(\vec{x}^{\,\prime})\approx\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle^{1/2}\approx\sigma_{1}$.
We find

$\displaystyle|\vec{x}^{\,\prime}|\approx\frac{\sigma_{1}}{\sigma_{2}}\left(\frac{1}{\lambda_{i}/\sigma_{2}}\right)\ll\frac{\sigma_{1}}{\sigma_{2}}=\frac{e^{-3v^{2}}}{k_{\star}}\,.$
(100)

As we will see in a moment (see also appendix F), we have that $1/k_{\star}$
is typically of the order of the comoving horizon length $1/aH$ at the time
when the perturbations re-enter the horizon and become causally connected.
Therefore, we find $|\vec{x}^{\,\prime}|\ll 1/aH$. This means that high peaks
of the overdensity field lie “close” (that is within an Hubble radius) to
local maxima of the curvature perturbation.

Figure 7:  Numerical simulation of the random field $\mathcal{R}$ together
with its first and second derivatives. At each point in space (discretized in
steps $\Delta x=5$ and $\Delta y=5$) we associate a vector of values
$\\{\mathcal{R},\mathcal{R}_{x},\mathcal{R}_{y},\mathcal{R}_{xx},\mathcal{R}_{xy},\mathcal{R}_{yy}\\}$
randomly generated from eq. (70). We use the power spectrum in eq. (89) to
compute the correlation matrix. We set $v=0.7$, $k_{\star}=1.5\times 10^{14}$
Mpc-1 and $A_{g}=2.5\times 10^{-3}$. Left panel. We show the density plot of
the random variable $2s/\sigma_{2}=-\triangle\mathcal{R}/\sigma_{2}$ which is
related to the overdensity field. Right panel. Same as in the left panel but
zoomed in the region where the random variable
$-\triangle\mathcal{R}/\sigma_{2}$ has a pronounced peak. We superimpose (blue
arrows) a vector plot that keeps track of the gradient field
$\\{\mathcal{R}_{x},\mathcal{R}_{y}\\}$. The yellow star marks the position of
the local maximum of $\mathcal{R}$ that lies close to the peak of
$-\triangle\mathcal{R}/\sigma_{2}$.

We can validate the analytical approximation by means of a numerical check. To
this end, we use the full joint probability density distribution in eq. (70)
to generate a sample of random values that we distribute on a two-dimensional
grid (see caption of fig. 7). In other words, at each point on the spatial
grid corresponds a value of
${R}=(\mathcal{R},\mathcal{R}_{x},\mathcal{R}_{y},\mathcal{R}_{xx},\mathcal{R}_{xy},\mathcal{R}_{yy})^{\rm
T}$ randomly generated from eq. (70). In the left panel of fig. 7 we show the
spatial distribution of the random variable
$2s/\sigma_{2}=-\triangle\mathcal{R}/\sigma_{2}$. We focus on a region in
which $-\triangle\mathcal{R}$ takes a large value (in units of $\sigma_{2}$).
In the left panel of fig. 7, we indicate this region with a red contour. We
zoom in this part of the plot in the right panel of fig. 7. The analytical
argument explained before suggests that we should find a maximum close to the
point where the curvature field peaks. We look for this maximum numerically by
looking at the behavior of the gradient field
$\\{\mathcal{R}_{x},\mathcal{R}_{y}\\}$ that we plot using blue arrows. We
indeed find a local maximum that we mark with a yellow star. We checked
numerically that at the position of the yellow star where the gradient field
vanishes the conditions on the second derivatives that define a maximum are
verified. We note that at the position of the peak of the overdensity field
the gradient field $\\{\mathcal{R}_{x},\mathcal{R}_{y}\\}$ does not vanish but
we find that its magnitude is of order $O(1)$ (in units of $\sigma_{1}$) as we
argued in eq. (100). Furthermore, we checked that the local maximum lies
closer to the peak of $-\triangle\mathcal{R}/\sigma_{2}$ for increasing higher
values of the latter.

This numerical result corroborates the validity of the analytical argument in
eq. (100), and we conclude, therefore, that eq. (87) is the right distribution
to consider: We count the peaks of $\delta$ by looking at the maxima of
$\mathcal{R}$ with large curvature.

Next, we ask if there exists some relation between the curvature of a local
maximum of $\mathcal{R}$ and the value of $\mathcal{R}$ at the
maximum.222222The result of this exercise will be useful later to discuss the
impact of non-linearities in eq. (5), see appendix G. We can use the
probability density distribution $\bar{n}_{\rm max}(\mathcal{R},s)$ defined in
eq. (87) to generate numerically, in position space, a sample of maxima by
extracting randomly the value of $\mathcal{R}$ and curvature
$2s=-\triangle\mathcal{R}$. We can then use eq. (5) to extract from the
distribution of $s$ the distribution of the overdensity field (of course, with
$\mathcal{R}$ instead of $h$ since we are considering here the gaussian case).
The outcome of this exercise is shown in fig. 8, fig. 9 and fig. 10 (see
captions for details).

Let us consider the case in which we take $v=0.1$. We remind that this choice
corresponds to a very narrow power spectrum. We have $\gamma\simeq 0.98$. As
noticed before, in this case we expect a strong correlation between
$\mathcal{R}$ and $s$. This means that regions with large
$\mathcal{R}/\sigma_{0}$ are likely to be also regions with large
$2s/\sigma_{2}$ (and, consequently, regions where the overdensity field
peaks).

Figure 8:  Numerical simulation of maxima of $\mathcal{R}$ in two spatial
dimensions. At each point $(x,y)$ in space (discretized in steps $\Delta x=10$
and $\Delta y=10$), we extract randomly from the density distribution
$\bar{n}_{\rm max}(\mathcal{R},s)$ (defined in eq. (87)) the value of
$\mathcal{R}$ and $s=-\triangle\mathcal{R}/2$. The former are shown in the
left panel. We remark that, by construction, all point generated by means of
$\bar{n}_{\rm max}(\mathcal{R},s)$ are maxima of $\mathcal{R}$. As far as the
values of $s$ are concerned, we plot on the right panel the corresponding
values of $\delta=(8/9)(1/aH)^{2}s$. We use the power spectrum in eq. (89),
and we set $v=0.1$, $k_{\star}=1.5\times 10^{14}$ Mpc-1 and $A_{g}=2.5\times
10^{-3}$. Furthermore, we take $(k_{\star}/aH)^{2}\simeq 0.36$ as suggested by
numerical simulations (see appendix F). For illustrative purposes, we set the
threshold $\delta_{c}=0.675$ (which is significantly smaller compared to the
value expected from numerical simulation of gravitational collapse into black
holes that is $\delta_{c}\simeq 1.19$, see appendix F; we use a smaller value
of $\delta_{c}$ otherwise events over the threshold would be too rare to be
simulated in our simplified numerical analysis). Points in the simulation with
$\delta>\delta_{c}$ are marked with a black dot.

This is evident in fig. 8, fig. 9 and fig. 10 (left panel).

Figure 9:  Same as in fig. 8 but zoomed in the region delimited by dashed
lines and displayed in the form of a tri-dimensional density plot. We see that
spiky maxima with large values of $\mathcal{R}$ (left panel) coincide with
peaks of $\delta$ (right panel).

In particular, in fig. 9 we see that maxima with large values of $\mathcal{R}$
(left panel) maps precisely regions where the overdensity field $\delta$ peaks
(right panel). As stated before, this is a consequence of the fact that
$\mathcal{R}$ and $s$ (hence $\delta$) are highly correlated.

Figure 10:  Left panel. Simulated points in figs. 8, 9 shown in the plane
$\\{\mathcal{R}/\sigma_{0},-\triangle\mathcal{R}/2\sigma_{2}\\}$. The
horizontal dashed line corresponds to the threshold value $\delta_{c}$ (see
caption of fig. 8). Points above threshold collapse into black holes (marked
with black dots here and in figs. 8, 9). This simulation corresponds to
$v=0.1$, and we see that $\mathcal{R}/\sigma_{0}$ and
$-\triangle\mathcal{R}/2\sigma_{2}$ are strongly correlated. Right panel.
Comparison between two numerical simulations with $v=0.1$ (green) and $v=0.7$
(red). The two ellipses represent the 95% confidence regions.

Furthermore, in the left panel of fig. 10 the symmetry of the randomly
generated points under the exchange $\mathcal{R}/\sigma_{0}\leftrightarrow
2s/\sigma_{2}$ is evident.

It is instructive to consider what happens if we take a different value of
$v$. Consider, for instance, the case with $v=0.7$. We have $\gamma\simeq
0.37$, and $\mathcal{R}$ and $s$ are now much less correlated compared to the
case with $v=0.1$. This is shown in the right panel of fig. 10 (see caption
for details). For $v=0.7$, spiky maxima (that is maxima with large curvature)
are seldom characterized also by a large value of $\mathcal{R}$. This means
that if we aim at deriving a generic formula for the number density of spiky
maxima we can not restrict eq. (86) to special values of $\mathcal{R}$.

Let us summarize our findings so far. We have that the number density of local
maxima of $\mathcal{R}$ with large curvature gives the number density of
regions where the overdensity field peaks. Physically, we are only interested
in local maxima with large curvature irrespectively on their value of
$\mathcal{R}$. Since there is no special relation—for generic values of $v$
(hence $\gamma$)—between $s$ and $\mathcal{R}$, eq. (86) must be integrated
over the entire range of variability of $\mathcal{R}$.

If we consider the case in which $s_{\rm min}$ does not depend on
$\mathcal{R}$, we can integrate eq. (86) analytically. We find the number
density

$\displaystyle\mathcal{N}_{\rm max}(s_{\rm
min})\equiv\int_{-\infty}^{\infty}d\mathcal{R}\int_{s_{\rm
min}}^{\infty}ds\,\bar{n}_{\rm
max}(\mathcal{R},s)=\frac{\sigma_{2}^{2}}{2\sqrt{2}\pi^{3/2}\sigma_{1}^{2}}\left[\frac{s_{\rm
min}}{\sigma_{2}}\exp\left(-\frac{2s_{\rm
min}^{2}}{\sigma_{2}^{2}}\right)+\frac{1}{2}\sqrt{\frac{\pi}{6}}{\rm
Erfc}\left(\frac{\sqrt{6}s_{\rm min}}{\sigma_{2}}\right)\right]\,.$ (101)

We can define the dimensionful quantity
$R_{*}\equiv\sqrt{d}\,\sigma_{1}/\sigma_{2}$, where $d$ is the number of
spatial dimensions. From eq. (90) we find

$\displaystyle
R_{*}=\frac{\sqrt{2}\sigma_{1}}{\sigma_{2}}=\frac{\sqrt{2}}{k_{\star}}e^{-3v^{2}}\,.$
(102)

In this simplified two-dimensional set-up, we can estimate the mass fraction
of black holes by means of the dimensionless quantity (for more details, see
appendix F)

$\displaystyle R_{*}^{2}\mathcal{N}_{\rm max}(s_{\rm
min})=\frac{1}{\sqrt{2}\pi^{3/2}}\left[\frac{s_{\rm
min}}{\sigma_{2}}\exp\left(-\frac{2s_{\rm
min}^{2}}{\sigma_{2}^{2}}\right)+\frac{1}{2}\sqrt{\frac{\pi}{6}}{\rm
Erfc}\left(\frac{\sqrt{6}s_{\rm
min}}{\sigma_{2}}\right)\right]\simeq\frac{1}{\sqrt{2}\pi^{3/2}}\left(\frac{s_{\rm
min}}{\sigma_{2}}\right)\exp\left(-\frac{2s_{\rm
min}^{2}}{\sigma_{2}^{2}}\right)\,,$ (103)

where the last approximation is valid if $s_{\rm min}/\sigma_{2}\gtrsim 1$.
The mass fraction is controlled by an exponential decaying function with
argument

$\displaystyle\frac{s_{\rm
min}}{\sigma_{2}}=\frac{9}{8}\frac{(aH)^{2}}{\sigma_{2}}\delta_{c}=\frac{9e^{-4v^{2}}}{8A_{g}^{1/2}}\left(\frac{aH}{k_{\star}}\right)^{2}\delta_{c}\,.$
(104)

The value of $s_{\rm min}/\sigma_{2}$, which is crucial for the determination
of the correct order-of-magnitude of the mass fraction, depends on the
properties of the power spectrum via the factor
$e^{-4v^{2}}/A_{g}^{1/2}k_{\star}^{2}$ and the details of the gravitational
collapse that leads to black hole formation via the factor
$(aH)^{2}\delta_{c}$. The comoving horizon length $1/aH$ depends on time. The
computation of the threshold for black hole production introduces the time
$t_{m}$ that is defined by the time when the curvature perturbations cross the
horizon and become causally connected. This is the time at which eq. (104) has
to be computed.232323More precisely, the process of black hole formation
involves three different times. First of all, the number density of peaks that
eventually form black holes has to be calculated at some initial time when
perturbations are still super-horizon. All the analysis done so far is based
on this time (even though we do not specify it explicitly) since we are
considering comoving curvature perturbations which are constant (see appendix
A). Second, we have the time $t_{m}$ when the perturbations re-enter the
horizon and become causally connected. Finally, we have the time $t_{f}$ at
which black holes form. In general $t_{f}\neq t_{m}$, and from simulations of
gravitational collapse in numerical relativity we have
$(t_{f}/t_{m})^{1/2}\simeq 3$ Musco:2018rwt ; eq. (6) is usually evaluated at
time $t_{m}$ while the factor $(t_{f}/t_{m})^{1/2}\simeq 3$ can be included in
the computation of the black hole abundance (see, e.g., discussion in ref.
Kalaja:2019uju ). The time $t_{m}$ defines implicitly, by means of the
condition $a(t_{m})H(t_{m})r_{m}=1$, the length scale $r_{m}$. This length
scale enters in the numerical evaluation of eq. (104), and one typically gets
$(a_{m}H_{m}/k_{\star})^{2}\delta_{c}=O(1)$. A precise numerical evaluation of
this factor is needed in order to set the size of the exponential suppression
in eq. (103). We postpone to appendix F a more detailed discussion about this
point.

It is important to check the validity of eq. (103) against the standard
result. In ref. Bardeen:1985tr , the number density of peaks of the
overdensity field is controlled by the exponential function
$\exp(-\delta_{c}^{2}/2\sigma_{\delta}^{2})$ where $\sigma_{\delta}^{2}$
refers to the variance of the overdensity field. By means of eq. (5)
(linearized, and with $\mathcal{R}$ instead of $h$), one finds
$\sigma_{\delta}^{2}=(16/81)(1/aH)^{4}\sigma_{2}^{2}$. Consequently,
$\exp(-\delta_{c}^{2}/2\sigma_{\delta}^{2})$ matches precisely the exponential
function in eq. (103). This is another indication that computing the number
density of peaks of the overdensity field via the number density of maxima of
$\mathcal{R}$ that are spiky enough leads to sensible results; importantly,
this is in spite of the fact that we derived eq. (103) in two spatial
dimensions.

After this long and detailed discussion about the gaussian case, we are ready
to move to the more interesting situation in which local non-gaussianities are
present. Actually, we are already in the position to make an interesting
comment. Suppose that we compute the analogue of eq. (87) with local non-
gaussianities with of course now $h$ and $-\triangle h$ instead of
$\mathcal{R}$ and $-\triangle\mathcal{R}$. It is crucial to answer the same
question that we asked in the gaussian case: Is the number density of maxima
of $h$ with large curvature $-\triangle h$ a good proxy for the number density
of peaks of the overdensity field? The same analytical argument discussed in
eqs. (91-100) can be repeated with $h$ and $-\triangle h$ instead of
$\mathcal{R}$ and $-\triangle\mathcal{R}$. The estimate of the eigenvalues
$\lambda_{i=1,2}$ in eq. (95) remains the same. In eq. (100), the only
difference is that local non-gaussianities alter the entries of the covariance
matrix that we used to estimate the magnitude of the gradient field. However,
we anticipate that we find (see eq. (204)) $\langle h_{x}h_{x}\rangle=\langle
h_{y}h_{y}\rangle=\sigma_{1}^{2}(1+4\alpha^{2}\sigma_{0}^{2})/2$ and $\langle
h_{x}h_{y}\rangle=0$; we conclude that the non-gaussian correction in this
case is negligible since we have $\alpha^{2}\sigma_{0}^{2}=\alpha^{2}A_{g}\ll
1$. Remember indeed that the size of $\sigma_{0}^{2}$ is controlled by the
amplitude of the power spectrum which is $\sigma_{0}^{2}=A_{g}\ll 1$. This is
particularly clear for the simple choice of $\mathcal{P}_{\mathcal{R}}$ in eq.
(89) but remains true in general.

We conclude that the same argument discussed in eqs. (91-100) is valid also in
the presence of local non-gaussianities. In appendix C and appendix D we will,
therefore, move to compute the number density of maxima of $h$ with curvature
above the threshold for black hole production that is needed to generalize eq.
(103) to the case $\alpha\neq 0$. We consider two different approaches.

  * $\ast$

In appendix C we come back to the formulation of the problem in three spatial
dimensions. We will derive an expression analogue to eq. (103) but valid in
the case $\alpha\neq 0$ and three spatial dimensions. We dub the result of
appendix C “exact” because no approximations will be used throughout the
computation (apart from the linearization in eq. (5); the case in which also
non-linearities are included will be discussed in appendix G).

  * $\ast$

In appendix D we follow a different route. We will consider an expansion in
cumulants around the gaussian probability density distribution. To make this
approach more transparent from the analytic point of view, we will work again
in two spatial dimensions.

##  Appendix C Peak statistics with local non-gaussianities

In appendix B we have shown that the number density of peaks of the
overdensity field can be approximated with good accuracy with the number
density of maxima of the comoving density perturbation which are spiky enough,
i.e. with a Laplacian smaller than a threshold. The analysis has been
performed assuming two spatial dimensions, to a have a better control of the
analytic expressions, and for ease of visualization of our simulations.
Nevertheless, the same conclusion is also valid in three dimensions, which is
the case of physical interest that we are going to consider in this section.
We shall now use these results to compute the number density of peaks of the
overdensity field in presence of local non-gaussianities. Our starting point
is the analogous of eq. (310) for the non-gaussian random field $h(\vec{x})$.
Furthermore, as explained above, we are interested into maxima of $h$ which
are spiky enough. Explicitly, we have

$\displaystyle n_{\rm pk}(h)\,dh\,d^{3}x=dh\,d^{3}x\,\int_{\rm
spiky\,max}d^{3}h_{i}\,d^{6}h_{ij}\,P_{\rm
NG}(h,h_{i},h_{ij})\,\delta^{3}(h_{i})\,|{\rm det}(h_{ij})|\,,$ (105)

where $P_{\rm NG}(h,h_{i},h_{ij})$ is the joint probability distribution
function of the non-gaussian field $h,$ the field gradient $h_{i}$ and the
second derivatives $h_{ij}.$ The Dirac delta $\delta^{3}(h_{i})$ enforces the
condition that the point is stationary. We will show in a moment how to
restrict the integration volume to the field configurations which represent
spiky maxima of $h$. The comoving number density of peaks is obtained
integrating over all the heights of the peaks

$\displaystyle\mathcal{N}_{\rm pk}$ $\displaystyle=\int_{h_{\rm
min}}^{\infty}dh\,n_{\rm pk}(h)=\int_{\rm
spiky\,max}dh\,d^{3}h_{i}\,d^{6}h_{ij}\,P_{\rm
NG}(h,h_{i},h_{ij})\,\delta^{3}(h_{i})\,|{\rm det}(h_{ij})|\,.$ (106)

Notice that for $\alpha>0$, which is the case relevant for the PBH production
(see sec. A), and from the relation
$h(\vec{x})=\mathcal{R}(\vec{x})+\alpha\left(\mathcal{R}(\vec{x})^{2}-\sigma_{0}^{2}\right),$
one realizes that $h(x)$ attains a minimum $h_{\rm
min}=-(1+4\alpha^{2}\sigma_{0}^{2})/4\alpha$ for $\mathcal{R}=-1/2\alpha.$ To
proceed, it turns out to be useful to consider the gaussian variables
$\mathcal{R},$ $\mathcal{R}_{i}$ and $\mathcal{R}_{ij}$ instead of the non-
gaussian fields $h$, $h_{i}$ and $h_{ij}.$ The reason is that the joint
probability distribution function of the former variables, that we denote with
$P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})$, is known and it has a
simple analytic expression, see sec. B for the explicit formula in the case of
two dimensions. Using the conservation of the probability in a differential
volume $P_{\rm
NG}(h,h_{i},h_{ij})\,dh\,d^{3}h_{i}\,d^{6}h_{ij}=P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})\,d\mathcal{R}\,d^{3}\mathcal{R}_{i}\,d^{6}\mathcal{R}_{ij},$
we can perform a change a variable and write

$\displaystyle\mathcal{N}_{\rm pk}$ $\displaystyle=\int_{\rm
spiky\,max}d\mathcal{R}\,d^{3}\mathcal{R}_{i}\,d^{6}\mathcal{R}_{ij}\,P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})\,\delta^{3}\left[h_{i}\left(\mathcal{R},\mathcal{R}_{i},\alpha\right)\right]\,\left|{\rm
det}\left[h_{ij}\left(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij},\alpha\right)\right]\right|\,.$
(107)

where the field gradient and the second derivatives of the non gaussian
variables are written in terms of the gaussian fields as

$\displaystyle
h_{i}\left(\mathcal{R},\mathcal{R}_{i},\alpha\right)=\mathcal{R}_{i}\left(1+2\alpha\mathcal{R}\right)\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}h_{ij}\left(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij},\alpha\right)$
$\displaystyle=\mathcal{R}_{ij}(1+2\alpha
R)+2\alpha\mathcal{R}_{i}\mathcal{R}_{j}\,.$ (108)

As noticed in sec. II, stationary points of $\mathcal{R}$ are also stationary
points of $h.$ Moreover, from the expressions above, one can see that also the
configuration $(1+2\alpha R)=0$ is a stationary point of $h.$ However, as
mentioned above, it is a minimum of $h$ (for $\alpha>0$), therefore we need to
focus only on the configurations with $\mathcal{R}_{i}=0.$ This means that one
can write eq. (107) as

$\displaystyle\mathcal{N}_{\rm pk}$ $\displaystyle=\int_{\rm
spiky\,max}d\mathcal{R}\,d^{6}\mathcal{R}_{ij}\,P(\mathcal{R},\mathcal{R}_{i}=0,\mathcal{R}_{ij})\,\left|{\rm
det}\left(h_{ij}\left(\mathcal{R},\mathcal{R}_{i}=0,\mathcal{R}_{ij},\alpha\right)\right)\right|\,\frac{1}{\left|1+2\alpha\mathcal{R}\right|^{3}}$
$\displaystyle=\int_{\rm
spiky\,max}d\mathcal{R}\,d^{6}\mathcal{R}_{ij}\,P(\mathcal{R},\mathcal{R}_{i}=0,\mathcal{R}_{ij})\,\left|{\rm
det}\left(\mathcal{R}_{ij}\right)\right|\,.$ (109)

In the equation above all the information about the non-gaussianities is
confined in ${}_{\rm spiky\,max},$ i.e. in the condition to impose on the
integration volume in order to restrict on maxima of $h$ spiky enough. Let us
see how these constraints can be implemented explicitly. The Hessian matrix,
for $\mathcal{R}_{i}=0$, is

$\displaystyle h_{ij}=\mathcal{R}_{ij}\left(1+2\alpha\mathcal{R}\right)\,.$
(110)

We should, therefore, select maxima (minima) of $\mathcal{R}$ for positive
(negative) values of $\left(1+2\alpha\mathcal{R}\right)$. The calculation can
be simplified aligning the coordinate axes along the eigenvectors of the
matrix $-\mathcal{R}_{ij}.$ The corresponding eigenvalues are denoted as
$\lambda_{i}$ with $i=1,2,3$ (the same strategy has been adopted in sec. B in
the case of two spatial dimensions). The three eigenvectors, and the three
Euler angles defining their orientation, can be used to parametrize the six
independent random fields of $\mathcal{R}_{ij}.$ Let us also define the
following variables

$\displaystyle\bar{\nu}=\mathcal{R}/\sigma_{0}\,,~{}~{}~{}~{}~{}~{}\left\\{\begin{array}[]{ccc}\sigma_{2}\,x&=&-\triangle\mathcal{R}=\lambda_{1}+\lambda_{2}+\lambda_{3}\\\
\sigma_{2}\,y&=&(\lambda_{1}-\lambda_{3})/2\\\
\sigma_{2}\,z&=&(\lambda_{1}-2\lambda_{2}+\lambda_{3})/2\end{array}\right.$
(114)

Changing variables and integrating over the Euler angles, from eq. (C) one
obtains Bardeen:1985tr

$\displaystyle\mathcal{N}_{\rm pk}=\int_{\rm spiky\,max}\bar{n}_{\rm
pk}(\bar{\nu},x,y,z)d\bar{\nu}\,dx\,dy\,dz\,~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}~{}~{}\bar{n}_{\rm
pk}(\bar{\nu},x,y,z)=A\,e^{-Q}\,|F(x,y,z)|\,,$ (115)

and

$\displaystyle A$
$\displaystyle\equiv\sqrt{\frac{5}{3}}\,\frac{\sigma_{2}^{3}}{\sigma_{1}^{3}}\frac{25}{16\pi^{2}\,\sqrt{1-\gamma^{2}}}\,,$
(116) $\displaystyle Q$
$\displaystyle\equiv\frac{\bar{\nu}^{2}}{2}+\frac{(x-x_{*})^{2}}{2(1-\gamma^{2})}+\frac{5}{2}\left(3y^{2}+z^{2}\right)\,,$
(117) $\displaystyle F(x,y,z)$ $\displaystyle\equiv
y\left(y^{2}-z^{2}\right)\left[\left(x+z\right)^{2}-9y^{2}\right]\left(x-2z\right)\,,$
(118) $\displaystyle x_{*}$ $\displaystyle\equiv\gamma\bar{\nu}\,.$ (119)

We can divide the integration on $\bar{\nu}$ in two parts.

  * $\circ$

$1+2\alpha\mathcal{R}>0$.
This implies that maxima of $\mathcal{R}$ are also maxima of $h$, and that we
are considering $\bar{\nu}>-1/2\sigma_{0}\alpha.$ One can choose an ordering
for the eigenvalues of the matrix $-\mathcal{R}_{ij}$:
$\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}.$ Therefore, we can select maxima
of $\mathcal{R}$ requiring that $\lambda_{3}>0.$ Under these conditions, and
using eqs. (114), the domain of integration for the variables $y$ and $z$
reads:

$\displaystyle\int_{0}^{x/4}dy\,\int_{-y}^{y}dz+\int_{x/4}^{x/2}dy\,\int_{3y-x}^{y}dz\,.$
(120)

Then, working within the linear approximation in eq. (5), eq. (6) implies that
only spiky maxima should be selected

$\displaystyle
x>\frac{9(a_{m}H_{m})^{2}}{4\sigma_{2}}\frac{\delta_{c}}{1+2\alpha\sigma_{0}\bar{\nu}}\equiv
x_{\delta}(\bar{\nu})\,.$ (121)

In the equation above, the horizon scale $1/a\,H$ has been evaluated at the
time $t_{m}$ when the curvature perturbations cross the horizon. Having
specified the appropriate domain of integration for all the variables, we can
now integrate eq. (115) and multiply by a factor 6 to take into account all
the possible orderings of the eigenvalues $\lambda_{i}$. We find

$\displaystyle\mathcal{N}_{\rm pk}^{({\rm I})}$
$\displaystyle=\int_{-\frac{1}{2\alpha\sigma_{0}}}^{\infty}d\bar{\nu}\,\int_{x_{\delta}(\bar{\nu})}^{\infty}dx\,\bar{n}_{\rm
pk}(\bar{\nu},x)\,,$ (122)

where

$\displaystyle\bar{n}_{\rm pk}(\bar{\nu},x)$
$\displaystyle=\frac{e^{-\bar{\nu}^{2}/2}}{(2\pi)^{2}R_{*}^{3}}\,f(x)\frac{e^{-\frac{(x-x_{*})^{2}}{2(1-\gamma^{2})}}}{\sqrt{2\pi(1-\gamma^{2})}}\,,$
(123) $\displaystyle f(x)=\frac{x^{3}-3x}{2}\left[{\rm
Erf}\left(\sqrt{\frac{5}{2}}\,x\right)+{\rm
Erf}\left(\sqrt{\frac{5}{2}}\,\frac{x}{2}\right)\right]+\sqrt{\frac{2}{5\pi}}\left[\left(\frac{31x^{2}}{4}+\frac{8}{5}\right)e^{-\frac{5}{8}x^{2}}+\left(\frac{x^{2}}{2}-\frac{8}{5}\right)e^{-\frac{5}{2}x^{2}}\right]\,,$
(124)

and we remind that $R_{*}\equiv\sqrt{3}\sigma_{1}/\sigma_{2}.$

  * $\circ$

$1+2\alpha\mathcal{R}<0$.
In this case minima of $\mathcal{R}$ are maxima of $h.$ Proceeding analogously
as before, we have:

$\displaystyle\mathcal{N}_{\rm pk}^{({\rm II})}$
$\displaystyle=\int_{-\infty}^{-\frac{1}{2\alpha\sigma_{0}}}d\bar{\nu}\,\int_{-\infty}^{x_{\delta}(\bar{\nu})}dx\,\bar{n}_{\rm
pk}(\bar{\nu},x)\,.$ (125)

The comoving number density of peaks of the overdensity field is approximated
as the sum of the two terms above: $\mathcal{N}_{\rm pk}=\mathcal{N}_{\rm
pk}^{({\rm I})}+\mathcal{N}_{\rm pk}^{({\rm II})}$. One can interpret these
two contributions along the lines exposed in section II. The term
$\mathcal{N}_{\rm pk}^{({\rm II})}$ counts then minima of $\mathcal{R}$ which
are maxima of $h.$ Instead $\mathcal{N}_{\rm pk}^{({\rm I})}$ corresponds to
maxima of $\mathcal{R}$, and the restriction in the integration range of
$\bar{\nu}$ (the lower limit) is designed to subtract those maxima of
$\mathcal{R}$ which are minima of $h$. Obviously, the gaussian result is
recovered for $\alpha\to 0$. In this limit we can compare with the standard
expression in ref. Bardeen:1985tr . Numerically, and for the power spectra
under consideration, we found that the two calculations agree within a factor
$\simeq 2$. As already mentioned in section B for the case of two spatial
dimensions, this confirms that peaks of the overdensity field are well
approximated by peaks of the comoving density perturbation which are spiky
enough. In appendix F we will discuss how to translate $\mathcal{N}_{\rm pk}$
into the primordial black hole abundance in eq. (8).

##  Appendix D Peak theory with local non-gaussianities: a perturbative
approach in 2D

The idea is to derive the joint probability density distribution for the
variables $h$, $h_{x}$, $h_{y}$, $h_{xx}$, $h_{xy}$, $h_{yy}$ (considering
again, for simplicity, the two-dimensional case). Once we get this probability
density distribution, we set $h_{x}=h_{y}=0$ and we integrate over $h_{xx}$,
$h_{xy}$ and $h_{yy}$ in the domain defining maxima. This strategy was simple
to implement in the case of the gaussian variable $\mathcal{R}$ because the
joint probability density distribution was a multivariate normal distribution.
The non-gaussian case is more complicated.

In order to tackle the problem, we shall use the approach based on the
characteristic function. In full generality, for a set of $N$ correlated
random variables $\xi_{i}$ the characteristic function is the Fourier
transform of their joint probability density distribution

$\displaystyle\chi(\lambda_{1},\dots,\lambda_{N})$ $\displaystyle\equiv\int
d\xi_{1}\dots
d\xi_{N}P(\xi_{1},\dots,\xi_{N})\exp\left[i(\xi_{1}\lambda_{1}+\dots+\xi_{N}\lambda_{N})\right]$
$\displaystyle=\int d\xi_{1}\dots
d\xi_{N}P(\xi_{1},\dots,\xi_{N})\left[1+i(\xi_{1}\lambda_{1}+\dots+\xi_{N}\lambda_{N})+\frac{i^{2}}{2!}(\xi_{1}\lambda_{1}+\dots+\xi_{N}\lambda_{N})^{2}+\dots\right]$
$\displaystyle=1+i\sum_{j}\langle\xi_{j}\rangle\lambda_{j}+\frac{i^{2}}{2!}\sum_{j_{1},j_{2}}\langle\xi_{j_{1}}\xi_{j_{2}}\rangle\lambda_{j_{1}}\lambda_{j_{2}}+\frac{i^{3}}{3!}\sum_{j_{1},j_{2},j_{3}}\langle\xi_{j_{1}}\xi_{j_{2}}\xi_{j_{3}}\rangle\lambda_{j_{1}}\lambda_{j_{2}}\lambda_{j_{3}}+\dots$
(126)

where, after a Taylor expansion, we introduced the moments of the joint
distribution by means of the integrals

$\displaystyle\langle\xi_{j_{1}}\dots\xi_{j_{k}}\rangle=\int d\xi_{1}\dots
d\xi_{N}P(\xi_{1},\dots,\xi_{N})\xi_{j_{1}}\dots\xi_{j_{k}}\,.$ (127)

If we take the natural log of the characteristic function and Taylor expand,
we define the cumulants

$\displaystyle\log\chi(\lambda_{1},\dots,\lambda_{N})\equiv
i\sum_{j}C_{1}(\xi_{j})\lambda_{j}+\frac{i^{2}}{2!}\sum_{j_{1},j_{2}}C_{2}(\xi_{j_{1}},\xi_{j_{2}})\lambda_{j_{1}}\lambda_{j_{2}}+\dots\,.$
(128)

The relation between moments and cumulants follows from the comparison of the
two Taylor expansions. We find the well-known relations

$\displaystyle C_{1}(\xi_{j})$ $\displaystyle=\langle\xi_{j}\rangle\,,$ (129)
$\displaystyle C_{2}(\xi_{j_{1}},\xi_{j_{2}})$
$\displaystyle=\langle\xi_{j_{1}}\xi_{j_{2}}\rangle-\langle\xi_{j_{1}}\rangle\langle\xi_{j_{2}}\rangle\,,$
(130) $\displaystyle C_{3}(\xi_{j_{1}},\xi_{j_{2}},\xi_{j_{3}})$
$\displaystyle=\langle\xi_{j_{1}}\xi_{j_{2}}\xi_{j_{3}}\rangle-\langle\xi_{j_{1}}\rangle\langle\xi_{j_{2}}\xi_{j_{3}}\rangle-\langle\xi_{j_{2}}\rangle\langle\xi_{j_{1}}\xi_{j_{3}}\rangle-\langle\xi_{j_{3}}\rangle\langle\xi_{j_{1}}\xi_{j_{2}}\rangle+2\langle\xi_{j_{1}}\rangle\langle\xi_{j_{2}}\rangle\langle\xi_{j_{3}}\rangle\,,$
(131) $\displaystyle C_{4}(\xi_{j_{1}},\xi_{j_{2}},\xi_{j_{3}},\xi_{j_{4}})$
$\displaystyle=\dots$ (132)

and so on. Working with cumulants instead of moments is more efficient. This
is particularly true in the gaussian case. For a set of correlated gaussian
variables, all cumulants of order higher than two vanish. In the gaussian
case, eq. (128) gives

$\displaystyle\chi(\lambda_{1},\dots,\lambda_{N})=\exp\left[i\sum_{j}\langle\xi_{j}\rangle\lambda_{j}-\frac{1}{2}\sum_{j_{1},j_{2}}C_{2}(\xi_{j_{1}},\xi_{j_{2}})\lambda_{j_{1}}\lambda_{j_{2}}\right]\,,$
(133)

and the inverse Fourier transform that gives the probability density
distribution reads

$\displaystyle
P(\xi_{1},\dots,\xi_{N})=\int\frac{d\lambda_{1}}{(2\pi)}\dots\frac{d\lambda_{N}}{(2\pi)}\exp\left[i\sum_{j}\left(\langle\xi_{j}\rangle-\xi_{j}\right)\lambda_{j}-\frac{1}{2}\sum_{j_{1},j_{2}}C_{2}(\xi_{j_{1}},\xi_{j_{2}})\lambda_{j_{1}}\lambda_{j_{2}}\right]\,.$
(134)

If we complete the square inside the integrand and compute the resulting
multivariate gaussian integral, we find precisely eq. (70) where the second-
order cumulants reconstruct the covariance matrix elements in eq. (71). This
was precisely the strategy that we followed in the previous section: we
computed the elements of the covariance matrix (that are the second-order
cumulants) and we (implicitly) performed an inverse Fourier transform to get
back the probability density distribution.

In the non-gaussian case, we can try to apply the same logic. First, we
compute the cumulants; second, we reconstruct the probability density
distribution by means of an inverse Fourier transform.

The computation of the cumulants require some mathematical tricks that we
shall explain in the following. In full generality, we need to compute the
$n^{\rm th}$-order cumulant $C_{n}(\partial_{j}h,\dots,\partial_{l}h)$ where
each $\partial_{j}h$ represents a certain number of spatial derivatives (zero,
one or two) acting on $h$. Remember also that the random field $h$ (and its
derivatives) is computed at a specific spatial position (say, $\vec{x}$) that
will be later identified with a stationary point. We can write

$\displaystyle C_{n}(\partial_{j}h,\dots,\partial_{l}h)$
$\displaystyle=C_{n}[\partial_{j}h(\vec{x}),\dots,\partial_{l}h(\vec{x})]$
$\displaystyle=\left.C_{n}[\partial_{j_{1}}h(\vec{x}_{1}),\dots,\partial_{l_{n}}h(\vec{x}_{n})]\right|_{\vec{x}_{1}=\dots=\vec{x}_{n}=\vec{x}}=\partial_{j_{1}}\dots\partial_{l_{n}}\left.C_{n}[h(\vec{x}_{1}),\dots,h(\vec{x}_{n})]\right|_{\vec{x}_{1}=\dots=\vec{x}_{n}=\vec{x}}\,.$
(135)

In the first step of eq. (135) we consider each $\partial_{j_{k}}h$ to act at
a different point $\vec{x}_{k}$, and later we set all points equal again (we
already used a similar trick in eq. (77)). The advantage of this step is that
since now each derivative acts at a different spatial point, we can bring them
outside the cumulant (as done in the last step of eq. (135)). It is a very
simple exercise to check explicitly (for instance, by computing second-order
cumulants of a gaussian variable with its derivatives) that this procedure is
completely legitimate. Using our notation $h(\vec{x}_{n})=h_{n}$, the problem
reduces to the computation of $C_{n}(h_{1},\dots,h_{n})$ where now the random
fields are evaluated at different spatial points. After computing
$C_{n}(h_{1},\dots,h_{n})$, we will take the spatial derivatives and finally
set all points equal according to the prescription in eq. (135).

To compute $C_{n}(h_{1},\dots,h_{n})$, we consider the relation
$h=\mathcal{R}+\alpha\mathcal{R}^{2}$. We write in general
$h=\mathcal{R}+f_{\rm NL}(\mathcal{R})$ where $f_{\rm NL}(\mathcal{R})$ can be
a generic non-linear function controlled by the parameter $\alpha$. Let us
consider the term $\alpha\mathcal{R}^{2}$ as a perturbation and expand in
$\alpha$. We have

$\displaystyle
C_{n}(h_{1},\dots,h_{n})=C_{n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$
$\displaystyle+C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]+\dots+C_{n}[\mathcal{R}_{1},\dots,f_{\rm
NL}(\mathcal{R}_{n})]$ (136) $\displaystyle+\left\\{C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2}),\dots,\mathcal{R}_{n}]+\dots\right\\}+O(\alpha^{3})\,.$
(137)

This deconstruction makes clear the fact that we can compute the cumulants
$C_{n}(h_{1},\dots,h_{n})$ based on the cumulants computed for the gaussian
random field $\mathcal{R}$.

The problem is now the computation of $C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]$ where we remind that
$\mathcal{R}_{n}=\mathcal{R}(\vec{x}_{n})$.

In order to compute $C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]$, we need to work out an
intermediate result. First, remember the definition of random field that we
gave at the beginning of appendix B: The scalar random field
$\mathcal{R}(\vec{x})$ is a set of random variables, one for each point
$\vec{x}$ in space, equipped with a joint probability density distribution
$p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$. When the random field is gaussian,
$p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$ is a multivariate Gaussian
distribution and we can write

$\displaystyle
p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})=\frac{1}{(2\pi)^{n/2}\sqrt{{\rm
det}\sigma}}\exp\left[-\frac{1}{2}\sum_{i,j}(\sigma^{-1})_{ij}\mathcal{R}_{i}\mathcal{R}_{j}\right]\,,$
(138)

where $\sigma_{ij}=\langle\mathcal{R}_{i}\mathcal{R}_{j}\rangle$ (as before,
we consider zero-mean random field).

Consider now a new set of random variables
$\\{\mathcal{Y}_{1},\dots,\mathcal{Y}_{n}\\}$ related to the previous one via
the transformations
$\mathcal{Y}_{1}=g_{1}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$, $\dots$,
$\mathcal{Y}_{n}=g_{n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$ with inverse
$\mathcal{R}_{1}=g_{1}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$, $\dots$,
$\mathcal{R}_{n}=g_{n}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$. The joint
probability density distribution for the transformed variables is given by

$\displaystyle
p_{\mathcal{Y}}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})=p[g_{1}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n}),\dots,g_{n}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})]\left|{\rm
det}J(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\right|\,,$ (139)

where the Jacobian matrix is

$\displaystyle
J(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\equiv\left(\begin{array}[]{ccc}\frac{\partial\mathcal{R}_{1}}{\partial\mathcal{Y}_{1}}&\dots&\frac{\partial\mathcal{R}_{1}}{\partial\mathcal{Y}_{n}}\\\
\vdots&\ddots&\vdots\\\
\frac{\partial\mathcal{R}_{n}}{\partial\mathcal{Y}_{1}}&\dots&\frac{\partial\mathcal{R}_{n}}{\partial\mathcal{Y}_{n}}\end{array}\right)\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}~{}~{}\frac{\partial\mathcal{R}_{j}}{\partial\mathcal{Y}_{k}}=\frac{\partial}{\partial\mathcal{Y}_{k}}\left[g_{j}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\right]\,.$
(143)

Eq. (139) is known as Jabobi’s multivariate theorem.

We note that in eq. (139) we can use the relation $\left|{\rm
det}J(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\right|=1/\left|{\rm
det}J(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\right|$ where
$J(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$ is the Jacobian matrix of the
inverse transformation with elements
$\partial\mathcal{Y}_{j}/\partial\mathcal{R}_{k}=\partial
g_{j}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})/\partial\mathcal{R}_{k}$. This
implies that we can rewrite

$\displaystyle p_{\mathcal{Y}}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})=\int
d\mathcal{R}_{1}\dots
d\mathcal{R}_{n}p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\delta[\mathcal{Y}_{1}-g_{1}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})]\dots\delta[\mathcal{Y}_{n}-g_{n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})]\,,$
(144)

and eq. (139) follows from eq. (144) if one applies the transformation
property of the multi-dimensional delta function

$\displaystyle\delta[\mathcal{Y}_{1}-g_{1}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})]\dots\delta[\mathcal{Y}_{n}-g_{n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})]=\frac{\delta[\mathcal{R}_{1}-g_{1}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})]\dots\delta[\mathcal{R}_{n}-g_{n}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})]}{\left|{\rm
det}J(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\right|}\,,$ (145)

and integrates over $\mathcal{R}_{i}$. Eq. (144) is a very useful formula.
Consider the characteristic function of
$p_{\mathcal{Y}}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$ defined by the
Fourier transform

$\displaystyle\chi_{\mathcal{Y}}(\lambda_{1},\dots,\lambda_{n})=\int
d\mathcal{Y}_{1}\dots
d\mathcal{Y}_{n}p_{\mathcal{Y}}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\exp\left[i(\mathcal{Y}_{1}\lambda_{1}+\dots+\mathcal{Y}_{n}\lambda_{n})\right]\,.$
(146)

Using eq. (144) and integrating over $\mathcal{Y}_{i}$ thanks to the delta
functions, we find

$\displaystyle\chi_{\mathcal{Y}}(\lambda_{1},\dots,\lambda_{n})=\int
d\mathcal{R}_{1}\dots
d\mathcal{R}_{n}e^{i[g_{1}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\lambda_{1}+\dots+g_{n}(\mathcal{R}_{1},\,\dots\,,\mathcal{R}_{n})\lambda_{n}]}\,p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\,.$
(147)

From the definition of cumulants with respect to the characteristic function
$\chi_{\mathcal{Y}}(\lambda_{1},\dots,\lambda_{n})$, we finally find (see eq.
(128))

$\displaystyle C_{n}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$
$\displaystyle=\left.(-i)^{n}\frac{\partial}{\partial\lambda_{1}}\dots\frac{\partial}{\partial\lambda_{n}}\log\chi_{\mathcal{Y}}(\lambda_{1},\dots,\lambda_{n})\right|_{\lambda_{1}=\,\dots\,=\lambda_{n}=0}$
$\displaystyle=\left.(-i)^{n}\frac{\partial}{\partial\lambda_{1}}\dots\frac{\partial}{\partial\lambda_{n}}\log\int
d\mathcal{R}_{1}\dots
d\mathcal{R}_{n}e^{i[g_{1}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\lambda_{1}+\dots+g_{n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\lambda_{n}]}\,p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\right|_{\lambda_{1}=\,\dots\,=\lambda_{n}=0}\,.$
(148)

This is a generic formula that can be applied to our case.

Consider the linear order in $\alpha$. From eq. (148), the cumulant
$C_{n}[f_{\rm NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]$ is given by

$\displaystyle C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]=\left.(-i)^{n}\frac{\partial}{\partial\lambda_{1}}\dots\frac{\partial}{\partial\lambda_{n}}\log\int
d\mathcal{R}_{1}\dots d\mathcal{R}_{n}e^{i[f_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}+\mathcal{R}_{2}\lambda_{2}+\dots+\mathcal{R}_{n}\lambda_{n}]}\,p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\right|_{\lambda_{1}=\,\dots\,=\lambda_{n}=0}\,.$
(149)

A similar formula is applicable to all remaining first-order terms in eq.
(136). Similarly, at order $O(\alpha^{2})$ one needs to compute, for instance,
cumulants like

$\displaystyle C_{n}[f_{\rm NL}$ $\displaystyle(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2}),\dots,\mathcal{R}_{n}]=$
$\displaystyle\left.(-i)^{n}\frac{\partial}{\partial\lambda_{1}}\dots\frac{\partial}{\partial\lambda_{n}}\log\int
d\mathcal{R}_{1}\dots d\mathcal{R}_{n}e^{i[f_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}+f_{\rm
NL}(\mathcal{R}_{2})\lambda_{2}+\dots+\mathcal{R}_{n}\lambda_{n}]}\,p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\right|_{\lambda_{1}=\,\dots\,=\lambda_{n}=0}\,.$
(150)

All we need to do is computing the integral and taking derivatives with
respect to $\lambda_{i}$. The computation of the integrals is simplified by
the fact that the variables $\mathcal{R}_{i}$ are gaussian with joint
distribution given by eq. (138).

Before proceeding, an important remark is in order. For the sake of
simplicity, we derived eq. (139) assuming that the functions
$\mathcal{Y}_{i=1,\dots,n}=g_{i=1,\dots,n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$
define one-to-one mappings. In this case, there exist unique inverse functions
so that
$\mathcal{R}_{i=1,\dots,n}=g_{i=1,\dots,n}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$.
The case we are interested in, however, is not exactly of this form. If we
solve $h=\mathcal{R}+\alpha\mathcal{R}^{2}$ for $\mathcal{R}$, we indeed find
$\mathcal{R}=(-1\pm\sqrt{1+4\alpha h})/2\alpha$ which has two roots.

However, this is not an insurmountable problem since Jabobi’s multivariate
theorem in eq. (139) can be generalized to the case in which the system
$\mathcal{Y}_{i=1,\dots,n}=g_{i=1,\dots,n}(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$
admits at most a countable number of roots. Let us indicate these
$q=1,\dots,Q$ roots in the form
$\mathcal{R}_{q,i=1,\dots,n}=g_{q,i=1,\dots,n}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$.
Eq. (139) becomes

$\displaystyle
p_{\mathcal{Y}}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})=\sum_{q=1}^{Q}p[g_{q,1}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n}),\dots,g_{q,n}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})]\left|{\rm
det}J_{q}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\right|\,,$ (151)

where now $J_{q}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})$ is the Jacobian
corresponding to the $q^{\rm th}$ root

$\displaystyle
J_{q}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\equiv\left(\begin{array}[]{ccc}\frac{\partial\mathcal{R}_{q,1}}{\partial\mathcal{Y}_{1}}&\dots&\frac{\partial\mathcal{R}_{q,1}}{\partial\mathcal{Y}_{n}}\\\
\vdots&\ddots&\vdots\\\
\frac{\partial\mathcal{R}_{q,n}}{\partial\mathcal{Y}_{1}}&\dots&\frac{\partial\mathcal{R}_{q,n}}{\partial\mathcal{Y}_{n}}\end{array}\right)\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}~{}~{}\frac{\partial\mathcal{R}_{q,j}}{\partial\mathcal{Y}_{k}}=\frac{\partial}{\partial\mathcal{Y}_{k}}\left[g_{q,j}^{-1}(\mathcal{Y}_{1},\dots,\mathcal{Y}_{n})\right]\,.$
(155)

The multi-dimensional delta-function identity in eq. (145) changes
accordingly, and the final result in eq. (148) remains unaltered.

After this digression, we are ready to compute the integral in eq. (148).
Instead of looking for a generic expression, let us consider specific cases
organized for increasing level of difficulty.

  * $\circ$

The simplest possibility is to truncate the analysis at the linear order in
$\alpha$, eq. (149). Consider the integral in eq. (149) which we rewrite as

$\displaystyle\mathcal{I}(\lambda_{1},\dots,\lambda_{n})\equiv\int
d\mathcal{R}_{1}e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}}\bigg{[}\underbrace{\int d\mathcal{R}_{2}\dots
d\mathcal{R}_{n}e^{i[\mathcal{R}_{2}\lambda_{2}+\,\dots\,+\mathcal{R}_{n}\lambda_{n}]}\,p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})}_{\equiv\,\mathcal{I}^{\prime}(\mathcal{R}_{1},\lambda_{2},\dots,\lambda_{n})}\bigg{]}\,.$
(156)

The key observation is that the integral inside the square brackets would
precisely match the definition of the characteristic function
$\chi(\lambda_{1},\dots,\lambda_{n})$ of
$p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})$ if it were completed by an
additional integration over $\mathcal{R}_{1}$ (see the definition in the first
line of eq. (D)). In such case, we could simply use (see eq. (133))

$\displaystyle\int
d\mathcal{R}_{1}e^{i\lambda_{1}\mathcal{R}_{1}}\mathcal{I}^{\prime}(\mathcal{R}_{1},\lambda_{2},\dots,\lambda_{n})=\chi(\lambda_{1},\dots,\lambda_{n})=\exp\bigg{[}-\frac{1}{2}\sum_{ij}\sigma_{ij}\lambda_{i}\lambda_{j}\bigg{]}\,,$
(157)

since $\\{\mathcal{R}_{1},\dots,\mathcal{R}_{n}\\}$ are gaussian. However, the
fact that integration over $\mathcal{R}_{1}$ is missing implies that the
integral $\mathcal{I}^{\prime}(\mathcal{R}_{1},\lambda_{2},\dots,\lambda_{n})$
is actually equal to the inverse Fourier transform of
$\chi(\lambda_{1},\dots,\lambda_{n})$ with respect to $\lambda_{1}$. In
formulas, eq. (156) becomes

$\displaystyle\mathcal{I}(\lambda_{1},\dots,\lambda_{n})$ $\displaystyle=\int
d\mathcal{R}_{1}e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}}\int\frac{d\lambda_{1}}{(2\pi)}e^{-i\lambda_{1}\mathcal{R}_{1}}\exp\bigg{[}-\frac{1}{2}\sum_{ij}\sigma_{ij}\lambda_{i}\lambda_{j}\bigg{]}$
$\displaystyle=\int d\mathcal{R}_{1}e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}}\exp\bigg{[}-\frac{1}{2}\sum_{ij\geqslant
2}\sigma_{ij}\lambda_{i}\lambda_{j}\bigg{]}\int\frac{d\lambda_{1}}{(2\pi)}\exp\bigg{[}-\frac{1}{2}\sigma_{11}\lambda_{1}^{2}-\bigg{(}i\mathcal{R}_{1}+\sum_{j\geqslant
2}\sigma_{1j}\lambda_{j}\bigg{)}\lambda_{1}\bigg{]}$ $\displaystyle=\int
d\mathcal{R}_{1}e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}}\frac{1}{\sqrt{2\pi\sigma_{11}}}\exp\bigg{[}-\frac{1}{2}\sum_{ij\geqslant
2}\sigma_{ij}\lambda_{i}\lambda_{j}+\frac{1}{2\sigma_{11}}\bigg{(}i\mathcal{R}_{1}+\sum_{j\geqslant
2}\sigma_{1j}\lambda_{j}\bigg{)}^{2}\bigg{]}\,,$ (158)

where we just reorganized terms before integrating over $\lambda_{1}$ in the
last step. If we take the natural log, we find

$\displaystyle\log[\mathcal{I}(\lambda_{1},\dots,\lambda_{n})]=-\frac{1}{2}\sum_{ij\geqslant
2}\sigma_{ij}\lambda_{i}\lambda_{j}+\log\int d\mathcal{R}_{1}e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}}\frac{1}{\sqrt{2\pi\sigma_{11}}}\exp\bigg{[}\frac{1}{2\sigma_{11}}\bigg{(}i\mathcal{R}_{1}+\sum_{j\geqslant
2}\sigma_{1j}\lambda_{j}\bigg{)}^{2}\bigg{]}\,.$ (159)

According to eq. (149), we now need to compute derivatives with respect to
$\lambda_{i=1,\,\dots\,,n}$ and finally set $\lambda_{i=1,\,\dots\,,n}=0$. We
find

$\displaystyle C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]=\frac{1}{\sqrt{2\pi\sigma_{11}}}\bigg{(}\prod_{k=2}^{n}\sigma_{1k}\bigg{)}\int
d\mathcal{R}_{1}f_{\rm
NL}(\mathcal{R}_{1})g_{n}(\mathcal{R}_{1}^{n-1})e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}\,,$
(160)

where $g_{n}(\mathcal{R}_{1}^{n-1})$ is a polynomial of order
$\mathcal{R}_{1}^{n-1}$ whose explicit expression is given by242424We can
notice that the $g_{n}$ polynomials are nothing but the Hermite polynomials,
except for a pre-factor. In particular:
$g_{n}(\mathcal{R}_{1}^{n-1})=\frac{1}{(2\sigma_{11})^{n/2}}H_{n-1}\left(\frac{\mathcal{R}_{1}}{\sqrt{2\sigma_{11}}}\right)$
(161)

$\displaystyle
g_{n}(\mathcal{R}_{1}^{n-1})=(-1)^{n-1}e^{\mathcal{R}_{1}^{2}/2\sigma_{11}}\frac{d^{n-1}}{d\mathcal{R}_{1}^{n-1}}e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}\,.$
(162)

This result implies that we can just integrate by parts $n-1$ times in eq.
(160). We find

$\displaystyle C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]=\bigg{(}\prod_{k=2}^{n}\sigma_{1k}\bigg{)}\underbrace{\frac{1}{\sqrt{2\pi\sigma_{11}}}\int
d\mathcal{R}_{1}f_{\rm
NL}^{(n-1)}(\mathcal{R}_{1})e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}}_{=\,\langle
f_{\rm NL}^{(n-1)}(\mathcal{R}_{1})\rangle}\,.$ (163)

The leftover integral is nothing but the expectation value of $f_{\rm
NL}^{(n-1)}(\mathcal{R}_{1})$. All in all, we find (remember that
$\sigma_{ij}=\langle\mathcal{R}_{i}\mathcal{R}_{j}\rangle$)

$\displaystyle C_{n}[f_{\rm
NL}(\mathcal{R}_{1}),\dots,\mathcal{R}_{n}]=\langle f_{\rm
NL}^{(n-1)}(\mathcal{R}_{1})\rangle\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle\dots\langle\mathcal{R}_{1}\mathcal{R}_{n}\rangle$
(164)

  * $\circ$

At order $O(\alpha^{2})$, consider eq. (150) and the integral

$\displaystyle\mathcal{J}(\lambda_{1},\dots,\lambda_{n})\equiv\int
d\mathcal{R}_{1}d\mathcal{R}_{2}e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}+if_{\rm
NL}(\mathcal{R}_{2})\lambda_{2}}\bigg{[}\int d\mathcal{R}_{3}\dots
d\mathcal{R}_{n}e^{i[\mathcal{R}_{3}\lambda_{3}+\,\dots\,+\mathcal{R}_{n}\lambda_{n}]}\,p(\mathcal{R}_{1},\dots,\mathcal{R}_{n})\bigg{]}\,.$
(165)

The computation goes as before with the only difference that now we need to
consider the inverse Fourier transform of the characteristic function with
respect to both $\lambda_{1}$ and $\lambda_{2}$. For the natural log of
$\mathcal{J}(\lambda_{1},\dots,\lambda_{n})$ we find

$\displaystyle\log[\mathcal{J}$
$\displaystyle(\lambda_{1},\dots,\lambda_{n})]=-\frac{1}{2}\sum_{ij\geqslant
3}\sigma_{ij}\lambda_{i}\lambda_{j}$ (166) $\displaystyle+\log\int
d\mathcal{R}_{1}d\mathcal{R}_{2}\frac{e^{if_{\rm
NL}(\mathcal{R}_{1})\lambda_{1}+if_{\rm
NL}(\mathcal{R}_{2})\lambda_{2}}}{2\pi\sqrt{\sigma_{11}\sigma_{22}-\sigma_{12}^{2}}}\exp\bigg{[}\frac{1}{2(\sigma_{11}\sigma_{22}-\sigma_{12}^{2})}\bigg{(}A_{1}^{2}\sigma_{22}+A_{2}^{2}\sigma_{11}-2A_{1}A_{2}\sigma_{12}\bigg{)}\bigg{]}\,,$

where $A_{1}\equiv i\mathcal{R}_{1}+\sum_{j\geqslant 3}\sigma_{1j}\lambda_{j}$
and $A_{2}\equiv i\mathcal{R}_{2}+\sum_{j\geqslant 3}\sigma_{2j}\lambda_{j}$.
This equation is the analogue of eq. (159) at order $O(\alpha^{2})$. As
before, we need to to compute derivatives with respect to
$\lambda_{i=1,\,\dots\,,n}$ and finally set $\lambda_{i=1,\,\dots\,,n}=0$. Let
us consider first the derivatives with respect to $\lambda_{1}$ and
$\lambda_{2}$. We find

$\displaystyle(-i)^{2}$
$\displaystyle\frac{\partial}{\partial\lambda_{2}}\frac{\partial}{\partial\lambda_{1}}\log[\mathcal{J}(\lambda_{1},\dots,\lambda_{n})]\bigg{|}_{\lambda_{1}=\lambda_{2}=0}=$
(167)
$\displaystyle+\underbrace{\frac{1}{2\pi\sqrt{\sigma_{11}\sigma_{22}-\sigma_{12}^{2}}}}_{2\pi\sqrt{{\rm
det}\,\tilde{\sigma}}}\int d\mathcal{R}_{1}d\mathcal{R}_{2}f_{\rm
NL}(\mathcal{R}_{1})f_{\rm
NL}(\mathcal{R}_{2})\exp\bigg{[}\underbrace{\frac{1}{2(\sigma_{11}\sigma_{22}-\sigma_{12}^{2})}\bigg{(}A_{1}^{2}\sigma_{22}+A_{2}^{2}\sigma_{11}-2A_{1}A_{2}\sigma_{12}\bigg{)}}_{=\,(1/2)\sum_{i,j=1}^{2}(\tilde{\sigma}^{-1})_{ij}A_{i}A_{j}}\bigg{]}$
$\displaystyle-\bigg{[}\frac{1}{\sqrt{2\pi\sigma_{11}}}\int
d\mathcal{R}_{1}f_{\rm
NL}(\mathcal{R}_{1})e^{A_{1}^{2}/2\sigma_{11}}\bigg{]}\bigg{[}\frac{1}{\sqrt{2\pi\sigma_{22}}}\int
d\mathcal{R}_{2}f_{\rm
NL}(\mathcal{R}_{2})e^{A_{2}^{2}/2\sigma_{22}}\bigg{]}\,,$

where $\tilde{\sigma}$ is the two-by-two sub-matrix of $\sigma$ formed by the
first two rows and columns. Before proceeding, let us consider a simple
example. If we are interested to the computation of the second-order cumulant,
these two derivatives are the only ones that we need to take. Furthermore, in
this case we have $A_{j=1,2}=i\mathcal{R}_{j=1,2}$. We find the simple result

$\displaystyle C_{2}[f_{\rm NL}(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2})]=\langle f_{\rm NL}(\mathcal{R}_{1})f_{\rm
NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm NL}(\mathcal{R}_{1})\rangle\langle
f_{\rm NL}(\mathcal{R}_{2})\rangle$ (168)

where the expectation values $\langle f_{\rm NL}(\mathcal{R}_{1})\rangle$ and
$\langle f_{\rm NL}(\mathcal{R}_{2})\rangle$ are defined as in eq. (163) while
$\langle f_{\rm NL}(\mathcal{R}_{1})f_{\rm NL}(\mathcal{R}_{2})\rangle$ is
computed by means of the joint probability distribution of $\mathcal{R}_{1}$
and $\mathcal{R}_{2}$; in formulas, we have

$\displaystyle\langle f_{\rm NL}(\mathcal{R}_{1})f_{\rm
NL}(\mathcal{R}_{2})\rangle=\frac{1}{2\pi\sqrt{{\rm det}\,\tilde{\sigma}}}\int
d\mathcal{R}_{1}d\mathcal{R}_{2}f_{\rm NL}(\mathcal{R}_{1})f_{\rm
NL}(\mathcal{R}_{2})\exp\bigg{[}-\frac{1}{2}\sum_{i,j=1}^{2}(\tilde{\sigma}^{-1})_{ij}\mathcal{R}_{i}\mathcal{R}_{j}\bigg{]}\,.$
(169)

For $n\geqslant 3$, we need to compute in eq. (167) derivatives with respect
to $\lambda_{i=3,\,\dots\,,n}$ and finally set $\lambda_{i=3,\,\dots\,,n}=0$.
These derivatives act on the exponential functions in the integrands of the
two terms on the right-hand side of eq. (167). Let us start from the second
term; if we compute the $\lambda$-derivatives on the exponential function, we
find the following result

$\displaystyle(-i$
$\displaystyle)^{n-2}\frac{\partial}{\partial\lambda_{n}}\dots\frac{\partial}{\partial\lambda_{3}}\bigg{(}e^{A_{1}^{2}/2\sigma_{11}}e^{A_{2}^{2}/2\sigma_{22}}\bigg{)}\bigg{|}_{\lambda_{3}=\,\dots\,\lambda_{n}=0}=$
(170)
$\displaystyle(-1)^{n-2}\bigg{(}\prod_{k=3}^{n}\sigma_{1k}\bigg{)}\bigg{(}\frac{d^{n-2}}{d\mathcal{R}_{1}^{n-2}}e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}\bigg{)}e^{-\mathcal{R}_{2}^{2}/2\sigma_{22}}+(-1)^{n-2}\bigg{(}\prod_{k=3}^{n}\sigma_{2k}\bigg{)}\bigg{(}\frac{d^{n-2}}{d\mathcal{R}_{2}^{n-2}}e^{-\mathcal{R}_{2}^{2}/2\sigma_{22}}\bigg{)}e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}$
$\displaystyle+(-1)^{n-2}\bigg{[}\sum_{a,b>0}^{a+b=n-2}\bigg{(}\prod_{k=3+b}^{n}\sigma_{1k}\bigg{)}\bigg{(}\prod_{j=3}^{n-a}\sigma_{2j}\bigg{)}\bigg{(}\frac{d^{a}}{d\mathcal{R}_{1}^{a}}e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}\bigg{)}\bigg{(}\frac{d^{b}}{d\mathcal{R}_{2}^{b}}e^{-\mathcal{R}_{2}^{2}/2\sigma_{22}}\bigg{)}+{\rm
perms\,of}\,(3,\dots,n)\bigg{]}\,.$

To fully understand the content of this equation, few comments are in order.
Inside the square brackets, the sum $\sum_{a,b}$ runs over all possible
combinations of $a>0$ and $b>0$ such that $a+b=n-2$. For instance, if we take
$n=4$ we have only $\\{a=1,b=1\\}$ while for $n=5$ we have two combinations,
namely $\\{a=1,b=2\\}$ and $\\{a=2,b=1\\}$ (notice that if $n=3$, there are no
combinations that are allowed, and the term in the square brackets does not
contribute to the third-order cumulant). Furthermore, after computing the sum
over $a$ and $b$, a further sum over all permutations of the indices
$(3,\dots,n)$ that give distinct results is required. Let us clarify this
point with an explicit example. Consider $n=4$. The sum over $\\{a=1,b=1\\}$
gives
$(e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}e^{-\mathcal{R}_{2}^{2}/2\sigma_{22}}\mathcal{R}_{1}\mathcal{R}_{2}/\sigma_{11}\sigma_{22})\sigma_{14}\sigma_{23}$;
we now need include the term with $(3,4)$ exchanged so that the final result
is
$(e^{-\mathcal{R}_{1}^{2}/2\sigma_{11}}e^{-\mathcal{R}_{2}^{2}/2\sigma_{22}}\mathcal{R}_{1}\mathcal{R}_{2}/\sigma_{11}\sigma_{22})(\sigma_{14}\sigma_{23}+\sigma_{13}\sigma_{24})$.
Only permutations that give distinct results need to be included. For
instance, take for $n=5$ the term proportional to
$\sigma_{14}\sigma_{15}\sigma_{23}$; in this case, the sum over permutations
of $(3,4,5)$ gives only three distinct terms,
$\sigma_{14}\sigma_{15}\sigma_{23}\to\sigma_{14}\sigma_{15}\sigma_{23}+\sigma_{13}\sigma_{15}\sigma_{24}+\sigma_{14}\sigma_{13}\sigma_{25}$
(instead of six—that is the number of permutations of three objects—since the
remaining three do not give distinct results).

As far as the first term on the right-hand side of eq. (167) is concerned, we
find

$\displaystyle(-i)^{n-2}\frac{\partial}{\partial\lambda_{n}}\dots\frac{\partial}{\partial\lambda_{3}}\bigg{[}e^{\frac{A_{1}^{2}\sigma_{22}+A_{2}^{2}\sigma_{11}-2A_{1}A_{2}\sigma_{12}}{2(\sigma_{11}\sigma_{22}-\sigma_{12}^{2})}}\bigg{]}\bigg{|}_{\lambda_{3}=\,\dots\,\lambda_{n}=0}=$
(171)
$\displaystyle(-1)^{n-2}\bigg{(}\prod_{k=3}^{n}\sigma_{1k}\bigg{)}\bigg{[}\frac{d^{n-2}}{d\mathcal{R}_{1}^{n-2}}e^{-\frac{\mathcal{R}_{1}^{2}\sigma_{22}+\mathcal{R}_{2}^{2}\sigma_{11}-2\mathcal{R}_{1}\mathcal{R}_{2}\sigma_{12}}{2(\sigma_{11}\sigma_{22}-\sigma_{12}^{2})}}\bigg{]}+(-1)^{n-2}\bigg{(}\prod_{k=3}^{n}\sigma_{2k}\bigg{)}\bigg{[}\frac{d^{n-2}}{d\mathcal{R}_{2}^{n-2}}e^{-\frac{\mathcal{R}_{1}^{2}\sigma_{22}+\mathcal{R}_{2}^{2}\sigma_{11}-2\mathcal{R}_{1}\mathcal{R}_{2}\sigma_{12}}{2(\sigma_{11}\sigma_{22}-\sigma_{12}^{2})}}\bigg{]}$
$\displaystyle+(-1)^{n-2}\bigg{\\{}\sum_{a,b>0}^{a+b=n-2}\bigg{(}\prod_{k=3+b}^{n}\sigma_{1k}\bigg{)}\bigg{(}\prod_{j=3}^{n-a}\sigma_{2j}\bigg{)}\bigg{[}\frac{\partial^{n-2}}{\partial\mathcal{R}_{1}^{a}\partial\mathcal{R}_{2}^{b}}e^{-\frac{\mathcal{R}_{1}^{2}\sigma_{22}+\mathcal{R}_{2}^{2}\sigma_{11}-2\mathcal{R}_{1}\mathcal{R}_{2}\sigma_{12}}{2(\sigma_{11}\sigma_{22}-\sigma_{12}^{2})}}\bigg{]}+{\rm
perms\,of}\,(3,\dots,n)\bigg{\\}}\,,$

where the sum $\sum_{a,b}$ and the subsequent one over different permutations
have the same meaning explained before.

Eqs. (170, 171) allow to compute the integrals in eq. (167) by means of $n-2$
integration by parts. We find the final result (for $n\geqslant 3$)

$\displaystyle C_{n}[f_{\rm NL}(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2}),\dots,\mathcal{R}_{n}]=$ (172)
$\displaystyle~{}~{}\left[\langle f_{\rm NL}^{(n-2)}(\mathcal{R}_{1})f_{\rm
NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm
NL}^{(n-2)}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}(\mathcal{R}_{2})\rangle\right]\langle\mathcal{R}_{1}\mathcal{R}_{3}\rangle\dots\langle\mathcal{R}_{1}\mathcal{R}_{n}\rangle~{}~{}+$
$\displaystyle~{}~{}\left[\langle f_{\rm NL}(\mathcal{R}_{1})f^{(n-2)}_{\rm
NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm NL}(\mathcal{R}_{1})\rangle\langle
f_{\rm
NL}^{(n-2)}(\mathcal{R}_{2})\rangle\right]\langle\mathcal{R}_{2}\mathcal{R}_{3}\rangle\dots\langle\mathcal{R}_{2}\mathcal{R}_{n}\rangle~{}~{}+$
$\displaystyle~{}~{}\bigg{\\{}\sum_{a,b>0}^{a+b=n-2}\left[\langle f^{(a)}_{\rm
NL}(\mathcal{R}_{1})f^{(b)}_{\rm NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm
NL}^{(a)}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}^{(b)}(\mathcal{R}_{2})\rangle\right]\prod_{k=3+b}^{n}\langle\mathcal{R}_{1}\mathcal{R}_{k}\rangle\prod_{j=3}^{n-a}\langle\mathcal{R}_{2}\mathcal{R}_{j}\rangle+{\rm
perms\,of}\,(3,\dots,n)\bigg{\\}}$

which generalizes eq. (164) at order $O(\alpha^{2})$. We remark that eqs.
(164, 172) are completely generic and do not depend on the functional form of
$f_{\rm NL}(\mathcal{R})$.

  * $\circ$

One can in principle continue the computation and include higher-order
$\alpha$-corrections to the cumulants. Needless to say, the corresponding
expressions quickly become quite unwieldy (even though the computation remains
conceptually simple). One possible strategy is to include only $O(\alpha)$
corrections, and check that $O(\alpha^{2})$ terms remain sub-leading. In order
for this criterium to be satisfactory, one should correctly identify the
expansion parameter. As we shall see, in the case with no derivatives involved
in eq. (135) the latter turns out to be $\alpha^{2}\sigma_{0}^{2}$.

  * $\circ$

We can conclude this part with few checks of eqs. (164, 172). For simplicity,
let us consider the case in which there are no derivatives in eq. (135).

$n=2$

We compute $C_{2}(h_{1},h_{2})$. At order $O(\alpha)$, we have

$\displaystyle C_{2}(h_{1},h_{2})$
$\displaystyle=C_{2}(\mathcal{R}_{1},\mathcal{R}_{2})+C_{2}[f_{\rm
NL}(\mathcal{R}_{1}),\mathcal{R}_{2}]+C_{2}[\mathcal{R}_{1},f_{\rm
NL}(\mathcal{R}_{2})]$ (173)
$\displaystyle=\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle+\langle f_{\rm
NL}^{(1)}(\mathcal{R}_{1})\rangle\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle+\langle
f_{\rm
NL}^{(1)}(\mathcal{R}_{2})\rangle\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle\,,$
(174)

where we use eq. (130) to rewrite the first term (which is the gaussian one)
and eq. (164) to rewrite the last two ones. Since we have $f_{\rm
NL}(\mathcal{R})=\alpha\mathcal{R}^{2}$, we have $\langle f_{\rm
NL}^{(1)}(\mathcal{R}_{i=1,2})\rangle=0$ since we are considering zero-mean
random fields. This means that there are no corrections at order $O(\alpha)$.
We now move to consider corrections at order $O(\alpha^{2})$. We use eq.
(168). At order $O(\alpha^{2})$, we have

$\displaystyle
C_{2}(h_{1},h_{2})=\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle+\langle f_{\rm
NL}(\mathcal{R}_{1})f_{\rm NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm
NL}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}(\mathcal{R}_{2})\rangle=\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle+\alpha^{2}\langle\mathcal{R}_{1}^{2}\mathcal{R}_{2}^{2}\rangle-\alpha^{2}\langle\mathcal{R}_{1}^{2}\rangle\langle\mathcal{R}_{2}^{2}\rangle\,.$
(175)

We now set the two points equal (since we are not considering derivatives of
cumulants, there is no need of considering different spatial point). We have
$\langle\mathcal{R}^{2}\rangle=\sigma_{0}^{2}$ and
$\langle\mathcal{R}^{4}\rangle=3\sigma_{0}^{4}$. In conclusion, we find

$\displaystyle
C_{2}(h,h)=\sigma_{0}^{2}\left(1+2\alpha^{2}\sigma_{0}^{2}\right)\,.$ (176)

This is an exact result since corrections of order $O(\alpha^{3})$ enter only
at higher-orders. As a rule of thumb, for a cumulant $C_{n}$ of order $n$,
only corrections up to order $O(\alpha^{n})$ are possible. Eq. (176) shows
that the correct expansion parameter is $\alpha^{2}\sigma_{0}^{2}$. In
realistic models of inflation we expect $\alpha=O(1)$ while we have
$\sigma_{0}^{2}=A_{g}$ for the specific power spectrum in eq. (89). Since one
typically has $A_{g}=O(10^{-3})$, we expect $\alpha^{2}\sigma_{0}^{2}\ll 1$.
Let us also notice that in applications of cosmological interest it is
customary to use $f_{\rm
NL}(\mathcal{R})=\alpha(\mathcal{R}^{2}-\langle\mathcal{R}^{2}\rangle)$ since
in this way one has a zero-mean non-gaussian field $h$, $\langle h\rangle=0$.
It is simple to see that the constant shift in $f_{\rm NL}(\mathcal{R})$ that
is implied by this particular choice does not affect the computation of
cumulants of order equal or higher than two. However, it does change the first
order cumulant since $C_{1}(h)=\langle h\rangle$ (see eq. (129)). The choice
$f_{\rm
NL}(\mathcal{R})=\alpha(\mathcal{R}^{2}-\langle\mathcal{R}^{2}\rangle)$,
therefore, leads to $C_{1}(h)=0$.

$n=3$

Consider now $C_{3}(h_{1},h_{2},h_{3})$. This cumulant vanishes in the
gaussian limit. At order $O(\alpha)$, we use eq. (164) and compute

$\displaystyle C_{3}(h_{1},h_{2},h_{3})$ $\displaystyle=C_{3}[f_{\rm
NL}(\mathcal{R}_{1}),\mathcal{R}_{2},\mathcal{R}_{3}]+C_{3}[\mathcal{R}_{1},f_{\rm
NL}(\mathcal{R}_{2}),\mathcal{R}_{3}]+C_{3}[\mathcal{R}_{1},\mathcal{R}_{2},f_{\rm
NL}(\mathcal{R}_{3})]$ (177)
$\displaystyle=2\alpha\left(\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle\langle\mathcal{R}_{1}\mathcal{R}_{3}\rangle+\langle\mathcal{R}_{2}\mathcal{R}_{1}\rangle\langle\mathcal{R}_{2}\mathcal{R}_{3}\rangle+\langle\mathcal{R}_{3}\mathcal{R}_{1}\rangle\langle\mathcal{R}_{3}\mathcal{R}_{2}\rangle\right)\,.$
(178)

If we set all points equal, we find

$\displaystyle C_{3}(h,h,h)=6\alpha\sigma_{0}^{4}\,.$ (179)

We now move to consider the possible presence of a correction at order
$O(\alpha^{2})$ that can be computed by means of eq. (172). As already
noticed, for $n=3$ the sum $\sum_{a,b}$ does not contribute to the final
result. We are left with the two terms

$\displaystyle\left[\langle f_{\rm NL}^{(1)}(\mathcal{R}_{1})f_{\rm
NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm
NL}^{(1)}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}(\mathcal{R}_{2})\rangle\right]\langle\mathcal{R}_{1}\mathcal{R}_{3}\rangle+\left[\langle
f_{\rm NL}(\mathcal{R}_{1})f_{\rm NL}^{(1)}(\mathcal{R}_{2})\rangle-\langle
f_{\rm NL}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}^{(1)}(\mathcal{R}_{2})\rangle\right]\langle\mathcal{R}_{2}\mathcal{R}_{3}\rangle\,.$
(180)

If we set the three points equal and consider explicitly $f_{\rm
NL}(\mathcal{R})=\alpha\mathcal{R}^{2}$, we find that eq. (180) reduces to
$4\alpha^{2}\langle\mathcal{R}^{3}\rangle\sigma_{0}^{2}$ which is zero because
$\langle\mathcal{R}^{3}\rangle=0$. Eq. (179), however, is not an exact result
since for $n=3$ corrections of order $O(\alpha^{3})$ are possible. Mimicking
eq. (176), we expect

$\displaystyle
C_{3}(h,h,h)=6\alpha\sigma_{0}^{4}[1+O(\alpha^{2}\sigma_{0}^{2})]\,.$ (181)

The correction of order $O(\alpha^{3})$ can not be computed by means of eq.
(172) but it is expected to be sub-leading since we work under the assumption
that $\alpha^{2}\sigma_{0}^{2}\ll 1$.

$n=4$

Consider the fourth-order cumulant $C_{4}(h_{1},h_{2},h_{3},h_{4})$. If we
limit the analysis to quadratic non-gaussianities as in $f_{\rm
NL}(\mathcal{R})=\alpha\mathcal{R}^{2}$, it is clear that there are no
corrections of order $O(\alpha)$ since $f_{\rm NL}^{(3)}(\mathcal{R})=0$. The
first non-trivial correction arises at order $O(\alpha^{2})$

$\displaystyle C_{4}(h_{1},h_{2},h_{3},h_{4})=C_{4}[f_{\rm
NL}(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2}),\mathcal{R}_{3},\mathcal{R}_{4}]+5\,{\rm
combinations}\,.$ (182)

If we take eq. (172) it is simple to see that the first two terms on the
right-hand side give vanishing contribution. In the last line, we have one
contribution corresponding to $a=b=1$ that gives two terms since one has to
sum $3\leftrightarrow 4$ exchange. We have

$\displaystyle C_{4}[f_{\rm NL}(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2}),\mathcal{R}_{3},\mathcal{R}_{4}]$
$\displaystyle=\left[\langle f_{\rm NL}^{(1)}(\mathcal{R}_{1})f_{\rm
NL}^{(1)}(\mathcal{R}_{2})\rangle-\langle f_{\rm
NL}^{(1)}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}^{(1)}(\mathcal{R}_{2})\rangle\right](\langle\mathcal{R}_{1}\mathcal{R}_{4}\rangle\langle\mathcal{R}_{2}\mathcal{R}_{3}\rangle+3\leftrightarrow
4\rangle)$
$\displaystyle=4\alpha^{2}\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle(\langle\mathcal{R}_{1}\mathcal{R}_{4}\rangle\langle\mathcal{R}_{2}\mathcal{R}_{3}\rangle+\langle\mathcal{R}_{1}\mathcal{R}_{3}\rangle\langle\mathcal{R}_{2}\mathcal{R}_{4}\rangle\rangle)\,.$
(183)

If we set all points equal and multiply by six because of eq. (182), we find

$\displaystyle
C_{4}(h,h,h,h)=48\alpha^{2}\sigma_{0}^{6}[1+O(\alpha^{2}\sigma_{0}^{2})]\,,$
(184)

where we introduced again a sub-leading correction of order $O(\alpha^{4})$.

We can check the validity of eqs. (176, 181, 184) explicitly. The reason is
that if we limit the analysis—as done here—to random variables without
derivatives, the computation of the cumulants of $h$ admits an exact
analytical solution. We can indeed extract the probability density
distribution of the random variable $h$—at a given spatial point, that is what
we need in order to compare with eqs. (176, 181, 184)—by means of the Jacobi’s
multivariate theorem in eq. (151) using the fact that $\mathcal{R}$ is a
gaussian variable with zero mean and variance $\sigma_{0}^{2}$. We find

$\displaystyle
p(h)=\frac{e^{-\mathcal{R}_{+}^{2}/2\sigma_{0}^{2}}+e^{-\mathcal{R}_{-}^{2}/2\sigma_{0}^{2}}}{\sqrt{2\pi}\sigma_{0}\sqrt{1+4\alpha
h}}\,,~{}~{}~{}~{}~{}{\rm
where}~{}~{}\mathcal{R}_{\pm}=\frac{-1\pm\sqrt{1+4\alpha h}}{2\alpha}\,,$
(185)

with $h\in[-1/4\alpha,\infty)$. By means of the Faà di Bruno’s formula, we
find the generic $n^{\rm th}$-order cumulant252525From the probability density
distribution we compute the moments
$\mu_{n}=\int_{-1/4\alpha}^{\infty}dh\,h^{n}\,p(h)$. The explicit expression
for the $n^{\rm th}$ cumulant in terms of the first $n$ moments can be
obtained by using Faà di Bruno’s formula for higher derivatives of composite
functions; explicitly, for $n\geqslant 2$, we have
$C_{n}=\sum_{k=1}^{n}(-1)^{k-1}(k-1)!B_{n,k}(0,\mu_{2},\dots,\mu_{n-k+1})$,
where $B_{n,k}$ are the incomplete Bell polynomials.

$\displaystyle
C_{n}(h,\,\dots\,,h)=2^{n-3}\sigma_{0}^{2}(\alpha\sigma_{0}^{2})^{n-2}(n+4\alpha^{2}\sigma_{0}^{2})(n-1)!\,,$
(186)

One can check that eqs. (176, 181, 184) are correctly reproduced in the
appropriate limits.

  * $\circ$

Since in the case with no derivatives the probability density distribution can
be obtained analytically (see eq. (185)), it is instructive to do the
following exercise. First of all, let us rewrite eq. (185) for the case with
$f_{\rm
NL}(\mathcal{R})=\alpha(\mathcal{R}^{2}-\langle\mathcal{R}^{2}\rangle)$ with
$\langle\mathcal{R}^{2}\rangle=\sigma_{0}^{2}$. We find

$\displaystyle
p(h)=\frac{e^{-\mathcal{R}_{+}^{2}/2\sigma_{0}^{2}}+e^{-\mathcal{R}_{-}^{2}/2\sigma_{0}^{2}}}{\sqrt{2\pi}\sigma_{0}\sqrt{1+4\alpha(h+\alpha\sigma_{0}^{2})}}\,,~{}~{}~{}~{}~{}{\rm
where}~{}~{}\mathcal{R}_{\pm}=\frac{-1\pm\sqrt{1+4\alpha(h+\alpha\sigma_{0}^{2})}}{2\alpha}\,.$
(187)

In this case we have $C_{1}(h)=0$ while $C_{n}(h,\,\dots\,,h)\equiv C_{n}(h)$
are still given by eq. (186). We can compute analytically the characteristic
function $\log\chi(\lambda)=\sum_{n=2}^{\infty}(i^{n}/n!)C_{n}(h)\lambda^{n}$.
We find

$\displaystyle\chi(\lambda)=\exp\bigg{[}-\frac{\lambda^{2}\sigma_{0}^{2}}{2(1-2i\alpha\lambda\sigma_{0}^{2})}-i\alpha\lambda\sigma_{0}^{2}\bigg{]}\frac{1}{\sqrt{1-2i\alpha\lambda\sigma_{0}^{2}}}\,,$
(188)

and one can check for consistency that $p(h)$ in eq. (187) is given by the
inverse Fourier transform

$\displaystyle
p(h)=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\lambda\,\chi(\lambda)\,e^{-ih\lambda}\,.$
(189)

Once we know $p(h)$, we can compute the integral
$\bar{p}(h_{c})\equiv\int_{h_{c}}^{\infty}dh\,p(h)$, which is the probability
to find $h$ above the threshold $h_{c}$. We find

$\displaystyle\bar{p}_{\rm NG}(h_{c})=\frac{1}{2}\bigg{[}2-{\rm
Erf}\bigg{(}\frac{1+\sqrt{1+4h_{c}\alpha+4\alpha^{2}\sigma_{0}^{2}}}{2\sqrt{2}\alpha\sigma_{0}}\bigg{)}-{\rm
Erf}\bigg{(}\frac{-1+\sqrt{1+4h_{c}\alpha+4\alpha^{2}\sigma_{0}^{2}}}{2\sqrt{2}\alpha\sigma_{0}}\bigg{)}\bigg{]}\,,$
(190)

with ${\rm Erf}(x)$ the error function. This simple case, therefore, can be
solved exactly. This means that we can use it as a playground to test the
following approximation. Let us organize the sum over cumulants in a power-
series expansion in $\alpha$. We have

$\displaystyle\sum_{n=2}^{\infty}\frac{i^{n}}{n!}C_{n}(h)\lambda^{n}=-\frac{\sigma_{0}^{2}\lambda^{2}}{2}-i\alpha\lambda^{3}\sigma_{0}^{4}+\alpha^{2}\sigma_{0}^{4}(-\lambda^{2}+2\lambda^{4}\sigma_{0}^{2})+O(\alpha^{3})\,,$
(191)

where the first term on the right-hand side corresponds to the quadratic
cumulant that reproduces the gaussian limit. We can truncate eq. (191) at some
order in $\alpha$, compute the corresponding $\bar{p}(h_{c})$ and compare with
eq. (190). At order $\alpha^{0}$, we find the gaussian result

$\displaystyle\bar{p}_{\rm G}(h_{c})=\frac{1}{2}{\rm
Erf}\left(\frac{h_{c}}{\sqrt{2}\sigma_{0}}\right)\simeq\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)\,,$
(192)

where we define $v_{c}\equiv h_{c}/\sigma_{0}$, and the last approximation
corresponds to $v_{c}\gg 1$ which is the limit that is relevant for our
analysis. Consider now the order $O(\alpha)$. Let us write eq. (190) at order
$O(\alpha)$ as
$\bar{p}_{\alpha}(h_{c})=\int_{h_{c}}^{\infty}dh\,p_{\alpha}(h)$ with

$\displaystyle p_{\alpha}(h)$
$\displaystyle=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\lambda\,\exp\bigg{[}-\frac{\sigma_{0}^{2}\lambda^{2}}{2}-i(\alpha\sigma_{0})\lambda^{3}\sigma_{0}^{3}\bigg{]}\,e^{-ih\lambda}=\frac{1}{2\pi}\int_{-\infty}^{\infty}d\lambda\,e^{-\sigma_{0}^{2}\lambda^{2}/2}\,\bigg{[}\sum_{m=0}^{\infty}\frac{(-i\alpha\sigma_{0}^{4})^{m}}{m!}\bigg{]}\,\underbrace{\lambda^{3m}\,e^{-ih\lambda}}_{i^{3m}\frac{\partial^{3m}e^{-ih\lambda}}{\partial
h^{3m}}}$
$\displaystyle=\sum_{m=0}^{\infty}\frac{(-i\alpha\sigma_{0}^{4})^{m}}{m!}i^{3m}\frac{\partial^{3m}}{\partial
h^{3m}}\underbrace{\frac{1}{2\pi}\int_{-\infty}^{\infty}d\lambda\,e^{-\sigma_{0}^{2}\lambda^{2}/2}\,e^{-ih\lambda}}_{{\rm
gaussian\,integral}\,\frac{1}{\sqrt{2\pi}\sigma_{0}}e^{-h^{2}/2\sigma_{0}^{2}}}\equiv\sum_{m=0}^{\infty}p_{\alpha}^{(m)}(h)\,.$
(193)

In practice, we can approximate $\bar{p}_{\alpha}(h_{c})$ as a sum of terms,
$\bar{p}_{\alpha}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha}(h_{c})$
each one of them defined by means of eq. ($\circ$ ‣ D), that is
$\bar{p}^{(m)}_{\alpha}(h_{c})=\int_{h_{c}}^{\infty}dh\,p^{(m)}_{\alpha}(h)$.
A similar decomposition can be defined at higher orders in $\alpha$. At order
$O(\alpha)$, we find

$\displaystyle\bar{p}^{(1)}_{\alpha}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})v_{c}^{3}\bigg{(}1-\frac{1}{v_{c}^{2}}\bigg{)}\,,$
(194) $\displaystyle\bar{p}^{(2)}_{\alpha}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})^{2}\frac{v_{c}^{6}}{2}\bigg{(}1-\frac{10}{v_{c}^{2}}+\frac{15}{v_{c}^{4}}\bigg{)}\,,$
(195) $\displaystyle\bar{p}^{(3)}_{\alpha}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})^{3}\frac{v_{c}^{9}}{6}\bigg{(}1-\frac{28}{v_{c}^{2}}+\frac{210}{v_{c}^{4}}-\frac{420}{v_{c}^{6}}+\frac{105}{v_{c}^{8}}\bigg{)}\,$
(196) $\displaystyle\bar{p}^{(4)}_{\alpha}(h_{c})$ $\displaystyle=\dots$

Since $v_{c}\gg 1$, we are tempted to approximate (by neglecting all sub-
leading terms in each of the round brackets)

$\displaystyle\bar{p}^{(m)}_{\alpha}(h_{c})=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)\frac{(\alpha\sigma_{0}v_{c}^{3})^{m}}{m!}~{}~{}~{}~{}~{}~{}\Longrightarrow~{}~{}~{}~{}~{}~{}\bar{p}_{\alpha}(h_{c})=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}+\alpha\sigma_{0}v_{c}^{3}\right)\,.$
(197)

We are now in the position to compare i) the gaussian result in eq. (192), ii)
the exact non-gaussian result in eq. (190), iii) the order $O(\alpha)$
approximation
$\bar{p}_{\alpha}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha}(h_{c})$
obtained by truncating the series for increasing values of $m$ and iv) the
order $O(\alpha)$ approximation resummed as in eq. (197). This comparison is
shown in the left panel of fig. 11.

We can now consider the order $O(\alpha^{2})$. We find

$\displaystyle\bar{p}^{(1)}_{\alpha^{2}}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})\bigg{[}v_{c}^{3}\bigg{(}1-\frac{1}{v_{c}^{2}}\bigg{)}+2v_{c}^{4}\alpha\sigma_{0}\bigg{(}1-\frac{5}{2v_{c}^{2}}\bigg{)}\bigg{]}\,,$
(198) $\displaystyle\bar{p}^{(2)}_{\alpha^{2}}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})^{2}\times$
$\displaystyle\bigg{[}\frac{v_{c}^{6}}{2}\bigg{(}1-\frac{10}{v_{c}^{2}}+\frac{15}{v_{c}^{4}}\bigg{)}+2\alpha\sigma_{0}v_{c}^{7}\bigg{(}1-\frac{29}{2v_{c}^{2}}+\frac{42}{v_{c}^{4}}-\frac{27}{2v_{c}^{6}}\bigg{)}+2(\alpha\sigma_{0})^{2}v_{c}^{8}\bigg{(}1-\frac{20}{v_{c}^{2}}+\frac{381}{4v_{c}^{4}}-\frac{363}{4v_{c}^{6}}\bigg{)}\bigg{]}\,,$
(199) $\displaystyle\bar{p}^{(3)}_{\alpha^{2}}(h_{c})$ $\displaystyle=\dots$

We are again tempted to consider the limit $v_{c}\gg 1$ in which we keep only
the leading term inside the round brackets, meaning that we write

$\displaystyle\bar{p}^{(1)}_{\alpha^{2}}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})\big{(}v_{c}^{3}+2v_{c}^{4}\alpha\sigma_{0}\big{)}\,,$
(200) $\displaystyle\bar{p}^{(2)}_{\alpha^{2}}(h_{c})$
$\displaystyle=\frac{1}{\sqrt{2\pi}v_{c}}\exp\left(-\frac{v_{c}^{2}}{2}\right)(\alpha\sigma_{0})^{2}\bigg{[}\frac{v_{c}^{6}}{2}+2\alpha\sigma_{0}v_{c}^{7}+2(\alpha\sigma_{0})^{2}v_{c}^{8}\bigg{]}\,,$
(201) $\displaystyle\bar{p}^{(3)}_{\alpha^{2}}(h_{c})$ $\displaystyle=\dots$

In this case, we find that the series
$\bar{p}_{\alpha^{2}}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha^{2}}(h_{c})$
can be resummed, and we get the neat expression

$\displaystyle\bar{p}_{\alpha^{2}}(h_{c})\simeq\frac{1}{\sqrt{2\pi}v_{c}}\exp\left[-\frac{v_{c}^{2}}{2}+\alpha\sigma_{0}v_{c}^{3}(1+2\alpha\sigma_{0}v_{c})\right]\,.$
(202)

We can do the same comparison we did before, and compare i) the gaussian
result in eq. (192), ii) the exact non-gaussian result in eq. (190), iii) the
order $O(\alpha^{2})$ approximation
$\bar{p}_{\alpha^{2}}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha^{2}}(h_{c})$
obtained by truncating the series at increasing values of $m$ and iv) the
order $O(\alpha^{2})$ approximation resummed as in eq. (202). This comparison
is shown in the right panel of fig. 11.

Figure 11:  Left panel. We compare at order $O(\alpha)$ the exact non-gaussian
expression for $\bar{p}(h_{c})$ (defined below eq. (189)) given in eq. (190)
with: i) its gaussian limit (eq. (192)), ii) its power-series expansion
$\bar{p}_{\alpha}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha}(h_{c})$ for
increasing values of $m$ (green region) and iii) the resummed expression
(labelled “exp approx”) given in eq. (197). For illustration, we take
$\sigma_{0}=0.1$ and $\alpha=0.2$. Right panel. Same as in the left panel but
at order $O(\alpha^{2})$.

The simplified setup studied in this exercise is conceptually different
compared to the actual problem we are facing (since we are only considering
here the non-gaussian random field $h$ without any information about its
derivatives) but, as we shall discuss later, we will use an approximation
scheme very similar to the one adopted above. It is, therefore, instructive to
draw some conclusions (this is particularly important in light of the fact
that in the case without derivatives we know the exact form of the probability
density distribution). Consider first the $O(\alpha)$ approximation
illustrated in the left panel of fig. 11. The sum
$\bar{p}_{\alpha}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha}(h_{c})$
(truncated at some finite $m$; in the plot we show the cases with $m\leqslant
9$) already gives a good—even though not optimal—approximation of
$\bar{p}_{\rm NG}(h_{c})$. On the contrary, the resummed expression in eq.
(197), more simple to handle, tends to overestimate the true result. Consider
now the $O(\alpha^{2})$ approximation illustrated in the right panel of fig.
11. The sum
$\bar{p}_{\alpha^{2}}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha^{2}}(h_{c})$
(truncated at some finite $m$; in the plot we show the cases with $m\leqslant
9$) gives an excellent approximation of $\bar{p}_{\rm NG}(h_{c})$. On the
contrary, it is evident that the resummed expression in eq. (202) quickly
diverges from the true result and can not be trusted.

The reason why the exponential approximation deviates from the exact result
can be identified already at order $O(\alpha)$. Consider (in units of
$\bar{p}_{\rm G}(h_{c})$) the first term in eq. (194) that we kept in our
expansion, that is $(\alpha\sigma_{0})v_{c}^{3}$, and compare it with one of
the terms that we neglected in eq. (196), for instance the fourth in the round
brackets $-70(\alpha\sigma_{0})^{3}v_{c}^{3}$. These two terms have the same
power of $v_{c}$ but the latter is parametrically suppressed by
$(\alpha\sigma_{0})^{2}$. It is true that we expect
$(\alpha\sigma_{0})^{2}\sim 10^{-2}\ll 1$ but we also have an extra numerical
factor $-70$ that partially compensate the suppression. On similar ground,
neglecting the term $-(14/3)(\alpha\sigma_{0})^{3}v_{c}^{7}$ in eq. (196)
compared to the terms that we kept in eq. (194) and eq. (195) seems not
justifiable.

We conclude that the exponential approximations in eq. (197) and eq. (202) can
not be considered as acceptably good proxies for the exact result since they
may overestimate it by many orders of magnitude. However, in order to have a
good approximation of the exact result it is more proper to decompose
$\bar{p}_{\alpha}(h_{c})=\sum_{m=0}^{\infty}\bar{p}^{(m)}_{\alpha}(h_{c})$ (or
even better at order $O(\alpha^{2})$) and sum over $m$ up to some finite value
according to the desired accuracy.

After this digression, we are finally ready to compute the non-gaussian
cumulants for the random variables $h$, $h_{i}$, $h_{ij}$. For simplicity, we
start again from the two-dimensional case, and we have six random variables
$\\{h,h_{x},h_{y},h_{xx},h_{xy},h_{yy}\\}$. We consider the non-gaussian
function $f_{\rm
NL}(\mathcal{R})=\alpha(\mathcal{R}^{2}-\langle\mathcal{R}^{2}\rangle)$ since
in this case we have $C_{1}(h)=\langle h\rangle=0$. Furthermore, in order to
fully exploit the properties of homogeneity and isotropy (see discussion below
eq. (73)), it is useful to change basis to the complex conjugated random
variables $\\{h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}}\\}$.262626
The following relations turn out to be useful. For a generic function
$f(x,y):\mathbb{R}^{2}\to\mathbb{R}$ we have $f_{z}=(f_{x}-if_{y})/2$,
$f_{z^{*}}=(f_{x}+if_{y})/2$ and  $\displaystyle
f_{zz^{*}}=\frac{1}{4}\left(f_{xx}+f_{yy}\right)\,,\hskip
14.22636ptf_{zz}=\frac{1}{4}\left(f_{xx}-f_{yy}-2if_{xy}\right)\,,\hskip
14.22636ptf_{z^{*}z^{*}}=\frac{1}{4}\left(f_{xx}-f_{yy}+2if_{xy}\right)\,.$
(203) The stationary point condition $f_{x}=f_{y}=0$ becomes $f_{z}=0$ (hence
$f_{z^{*}}=0$). The condition $f_{xx}f_{yy}-f_{xy}^{2}>0$ that separates
extrema from saddle points becomes
$f_{zz^{*}}^{2}>f_{zz}f_{z^{*}z^{*}}=|f_{zz}|^{2}$ that is
$f_{zz^{*}}<-|f_{zz}|\lor f_{zz^{*}}>|f_{zz}|$ (equivalently,
$|f_{zz^{*}}|>|f_{zz}|$). Notice that $f_{zz^{*}}$ is real. Maxima correspond
to the condition $f_{zz^{*}}<0$ while minima are identified by $f_{zz^{*}}>0$.

All first-order cumulants vanish. This is true for $C_{1}(h)=\langle
h\rangle=0$ as discussed before. As far as the other first-order cumulants are
concerned, we have for instance $C_{1}(h_{z})=\langle
h_{z}\rangle=\partial_{z}\langle h\rangle=0$ which vanishes (even without
imposing $C_{1}(h)=0$) because $\langle h\rangle$ does not depend, as a
consequence of homogeneity, on the specific spatial point at which it is
computed. Similarly, we have $C_{1}$ also for the remaining random variables.

As far as the second-order cumulants are concerned, we find the following non-
zero entries

$\bordermatrix{&{\scriptstyle h}&{\scriptstyle h_{z}}&{\scriptstyle
h_{z^{*}}}&{\scriptstyle h_{zz}}&{\scriptstyle h_{zz^{*}}}&{\scriptstyle
h_{z^{*}z^{*}}}\cr{\scriptstyle
h}&C_{2}(h,h)&0&0&0&C_{2}(h,h_{zz^{*}})&0\cr{\scriptstyle
h_{z}}&0&0&C_{2}(h_{z},h_{z^{*}})&0&0&0\cr{\scriptstyle
h_{z^{*}}}&0&C_{2}(h_{z},h_{z^{*}})&0&0&0&0\cr{\scriptstyle
h_{zz}}&0&0&0&0&0&C_{2}(h_{zz},h_{z^{*}z^{*}})\cr{\scriptstyle
h_{zz^{*}}}&C_{2}(h,h_{zz^{*}})&0&0&0&C_{2}(h_{zz^{*}},h_{zz^{*}})&0\cr{\scriptstyle
h_{z^{*}z^{*}}}&0&0&0&C_{2}(h_{zz},h_{z^{*}z^{*}})&0&0}\,,$ (204)

with

$\displaystyle C_{2}(h,h)$ $\displaystyle=$
$\displaystyle\sigma_{0}^{2}(1+2\alpha^{2}\sigma_{0}^{2})\,,$ (205)
$\displaystyle C_{2}(h_{z},h_{z^{*}})$ $\displaystyle=$
$\displaystyle\frac{\sigma_{1}^{2}}{4}(1+4\alpha^{2}\sigma_{0}^{2})\,,$ (206)
$\displaystyle C_{2}(h,h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-\frac{\sigma_{1}^{2}}{4}(1+4\alpha^{2}\sigma_{0}^{2})\,,$ (207)
$\displaystyle C_{2}(h_{zz},h_{z^{*}z^{*}})=C_{2}(h_{zz^{*}},h_{zz^{*}})$
$\displaystyle=$
$\displaystyle\frac{\sigma_{2}^{2}}{16}\left[1+4\alpha^{2}\left(\frac{2\sigma_{1}^{4}}{\sigma_{2}^{2}}+\sigma_{0}^{2}\right)\right]\,.$
(208)

At order $O(\alpha)$, all non-gaussian corrections vanish. We already computed
$C_{2}(h,h)$ in eq. (176). It is instructive to consider explicitly few more
cases. The relevant equations are eq. (135), eq. (164), eq. (168) and eq.
(172).

Consider for instance the computation of $C_{2}(h_{z},h_{z^{*}})$. We have

$\displaystyle C_{2}(h_{z},h_{z^{*}})\overset{\makebox[0.0pt]{\mbox{\tiny
eq.\,(\ref{eq:Cumutrick})}}}{~{}=~{}}\left.\partial_{z_{1}}\partial_{z_{2}^{*}}C_{2}(h_{1},h_{2})\right|_{\vec{x}_{1}=\vec{x}_{2}}$
(209)
$\displaystyle=\left.\partial_{z_{1}}\partial_{z_{2}^{*}}\left\\{C_{2}(\mathcal{R}_{1},\mathcal{R}_{2})+C_{2}[f_{\rm
NL}(\mathcal{R}_{1}),\mathcal{R}_{2}]+C_{2}[\mathcal{R}_{1},f_{\rm
NL}(\mathcal{R}_{2})]+C_{2}[f_{\rm NL}(\mathcal{R}_{1}),f_{\rm
NL}(\mathcal{R}_{2})]\right\\}\right|_{\vec{x}_{1}=\vec{x}_{2}}$
$\displaystyle=\partial_{z_{1}}\partial_{z_{2}^{*}}\big{[}\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle+\underbrace{\langle
f_{\rm
NL}^{(1)}(\mathcal{R}_{1})\rangle\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle+\langle
f_{\rm
NL}^{(1)}(\mathcal{R}_{2})\rangle\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle}_{\tiny{O(\alpha),\,\,\rm
eq.\,(\ref{eq:CumufNL1})}}+\underbrace{\langle f_{\rm
NL}(\mathcal{R}_{1})f_{\rm NL}(\mathcal{R}_{2})\rangle-\langle f_{\rm
NL}(\mathcal{R}_{1})\rangle\langle f_{\rm
NL}(\mathcal{R}_{2})\rangle}_{\tiny{O(\alpha^{2}),\,\,\rm
eq.\,(\ref{eq:CCC})}}\big{]}\big{|}_{\vec{x}_{1}=\vec{x}_{2}}$
$\displaystyle=\langle\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle+2\langle
f_{\rm
NL}^{(1)}(\mathcal{R})\rangle\langle\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle+\alpha^{2}\partial_{z_{1}}\partial_{z_{2}^{*}}\langle\mathcal{R}_{1}^{2}\mathcal{R}_{2}^{2}\rangle\big{|}_{\vec{x}_{1}=\vec{x}_{2}}$
(210)
$\displaystyle=\langle\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle+4\alpha^{2}\langle\mathcal{R}^{2}\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle\,.$
(211)

Eq. (210) follows from the fact that the expectation values $\langle f_{\rm
NL}(\mathcal{R})\rangle$ and $\langle f_{\rm NL}^{(1)}(\mathcal{R})\rangle$ do
not depend, because of homogeneity, on the specific spatial point at which
they are evaluated and, therefore, their spatial derivatives vanish. Eq. (211)
follows from the fact that $\langle f_{\rm
NL}^{(1)}(\mathcal{R})\rangle=2\alpha\langle\mathcal{R}\rangle=0$; at order
$O(\alpha^{2})$, we moved the derivatives inside the statistical average and
set $\vec{x}_{1}=\vec{x}_{2}$. In eq. (211) we have
$\langle\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle=(\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle+\langle\mathcal{R}_{y}\mathcal{R}_{y}\rangle)/4=\sigma_{1}^{2}/4$
and
$\langle\mathcal{R}^{2}\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle=\langle\mathcal{R}^{2}(\mathcal{R}_{x}^{2}+\mathcal{R}_{y}^{2})\rangle/4$.
This statistical average can be computed by means of the multivariate normal
distribution in eq. (78)

$\displaystyle\langle\mathcal{R}^{2}(\mathcal{R}_{x}^{2}+\mathcal{R}_{y}^{2})\rangle$
$\displaystyle=\int
d\mathcal{R}d\mathcal{R}_{x}d\mathcal{R}_{y}d\mathcal{R}_{xx}d\mathcal{R}_{xy}d\mathcal{R}_{yy}\,\mathcal{R}^{2}(\mathcal{R}_{x}^{2}+\mathcal{R}_{y}^{2})\,P(\mathcal{R}_{x})P(\mathcal{R}_{y})P(\mathcal{R}_{xy})P(\mathcal{R},\mathcal{R}_{xx},\mathcal{R}_{yy})$
(212) $\displaystyle=\underbrace{\int
d\mathcal{R}_{x}d\mathcal{R}_{y}(\mathcal{R}_{x}^{2}+\mathcal{R}_{y}^{2})P(\mathcal{R}_{x})P(\mathcal{R}_{y})}_{=\sigma_{1}^{2}}\underbrace{\int
d\mathcal{R}_{xy}P(\mathcal{R}_{xy})}_{=\,1}\int
d\mathcal{R}\mathcal{R}^{2}\underbrace{\int
d\mathcal{R}_{xx}d\mathcal{R}_{yy}P(\mathcal{R},\mathcal{R}_{xx},\mathcal{R}_{yy})}_{=\,\frac{1}{\sqrt{2\pi\sigma_{0}^{2}}}\exp(-\mathcal{R}^{2}/2\sigma_{0}^{2})}\,,$

and, from the explicit computation of the integrals, we get
$\langle\mathcal{R}^{2}(\mathcal{R}_{x}^{2}+\mathcal{R}_{y}^{2})\rangle=\sigma_{0}^{2}\sigma_{1}^{2}$.
All in all, we find
$C_{2}(h_{z},h_{z^{*}})=\sigma_{1}^{2}(1+4\alpha^{2}\sigma_{0}^{2})/4$. All
the remaining non-zero entries in eq. (204) can be derived in the same way.

In eq. (204) the zero entries vanish as a consequence of isotropy. We can
introduce again the parameter $\kappa\equiv(\\#\,z^{*}\,{\rm
derivatives\,in\,}C_{2})-(\\#\,z\,{\rm derivatives\,in\,}C_{2})$ and set to
zero all cumulants with $\kappa\neq 0$ since not invariant under spatial
rotations. The reasoning that we followed in the gaussian case for the
computations of the two-point correlators (see discussion below eq. (73))
applies also in the case of generic $n^{\rm th}$-order cumulants. This is
because cumulants are functions of moments (see eqs. (129-132)). Furthermore,
the relation between $C_{2}(h_{z},h_{z^{*}})$ and $C_{2}(h,h_{zz^{*}})$ can be
understood as a consequence of homogeneity.

We now move to consider the third-order cumulants. Based on isotropy, we
expect only a handful of non-zero cumulants (that are those with $\kappa=0$)
that we list in table 2.

No derivatives 2 derivatives 4 derivatives 6 derivatives $C_{3}(h,h,h)$
$C_{3}(h,h,h_{zz^{*}})$ $C_{3}(h,h_{zz^{*}},h_{zz^{*}})$
$C_{3}(h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})$ $\times 1$ $\times 3$ $\times 3$
$\times 1$ $C_{3}(h,h_{z},h_{z^{*}})$ $C_{3}(h_{z},h_{z^{*}},h_{zz^{*}})$
$C_{3}(h_{zz},h_{z^{*}z^{*}},h_{zz^{*}})$ $\times 6$ $\times 6$ $\times 6$
$C_{3}(h,h_{zz},h_{z^{*}z^{*}})$ $\times 6$
$C_{3}(h_{z},h_{z},h_{z^{*}z^{*}})$ $\times 3$
$C_{3}(h_{z^{*}},h_{z*},h_{zz})$ $\times 3$

Table 2: Third-order cumulants for the random variables
$\\{h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}}\\}$ that are non-zero
based on isotropy. For these cumulants, we have
$\kappa\equiv(\\#\,z^{*}\,{\rm{\it derivatives\,in\,}}C_{3})-(\\#\,z\,{\rm{\it
derivatives\,in\,}}C_{3})=0$. We gather together in each column cumulants with
the same number of spatial derivatives (cumulants without derivatives in the
first column, with two derivatives in the second, four in the third and so
on). For each entry, the tiny numbers in the second row indicate the
multiplicity of the corresponding cumulant due to distinct permutations of its
arguments.

We already computed $C_{3}(h,h,h)=6\alpha\sigma_{0}^{4}$ in eq. (181). As far
as the remaining cumulants are concerned, we find

$\displaystyle C_{3}(h,h_{z},h_{z^{*}})$ $\displaystyle=$
$\displaystyle\alpha\sigma_{0}^{2}\sigma_{1}^{2}\,,$ (213) $\displaystyle
C_{3}(h,h,h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-2\alpha\sigma_{0}^{2}\sigma_{1}^{2}\,,$ (214) $\displaystyle
C_{3}(h,h_{zz^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle\frac{\alpha}{8}\left(3\sigma_{1}^{4}+2\sigma_{0}^{2}\sigma_{2}^{2}\right)\,,$
(215) $\displaystyle C_{3}(h,h_{zz},h_{z^{*}z^{*}})$ $\displaystyle=$
$\displaystyle\frac{\alpha}{4}\sigma_{0}^{2}\sigma_{2}^{2}\,,$ (216)
$\displaystyle C_{3}(h_{z},h_{z^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-\frac{\alpha}{8}\sigma_{1}^{4}\,,$ (217) $\displaystyle
C_{3}(h_{z},h_{z},h_{z^{*}z^{*}})=C_{3}(h_{z^{*}},h_{z^{*}},h_{zz})$
$\displaystyle=$ $\displaystyle\frac{\alpha}{4}\sigma_{1}^{4}\,,$ (218)
$\displaystyle C_{3}(h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-\frac{3\alpha}{16}\sigma_{1}^{2}\sigma_{2}^{2}\,,$ (219)
$\displaystyle C_{3}(h_{zz},h_{z^{*}z^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-\frac{\alpha}{16}\sigma_{1}^{2}\sigma_{2}^{2}\,.$ (220)

We find that corrections at order $O(\alpha^{2})$ vanish for all third-order
cumulants listed before272727This happens because the $O(\alpha^{2})$
contribute to $C_{3}(h_{1},h_{2},h_{3})$ is:  $\Delta
C_{3}(h_{1},h_{2},h_{3})=2\alpha^{2}\left[\langle\mathcal{R}_{1}\mathcal{R}_{2}\rangle\langle\mathcal{R}_{3}^{2}(\mathcal{R}_{1}+\mathcal{R}_{2})\rangle+\langle\mathcal{R}_{2}\mathcal{R}_{3}\rangle\langle\mathcal{R}_{1}^{2}(\mathcal{R}_{2}+\mathcal{R}_{3})\rangle+\langle\mathcal{R}_{3}\mathcal{R}_{1}\rangle\langle\mathcal{R}_{2}^{2}(\mathcal{R}_{3}+\mathcal{R}_{1})\rangle\right]$
Thus, when we derive and then we set $\vec{x_{1}}=\vec{x_{2}}=\vec{x_{3}}$,
following the prescription (135) and we return to the $\\{x,y\\}$ basis, we
see that all the terms have the form
$\langle\xi_{j_{1}}\xi_{j_{2}}\xi_{j_{3}}\rangle$, where
$\xi_{j_{k}}\in\\{\mathcal{R},\mathcal{R}_{x},\mathcal{R}_{y},\mathcal{R}_{xx},\mathcal{R}_{xy},\mathcal{R}_{yy}\\}$.
These are third order moments of a gaussian multivariate distribution, which
are vanishing.. However, as already discussed in the case of eq. (181), eqs.
(213-220) are not exact because we expect the presence of $O(\alpha^{3})$
corrections we did not include. It is simple to check that properties of
homogeneity are respected by the explicit expressions in eqs. (213-220). For
instance, from $\partial_{z^{*}}C_{3}(h_{z},h_{z},h_{z^{*}})=0$ we find
$2C_{3}(h_{z},h_{z^{*}},h_{zz^{*}})=-C_{3}(h_{z},h_{z},h_{z^{*}z^{*}})$ that
is indeed verified by eq. (217) and eq. (218). Similarly, from
$\partial_{z^{*}}C_{3}(h,h,h_{z})=0$ we have
$C_{3}(h,h,h_{zz^{*}})=-2C_{3}(h,h_{z},h_{z^{*}})$ that is indeed verified by
eq. (213) and eq. (214).

After computing the cumulants at the desired order in $\alpha$, we are finally
ready to take the last step in our computation. From the cumulants, we will
reconstruct the characteristic function and, via an inverse Fourier transform,
the probability density distribution. Few technical remarks are in order
because it is crucial to properly identify the variables participating to the
Fourier transforms. The inverse Fourier transform of the first line in eq. (D)
is given by

$\displaystyle
P(\xi_{1},\dots,\xi_{N})=\int\frac{d\lambda_{1}}{(2\pi)}\dots\frac{d\lambda_{N}}{(2\pi)}\,\chi(\lambda_{1},\dots,\lambda_{N})\exp\left[-i(\xi_{1}\lambda_{1}+\dots+\xi_{N}\lambda_{N})\right]\,.$
(221)

Both eq. (D) and eq. (221) are valid for real variables. The simplest
identification, therefore, would be
$\\{\xi_{i=1,\dots,6}\\}=\\{h,h_{x},h_{y},h_{xx},h_{xy},h_{yy}\\}$. However,
we found more efficient to work with
$\\{h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}}\\}$ since properties
like isotropy become more transparent. Among these variables, $h$ and
$h_{zz^{*}}$ are real while $h_{z}$ and $h_{zz}$ are complex (with conjugated
variables given by $h_{z^{*}}$ and $h_{z^{*}z^{*}}$, respectively). In this
case a suitable choice of real variables, therefore, is
$\\{\xi_{i=1,\dots,6}\\}=\\{h,{\rm Re}h_{z},{\rm Im}h_{z},{\rm
Re}h_{zz},h_{zz^{*}},{\rm Im}h_{zz}\\}$ with the obvious relations

$\displaystyle{\rm Re}h_{z}=(h_{z}+h_{z^{*}})/2\,,~{}~{}~{}~{}~{}{\rm
Im}h_{z}=-i(h_{z}-h_{z^{*}})/2\,,~{}~{}~{}~{}~{}{\rm
Re}h_{zz}=(h_{zz}+h_{z^{*}z^{*}})/2\,,~{}~{}~{}~{}~{}{\rm
Im}h_{zz}=-i(h_{zz}-h_{z^{*}z^{*}})/2\,.$ (222)

Consequently, we introduce their Fourier counterparts
$\\{\lambda_{i=1,\dots,6}\\}\equiv\\{\lambda,\lambda_{{\rm
Re}h_{z}},\lambda_{{\rm Im}h_{z}},\lambda_{{\rm
Re}h_{zz}},\lambda_{zz^{*}},\lambda_{{\rm Im}h_{zz}}\\}$ such that the
characteristic function in eq. (D) is given by the integral

$\displaystyle\int dh\,d{\rm Re}h_{z}\,d{\rm Im}h_{z}\,d{\rm
Re}h_{zz}\,dh_{zz^{*}}\,d{\rm
Im}h_{zz}P(h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}})\exp[i(h\lambda+{\rm
Re}h_{z}\lambda_{{\rm Re}h_{z}}+{\rm Im}h_{z}\lambda_{{\rm
Im}h_{z}}+\dots)]\,.$ (223)

Consider, for instance, the two terms ${\rm Re}h_{z}\lambda_{{\rm
Re}h_{z}}+{\rm Im}h_{z}\lambda_{{\rm Im}h_{z}}$ in the argument of the
exponential. By means of eqs. (222) we can simply rewrite ${\rm
Re}h_{z}\lambda_{{\rm Re}h_{z}}+{\rm Im}h_{z}\lambda_{{\rm
Im}h_{z}}=h_{z}\lambda_{z}^{*}+h_{z^{*}}\lambda_{z}$ if we define
$\lambda_{z}\equiv(\lambda_{{\rm Re}h_{z}}+i\lambda_{{\rm Im}h_{z}})/2$.
Notice that in this case we have $\lambda_{{\rm Re}h_{z}}=2{\rm
Re}\lambda_{z}$ and $\lambda_{{\rm Im}h_{z}}=2{\rm Im}\lambda_{z}$. We can,
therefore, use $\lambda_{z}$ as the complex Fourier variable associated to
$h_{z^{*}}$ (and, correspondingly, its conjugated $\lambda_{z}^{*}$ associated
to $h_{z}$). Similarly, we introduce $\lambda_{zz}\equiv(\lambda_{{\rm
Re}h_{zz}}+i\lambda_{{\rm Im}h_{zz}})/2$, and the exponential function in eq.
(223) reads
$\exp\\{i[h\lambda+h_{zz^{*}}\lambda_{zz^{*}}+(h_{z}\lambda_{z}^{*}+h_{zz}\lambda_{zz}^{*}+c.c.)]\\}$.
We are now in the position to consider explicitly the inverse Fourier
transform of eq. (223). As integration variables, instead of $\lambda_{{\rm
Re}h_{z}},\lambda_{{\rm Im}h_{z}},\lambda_{{\rm Re}h_{zz}},\lambda_{{\rm
Im}h_{zz}}$ we use $\lambda_{z}$ and $\lambda_{zz}$ introduced before by means
of the relations $\lambda_{{\rm Re}h_{z}}=2{\rm Re}\lambda_{z}$,
$\lambda_{{\rm Im}h_{z}}=2{\rm Im}\lambda_{z}$, $\lambda_{{\rm
Re}h_{zz}}=2{\rm Re}\lambda_{zz}$ and $\lambda_{{\rm Im}h_{zz}}=2{\rm
Im}\lambda_{zz}$. We have

$\displaystyle P(h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}})=$ (224)
$\displaystyle\int\frac{d\lambda}{(2\pi)}\frac{d{\rm
Re}\lambda_{z}}{\pi}\frac{d{\rm
Im}\lambda_{z}}{\pi}\frac{d\lambda_{zz^{*}}}{(2\pi)}\frac{d{\rm
Re}\lambda_{zz}}{\pi}\frac{d{\rm
Im}\lambda_{zz}}{\pi}\,\chi(\lambda,\lambda_{z},\lambda_{zz},\lambda_{zz^{*}})\,\exp\left\\{-i\left[h\lambda+h_{zz^{*}}\lambda_{zz^{*}}+\left(h_{z}\lambda_{z}^{*}+h_{zz}\lambda_{zz}^{*}+c.c.\right)\right]\right\\}\,,$

with of course $\lambda_{z}={\rm Re}\lambda_{z}+i{\rm Im}\lambda_{z}$ and
$\lambda_{zz}={\rm Re}\lambda_{zz}+i{\rm Im}\lambda_{zz}$. The natural log of
the characteristic function
$\log\chi(\lambda,\lambda_{z},\lambda_{zz},\lambda_{zz^{*}})$ is written,
according to eq. (128), as a series expansion in terms of the cumulants with
respect to the Fourier variables
$\\{\lambda,\lambda_{z},\lambda_{zz},\lambda_{zz^{*}}\\}$ (remember that
$\lambda_{z}$ and $\lambda_{zz}$ are complex variables while $\lambda$ and
$\lambda_{zz^{*}}$ are real). The first non-gaussian correction arises at
order $O(\alpha)$ and reads

$\displaystyle\log\chi(\lambda,$
$\displaystyle\lambda_{z},\lambda_{zz},\lambda_{zz^{*}})=$ (225)
$\displaystyle-\frac{1}{2}C_{2}(h,h)\lambda^{2}-C_{2}(h_{z},h_{z^{*}})|\lambda_{z}|^{2}-C_{2}(h,h_{zz^{*}})\lambda\lambda_{zz^{*}}-C_{2}(h_{zz},h_{z^{*}z^{*}})|\lambda_{zz}|^{2}-\frac{1}{2}C_{2}(h_{zz^{*}},h_{zz^{*}})\lambda_{zz^{*}}^{2}$
$\displaystyle-\frac{i}{6}C_{3}(h,h,h)\lambda^{3}-\frac{i}{2}C_{3}(h,h,h_{zz^{*}})\lambda^{2}\lambda_{zz^{*}}-iC_{3}(h,h_{z},h_{z^{*}})\lambda|\lambda_{z}|^{2}-\frac{i}{2}C_{3}(h,h_{zz^{*}},h_{zz^{*}})\lambda\lambda_{zz^{*}}^{2}$
$\displaystyle-
iC_{3}(h,h_{zz},h_{z^{*}z^{*}})\lambda|\lambda_{zz}|^{2}-iC_{3}(h_{z},h_{z^{*}},h_{zz^{*}})|\lambda_{z}|^{2}\lambda_{zz^{*}}-\frac{i}{2}C_{3}(h_{z},h_{z},h_{z^{*}z^{*}})(\lambda_{z}^{*})^{2}\lambda_{zz}$
$\displaystyle-\frac{i}{2}C_{3}(h_{z^{*}},h_{z^{*}},h_{zz})\lambda_{z}^{2}\lambda_{zz}^{*}-\frac{i}{6}C_{3}(h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})\lambda_{zz^{*}}^{3}-iC_{3}(h_{zz},h_{z^{*}z^{*}},h_{zz^{*}})|\lambda_{zz}|^{2}\lambda_{zz^{*}}\,,$

where the second-order cumulants given in eq. (204) are taken at order
$O(\alpha)$ and, therefore, coincide with their gaussian limit. Notice that
each term in the series expansion enters with a coefficient that counts its
multiplicity (explicitly written in table 2). This follows from the fact that
in the last line in eq. (D) there are different terms that give the same
contribution. The integration in eq. (224) gives the probability density
distribution $P(h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}})$. After
computing $P(h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}})$, it is also
possible to reconstruct $P(h,h_{x},h_{y},h_{xx},h_{xy},h_{yy})$ by means of
the Jacobi’s multivariate theorem in eq. (139) using the transformation in eq.
(203).

It is possible to check this procedure in the gaussian limit $\alpha\to 0$. In
the gaussian limit only the second-order cumulants survive (with $\alpha\to 0$
taken in the corresponding expressions given in eq. (204)). The inverse
Fourier transform in eq. (224) can be computed analytically in the gaussian
limit, and one gets
$P(\mathcal{R},\mathcal{R}_{z},\mathcal{R}_{z^{*}},\mathcal{R}_{zz},\mathcal{R}_{zz^{*}},\mathcal{R}_{z^{*}z^{*}})$
where we used $\mathcal{R}$ instead of $h$ since $\alpha\to 0$. The
determinant of the Jacobian matrix is $\left|{\rm
det}J(\mathcal{R},\mathcal{R}_{x},\mathcal{R}_{y},\mathcal{R}_{xx},\mathcal{R}_{xy},\mathcal{R}_{yy})\right|=1/64$,
and Jacobi’s multivariate theorem gives precisely eq. (78). This is a non-
trivial check that our procedure is technically correct.

In the presence of local non-gaussianities, the computation of eq. (224) is
much more complicated because
$\chi(\lambda,\lambda_{z},\lambda_{zz},\lambda_{zz^{*}})$ in eq. (225) is an
exponential function whose argument is a polynomial of cubic order. A possible
solution strategy is the following. Let us write eq. (224) in the schematic
form (we use for illustration a generic set of variables
$\\{\lambda_{i=1,\dots,N}\\}$ that in the actual computation must be replaced
with $\\{\lambda,\lambda_{z},\lambda_{zz},\lambda_{zz^{*}}\\}$)

$\displaystyle\chi(\lambda_{i})=\exp[p_{2}(\lambda_{i})+p_{3}(\lambda_{i})]=\exp[p_{2}(\lambda_{i})]\exp[p_{3}(\lambda_{i})]=\exp[p_{2}(\lambda_{i})]\sum_{m=0}^{\infty}\frac{1}{m!}p_{3}(\lambda_{i})^{m}\,,$
(226)

where $p_{2}(\lambda_{i})$ is the quadratic polynomial in $\\{\lambda_{i}\\}$
that corresponds to the gaussian limit $\alpha\to 0$ while
$p_{3}(\lambda_{i})$ is the cubic polynomial that contains $O(\alpha)$
deviations. At a given order $m$ in the series expansion defined by eq. (226),
the inverse Fourier integral in eq. (224) can be computed analytically. This
is because we can make use of the following identity (the sum over $k$ is
understood)

$\displaystyle\int
d\lambda_{1}\dots\lambda_{N}(\lambda_{i}^{n}\dots\lambda_{j}^{m})\exp[p_{2}(\lambda_{i})]\exp(-i\lambda_{k}h_{k})=(i)^{n}\dots(i)^{m}\frac{\partial^{n}}{\partial
h_{i}^{n}}\dots\frac{\partial^{m}}{\partial h_{j}^{m}}\int
d\lambda_{1}\dots\lambda_{N}\exp[p_{2}(\lambda_{i})]\exp(-i\lambda_{k}h_{k})\,,$
(227)

where the last integral is nothing but the gaussian integral that can be
easily computed. We go along with this strategy, that we summarize for ease of
reading in the following.

  * i)

For fixed $m$, we compute the probability density distribution in eq. (224)
performing the inverse Fourier transform by means of the trick in eq. (226)
and eq. (227).

  * ii)

We transform $P(h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}})$ into
$P(h,h_{x},h_{y},h_{xx},h_{xy},h_{yy})$ by means of Jacobi’s multivariate
theorem (as illustrated before for the gaussian case).

  * iii)

We set $h_{x}=h_{y}=0$, and compute the number density of local maxima $n_{\rm
max}(h)$ according to the general definition in eq. (69) (with of course $h$
instead of $\mathcal{R}$, and working in two spatial dimensions).

  * iv)

We use the same change of variables proposed in eq. (82) (with again $h$
instead of $\mathcal{R}$), and we integrate over $\theta$, $r$, $h$ and $s$ in
the domain defining maxima, as discussed in the gaussian case; similarly, we
implement the same lower bound of integration over $s$ given by the threshold
condition for black hole formation in eq. (104). Notice that in the following
we will set the lower limit of integration over $h$ to be $h_{\rm
min}=-\infty$ despite the fact that we have $h_{\rm
min}=-(1+4\alpha^{2}\sigma_{0}^{2})/4\alpha$ for
$h=\mathcal{R}+\alpha(\mathcal{R}^{2}-\sigma_{0}^{2})$ and
$\mathcal{R}\in(-\infty,+\infty)$. Approximating $h_{\rm min}$ with its
gaussian value allows to obtain close analytical formulas. We shall justify
the validity of this approximation at the end of this section.

  * v)

Finally, we define the mass fraction of black holes as in eq. (103).

We now illustrate the result of this procedure. Needless to say, the term with
$m=0$ reproduces the gaussian result in eq. (103). The first correction is
given by the term with $m=1$. Let us introduce for the sake of simplicity the
notation $\beta=\beta_{\rm G}+\sum_{m=1}^{\infty}\beta_{{\rm
NG}}^{(\alpha,m)}$ with $\beta_{\rm G}$ given by eq. (103). We find

$\displaystyle\beta_{\rm
NG}^{(\alpha,1)}=\left(\frac{\alpha\sigma_{1}^{2}}{\sigma_{2}}\right)\frac{2\sqrt{2}}{\pi^{3/2}}e^{-6s_{\rm
min}^{2}/\sigma_{2}^{2}}\left[\frac{s_{\rm
min}^{2}}{\sigma_{2}^{2}}+e^{4s_{\rm
min}^{2}/\sigma_{2}^{2}}\left(-\frac{s_{\rm
min}^{2}}{\sigma_{2}^{2}}+4\frac{s_{\rm
min}^{4}}{\sigma_{2}^{4}}\right)\right]\,.$ (228)

For $m=2$, we find

$\displaystyle\beta_{\rm
NG}^{(\alpha,2)}=\left(\frac{\alpha\sigma_{1}^{2}}{\sigma_{2}}\right)^{2}\frac{e^{-6s_{\rm
min}^{2}/\sigma_{2}^{2}}}{3\pi^{3/2}}\bigg{\\{}$ $\displaystyle
2\sqrt{3\pi}e^{6s_{\rm min}^{2}/\sigma_{2}^{2}}{\rm
Erfc}\bigg{(}\frac{\sqrt{6}s_{\rm
min}}{\sigma_{2}}\bigg{)}+3\sqrt{2}\frac{s_{\rm
min}}{\sigma_{2}}\bigg{[}7-40\frac{s_{\rm
min}^{2}}{\sigma_{2}^{2}}+48\frac{s_{\rm min}^{4}}{\sigma_{2}^{4}}$
$\displaystyle+e^{4s_{\rm min}^{2}/\sigma_{2}^{2}}\bigg{(}-3+68\frac{s_{\rm
min}^{2}}{\sigma_{2}^{2}}-176\frac{s_{\rm
min}^{4}}{\sigma_{2}^{4}}+64\frac{s_{\rm
min}^{6}}{\sigma_{2}^{6}}\bigg{)}\bigg{]}\bigg{\\}}\,.$ (229)

For increasing values of $m$ the corresponding expressions of $\beta_{\rm
NG}^{(m)}$ become more lengthy and less transparent, and we do not report them
explicitly. On the contrary, we make use of the same expansion introduced in
eq. (103) for $s_{\rm min}/\sigma_{2}\gtrsim 1$. Interestingly, in this limit
we find that the generic term $\beta_{\rm NG}^{(\alpha,m)}$ can be written as

$\displaystyle\beta_{\rm NG}^{(\alpha,m)}\simeq\bigg{[}\bigg{(}\frac{s_{\rm
min}^{3}}{\sigma_{2}^{3}}\bigg{)}\frac{\alpha\sigma_{1}^{2}}{\sigma_{2}}\bigg{]}^{m}\frac{2^{4m-1/2}}{\pi^{3/2}m!}\bigg{(}\frac{s_{\rm
min}}{\sigma_{2}}\bigg{)}e^{-2s_{\rm min}^{2}/\sigma_{2}^{2}}\,,$ (230)

and the sum $\sum_{m=1}^{\infty}\beta_{\rm NG}^{(\alpha,m)}$ admits the
following analytical form

$\displaystyle\sum_{m=1}^{\infty}\beta_{\rm
NG}^{(\alpha,m)}\simeq\frac{1}{\sqrt{2}\pi^{3/2}}\bigg{(}\frac{s_{\rm
min}}{\sigma_{2}}\bigg{)}e^{-2s_{\rm
min}^{2}/\sigma_{2}^{2}}\bigg{\\{}-1+\exp\bigg{[}16\left(\frac{s_{\rm
min}^{3}}{\sigma_{2}^{3}}\right)\frac{\alpha\sigma_{1}^{2}}{\sigma_{2}}\bigg{]}\bigg{\\}}\,.$
(231)

The sum $\beta=\beta_{\rm G}+\sum_{m=1}^{\infty}\beta_{\rm NG}^{(\alpha,m)}$
gives

$\displaystyle\beta\simeq\frac{1}{\sqrt{2}\pi^{3/2}}\left(\frac{s_{\rm
min}}{\sigma_{2}}\right)\exp\left[-\frac{2s_{\rm
min}^{2}}{\sigma_{2}^{2}}+16\left(\frac{s_{\rm
min}^{3}}{\sigma_{2}^{3}}\right)\frac{\alpha\sigma_{1}^{2}}{\sigma_{2}}\right]\,.$
(232)

Figure 12:  Left panel. We compare the gaussian approximation for $\beta$
given in eq. (103) with i) its power-series expansion $\beta=\beta_{\rm
G}+\sum_{m=1}^{\infty}\beta_{\rm NG}^{(\alpha,m)}$ computed at order
$O(\alpha)$ for increasing values of $m$ (green region, with the thin solid
green lines that label the case $m=1,\dots,6$) and ii) the resummed expression
(labelled “exp approx”) given in eq. (232). Right panel. Same as in the left
panel but at order $O(\alpha^{2})$ and up to $m=4$. In both panels we consider
the log-normal model for the power spectrum, and we use the benchmark values
$v=0.6$ and $\mathcal{P}_{\mathcal{R}}(k_{\star})=5\times 10^{-3}$. We take
$\alpha=0.2$.

Notice that the structure of this equation is analogue to the one that we
found in eq. (197).

Before proceeding, let us come back to the issue of the lower limit of
integration over $h$. As anticipated, from the very same definition
$h=\mathcal{R}+\alpha(\mathcal{R}^{2}-\sigma_{0}^{2})$ it follows that
$h>h_{\rm min}=-(1+\alpha^{2}\sigma_{0}^{2})/4\alpha$. To be precise,
therefore, this lower limit of integration should be implemented when
computing the non-gaussian contributions $\beta_{\rm NG}^{(\alpha,m)}$.
However, implementing the gaussian lower limit of integration $h_{\rm
min}=-\infty$ does not change qualitatively our results, and the reason is the
following. When computing $\beta$, we only focus on maxima with large
curvature $-\triangle h/\sigma_{2}$ since they correspond to peak of the
overdensity field. The value of $-\triangle h/\sigma_{2}$ is positively
correlated with the value of $h/\sigma_{0}$, and the amount of correlation is
controlled by the parameter $\gamma$. This fact was illustrated in appendix B
in the gaussian case, and it remains valid also in the presence of non-
gaussianities (for the values of $\alpha$ that are relevant for our analysis).
When $\gamma\to 1$, the correlation is maximal and regions with large
$-\triangle h/\sigma_{2}$ are also regions with positive $h/\sigma_{0}\gtrsim
1$. This means that regions where $h$ takes negative values are completely
irrelevant for our purposes. In the cases that are relevant for our analysis,
we typically have $\gamma\sim 0.6$; we checked numerically (using for
simplicity the log-normal power spectrum, and integrating numerically over
$s$) that the value of $\beta_{\rm NG}^{(\alpha,m)}$ computed with the proper
lower bound of integration over $h$ remain equal to those obtained with the
gaussian limit $h_{\rm min}=-\infty$. In order to see some deviation, one
should consider the opposite limit $\gamma\to 0$. In this limit $-\triangle
h/\sigma_{2}$ and $h/\sigma_{0}$ are only weakly correlated and it is
therefore possible to find regions with large $-\triangle h/\sigma_{2}$ in
which $h$ takes negative values. In such situation, of course, implementing
the correct lower limit of integration over $h$ becomes important. However,
cases with $\gamma\to 0$ fall outside the class of models we are considering
in this paper.

We compare in the left panel of fig. 12 the series $\beta=\beta_{\rm
G}+\sum_{m=1}^{m_{\rm max}}\beta_{\rm NG}^{(\alpha,m)}$, truncated at some
finite $m_{\rm max}$, with the resummed expression in eq. (232). We show the
cases $m_{\rm max}=1,\dots,6$ and we use the log-normal power spectrum in eq.
(89) to compute the spectral parameters. The situation is very similar to what
already discussed in fig. 11, and we can exploit the knowledge that we gained
from that case to argue the following conclusions. The series
$\beta=\beta_{\rm G}+\sum_{m=1}^{m_{\rm max}}\beta_{\rm NG}^{(\alpha,m)}$
shows a convergence towards the order $O(\alpha)$ approximation (that is, we
recall, the approximation in which we only include the leading part of the
third-order cumulants) of the exact non-gaussian result while the exponential
approximation overestimates the abundance.

We move to consider corrections at order $O(\alpha^{2})$. We already computed
in eqs. (205-208) the second-order cumulants at order $O(\alpha^{2})$. In
addition, we need to compute the fourth-order cumulants at order
$O(\alpha^{2})$. We computed in eq. (184) the simplest one without
derivatives, $C_{4}(h,h,h,h)=48\alpha^{2}\sigma_{0}^{6}$.

No derivatives 2 derivatives 4 derivatives 6 derivatives 8 derivatives
$C_{4}(h,h,h,h)$ $C_{4}(h,h,h,h_{zz^{*}})$ $C_{4}(h,h,h_{zz},h_{z^{*}z^{*}})$
$C_{4}(h_{z},h_{z^{*}},h_{zz^{*}},h_{zz^{*}})$
$C_{4}(h_{zz},h_{zz},h_{z^{*}z^{*}},h_{z^{*}z^{*}})$ $\times 1$ $\times 4$
$\times 12$ $\times 12$ $\times 6$ $C_{4}(h,h,h_{z},h_{z^{*}})$
$C_{4}(h,h,h_{zz^{*}},h_{zz^{*}})$
$C_{4}(h_{z},h_{z},h_{z^{*}z^{*}},h_{zz^{*}})$
$C_{4}(h_{zz^{*}},h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})$ $\times 12$ $\times 6$
$\times 12$ $\times 1$ $C_{4}(h,h_{z},h_{z^{*}},h_{zz^{*}})$
$C_{4}(h_{z^{*}},h_{z^{*}},h_{zz},h_{zz^{*}})$
$C_{4}(h_{zz},h_{z^{*}z^{*}},h_{zz^{*}},h_{zz^{*}})$ $\times 24$ $\times 12$
$\times 12$ $C_{4}(h,h_{z},h_{z},h_{z^{*}z^{*}})$
$C_{4}(h,h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})$ $\times 12$ $\times 4$
$C_{4}(h,h_{z^{*}},h_{z^{*}},h_{zz})$
$C_{4}(h,h_{zz},h_{z^{*}z^{*}},h_{zz^{*}})$ $\times 12$ $\times 24$
$C_{4}(h_{z},h_{z},h_{z^{*}},h_{z^{*}})$
$C_{4}(h_{z},h_{z^{*}},h_{zz},h_{z^{*}z^{*}})$ $\times 6$ $\times 24$

Table 3: Fourth-order cumulants for the random variables
$\\{h,h_{z},h_{z^{*}},h_{zz},h_{zz^{*}},h_{z^{*}z^{*}}\\}$ that are non-zero
based on isotropy. For these cumulants, we have
$\kappa\equiv(\\#\,z^{*}\,{\rm{\it derivatives\,in\,}}C_{4})-(\\#\,z\,{\rm{\it
derivatives\,in\,}}C_{4})=0$. As in table 2, we gather together in each column
cumulants with the same number of spatial derivatives (cumulants without
derivatives in the first column, with two derivatives in the second, four in
the third and so on). For each entry, the tiny numbers in the second row
indicate the multiplicity of the corresponding cumulant due to distinct
permutations of its arguments.

The remaining non-zero cumulants are summarized in table 3. We find the
following explicit expressions

$\displaystyle C_{4}(h,h,h,h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-18\alpha^{2}\sigma_{0}^{4}\sigma_{1}^{2}\,,$ (233)
$\displaystyle C_{4}(h,h,h_{z},h_{z^{*}})$ $\displaystyle=$ $\displaystyle
6\alpha^{2}\sigma_{0}^{4}\sigma_{1}^{2}\,,$ (234) $\displaystyle
C_{4}(h,h,h_{zz},h_{z^{*}z^{*}})$ $\displaystyle=$
$\displaystyle\frac{3\alpha^{2}}{2}\sigma_{0}^{4}\sigma_{2}^{2}\,,$ (235)
$\displaystyle C_{4}(h,h,h_{zz^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{2}\sigma_{0}^{2}(10\sigma_{1}^{4}+3\sigma_{0}^{2}\sigma_{2}^{2})\,,$
(236) $\displaystyle C_{4}(h,h_{z},h_{z^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-\frac{3\alpha^{2}}{2}\sigma_{0}^{2}\sigma_{1}^{4}\,,$ (237)
$\displaystyle
C_{4}(h,h_{z},h_{z},h_{z^{*}z^{*}})=C_{4}(h,h_{z^{*}},h_{z^{*}},h_{zz})$
$\displaystyle=$ $\displaystyle\alpha^{2}\sigma_{0}^{2}\sigma_{1}^{4}\,,$
(238) $\displaystyle C_{4}(h_{z},h_{z},h_{z^{*}},h_{z^{*}})$ $\displaystyle=$
$\displaystyle 2\alpha^{2}\sigma_{0}^{2}\sigma_{1}^{4}\,,$ (239)
$\displaystyle C_{4}(h_{z},h_{z^{*}},h_{zz^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{4}\sigma_{1}^{2}(\sigma_{1}^{4}+\sigma_{0}^{2}\sigma_{2}^{2})\,,$
(240) $\displaystyle
C_{4}(h_{z},h_{z},h_{z^{*}z^{*}},h_{zz^{*}})=C_{4}(h_{z^{*}},h_{z^{*}},h_{zz},h_{zz^{*}})$
$\displaystyle=$ $\displaystyle 0\,,$ (241) $\displaystyle
C_{4}(h,h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})$ $\displaystyle=$
$\displaystyle-\frac{3\alpha^{2}}{4}\sigma_{1}^{2}(\sigma_{1}^{4}+2\sigma_{0}^{2}\sigma_{2}^{2})\,,$
(242) $\displaystyle C_{4}(h,h_{zz},h_{z^{*}z^{*}},h_{zz^{*}})$
$\displaystyle=$
$\displaystyle-\frac{\alpha^{2}}{2}\sigma_{0}^{2}\sigma_{1}^{2}\sigma_{2}^{2}\,,$
(243) $\displaystyle C_{4}(h_{z},h_{z^{*}},h_{zz},h_{z^{*}z^{*}})$
$\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{4}\sigma_{1}^{2}(\sigma_{1}^{4}+\sigma_{0}^{2}\sigma_{2}^{2})\,,$
(244) $\displaystyle C_{4}(h_{zz},h_{zz},h_{z^{*}z^{*}},h_{z^{*}z^{*}})$
$\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{8}\sigma_{0}^{2}\sigma_{2}^{4}\,,$ (245)
$\displaystyle C_{4}(h_{zz^{*}},h_{zz^{*}},h_{zz^{*}},h_{zz^{*}})$
$\displaystyle=$
$\displaystyle\frac{3\alpha^{2}}{16}\sigma_{2}^{2}(3\sigma_{1}^{4}+\sigma_{0}^{2}\sigma_{2}^{2})\,,$
(246) $\displaystyle C_{4}(h_{zz},h_{z^{*}z^{*}},h_{zz^{*}},h_{zz^{*}})$
$\displaystyle=$
$\displaystyle\frac{\alpha^{2}}{32}\sigma_{2}^{2}(3\sigma_{1}^{4}+2\sigma_{0}^{2}\sigma_{2}^{2})\,,$
(247)

that we derive as before using eq. (135) and eq. (172). Conceptually, the rest
of the computation follows again points i)-v) discussed below eq. (227). At
the technical level, instead of eq. (226) we now have

$\displaystyle\chi(\lambda_{i})=\exp[p_{2}(\lambda_{i})+p_{4}(\lambda_{i})]=\exp[p_{2}(\lambda_{i})]\exp[p_{4}(\lambda_{i})]=\exp[p_{2}(\lambda_{i})]\sum_{m=0}^{\infty}\frac{1}{m!}p_{4}(\lambda_{i})^{m}\,,$
(248)

where $p_{4}(\lambda_{i})$ is the quartic polynomial that contains $O(\alpha)$
and $O(\alpha^{2})$ deviations. The computation of the abundance follows the
same prescription discussed before, and we introduce the series expansion
$\beta=\beta_{\rm G}+\sum_{m=1}^{m_{\rm max}}\beta_{{\rm
NG}}^{(\alpha^{2},m)}$ with $\beta_{\rm G}$ given by eq. (103) and
$\beta_{{\rm NG}}^{(\alpha^{2},m)}$ which corresponds to the series in eq.
(248) truncated at some $m>0$. For instance, we find (with $\bar{s}_{\rm
min}\equiv s_{\rm min}/\sigma_{2}$ and
$\gamma=\sigma_{1}^{2}/\sigma_{2}\sigma_{0}$)

$\displaystyle\beta_{{\rm NG}}^{(\alpha^{2},1)}=$ (249)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}\frac{\sqrt{2}\alpha\sigma_{0}\bar{s}_{\rm
min}}{\pi^{3/2}}e^{-2\bar{s}_{\rm
min}^{2}}\left[\alpha(2+\gamma^{2})\sigma_{0}-2\gamma\bar{s}_{\rm
min}-4\alpha(5+8\gamma^{2})\sigma_{0}\bar{s}_{\rm min}^{2}+8\gamma\bar{s}_{\rm
min}^{3}+16\alpha(1+3\gamma^{2})\sigma_{0}\bar{s}_{\rm
min}^{4}\right]+O(e^{-6\bar{s}_{\rm min}^{2}})\,.$

The exact analytic expressions for $\beta_{{\rm NG}}^{(\alpha^{2},m)}$ are
quite lengthly and we do not report them here explicitly. The sum
$\sum_{m=1}^{\infty}\beta_{{\rm NG}}^{(\alpha^{2},m)}$ admits an analytical
expression only if we retain, in each one of the $\beta_{{\rm
NG}}^{(\alpha^{2},m)}$, the highest power of $\bar{s}_{\rm min}$ in the
polynomial that multiplies the leading exponential suppression. We find

$\displaystyle\beta\simeq\frac{1}{\sqrt{2}\pi^{3/2}}\left(\frac{s_{\rm
min}}{\sigma_{2}}\right)\exp\left\\{-\frac{2s_{\rm
min}^{2}}{\sigma_{2}^{2}}+16\left(\frac{s_{\rm
min}^{3}}{\sigma_{2}^{3}}\right)\frac{\alpha\sigma_{1}^{2}}{\sigma_{2}}+32\left(\frac{s_{\rm
min}^{4}}{\sigma_{2}^{4}}\right)\alpha^{2}\left[3\left(\frac{\sigma_{1}^{2}}{\sigma_{2}}\right)^{2}+\sigma_{0}^{2}\right]\right\\}\,.$
(250)

The structure of this expression is analogue to eq. (202). As shown in fig.
12, the exponential approximation contained in eqs. (232) and (250) gets worse
for increasing $s_{\rm min}/\sigma_{2}$ and diverges from the actual result,
overstimating the abundance. The failure of the exponential approximation is
analogue to the one discussed in section III.2 in the context of threshold
statistics.

##  Appendix E Towards an exact computation

The generic $n$-th order cumulant can be exactly computed if we specify the
explicit form of the non-gaussianity. In our case
$f_{\mathrm{NL}}(\mathcal{R})=\alpha(\mathcal{R}^{2}-\langle\mathcal{R}^{2}\rangle)$.

The translation only affects the first-order cumulant:

$\displaystyle
C_{1}[\mathcal{R}+\alpha(\mathcal{R}^{2}-\langle\mathcal{R}^{2}\rangle)]=C_{1}(\mathcal{R}+\alpha\mathcal{R}^{2})-\alpha\langle\mathcal{R}^{2}\rangle,$
(251) $\displaystyle