# Solving peak theory in the presence of local non-gaussianities

Flavio Riccardia,b<EMAIL_ADDRESS>Marco Taosoc
<EMAIL_ADDRESS>Alfredo Urbanod,e<EMAIL_ADDRESS>aSISSA,
via Bonomea 265, I-34132 Trieste, Italy  bI.N.F.N. sezione di Trieste, SISSA,
via Bonomea 265, I-34132 Trieste, Italy  cI.N.F.N. sezione di Torino, via P.
Giuria 1, I-10125 Torino, Italy  dDipartimento di Fisica, “Sapienza”
Università di Roma, Piazzale Aldo Moro 5, 00185, Roma, Italy  eI.F.P.U.,
Institute for Fundamental Physics of the Universe, via Beirut 2, I-34014
Trieste, Italy.

###### Abstract

We compute the probability density distribution of maxima for a scalar random
field in the presence of local non-gaussianities. The physics outcome of this
analysis is the following. If we focus on maxima whose curvature is larger
than a certain threshold for gravitational collapse, our calculations
illustrate how the fraction of the Universe’s mass in the form of primordial
black holes (PBHs) changes in the presence of local non-gaussianities. We find
that previous literature on the subject overestimates, by many orders of
magnitude, the impact of local non-gaussianities on the PBH abundance. We
explain the origin of this discrepancy, and conclude that, in realistic
single-field inflationary models with ultra slow-roll, one can obtain the same
abundance found with the gaussian approximation simply changing the peak
amplitude of the curvature power spectrum by no more than a factor of two. We
comment about the relevance of non-gaussianities for second-order
gravitational waves.

##  I Introduction

The possibility that the totality of dark matter in the Universe consists of
primordial black holes (PBHs) still holds the stage even though almost half-a-
century has passed after the pioneering proposal of Hawking and Carr
Hawking:1971ei ; Carr:1974nx . This is especially true in the mass range
$10^{18}\lesssim M_{\rm PBH}\,[{\rm g}]\lesssim 10^{21}$ in which black holes
are neither too light (otherwise they would have evaporated in the past
through Hawking radiation Carr:2009jm ) or too heavy (otherwise they would
distort space-time in a way that contradicts present bounds from lensing
experiments Niikura:2017zjd ; Katz:2018zrn ).

PBHs could have formed in the very early Universe during the radiation
dominated era.111It is also possible to have PBH formation during matter
domination Harada:2016mhb . The key ingredient that triggers the formation of
a PBH is the presence of an over-fluctuation in the density of the Universe
which, if large enough, gravitationally collapses dragging down any matter
within its horizon, that is the parcel of space around any point reachable at
the speed of light.

The theory of inflation provides an elegant mechanism that explains the origin
of density perturbations in the Universe. In the inflationary picture, space-
time fluctuates quantum mechanically around a background that is expanding
exponentially fast. After the end of inflation, these curvature fluctuations
are transferred to the radiation field, creating slightly overdense and under-
dense regions. It is, therefore, fascinating to ask whether the formation of
PBHs fits in the inflationary picture of structure formation.222This is not
the only option. The formation of PBHs may have been independent of
inflationary physics; PBHs may have been originated from topological defects
formed during symmetry breaking phase transition, for instance from the
collapse of string loops Vilenkin:1981iu ; Hawking:1990tx ; Fort:1993zb . To
answer this question, two (related) aspects need to be addressed.

  * i)

The inflaton dynamics should give rise to a peak in the power spectrum of
curvature perturbations.333 An exception, where no amplification of the power
spectrum is needed, are models where PBHs form from the collapse of domain
walls created during inflation Deng:2016vzb .  This translates into a large
variance for density perturbations that, in turn, enhances the chance to
create overdense regions above the threshold for gravitational collapse.

  * ii)

The abundance of such collapsing regions should be large enough to explain the
totality of dark matter.

In this paper we focus on simple single-field inflationary models. It is known
that in order to fulfil point i) slow-roll conditions must be violated. The
simplest option is to introduce, few $e$-folds before the end of inflation, an
approximate stationary inflection point in the inflaton potential (see refs.
Starobinsky:1992ts ; Ivanov:1994pa ; Saito:2008em  for the earliest proposal
in this direction). When the inflaton, during its classical dynamics, crosses
this region (during the so-called “ultra slow-roll phase”), curvature
perturbations, due to the presence of negative friction, get exponentially
enhanced Leach:2000yw ; Leach:2001zf ; Tsamis:2003px ; Kinney:2005vj ;
Kinney:1997ne .444Alternatively, a parametric amplification of curvature
perturbations could be caused by resonance with oscillations in the sound
speed of their propagation Cai:2018tuh . Another possibility is that, after
the inflationary phase, the inflaton begins to oscillate near the minimum of
the potential and fragments into oscillons which, in turn, lead to copious
production of PBHs Cotner:2018vug . Point ii) is more subtle. Due of their
intrinsic quantum-mechanical origin, the way in which quantum fluctuations
lead to a classical pattern of perturbations can be described only in a
probabilistic sense. Consequently, the computation of the abundance of
collapsing regions requires informations about the statistical distribution of
density perturbations. Most of the time, for simplicity, the gaussian
approximation is assumed. However, the very same fact that slow-roll
conditions are violated as a consequence of point i) suggests that non-
gaussianities may play a relevant role. Refs. Atal:2018neu ; Taoso:2021uvl
indeed find that during an ultra slow-roll phase sizable non-gaussianities of
local type are generated. In the rest of this paper we will dub these non-
gaussianities “primordial” to distinguish them from non-gaussianities that
arise from the non-linear relation between curvature and density
perturbations.

What is the impact of primordial non-gaussianities on the gaussian
approximation when computing the PBH abundance? Ref. Franciolini:2018vbk
addressed this question in the context of threshold statistics. The main
result of ref. Franciolini:2018vbk  is that the abundance of PBHs is
exponentially sensitive to primordial non-gaussianities 555More precisely this
means that in the context of threshold statistics the PBH abundance is given
by eq. (30), where $\mathcal{C}_{n}$ are the $n^{\rm th}$ normalized cumulants
of the non-gaussian distribution.. Based on this result, ref. Atal:2018neu
claims that, in the context of single-field inflationary models which feature
an approximate stationary inflection point, the gaussian approximation is
hardly trustable when computing the PBH abundance.

The goal of this work is to address the same question using a different
computational strategy inspired by peak theory Bardeen:1985tr . More
precisely, we associate regions where the overdensity field takes values above
the threshold for gravitational collapse with spiky local maxima of the
curvature perturbation field, and compute the number density of the latter
using peak theory that we extend to include local non-gaussianities.

Our main conclusion is that the impact of local non-gaussianities on the PBH
abundance is far less important compared to what previously thought. We
confirm that in models for PBH production (at least the class of models that
we are going to consider), local non-gaussianities are sizeable enough to
invalid the use of the the gaussian approximation to estimate their abundance.
However we find that their impact is modest when translated in terms of the
amplitude of curvature power spectrum, namely it is enough to change it by a
factor $\simeq 2$ or smaller to obtain the same PBH abundance predicted by the
gaussian calculation. This shift can be obtained by a small change of the
parameters of the inflationary model.

As a phenomenological application, we consider the impact of local non-
gaussianities on the computation of the amplitude of the induced second-order
gravitational-wave signal. Our main conclusion is that a careful treatment of
non-gaussianities is needed in order to provide a reliable comparison with the
expected experimental sensitivities of future gravitational-wave
interferometers. En route, we discuss the difference between threshold
statistics and peak theory, and we explain under which conditions (and why)
peak theory gives a PBH abundance which is larger than the one computed by
means of threshold statistics.

The structure of this paper is as follows.

  * $\ast$

In section II we introduce the problem and present our solution strategy. This
section is paired with appendix A where we explain in more detail the
cosmological interpretation of all quantities involved.

  * $\ast$

In section III we discuss our main results and we explain the discrepancy with
the previous literature. This section is paired with appendix B-G where we
collect all relevant technical details.

  * $\ast$

We conclude in section IV.

##  II Problem setup and solution strategy

Consider in position space

$\displaystyle
h(\vec{x})=\mathcal{R}(\vec{x})+\alpha\mathcal{R}(\vec{x})^{2}\,,$ (1)

where $\alpha$ is a constant, $\mathcal{R}(\vec{x})$ is a gaussian scalar
random field while $h(\vec{x})$ is non-gaussian because of the presence of the
non-linear term on the right-hand side. In this case, non-gaussianities are
called of local type because for a given $\vec{x}$ at which we evaluate $h$
the amount of non-gaussianity is localized at the same position. We briefly
discuss in appendix A the physical interpretation of eq. (1) and the
limitations of this parametrization of non-gaussianities.

The first observation is that

$\displaystyle\partial_{i}h(\vec{x})=[\partial_{i}\mathcal{R}(\vec{x})][1+2\alpha\mathcal{R}(\vec{x})]\,,$
(2)

meaning that stationary points of $\mathcal{R}$ are also stationary points of
$h$ ($\partial_{i}\mathcal{R}=0$ implies $\partial_{i}h=0$).

What is crucial, however, is that the nature (saddle points, maxima or minima)
of these “shared” stationary points depends on the sign of the factor
$1+2\alpha\mathcal{R}$. Consider the matrix of second derivatives evaluated at
a stationary point $\vec{x}_{\rm st}$ (we will further indicate a minimum with
$\vec{x}_{\rm m}$ and a maximum with $\vec{x}_{\rm M}$). One finds the Hessian
matrix $h_{ij}(\vec{x}_{\rm st})=[\mathcal{R}_{ij}(\vec{x}_{\rm
st})](1+2\alpha\mathcal{R}_{\rm st})$.666We use the short-hand notation
$f_{\rm k}\equiv f(\vec{x}_{\rm k})$, $f_{i}(\vec{x}_{\rm
k})\equiv\partial_{i}f(\vec{x})$ evaluated at $\vec{x}_{\rm k}$ and
$f_{ij}(\vec{x}_{\rm k})\equiv\partial_{ij}f(\vec{x})$ evaluated at
$\vec{x}_{\rm k}$ for some generic function $f$. The flat spatial Laplacian is
$\triangle f(\vec{x})\equiv\sum_{i}f_{ii}(\vec{x})$.

To fix ideas, consider the simple case of two spatial dimensions
$\vec{x}=\\{x,y\\}$ and the case in which the stationary point $\vec{x}_{\rm
m}$ is a minimum of $\mathcal{R}$.

Minima of $\mathcal{R}$ are identified by two conditions. The first one,
$\mathcal{R}_{xx}(\vec{x}_{\rm m})\mathcal{R}_{yy}(\vec{x}_{\rm
m})-\mathcal{R}_{xy}(\vec{x}_{\rm m})^{2}>0$ separates extrema from saddle
points. The second one, $\mathcal{R}_{xx}(\vec{x}_{\rm m})>0$ and
$\mathcal{R}_{yy}(\vec{x}_{\rm m})>0$, separates minima from maxima. Since we
have $h_{xx}(\vec{x}_{\rm m})h_{yy}(\vec{x}_{\rm m})-h_{xy}(\vec{x}_{\rm
m})^{2}=(1+2\alpha\mathcal{R}_{\rm m})^{2}[\mathcal{R}_{xx}(\vec{x}_{\rm
m})\mathcal{R}_{yy}(\vec{x}_{\rm m})-\mathcal{R}_{xy}(\vec{x}_{\rm m})^{2}]$
it is obvious that the condition $\mathcal{R}_{xx}(\vec{x}_{\rm
m})\mathcal{R}_{yy}(\vec{x}_{\rm m})-\mathcal{R}_{xy}(\vec{x}_{\rm m})^{2}>0$
is also satisfied by $h$.

On the contrary, since $h_{xx}(\vec{x}_{\rm m})=(1+2\alpha\mathcal{R}_{\rm
m})\mathcal{R}_{xx}(\vec{x}_{\rm m})$ and $h_{yy}(\vec{x}_{\rm
m})=(1+2\alpha\mathcal{R}_{\rm m})\mathcal{R}_{yy}(\vec{x}_{\rm m})$, it is
possible that a minimum of $\mathcal{R}$ becomes a maximum of $h$ if
$1+2\alpha\mathcal{R}_{\rm m}<0$. Viceversa, a maximum of $\mathcal{R}$ can
become a minimum of $h$.

The argument trivially generalizes to the more realistic case of three spatial
dimensions.

In peak theory, one computes the number density of maxima Bardeen:1985tr . We
are interested in the number density of maxima of the non-gaussian variable
$h$. As argued before, identifying this quantity with the number density of
maxima of $\mathcal{R}$ (based on the observation that $h$ and $\mathcal{R}$
have the same stationary points) is not correct. Let us give a quantitative
argument to support this claim. From the previous discussion, it is clear that
counting the maxima of $\mathcal{R}$ might be not enough. On the contrary, a
simple modification could be the following. One should

  * i)

Count the maxima of $\mathcal{R}$;

  * ii)

Add the minima of $\mathcal{R}$ that, depending on the value of
$(1+2\alpha\mathcal{R}_{\rm m})$, become maxima of $h$;

  * iii)

Subtract the maxima of $\mathcal{R}$ that, depending on the value of
$(1+2\alpha\mathcal{R}_{\rm M})$, become minima of $h$.

Ref. Yoo:2019pma  assumes i). However, the two operations ii) and iii) do not
balance between each others, and a sizable correction to i) will be introduced
if $\alpha$ is large enough. In fig. 1 we show how the number density of
maxima of the gaussian variable $\mathcal{R}$ changes (in percentage) as a
function of $\alpha$ when ii) and iii) are implemented. Schematically, we
compute

$\displaystyle\Delta n_{\rm
max}=\frac{({\rm\\#\,minima\,of\,}\mathcal{R}\to{\rm
maxima\,of\,}h)-({\rm\\#\,maxima\,of\,}\mathcal{R}\to{\rm
minima\,of\,}h)}{{\rm\\#\,maxima\,of\,}\mathcal{R}}\,.$ (3)

We obtain fig. 1 using gaussian peak theory (implementing the results of ref.
Bardeen:1985tr , see appendix B). If we take $\alpha\ll 1$, ii) and iii) do
not alter the estimate of i). However, for sizable $\alpha\gtrsim 0.2$ the
change in the number density of maxima of $\mathcal{R}$ becomes evident.

In situations of cosmological interest, the issue is further complicated by
the fact that we are not really interested in all maxima of $h$ but only in
those which are “spiky enough.” The reason is that the quantity which is
relevant is the density contrast $\delta(\vec{x},t)$ (also dubbed overdensity
field in the following) whose relation with $h(\vec{x})$ (assuming radiation
dominated epoch) reads Harada:2015yda

$\displaystyle\delta(\vec{x},t)=-\frac{4}{9}\left(\frac{1}{aH}\right)^{2}e^{-2h(\vec{x})}\bigg{[}\triangle
h(\vec{x})+\frac{1}{2}h_{i}(\vec{x})h_{i}(\vec{x})\bigg{]}\,,$ (4)

where the time dependence comes from the scale factor $a=a(t)$ and the Hubble
rate $H=H(t)$ while $h$ does not depend on time because eq. (4) assumes
perturbations to be on super-horizon scales. Eq. (4) can be thought as a
Poisson equation in which $h$ plays the role of gravitational potential while
the density contrast can be written more precisely as
$\delta(\vec{x},t)\equiv\delta\rho(\vec{x},t)/\rho_{b}(t)$ where $\rho_{b}(t)$
is the average background radiation energy density and
$\delta\rho(\vec{x},t)=\rho(\vec{x},t)-\rho_{b}(t)$ its perturbation. The
physics-case that is relevant for the present study is the one in which the
density contrast has a peak localized in some region of space that is high
enough to trigger the gravitational collapse into a black hole. If the number
of these peaks above threshold is large enough, these black holes can be part
of dark matter. In the range $10^{18}\lesssim M_{\rm PBH}\,[{\rm g}]\lesssim
10^{21}$, a population of PBHs may account for the totality of dark matter
observed in the Universe today.

Figure 1:  Percentage increase in the number density of maxima of
$\mathcal{R}$ when we ii) add the minima of $\mathcal{R}$ that become maxima
of $h$ and iii) subtract the maxima of $\mathcal{R}$ that become minima of $h$
(see eq. (3)). To make this (illustrative) plot we set $\sigma_{0}=1$ and
$\gamma=3/4$ (see appendix B for definitions).

Consider a peak of the overdensity field, located at some spatial point
$\vec{y}_{\rm pk}$.

$\displaystyle\delta(\vec{y}_{\rm
pk},t)=-\frac{4}{9}\left(\frac{1}{aH}\right)^{2}e^{-2h(\vec{y}_{\rm
pk})}\bigg{[}\triangle h(\vec{y}_{\rm pk})+\frac{1}{2}h_{i}(\vec{y}_{\rm
pk})h_{i}(\vec{y}_{\rm
pk})\bigg{]}\simeq-\frac{4}{9}\left(\frac{1}{aH}\right)^{2}\triangle
h(\vec{y}_{\rm pk})\,,$ (5)

where in the second step we linearized in $h$. We follow here the approach of
refs. Germani:2018jgr ; Musco:2018rwt  in which the linear approximation was
adopted. Since the peak amplitude of the overdensity must be larger than some
critical value $\delta_{c}$, we deduce the condition

$\displaystyle-\triangle h(\vec{y}_{\rm
pk})\gtrsim\frac{9}{4}(aH)^{2}\delta_{c}\,,$ (6)

on the curvature of $h$ at the peak of $\delta$.

If we assume that local maxima of $h$ coincide with peaks of $\delta$ (that is
$\vec{y}_{\rm pk}\simeq\vec{x}_{\rm M}$), then the condition in eq. (6) tells
that only maxima of $h$ which are “spiky enough” contribute to the formation
of black holes. Of course, the assumption that local maxima of $h$ coincide
with peaks of $\delta$ requires some care. Ref. DeLuca:2019qsy  argues, both
analytically and numerically, that this assumption is well justified. However,
ref. DeLuca:2019qsy  only considers the case in which $h$ is gaussian, that
is, $h=\mathcal{R}$ with $\alpha=0$ in our case (but they include the presence
of the non-linearities in eq. (5)). Since stationary points of $\mathcal{R}$
are also stationary points of $h$, we tend to believe that the same conclusion
holds true in the case with $\alpha\neq 0$ but of course this is an important
point that has to be checked explicitly.

All in all, the strategy we shall follow in the course of this work is the
following.

First, we will compute the number density of maxima of the non-gaussian random
field $h$ that are “spiky enough” according to the condition in eq. (6). This
requires a generalization of the work in ref. Bardeen:1985tr  such as to
implement local non-gaussianities. Second, we will check that these maxima are
also peaks of the overdensity field. If this last point will turn out to be
true, our computation of the number density of “spiky enough” maxima of $h$
will provide the abundance of peaks of the overdensity field that are large
enough to form black holes.

##  III Results and discussion

We present in this section the main results of our analysis. In section III.1
we discuss how primordial non-gaussianities of local type alter the abundance
of PBHs. In section III.2 we compare with the existing literature. In section
III.3 we (partially) include the effect of non-linearities in the relation
between curvature and density perturbations.

All technical details are collected in appendix A (where we discuss the origin
of eq. (1) from a cosmological viewpoint), appendix B (where we discuss the
gaussian limit), appendixes C and D (where we discuss how to construct the
non-gaussian part and the approximations that are involved), appendix E (where
we give formulas for computing cumulants of generic order), appendix F (where
we discuss how to compute the threshold value for collapse into black holes)
and appendix G (where we discuss non-linearities).

###  III.1 The abundance of PBHs in the presence of primordial non-
gaussianities of local type

The quantity of central interest is the fraction of the Universe’s mass in the
form of PBHs at the time of their formation. As customary in the literature,
we indicate this quantity with $\beta$. The present-day fractional abundance
of dark matter in the form of PBHs is given by (for a review, see ref.
Sasaki:2018dmp )

$\displaystyle\frac{\Omega_{\rm PBH}}{\Omega_{\rm
DM}}=O(1)\times\left(\frac{\beta}{10^{-16}}\right)\left[\frac{g_{*}(t_{f})}{106.75}\right]^{-1/4}\left(\frac{M_{\rm
PBH}}{10^{18}\,{\rm g}}\right)^{-1/2}\,,$ (7)

where $g_{*}(t_{f})$ is the number of relativistic degrees of freedom at the
time of black hole formation (that we normalize and set to its standard model
value). Eq. (7) is defined modulo an overall $O(1)$ factor whose precise value
depends on the detail of the gravitational collapse that leads to black hole
formation. In this paper we consider $M_{\rm PBH}\simeq 10^{18}$ g;
consequently, as an order-of-magnitude estimate, $\beta\gtrsim 10^{-16}$ is
excluded since it would imply overclosure of the present-day Universe,
$\Omega_{\rm PBH}>\Omega_{\rm DM}$.

We find the following formula

fraction of the Universe’s mass in PBH in the presence of primordial local
non-gaussianities
$\displaystyle\beta\simeq\frac{1}{4\sqrt{2\pi(1-\gamma^{2})}}\,\left[\int_{-\frac{1}{2\alpha\sigma_{0}}}^{\infty}d\bar{\nu}\,\int_{x_{\delta}(\bar{\nu})}^{\infty}dx\,e^{-\bar{\nu}^{2}/2}\,f(x)\,e^{-\frac{(x-x_{*})^{2}}{2(1-\gamma^{2})}}+\int_{-\infty}^{-\frac{1}{2\alpha\sigma_{0}}}d\bar{\nu}\,\int_{-\infty}^{x_{\delta}(\bar{\nu})}dx\,e^{-\bar{\nu}^{2}/2}\,f(x)\,e^{-\frac{(x-x_{*})^{2}}{2(1-\gamma^{2})}}\right]$
(8)

that we derive in detail in appendix B (as far as the gaussian limit is
concerned), appendix C (where we discuss how to construct the non-gaussian
part and the approximations that are involved) and appendix F (where we
discuss how to compute the threshold value for collapse into black holes). In
short:

  * $\ast$

The parameter $\alpha$, already defined in eq. (1), indicates the presence of
local non-gaussianities (of quadratic type). The limit $\alpha\to 0$
reproduces the gaussian result. A more physical interpretation of this
parameter is given in appendix A. In concrete models of inflation which
generate a sizable abundance of dark matter in the form of PBHs (see, for
instance, ref. Ballesteros:2020qam ), we expect $\alpha\simeq[0.24\div 0.61]$
Atal:2018neu ; Taoso:2021uvl .

We derive our result based on peak theory. More precisely, we associate
regions where the overdensity field takes large values with spiky local maxima
of the comoving curvature perturbation, and compute the number density of the
latter using peak theory that we extend to include local non-gaussianities.
Within this approach, eq. (8) represents an original result.

  * $\ast$

The spectral moments $\sigma_{j}^{2}$ are defined by (see eq. (73) and
discussion in appendix B)

$\sigma_{j}^{2}\equiv\int\frac{dk}{k}\mathcal{P}_{\mathcal{R}}(k)\,k^{2j}\,,$
(9)

where $\mathcal{P}_{\mathcal{R}}(k)$ is the dimensionless power spectrum of
the gaussian random field $\mathcal{R}$. The a-dimensional parameter $\gamma$
is defined as $\gamma=\sigma_{1}^{2}/\sigma_{2}\sigma_{0}$ and takes values
$0<\gamma<1$. In this paper we analyze two possible cases. In order to
elucidate some intermediate results of our computational strategy, we use in
appendix B a simple toy-model for the power spectrum given by the log-normal
function (see eq. (89) and related discussion)

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

since in this case the spectral moments can be computed analytically and they
are given by $\sigma_{j}^{2}=A_{g}k_{\star}^{2j}e^{2j^{2}v^{2}}$. The three
parameters $\\{k_{\star},A_{g},v\\}$ in eq. (10) control, respectively, the
position of the peak of the power spectrum, the peak amplitude of the power
spectrum and its width. However, we remark that in single-field inflationary
models the value of $\alpha$ that defines the amount of local non-
gaussianities and the shape of the power spectrum are intimately related, and
in general one can not take $\alpha$ as a free parameter and fix the power
spectrum to a specific functional form like the one introduced in eq. (10). A
more realistic example is the following. Consider the power spectrum defined
by the piecewise function

$\displaystyle{\rm
realistic\,power\,spectrum:}~{}~{}~{}~{}~{}~{}\mathcal{P}_{\mathcal{R}}(k)=\mathcal{P}_{\mathcal{R}}(k_{\star})\times\left\\{\begin{array}[]{ccc}\left(\frac{k}{k_{1}}\right)^{n_{1}}\exp\left[-\frac{\log^{2}(k_{1}/k_{\star})}{2v^{2}}\right]&{\rm
for}&k<k_{1}\\\ \exp\left[-\frac{\log^{2}(k/k_{\star})}{2v^{2}}\right]&{\rm
for}&k_{1}\leqslant k\leqslant k_{2}\\\
\left(\frac{k}{k_{2}}\right)^{n_{2}}\exp\left[-\frac{\log^{2}(k_{2}/k_{\star})}{2v^{2}}\right]&{\rm
for}&k>k_{2}\end{array}\right.$ (14)

with two power-law behaviors for $k\ll k_{\star}$ and $k\gg k_{\star}$ that,
for $k\approx k_{\star}$, are connected by the log-normal function in eq.
(10). In this case, it is possible to show that the spectral index of the
fall-off of the power spectrum after the peak at $k=k_{\star}$ is related to
$\alpha$ by the relation $n_{2}\approx-4\alpha$ that is twice the value of the
Hubble parameter $\eta$ after the end of the ultra slow-roll phase.777This can
be understood as follows. The modes that constitute the fall-off of the power
spectrum after the peak are those for which the horizon-crossing condition
$k=aH$ happens after the end of the ultra slow-roll phase Ballesteros:2020qam
; during this part of the dynamics the power spectrum can be approximated by
means of the conventional slow-roll relation
$\mathcal{P}_{\mathcal{R}}(k)=H^{2}/8\pi^{2}\epsilon$ where the Hubble
parameter $\epsilon$ evolves in time according to $\epsilon(N)\propto
e^{-2\eta_{0}N}$ where $\eta_{0}$ (which is a negative number, $\eta_{0}<0$,
see appendix A) is the value of the Hubble parameter $\eta$ after the end of
the ultra slow-roll phase. We neglect the contribution coming from the time-
evolution of $H$, which is sub-leading. This means that we have
$\mathcal{P}_{\mathcal{R}}(k)\propto e^{2\eta_{0}N}$. From $k=aH$ we have
$dk/k=dN$, and we can convert the $e$-fold time-dependence into a
$k$-dependence, $\mathcal{P}_{\mathcal{R}}(k)\propto
k^{2\eta_{0}}=k^{-4\alpha}$ where we used that $\alpha=-\eta_{0}/2$ (see
appendix A).  In our numerical analysis, therefore, we use the realistic power
spectrum in eq. (14) with the condition $n_{2}=-4\alpha$ for fixed $\alpha$.
This provides the above-mentioned relation between the amount of local non-
gaussianities and the shape of the power spectrum. We remark that this point
is often overlooked in the literature.

Figure 2:  Power spectra for the model in ref. Ballesteros:2020qam  (red, with
label $k^{-1.51}$), ref. Ballesteros:2017fsr  (blue, with label $k^{-2.45}$)
and ref. Dalianis:2018frf  (green, with label $k^{-0.88}$) computed
numerically by solving the Mukhanov-Sasaki equation (left panel). The slope of
the power-law falloff after the peak is $k^{2\eta_{0}}$ where $\eta_{0}$ is
the value of the Hubble parameter $\eta$ (whose evolution, as function of the
$e$-fold time $N$ with $N_{\rm in}$ the beginning of the ultra slow-roll
phase, is shown in the right panel) after the end of the ultra slow-roll
phase; $\eta_{0}$, in turn, is related to the parameter $\alpha$ that controls
the size of local non-gaussianities via $\alpha=-\eta_{0}/2$ (see appendix A).

As a numerical check, we show in fig. 2 (left panel) the power spectra of
comoving curvature perturbations computed numerically for three inflationary
models in which a ultra slow-roll phase takes place. All power spectra are
well described by the analytical ansatz in eq. (14). The slope of the power
spectra after the peak is related to the value of the Hubble parameter $\eta$
(shown in the right panel of the same figure) after the end of the ultra slow-
roll phase which, in turn, controls the size of local non-gaussianities (see
caption and appendix A for more details).

As far as the values of the other parameters in eq. (14) are concerned, we use
$k_{\star}=1.5\times 10^{14}$ Mpc-1, $n_{1}=3.4$, $k_{1}=k_{\star}/5$,
$k_{2}=3k_{\star}/2$ and $v=0.7$; the values of $k_{1,2}$, $v$ and $n_{1}$ are
motivated by a fit of eq. (14) done with respect to the numerical power
spectrum obtained in the context of the explicit models studied in ref.
Ballesteros:2020qam . In particular, notice that the spectral index $n_{1}$
describes the growth of the power spectrum that leads to the formation of the
peak; its value is related to the value of the Hubble parameter $\eta$ during
the ultra slow-roll phase and the duration of the latter. Semi-analytical
arguments (see ref. Byrnes:2018txb ) suggest that $n_{1}<4$.888However, see
ref. Ozsoy:2019lyy  for a special case (derived in the context of non-
attractor inflation) in which the growth of the power spectrum can be steeper,
although for a limited range of $k$. A steeper growth can also be attained in
multi-field inflationary scenarios Fumagalli:2020adf .  The value
$k_{\star}=1.5\times 10^{14}$ Mpc-1 is chosen because it implies $M_{\rm
PBH}\simeq 10^{18}$ (for which $\Omega_{\rm PBH}\simeq\Omega_{\rm DM}$ is
possible, see ref. Taoso:2021uvl ). We consider the peak amplitude of the
power spectrum $\mathcal{P}_{\mathcal{R}}(k_{\star})$ as a free parameter.
Notice that the power spectra that we consider lead to a narrow PBHs mass
function, see appendix F for details. In the following, to analyze the impact
of non-gaussianities, we will focus on the PBH mass around the peak of the
distribution.

  * $\ast$

The function $f(x)$ is given in eq. (124) and the quantity $x_{*}$ in eq.
(119). The function $x_{\delta}(\bar{\nu})$ is defined by the relation (see
eq. (121))

$\displaystyle(1+2\alpha\sigma_{0}\bar{\nu})\,x_{\delta}(\bar{\nu})=\frac{9(a_{m}H_{m})^{2}}{4\sigma_{2}}\delta_{c}\,,$
(15)

where $\delta_{c}=O(1)$ is a threshold value above which a peak of the
overdensity field collapses to form a black hole. The left hand side of eq.
(15) corresponds to the critical curvature of $h$ in eq. (6), see section C.
Formally, eq. (15) depends on time via the comoving Hubble radius $1/aH$. We
evaluate eq. (15) at the time $t_{m}$ when curvature perturbations re-enter
the horizon and become causally connected (see refs. Germani:2018jgr ;
Musco:2018rwt  and appendix F) 999Notice however that eq. (4) is valid on
super-horizon scales. This means that at horizon-crossing additional non
linear effects are present. Recently, these corrections have been considered
in ref. Musco:2020jjb .. In eq. (15) we use the short-hand notation
$a_{m}H_{m}\equiv a(t_{m})H(t_{m})$.

  * $\ast$

An important comment concerns the so-called smoothing procedure. Consider the
case in which one takes a very narrow power spectrum, like the toy-model
introduced in eq. (10) with a small value of $v$, say $v=0.1$. In this case it
is not strictly necessary to introduce a smoothing procedure because the power
spectrum is characterized by a well-defined scale in momentum space,
$k=k_{\star}$. The realistic case introduced in eq. (14), on the contrary,
requires more care. Although the power spectrum peaks at $k=k_{\star}$, the
peak is broadened by the relatively large value of $v$, and it also possesses
a pronounced power-law tail at large $k\gg k_{\star}$. In this situation we
can not blindly apply eq. (8) to compute the PBH abundance because the
spectral moment $\sigma_{2}^{2}$ is formally ultraviolet-divergent unless the
power spectrum decays fast enough, which is however not the case in the
realistic model.101010Of course, any power spectrum generated by the
inflationary dynamics has an intrinsic cut-off set by the smallest scale
(largest $k$) that exits the Hubble horizon before inflation ends. More
precisely, therefore, with the words “ultraviolet-divergent integral” we mean
that $\sigma_{2}^{2}$ is dominated by small scales. The solution to this issue
(discussed in appendix F) is to smooth-out small scales by introducing an
appropriate cut-off. At the operative level, we use, instead of eq. (14), the
power spectrum $\mathcal{P}_{\mathcal{R}}^{\rm
cut}(k)\equiv\mathcal{P}_{\mathcal{R}}(k)\exp(-k^{2}/k_{\rm cut}^{2})$ and we
choose $k_{\rm cut}$ such as to minimize the threshold value in the right-hand
side of eq. (15). Notice that this smoothing procedure is relevant for the
determination both of the threshold value and the spectral moments in eq. (9).

Physically, the fact that the power spectrum in eq. (14) does not possess a
well-defined scale means that the PBHs it generates will be characterized by a
relatively broad mass distribution (rather than sharply peaked at the value
associated to $k_{\star}$). The cut-off procedure described before selects the
scale (and, therefore, the value of the mass $M_{\rm PBH}$) at which the
abundance of PBH will be the largest.

  * $\ast$

We use the linear approximation in eq. (4). We discuss the role of non-
linearities in section III.3.

We now discuss the implications of eq. (8).

Figure 3:  Left panel: Fraction of the Universe’s mass in PBHs at the time of
their formation computed with (solid lines with colors) and without (dashed
black line) non-gaussianities as a function of the peak amplitude of the power
spectrum $\mathcal{P}_{\mathcal{R}}(k_{\star})$. We adopt the realistic model
for the power spectrum introduced in eq. (14) and the abundance is computed
using eq. (8). We show the impact of non-gaussianities for different benchmark
values of $\alpha$. In the hatched region we have $\beta>10^{-16}$ and the
Universe is overclosed (see eq. (7) and related discussion). Right panel: we
consider two spatial dimensions, and compare the exact computation of $\beta$
with the approximation obtained including only the third-order cumulants. We
also show the value of $\beta$ obtained by means of the exponential
approximation in eq. (232); the latter gives an estimate of the abundance off
by many orders of magnitude compared with the actual result.

In the left panel of fig. 3 we show the fraction of Universe’s mass in the
form of PBH computed according to eq. (8) for increasing values of the
parameter $\alpha$ starting from the gaussian case with $\alpha=0$. What
values of $\alpha$ are expected in concrete models? In popular single-field
models of inflation that generate a sizable abundance of dark matter in the
form of PBHs, we find $\alpha\simeq 0.38$ (ref. Ballesteros:2020qam ),
$\alpha\simeq 0.37$ (ref. Ozsoy:2018flq ), $\alpha\simeq 0.30$ (ref.
Cicoli:2018asa ), $\alpha\simeq 0.22$ (ref. Dalianis:2018frf ), $\alpha\simeq
0.61$ (ref. Ballesteros:2017fsr ). The plot shows that including local non-
gaussianities of primordial origin makes the formation of PBHs easier, and the
value of $\mathcal{P}_{\mathcal{R}}(k_{\star})$ required to reproduce the
benchmark abundance $\beta=10^{-16}$ turns out to be smaller than the gaussian
one. We find that the rescaling of the peak amplitude implied by the presence
of local non-gaussianities is modest, a factor a few in realistic models.
Moreover, let us mention that it can be obtained at a price of an even smaller
retuning of the parameters of the inflationary models. Similar results have
been obtained in Taoso:2021uvl  for the calculation of the PBH abundance with
threshold statistics.

In addition to the exact result presented in eq. (8), we also consider a
“perturbative” approach based on a power-series $\alpha$-expansion around the
gaussian distribution. To this end, we work in two (instead of three) spatial
dimensions. This simplifying assumption allows to perform most of the
computations in appendix B and appendix D analytically, and makes possible to
visualize and check numerically a number of intermediate results by means of
simple two-dimensional plots.

Let us first clarify the exact meaning of the word “perturbative.” From a
statistical viewpoint, the $\alpha$-expansion corresponds to an expansion in
cumulants of the joint non-gaussian probability distribution according to the
schematic summarized in table 1 (see appendix B, appendix D and appendix E for
details).

| Gaussian | $O(\alpha)$ | $O(\alpha^{2})$ | $O(\alpha^{3})$ | $O(\alpha^{4})$ | $O(\alpha^{5})$ | $\dots$
---|---|---|---|---|---|---|---
$\alpha=0$
second-order cumulants $C_{2}$ | ✓ | ✗ | ✓ | ✗ | ✗ | ✗ | $\dots$
third-order cumulants $C_{3}$ | ✗ | ✓ | ✗ | ✓ | ✗ | ✗ | $\dots$
fourth-order cumulants $C_{4}$ | ✗ | ✗ | ✓ | ✗ | ✓ | ✗ | $\dots$
fifth-order cumulants $C_{5}$ | ✗ | ✗ | ✗ | ✓ | ✗ | ✓ | $\dots$
$\dots$ | $\dots$ | $\dots$ | $\dots$ | $\dots$ | $\dots$ | $\dots$ | $\dots$
Table 1: We organize the cumulants $C_{n}$ of the joint six-dimensional
probability distribution $P(h,h_{x},h_{y},h_{xx},h_{xy},h_{yy})$ as a series
expansion in $\alpha$. In the case $\alpha=0$, only second-order cumulants are
non-vanishing (eqs. (205-208) with $\alpha=0$), and we reconstruct the
gaussian limit. At order $O(\alpha)$, the leading correction is given by
third-order cumulants (eqs. (213-220)). At order $O(\alpha^{2})$, we include
corrections to the second-order cumulants (eqs. (205-208)) and the leading
pieces in the expression of fourth-order cumulants (eqs. (233-247)).

In the gaussian approximation, all cumulants $C_{n\geqslant 3}$ vanish and the
second-order ones correspond to the entries of the covariance matrix of the
distribution. This result is valid in the limit $\alpha\to 0$. If $\alpha\neq
0$, all cumulants $C_{n\geqslant 3}$ are generated. However, the non-zero
cumulants can be organized in terms of an $\alpha$-expansion as shown in table
1 (see caption, and appendix D and appendix E for details). Crudely speaking,
we have

$\displaystyle C_{n\geqslant 2}\sim O(\alpha^{n-2})+O(\alpha^{n})\,.$ (16)

At order $O(\alpha)$, only the leading part of the third-order cumulants
$C_{3}$ appears.111111More precisely, for each cumulant the expansion is
controlled—considering for simplicity the log-normal power spectrum in eq.
(10)—by the dimensionless parameter $\alpha A_{g}\ll 1$. We can, therefore,
consider an expansion around the gaussian distribution including only the
leading part of the third-order cumulants $C_{3}$. The rationale for this
approximation is twofold.

  * i)

If we compare the perturbative approach with the exact computation (downgraded
in two spatial dimensions) we can estimate the validity of the approximation
in which only the leading part of the third-order cumulant $C_{3}$ is
included.

This exercise is useful for the following reason. In the approach based on
threshold statistics, one usually computes, using the tools of cosmological
perturbation theory, the cumulants in the form of (the connected part of)
correlators of the density perturbation field. In the presence of ultra slow-
roll, however, this computation is not simple (see refs. Atal:2018neu ;
Taoso:2021uvl ), and one typically includes only the leading term in the so-
called bispectrum (that is the three-point correlator). This approximation
precisely corresponds to the one in which only the leading part of the third-
order cumulant $C_{3}$ is included (that is the one proportional to $\alpha$
in eq. (16)). In the right panel of fig. 3 we show the comparison between the
exact computation of $\beta$ and the approximation in which only the leading
part of $C_{3}$ is included; as mentioned before, we work in two spatial
dimensions and we fix $\alpha=0.45$. The comparison shows that truncating the
expansion at the first order in the non-gaussian corrections does not fully
capture the impact of non-gaussianities on $\beta$.

Based on this result, we pose the attention on the fact that a similar
conclusion is likely to be valid also when computing correlators in
cosmological perturbation theory. The local non-gaussianity which is present
at the level of the three-point correlator (computed, for instance, in refs.
Atal:2018neu ; Taoso:2021uvl  and used in ref. Atal:2018neu  to estimate the
impact of non-gaussianities on $\beta$) induces corrections also at higher
orders. In terms of (cosmological) Feynman diagrams, the situation can be
sketched as follows

third-order cumulant

fourth-order cumulant from (cubic interactions)2

pure fourth-order cumulant

end of inflation

(17)

where the square of the third-order local interaction enters at the fourth-
order (as well as at higher ones).121212Of course, higher-order correlators
generated by pure higher-order interactions—for instance, a pure quartic
interaction in the violet vertex of the example above—could also be present
but, as discussed in appendix A, we do not consider them explicitly in this
paper. The comparison shown in the right panel of fig. 3 suggests that,
without including the contributions that third-order local interactions induce
in higher-order cumulants, a precise computation of $\beta$ cannot be claimed.

In addition to this observation there is a second, and by far more
problematic, issue that we shall discuss next.

  * ii)

When computing non-gaussian corrections in the form of a series expansion
around the gaussian distribution, some care must be taken. Seemingly harmless
approximations, indeed, may lead to erroneous conclusions. Consider the
simplified case discussed before in which only the leading part of the third-
order cumulants $C_{3}$ is included. A wrong approximation in the evaluation
of the non-gaussian probability distribution leads to the expression of the
abundance that we report in eq. (232). In this compact expression, the non-
gaussian correction exponentiates and alters the argument of the exponential
function in the gaussian distribution. For illustration, we add in the right
panel of fig. 3 the value of $\beta$ that one gets by means of eq. (232).
Clearly, eq. (232) overestimates the actual value of $\beta$ by many orders of
magnitude. In appendix D we explain why this approximation is wrong and what
is the correct procedure to follow.

The reader may wonder why we are wasting time discussing a wrong result. The
reason is that it rings a bell. Previous studies on the impact of primordial
non-gaussianities found, although in a slightly different statistical context,
that non-gaussian corrections alter exponentially the gaussian value of
$\beta$, analogously to what happens with eq. (232), an expression that we
just branded inaccurate. A closer look at this literature, therefore, is
mandatory. This will be the subject of the next section.

###  III.2 Comparison with the literature

The impact of primordial non-gaussianities on the computation of PBH abundance
was discussed in ref. Franciolini:2018vbk  in the context of the so-called
threshold statistics. Ref. Franciolini:2018vbk  finds that primordial non-
gaussianities may play a very relevant role because they alter the argument of
the exponential function that sets the value of the PBH abundance in the
gaussian case. Ref. Atal:2018neu  applies the results of ref.
Franciolini:2018vbk  to the case of single-field inflationary models which
feature the presence of an approximate stationary inflection point, and
concludes that the gaussian approximation does not give the correct estimate
of the PBHs abundance precisely because the non-gaussian correction
exponentiates and drastically changes the gaussian result.

However, this conclusion clashes with our result. In fact, as we are going to
show, the impact of non-gaussianities on the PBH abundance is much smaller
than what found in ref. Franciolini:2018vbk . To explain the discrepancy, let
us first re-derive the main result of ref. Franciolini:2018vbk  in a
simplified way.

In a nutshell, in the context of threshold statistics one computes, assuming a
gaussian distribution, the probability to find regions where the overdensity
field $\delta$ takes values above a given threshold. It is known that
threshold statistics gives a smaller PBH abundance if compared with peak
theory. Let us briefly explain the origin of this difference. This point is
not crucial to understand the discrepancy between our conclusions and ref.
Franciolini:2018vbk  but it will play an important role later when we will
discuss the relation with gravitational waves. In this section we work
directly with the density contrast $\delta$. Furthermore, we start considering
the gaussian limit. The gaussian probability density distribution of $\delta$
with variance $\sigma_{\delta}$ and spectral moments $\bar{\sigma}_{j}$ (where
$\bar{\sigma}_{0}=\sigma_{\delta}$) is given in the two cases by the
expressions

$\displaystyle P_{\rm threshold}(\delta)$
$\displaystyle=\frac{1}{\sqrt{2\pi}\sigma_{\delta}}e^{-\delta^{2}/2\sigma_{\delta}^{2}}\,,$
(18) $\displaystyle P_{\rm peak}(\delta)$
$\displaystyle=\frac{1}{(2\pi)^{2}\sigma_{\delta}R_{*}^{3}}\underbrace{\bigg{\\{}\int_{0}^{\infty}dx\frac{f(x)}{\sqrt{2\pi(1-\gamma^{2})}}\exp\left[-\frac{(x-\delta\gamma/\sigma_{\delta})^{2}}{2(1-\gamma^{2})}\right]\bigg{\\}}}_{\equiv\,G(\gamma,\delta)}e^{-\delta^{2}/2\sigma_{\delta}^{2}}\equiv\frac{1}{(2\pi)^{2}\sigma_{\delta}R_{*}^{3}}\int_{0}^{\infty}dxP_{\rm
peak}(x,\delta)\,,$ (19)

where the function $f(x)$ is given explicitly in eq. (124). The two quantities
$0\leqslant\gamma\equiv\bar{\sigma}_{1}^{2}/\bar{\sigma}_{2}\sigma_{\delta}<1$
and $R_{*}\equiv\sqrt{3}\bar{\sigma}_{1}/\bar{\sigma}_{2}$ are factors
depending on the power spectrum of $\delta$ (notice in particular that
$\gamma$ is the same defined below eq. (9) but now written in terms of the
spectral moments of $\delta$ instead of $\mathcal{R}$). The variable $x$ is
defined by $x\equiv-\triangle\delta/\bar{\sigma}_{2}$. The two expressions are
similar but there are two differences. First, $P_{\rm peak}(\delta)$ has
dimension of inverse spatial volume (because of the factor $1/R_{*}^{3}$).
This is because $P_{\rm peak}(\delta)$ is defined in peak theory as a number
density of maxima. This implies that the true comparison is between the
a-dimensional quantities $P_{\rm threshold}(\delta)$ and $R_{*}^{3}P_{\rm
peak}(\delta)$. Second, $P_{\rm peak}(\delta)$ contains an extra factor (the
one in curly brackets) with respect to $P_{\rm threshold}(\delta)$. This
factor arises because in peak theory one starts from the ten-dimensional
gaussian joint probability density distribution
$P(\delta,\delta_{i},\delta_{ij})$ of $\delta$ and its first ($\delta_{i}$)
and second ($\delta_{ij}$) spatial derivatives and imposes a number of
conditions that select among the stationary points those that are maxima
Bardeen:1985tr . In threshold statistics, on the contrary, one simply
integrates out all spatial informations by reducing
$P(\delta,\delta_{i},\delta_{ij})$ to the one-dimensional gaussian
distribution in eq. (18). The crucial aspect is that the function
$G(\gamma,\delta)$ in eq. (19) depends on $\delta$ in a way which is
proportional to the parameter $\gamma$. The latter controls the degree of
correlation between $\delta$ and $x$. When $\gamma=0$, the two variables are
completely uncorrelated, and we find that
$G(0,\delta)=(29-6\sqrt{6})/10\sqrt{10\pi}$. In this case, $G(0,\delta)$ does
not depend on $\delta$: peak theory and threshold statistics give the same
qualitative answer (in the sense that the only $\delta$-dependence is encoded
in the gaussian exponential which is in common between the two). When
$\gamma\to 1$, $\delta$ and $x$ are strongly correlated.

Figure 4:  Isocontours of constant probability density distribution $P_{\rm
peak}(x,\delta)$ defined in eq. (19). We use $\nu\equiv\delta/\sigma_{\delta}$
and we set (for illustrative purposes only) the threshold at $\nu_{c}=3$. The
probability above the threshold is obtained integrating $P_{\rm
peak}(x,\delta)$ for $x\in[0,\infty)$ and $\nu\in[\nu_{c},\infty)$. We show
the case with $\gamma=0$ (no correlation between $x$ and $\delta$, left panel)
and $\gamma=0.95$ (strong correlation between $x$ and $\delta$, right panel).

This fact deforms the shape of the probability density distribution of $x$ and
$\delta$ and gives more weight to the region in which $\delta$ crosses the
threshold for collapse. This is illustrated schematically in fig. 4 (using
$\nu\equiv\delta/\sigma_{\delta}$). In the case of strong correlation between
$\delta$ and $x$, therefore, it is well expected that after integrating above
the threshold $\nu>\nu_{c}$ peak theory gives a result which is larger than
threshold statistics.131313Notice that, for the sake of simplicity, we are
implicitly assuming in this example the same threshold value $\nu_{c}$ in the
case of peak theory and threshold statistics in order to highlight the main
difference between the two statistical approaches. In threshold statistics,
this information about the correlation between $\delta$ and $x$ is completely
lost since $\delta$ is treated independently from the spatial configuration.
We remark that the difference between threshold statistics and peak theory is
more and more relevant as we approach the limit $\gamma\to 1$. As already
pointed out, the value of $\gamma$ is dictated by the properties of the power
spectrum (in this case, strictly speaking, the power spectrum of $\delta$
which is however related to the power spectrum of comoving curvature
perturbations) with a very peaked power spectrum that corresponds to the limit
$\gamma\to 1$.

After this digression, we can back to the main point of this section. In the
context of threshold statistics, non-gaussianities are included in the form of
a Gram–Charlier A series around the gaussian ansatz bookStat

$\displaystyle P_{\rm NG}(\delta)$ $\displaystyle=P_{\rm
G}(\delta)\left[1+\sum_{n=3}^{\infty}\frac{1}{n!2^{n/2}}B_{n}(0,0,\mathcal{C}_{3},\dots,\mathcal{C}_{n})H_{n}\left(\frac{\nu}{\sqrt{2}}\right)\right]$
(20)
$\displaystyle=P_{\rm{G}}(\delta)\bigg{[}1+\frac{\mathcal{C}_{3}}{3!2^{3/2}}H_{3}\left(\frac{\nu}{\sqrt{2}}\right)+\frac{\mathcal{C}_{4}}{4!2^{4/2}}H_{4}\left(\frac{\nu}{\sqrt{2}}\right)+\frac{\mathcal{C}_{5}}{5!2^{5/2}}H_{5}\left(\frac{\nu}{\sqrt{2}}\right)+\frac{(10\mathcal{C}_{3}^{2}+\mathcal{C}_{6})}{6!2^{6/2}}H_{6}\left(\frac{\nu}{\sqrt{2}}\right)+\dots\bigg{]}\,,$

where
$P_{\rm{G}}(\delta)=(1/\sqrt{2\pi}\sigma_{\delta})\exp(-\delta^{2}/2\sigma_{\delta}^{2})$
is the gaussian probability density distribution of $\delta$ with variance
$\sigma_{\delta}$ (that is $P_{\rm threshold}(\delta)$ introduced in eq. (18))
and $H_{n}$ are the physicists’ Hermite polynomials.141414We remind that
$\displaystyle\frac{d^{n}}{dx^{n}}e^{-x^{2}/2\sigma^{2}}=\frac{(-1)^{n}}{2^{n/2}\sigma^{n}}e^{-x^{2}/2\sigma^{2}}H_{n}\left(\frac{x}{\sqrt{2}\,\sigma}\right)\,,~{}~{}~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}H_{n}(x)=(-1)^{n}e^{x^{2}}\frac{d^{n}}{dx^{n}}e^{-x^{2}}\,,$ (21)
so that the reader can immediately recognize in eq. (20) the structure of a
derivative expansion.  We define $\nu\equiv\delta/\sigma_{\delta}$ and we
introduce the $n$-th complete exponential Bell polynomial
$B_{n}(x_{1},\dots,x_{n})=\sum_{k=1}^{n}B_{n,k}(x_{1},\dots,x_{n-k+1})$ with
$B_{n,k}$ the partial exponential Bell polynomials. We indicate with
$\mathcal{C}_{n}$ the $n^{\rm th}$ normalized cumulant defined as the
connected part of the $n$-point correlator (evaluated at the same point) of
the overdensity field normalized by the $n^{\rm th}$ power of the standard
deviation

$\displaystyle\mathcal{C}_{n}\equiv\frac{\langle\overbrace{\delta(\vec{x})\dots\delta(\vec{x})}^{n\,{\rm
times}}\rangle_{\rm conn}}{\sigma_{\delta}^{n}}\,.$ (22)

The explicit computation of eq. (22) in the presence of an ultra-slow roll
phase is discussed, in the context of cosmological perturbation theory, in
ref. Taoso:2021uvl  (see also ref. Atal:2018neu ). Notice that the quantity
$\langle\delta(\vec{x})\dots\delta(\vec{x})\rangle$ is dimensionless, and so
is $\mathcal{C}_{n}$. We now integrate $P_{\rm NG}(\delta)$ over some
threshold $\delta_{c}=\nu_{c}\sigma_{\delta}$. We define the abundance
$\beta(\nu_{c})=\int_{\delta_{c}}^{\infty}d\delta\,P_{\rm NG}(\delta)$. We
find151515 We use the property
$\displaystyle\frac{1}{\sqrt{2\pi}\sigma}\int_{\delta_{c}}^{\infty}d\delta
e^{-\delta^{2}/2\sigma^{2}}H_{n}\left(\frac{\delta}{\sqrt{2}\sigma}\right)=\frac{1}{\sqrt{\pi}}e^{-\nu_{c}^{2}/2}H_{n-1}\left(\frac{\nu_{c}}{\sqrt{2}}\right)\,.$
(23)

$\displaystyle\beta(\nu_{c})$ $\displaystyle=\frac{1}{2}{\rm
Erfc\left(\frac{\nu_{c}}{\sqrt{2}}\right)}+\frac{1}{\sqrt{2\pi}\nu_{c}}e^{-\nu_{c}^{2}/2}\sum_{n=3}^{\infty}\frac{\sqrt{2}\nu_{c}}{n!2^{n/2}}B_{n}(0,0,\mathcal{C}_{3},\dots,\mathcal{C}_{n})H_{n-1}\left(\frac{\nu_{c}}{\sqrt{2}}\right)$
(24)
$\displaystyle\simeq\underbrace{\frac{1}{\sqrt{2\pi}\nu_{c}}e^{-\nu_{c}^{2}/2}}_{{\rm
gaussian\,approx}\,\beta_{\rm
G}(\nu_{c})}\underbrace{\left[1+\sum_{n=3}^{\infty}\frac{\sqrt{2}\nu_{c}}{n!2^{n/2}}B_{n}(0,0,\mathcal{C}_{3},\dots,\mathcal{C}_{n})H_{n-1}\left(\frac{\nu_{c}}{\sqrt{2}}\right)\right]}_{\rm
non\,gaussian\,correction}\,,$ (25)

where in the last step we use $(1/2){\rm
Erfc}(x/\sqrt{2})\simeq(1/\sqrt{2\pi}x)e^{-x^{2}/2}$ for $x\gg 1$.

We now adopt the approach suggested in ref. Franciolini:2018vbk  and
approximate

$\displaystyle
H_{n}\left(\frac{\nu_{c}}{\sqrt{2}}\right)=2^{n/2}\nu_{c}^{n}+O(\nu_{c}^{n-2})\,,$
(26)

which seems justified since $\nu_{c}\gg 1$. In this case eq. (25)
becomes161616 The complete exponential Bell polynomial is defined by the
exponential generating function Comtet
$\displaystyle\exp\left(\sum_{j=1}^{\infty}\frac{x_{j}t^{j}}{j!}\right)=\sum_{n=0}^{\infty}B_{n}(x_{1},\dots,x_{n})\frac{t^{n}}{n!}\,,$
(27) and the first few complete Bell polynomials are $B_{0}=1$,
$B_{1}(x_{1})=x_{1}$, $B_{2}(x_{1},x_{2})=x_{1}^{2}+x_{2}$,
$B_{3}(x_{1},x_{2},x_{3})=x_{1}^{3}+3x_{1}x_{2}+x_{3}$. In our case we have
$x_{1}=x_{2}=0$ and the previous definition takes the form
$\displaystyle\exp\left(\sum_{j=3}^{\infty}\frac{x_{j}t^{j}}{j!}\right)=1+\sum_{n=3}^{\infty}B_{n}(0,0,x_{3},\dots,x_{n})\frac{t^{n}}{n!}\,,$
(28) which applies to our case with $t=\nu_{c}$ and $x_{i\geqslant
3}=\mathcal{C}_{i\geqslant 3}$ since the application of eq. (26) to eq. (25)
gives
$\displaystyle\beta(\nu_{c})\simeq\frac{1}{\sqrt{2\pi}\nu_{c}}e^{-\nu_{c}^{2}/2}\left[1+\sum_{n=3}^{\infty}B_{n}(0,0,\mathcal{C}_{3},\dots,\mathcal{C}_{n})\frac{\nu_{c}^{n}}{n!}\right]\,.$
(29) Notice that, compared to our result, ref. Franciolini:2018vbk  finds an
additional factor $(-1)^{n}$ in eq. (30) which, however, does not appear in
our computation.

$\displaystyle\beta(\nu_{c})$
$\displaystyle\simeq\frac{1}{\sqrt{2\pi}\nu_{c}}\exp\left(-\frac{\nu_{c}^{2}}{2}+\sum_{n=3}^{\infty}\frac{\mathcal{C}_{n}\nu_{c}^{n}}{n!}\right)\,.$
(30)

We find, therefore, the main result of ref. Franciolini:2018vbk : The non-
gaussian correction exponentiates, and changes the argument of the exponential
in the gaussian distribution.

Consider the simplified case in which only the third-order normalized cumulant
(a.k.a. skewness) is non-vanishing, $\mathcal{C}_{3}\neq 0$ and
$\mathcal{C}_{n>3}=0$. This is the approximation studied in ref. Atal:2018neu
in which the third-order cumulant is computed in the case of local non-
gaussianities.171717Notice, however, that assuming $\mathcal{C}_{3}\neq 0$ and
$\mathcal{C}_{n>3}=0$ is not fully consistent with the local non-gaussianities
computed in ref. Atal:2018neu . Indeed, local non-gaussianities automatically
generate non-vanishing cumulants of any order. This is evident in our
approach, in which $\alpha\neq 0$ generates a whole tower of non-zero
cumulants. From a more mathematical viewpoint, consider the following theorem.
The sequence $\\{\mu_{n},n=0,1,2,\dots\\}$ corresponds to moments of a non-
negative probability density function if and only if the determinants
$\displaystyle D_{n+1}\equiv{\rm
det}\left(\begin{array}[]{ccccc}\mu_{0}&\mu_{1}&\mu_{2}&\cdots&\mu_{n}\\\
\mu_{1}&\mu_{2}&\mu_{3}&\cdots&\mu_{n+1}\\\
\mu_{2}&\mu_{3}&\mu_{4}&\cdots&\mu_{n+2}\\\
\vdots&\vdots&\vdots&\ddots&\vdots\\\
\mu_{n}&\mu_{n+1}&\mu_{n+2}&\cdots&\mu_{2n}\end{array}\right)\,,~{}~{}~{}~{}~{}~{}{n=0,1,2,\dots}\,,$
(36) are all non-negative moments . Consider the illustrative case with
$\mathcal{C}_{n}=0$ for $n\geqslant 4$. We have $\mu_{0}=1$, $\mu_{1}=0$,
$\mu_{2}=\sigma_{\delta}^{2}$, $\mu_{3}=\langle\delta^{3}\rangle$ and
$\mu_{n}=\sigma_{\delta}^{n}(n-1)!!$ for $n\geqslant 4$ even ($\mu_{n}=0$ for
$n\geqslant 4$ odd). The last condition simply means that cumulants of order
higher than three vanish (and the corresponding moments purely gaussian). The
first determinant $D_{3}>0$ already sets a non-trivial condition on the
skewness, that is $-\sqrt{2}<\mathcal{C}_{3}<\sqrt{2}$. However, there is more
than this. It is indeed trivial to see that higher-order conditions $D_{n}>0$
for $n>3$ impose increasingly strong bound on $\mathcal{C}_{3}$ so that only
$\mathcal{C}_{3}\to 0$ is allowed if all the infinite number of constraints
are implemented. In other words, the theorem above implies that skewness alone
can not consistently parameterize a non-gaussian probability density function.
This result resonates with what we discussed at point ii) in section III.1.
Computing only the three-point correlator of the overdensity field does not
fully describe non-gaussianities but one should at least include the
contributions that the three-point function generates at higher orders.  We
find

$\displaystyle\beta(\nu_{c})\simeq\frac{1}{\sqrt{2\pi}\nu_{c}}\exp\left(-\frac{\nu_{c}^{2}}{2}+\frac{\mathcal{C}_{3}\nu_{c}^{3}}{6}\right)\,,~{}~{}~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}\mathcal{C}_{3}\neq 0~{}~{}~{}{\rm
and}~{}~{}~{}\mathcal{C}_{n>3}=0\,.$ (37)

Ref. Franciolini:2018vbk  and ref. Atal:2018neu  used the above equations to
conclude that the gaussian estimate of the PBH abundance is hardly trustable.
This result is based on the approximation in eq. (26). However, as we shall
now discuss, the applicability of this approximation is not as straightforward
as one may think. Let us critically inspect the issue.

We define, for ease of reading, the quantity
$b_{n}(\nu_{c})\equiv(\sqrt{2}\nu_{c}/n!2^{n/2})B_{n}(0,0,\mathcal{C}_{3},\dots,\mathcal{C}_{n})H_{n-1}(\nu_{c}/\sqrt{2})$
which enters in the non-gaussian correction in eq. (25). As a consequence of
eq. (26), we have (slashed terms are neglected if we apply the approximation
in eq. (26))

$\displaystyle b_{3}$
$\displaystyle=\frac{\mathcal{C}_{3}\nu_{c}^{3}}{6}\bigg{(}1-\cancel{\frac{1}{\nu_{c}^{2}}}\bigg{)}\,,$
(38) $\displaystyle b_{6}$
$\displaystyle=\frac{\mathcal{C}_{3}^{2}\nu_{c}^{6}}{72}\bigg{(}1-\cancel{\frac{10}{\nu_{c}^{2}}+\frac{15}{\nu_{c}^{4}}}\bigg{)}\,,$
(39) $\displaystyle b_{9}$
$\displaystyle=\frac{\mathcal{C}_{3}^{3}\nu_{c}^{9}}{1296}\bigg{(}1-\cancel{\frac{28}{\nu_{c}^{2}}+\frac{210}{\nu_{c}^{4}}-\frac{420}{\nu_{c}^{6}}+\frac{105}{\nu_{c}^{8}}}\bigg{)}\,,$
(40) $\displaystyle b_{12}$ $\displaystyle=\dots$ (41)

and so on. We see, for instance, that we neglected the term
$-7\mathcal{C}_{3}^{3}\nu_{c}^{7}/324$ in eq. (40) but we kept the term
$\mathcal{C}_{3}\nu_{c}^{3}/6$ in eq. (38). This is hardly justifiable given
that $\nu_{c}\gg 1$ and we expect $\mathcal{C}_{3}\sim O(1)$ (see ref.
Taoso:2021uvl  for a careful computation of $\mathcal{C}_{3}$). This is enough
to question the validity of the result based on eq. (26). Furthermore, from
this simple comparison, we also see that the approximation in eq.
(26)—contrary to what naïvely expected—gets worse for larger $\nu_{c}$.

Figure 5:  Left panel. Numerical comparison between i) the gaussian result
$\beta_{\rm G}(\nu_{c})$ defined in eq. (25), ii) the exponential
approximation in eq. (37) and iii) the exact expansion of
$\beta(\nu_{c})=\beta_{\rm G}(\nu_{c})[1+\sum_{n=3}^{n_{\rm
max}}b_{n}(\nu_{c})]$ for increasing values of $n_{\rm max}$. We take
$\mathcal{C}_{3}=0.25$ and we show our results as function of $\nu_{c}$. Right
panel. Same as in the left panel but using the scaling in eq. (45) to draw
$\beta$ as function of $\mathcal{P}_{\mathcal{R}}(k_{\star})$. We take
$\delta_{c}=1.2$.

In particular, in this specific case, we can rewrite the sum (24) by using the
recurrence property of the Bell polynomials:

$B_{n+1}\left(x_{1},\dots,x_{n+1}\right)=\sum_{k=0}^{n}\binom{n}{k}B_{n-k}\left(x_{1},\dots,x_{n-k}\right)x_{k+1}.$
(42)

In our particular case $x_{k+1}=\mathcal{C}_{3}\,\delta_{k+1,3}$, therefore
$B_{3(k+1)}(0,0,\mathcal{C}_{3},0,\dots,0)=\frac{1}{2}(3k+1)(3k+2)\mathcal{C}_{3}B_{3k}(0,0,\mathcal{C}_{3},0,\dots,0)$.
This can be easily recasted in the form:

$B_{n}\left(0,0,\mathcal{C}_{3},0,\dots,0\right)=\begin{cases}&3\left(\frac{\mathcal{C}_{3}}{6}\right)^{k}\frac{(3k-1)!}{(k-1)!},\quad\textrm{if}\quad
n=3k,\\\ &0,\quad\textrm{if}\quad n\neq 3k.\end{cases}$ (43)

Using this identity inside the eq. (24) we get:

$\beta(\nu_{c})=\frac{1}{2}\mathrm{Erfc}\left(\frac{\nu_{c}}{\sqrt{2}}\right)+\frac{1}{\sqrt{\pi}}e^{-\nu_{c}^{2}/2}\sum_{k=1}^{\infty}\frac{1}{k!}\left(\frac{\mathcal{C}_{3}}{12\sqrt{2}}\right)^{k}H_{3k-1}\left(\frac{\nu_{c}}{\sqrt{2}}\right).$
(44)

To show more explicitly the error that one makes by taking the exponential
approximation in eq. (37) we plot in the left panel of fig. 5 the comparison
between i) the gaussian result $\beta_{\rm G}(\nu_{c})$ defined in eq. (25),
ii) the exponential approximation in eq. (37) and iii) the exact expansion of
$\beta(\nu_{c})=\beta_{\rm G}(\nu_{c})[1+\sum_{n=3}^{n_{\rm
max}}b_{n}(\nu_{c})]$ for increasing values of $n_{\rm max}$181818One can
argue that the power series $\sum_{k=1}^{\infty}\frac{z^{k}}{k!}H_{3k-1}(x)$
contained in eq. (44) is not convergent. Nevertheless it should be regarded as
an asymptotic series that gives a trustable value after an optimal truncation.
Moreover it is Borel summable, and its value can be computed numerically
evalutating the Borel sum
$\int_{0}^{\infty}e^{-t}\left[\sum_{k=1}^{\infty}\frac{(tz)^{k}}{(k!)^{2}}H_{3k-1}(x)\right]dt$.
This check has been done and the value obtained agrees with the one got by the
optimal truncation.. We take $\mathcal{C}_{3}=0.25$ and we set
$\mathcal{C}_{n>3}=0$. The approximation gives a wrong estimate of the actual
magnitude of the non-gaussianities. As expected, the exponential approximation
diverges from the actual result for larger values of $\nu_{c}$.

There is another point that is worth emphasizing. Taking $\mathcal{C}_{3}$
fixed and changing only $\nu_{c}$, as done in the left panel of fig. 5 is not
fully consistent since both $\mathcal{C}_{3}$ and $\nu_{c}$ are functions of
the power spectrum $\mathcal{P}_{\mathcal{R}}$. More explicitly, we have the
scaling $\nu_{c}\propto\mathcal{P}_{\mathcal{R}}(k_{\star})^{-1/2}$ and
$\mathcal{C}_{3}\propto\mathcal{P}_{\mathcal{R}}(k_{\star})^{1/2}$ (see
appendix A). This means that reducing $\mathcal{P}_{\mathcal{R}}(k_{\star})$
to decrease $\mathcal{C}_{3}$ does not improve on the applicability of the
exponential approximation since $\nu_{c}$, in turn, increases. To better
visualize this point let us consider the following benchmark scalings

$\displaystyle\sigma_{\delta}=0.35\left[\frac{\mathcal{P}_{\mathcal{R}}(k_{\star})}{0.05}\right]^{1/2}\,,~{}~{}~{}~{}~{}\mathcal{C}_{3}=0.8\left[\frac{\mathcal{P}_{\mathcal{R}}(k_{\star})}{0.05}\right]^{1/2}\,,$
(45)

and take $\delta_{c}=1.2$. We can now redo the comparison we did before but
now as function of $\mathcal{P}_{\mathcal{R}}(k_{\star})$. The result is shown
in the right panel of fig. 5. We conclude again that the exponential
approximation overestimates the impact of local non-gaussianities on the PBH
abundance by many orders of magnitude compared with the actual result. Using
the exponential approximation, one would wrongly conclude that, in order to
fit the reference value $\beta\simeq 10^{-16}$, an order-of-magnitude decrease
in the peak amplitude of the power spectrum is needed compared to the gaussian
result.

A precise determination of the value of $\mathcal{P}_{\mathcal{R}}(k_{\star})$
that is needed to get $\Omega_{\rm DM}\simeq\Omega_{\rm PBH}$ has important
phenomenological implications. Let us discuss a specific example. PBHs with
mass $M_{\rm PBH}\simeq 10^{18}$ g have approximatively the size of the atomic
nucleus (remember that $1\,M_{\odot}\simeq 1.48\,{\rm km}$ in the Planck unit
system). Testing the nature of such small objects in the form of dark matter
is extremely challenging. An interesting prospect is the following. The peak
in the power spectrum of scalar perturbations that is responsible for the
formation of PBHs also generates (as a second-order effect) a gravitational
wave signal. The position of the peak amplitude of the power spectrum of
curvature perturbations, the peak height in the mass distribution of PBHs and
the frequency of the peak of the induced gravitational wave signal
($f_{\star}$ in the following) are related by the approximate relation

$\displaystyle\left(\frac{M_{\rm PBH}}{10^{17}\,{\rm
g}}\right)^{-1/2}\simeq\frac{k_{\star}}{2\times 10^{14}\,{\rm
Mpc}^{-1}}\simeq\frac{f_{\star}}{0.3\,{\rm Hz}}\,,$ (46)

meaning that a peak amplitude around $k_{\star}=10^{14}$ Mpc-1 corresponds to
a frequency range detectable by future gravitational-wave interferometers like
LISA, DECIGO and MAGIS-100. This is, in principle, a powerful probe. By means
of the condition $\Omega_{\rm DM}\simeq\Omega_{\rm PBH}$, one can fix the
parameters of a given inflationary model that produces a sizable abundance of
dark matter in the form of PBHs; in turn, this “predicts” the induced signal
of gravity waves since the latter is completely determined once the former
step is taken. This is because the current energy density of gravitational
waves induced by scalar perturbations is given by

$\displaystyle\Omega_{\rm
GW}=\frac{c_{g}\Omega_{r}}{36}\int_{0}^{\frac{1}{\sqrt{3}}}dt\int_{\frac{1}{\sqrt{3}}}^{\infty}ds\left[\frac{(t^{2}-1/3)(s^{2}-1/3)}{t^{2}-s^{2}}\right]^{2}\left[\mathcal{I}_{c}(t,s)^{2}+\mathcal{I}_{s}(t,s)^{2}\right]\mathcal{P}_{\mathcal{R}}\left[\frac{k\sqrt{3}}{2}(s+t)\right]\mathcal{P}_{\mathcal{R}}\left[\frac{k\sqrt{3}}{2}(s-t)\right]\,,$
(47)

where $c_{g}\approx 0.4$, $\Omega_{r}$ is the current energy density of
radiation and $\mathcal{I}_{c}$ and $\mathcal{I}_{s}$ are two functions that
can be computed analytically (see, for instance, ref. Espinosa:2018eve  and
references therein); $\Omega_{\rm GW}$ is completely fixed once
$\mathcal{P}_{\mathcal{R}}(k)$ is given. In most of the analysis the power
spectrum $\mathcal{P}_{\mathcal{R}}(k)$ is fixed to a specific functional
form, and the abundance of PBHs computed assuming gaussian statistics. Typical
functional forms for $\mathcal{P}_{\mathcal{R}}(k)$ used in the literature are
the delta-function power spectrum (see, e.g., ref. Yuan:2019udt ), the log-
normal power spectrum (see, e.g., refs. Kapadia:2020pir ; Inomata:2018epa )
and the more realistic broken power-law in eq. (14) (see, e.g., ref.
Bhaumik:2020dor ); all these studies assume gaussian statistics for PBH
formation.

The first remark that we make is that, as already anticipated, taking a
specific functional form for $\mathcal{P}_{\mathcal{R}}(k)$ and assuming
gaussian statistics can be conceptually wrong; for instance, for single-field
inflationary models of phenomenological relevance, eq. (14) with a fixed value
of $n_{2}=-4\alpha$ implies the presence of local non-gaussianities, and the
value of $\beta$ should be computed accordingly. The second remark is that if
one takes the exponential approximation in eq. (37) (considering for example a
situation like the one given in eq. (45)) then one would incorrectly conclude
that primordial non-gaussianities imply a two orders-of-magnitude suppression
of the gravitational wave signal, since the latter is proportional to
$\mathcal{P}_{\mathcal{R}}(k_{\star})^{2}$. In reality (see fig. 5),
primordial non-gaussianities reduce the value of
$\mathcal{P}_{\mathcal{R}}(k_{\star})$ by no more than a factor of 2 compared
to the gaussian result. We will come back in more detail to this point in the
next section. Before that, let us mention some other previous analysis on the
impact of primordial non-gaussianities on the PBHs phenomenology. Ref.
Byrnes:2012yx  approximates the calculation of the PBH abundance working with
the comoving curvature perturbation instead of the density. This simplified
method is discussed in appendix D, see eq. (190). Ref. Shandera:2012ke
describes the non-gaussian probability density distribution using the
Edgeworth expansion, analogous to the Gram-Charlier series in eq. (20) but
with a different truncation and ordering of terms. The higher order cumulants
$\mathcal{C}_{n}$ are parametrized in terms of the third order one
$\mathcal{C}_{3}$ assuming scaling relations compatible with a certain ansatz
for non-gaussianties (one of the two cases considered corresponds to local
non-gaussianties). Ref. Young:2015cyn  considers several templates for the
bispectrum, which are then used to compute the non-gaussian smoothed density
field. The latter quantity is simulated numerically on a grid, and the PBHs
abundance is computed counting the number of grid points above the threshold
for collapse. Since the formation of the PBHs is a rare event, the limitation
of the computing power forces ref. Young:2015cyn  to focus an on large value
of the PBH abundance, $\beta=10^{-4},$ order of magnitudes above the realistic
ones. As evident, the methods and the statistical contexts of these analysis
are different from ours. For instance, all these works have been performed in
the context of threshold statistics, instead of peak theory. Nevertheless, at
a qualitative level, we find that their results are in agreement with ours,
namely primordial non-gaussianities as those considered here increases the PBH
abundance, but their effect can be compensated by a decrease of
$\mathcal{P}_{\mathcal{R}}(k_{\star})$ by a modest factor.

Finally, we shall mention that primordial non-gaussianities are constrained by
the non-observations of isocurvature modes at CMB scales Tada:2015noa ;
Young:2015kda . These constraints rely on the influence of long wavelength
fluctuations on the formation of PBHs. The role of such long modes is still
under discussion in the literature, e.g.Passaglia:2018ixg ; Suyama:2021adn ;
Matarrese:2020why .

###  III.3 Non gaussianities from the non-linear relation between density and
curvature perturbations.

In the previous sections we have assumed a linear relation between the density
and curvature perturbations. This is because our main focus was the study of
the non-gaussianities of primordial origin. Now we are going to include also
the effect of the non-linearities in eq. (4). In our approach, they enter in
two ways. First, they modify the relation between the overdensity field and
the curvature perturbation; this modification translates into a sort of
“renormalization” of the value of $\delta_{c}$ whose net effect is to make the
production of black holes harder.

Importantly, non-linearities also change the shape of the profile of the
overdensity peaks which eventually collapse into black holes. Since the
threshold value depends on the shape of the overdensity, non-linearities also
affect the way in which $\delta_{c}$ is computed. Ref. Kehagias:2019eil
concludes that $\delta_{c}$ is robust against non-linearities since its value
changes at the percent level. However, ref. Kehagias:2019eil  only considers
the idealized case of a $\delta$-function power spectrum. Deviations from the
$\delta$-function limit are considered in refs. Yoo:2018kvb ; Germani:2019zez
; Yoo:2020dkz ; Kawasaki:2019mbl  using different procedures whose consistency
with each other is not yet fully understood. We do not aim to clarify this
issue in the present work; however, we point out that refs. Yoo:2018kvb ;
Germani:2019zez ; Yoo:2020dkz ; Kawasaki:2019mbl  assume that the curvature
perturbation field is gaussian while, at least for the class of models which
are relevant for the present analysis, primordial non-gaussianities should be
included as well.

A comprehensive study of the interplay between non-gaussianities of primordial
origin and non-linearities in eq. (4) is left for future work. Here we shall
focus on the first effect above, and assume that the threshold $\delta_{c}$ is
not significantly affected by the presence of non-linearities. Under this
assumption, in appendix G we show how to compute the PBH abundance accounting
for primordial non-gaussianities of local type and non-linearities between
density and curvature perturbations. We find that one can still use eq. (8),
but replacing the threshold $x_{\delta}(\bar{\nu})$ with $x^{\rm
NL}_{\delta}(\bar{\nu}),$ defined as

$\displaystyle x_{\delta}^{\rm
NL}(\bar{\nu})\equiv\frac{9(a_{m}H_{m})^{2}}{4\sigma_{2}}\frac{\delta_{c}}{1+2\alpha\sigma_{0}\bar{\nu}}\exp\left\\{2\sigma_{0}[\bar{\nu}+\alpha\sigma_{0}(\bar{\nu}^{2}-1)]\right\\}\,.$
(48)

Figure 6:  Left panel. Same as in the left panel of fig. 3; we compare the
gaussian result (dashed black line) with i) the computation that includes non-
linearities between curvature and density perturbations (dot-dashed line,
label “NL”) and ii) the computation that includes both primordial non-
gaussianities and non-linearities (solid clored lines, label “PNG+NL”). Right
panel. Fraction of the energy density in gravitational waves relative to the
critical energy density of the Universe as a function of the frequency (see
eq. (47)). We compare the gaussian result (dashed black line) with the
computation that includes non-linearities between curvature and density
perturbations and the computation that includes both primordial non-
gaussianities and non-linearities. We also consider for illustration the
comparison with the exponential approximation (dot-dashed line, see the right
panel of fig. 3). We take $\alpha=0.45$. We superimpose the signal on the
expected sensitivity curves of the future gravitational wave detectors LISA
(assuming the C1 configuration Caprini:2015zlo ), DECIGO Yagi:2011wg ,
MAGIS-100 Coleman:2018ozp  and the Einstein Telescope Maggiore:2019uih . We
also show present bound LIGOScientific:2019vic  and future prospect fur  for
the aLIGO experiment. The diagonal gray band is the stochastic gravitational
wave background from binary black holes (BBH) and binary neutron stars (BNS),
which has been computed following ref. Chen:2018rzo .

  * $\ast$

Consider first the gaussian case $\alpha\to 0$. We have
$4\sigma_{2}x_{\delta}^{\rm
NL}(\bar{\nu})|_{\alpha=0}=9(a_{m}H_{m})^{2}\delta_{c}e^{2\sigma_{0}\bar{\nu}}$
to be compared with its linear limit
$4\sigma_{2}x_{\delta}(\bar{\nu})|_{\alpha=0}=9(a_{m}H_{m})^{2}\delta_{c}$.
The main difference, as already emphasized in ref. DeLuca:2019qsy , is that
non-linearities “renormalize” the value of $\delta_{c}$ by means of the
exponential factor $e^{2\sigma_{0}\bar{\nu}}$. This means that when
considering maxima of the curvature field their height, in addition to their
curvature, becomes important at the non-linear level. We remark that these two
quantities are statistically correlated, and the amount of correlation is
controlled by the power spectrum via the a-dimensional parameter $\gamma$ (see
appendix B). If the power spectrum is very narrow, maxima with large curvature
are likely to have also large $\bar{\nu}$; consequently, the threshold value
$\delta_{c}$ effectively increases because of the factor
$e^{2\sigma_{0}\bar{\nu}}$. This effect decreases the abundance with compared
to the gaussian case. However, if the power spectrum is not very narrow, the
correlation mentioned before is less accentuated, and maxima with large
curvature are not necessarily characterized by large $\bar{\nu}$; in this
case, $\delta_{c}$ increases by a smaller factor and the reduction of the
abundance turns out to be less important. Our realistic power spectrum in eq.
(14) belongs to the last case. In the left panel of fig. 6 we compare the
gaussian computation of $\beta$ in the case without (dashed black line) and
with (dot-dashed black line, label “NL”) non-linearities. As expected, $\beta$
decreases but the change in the peak amplitude of the power spectrum that is
needed to reproduce the reference value $\beta\simeq 10^{-16}$ is less than a
factor of 2 (while very narrow power spectra may increase the peak amplitude
of the power spectrum up to one order of magnitude DeLuca:2019qsy ).

  * $\ast$

We take $\alpha\neq 0$. As discussed in section III.1, primordial non-
gaussianities modify the statistics of the overdensity peaks. More peaks form,
and, at constants threshold for collapse, the abundance of PBHs increases.
Non-linearities, however, tend to enhance the value of $\delta_{c}$ in the
sense discussed in the previous paragraph, and this gives the opposite effect
of reducing the abundance. Moreover, when both non-linearities and primordial
non-gaussianities are present the renormalization of $\delta_{c}$ also depends
on the value of $\alpha$, as shown in eq. (48). We have, therefore, two
competing effects: More overdensity peaks form but their collapse becomes less
probable. We find that the result is a net decrease of the PBH abundance. This
is shown in the left panel of fig. 6 (solid color lines, label “PNG+NL”).
Quantitatively, in order to fit the reference value $\beta\simeq 10^{-16}$ the
combination of primordial non-gaussianities and non-linearities require a peak
amplitude of the power spectrum larger compared to its gaussian value.
However, we find that $\mathcal{P}_{\mathcal{R}}(k_{\star})$ only increases by
a factor which is at most 2 in all the realistic cases analyzed in this paper.

Before concluding, let us explain the reason why it is important to have
control on the computation of the PBH abundance at the level of the effects
discussed in the present work.

The only condition $\Omega_{\rm DM}\simeq\Omega_{\rm PBH}$ (which is
$\beta\simeq 10^{-16}$ under the assumptions stated at the beginning of this
section) is not enough to make the presence of non-gaussianities
“phenomenologically observable”. With this term in quotation marks we mean
that the presence (or absence) of non-gaussianities can be simply reabsorbed
in a re-tuning of $\mathcal{P}_{\mathcal{R}}(k_{\star})$ (which is equivalent
to a very modest re-tuning of model parameters, see ref. Taoso:2021uvl ) if
the only goal is to reproduce $\beta\simeq 10^{-16}$.

What changes the rules of the game is the fact that, once
$\mathcal{P}_{\mathcal{R}}(k_{\star})$ is fixed, the amplitude of the induced
gravity wave signal in eq. (47) is also predicted. Importantly, this implies
that partial or inaccurate treatments of non-gaussianities will result in
different values of $\mathcal{P}_{\mathcal{R}}(k_{\star})$ (required by
imposing the same condition $\beta\simeq 10^{-16}$) and, consequently,
different gravitational wave signals. The (inaccurate) exponential
approximation in eq. (30) provides a paradigmatic example. As discussed, if
taken at face value the exponential approximation in eq. (30) may imply a one
order-of-magnitude reduction in the value of
$\mathcal{P}_{\mathcal{R}}(k_{\star})$ which gives $\beta\simeq 10^{-16}$
(see, e.g., fig. 3 or fig. 5). Consequently, the amplitude of the induced
gravity wave signal in eq. (47) would be suppressed by two orders-of-
magnitude.

In the right panel of fig. 6 we compare the value of $h^{2}\Omega_{\rm GW}$
that one gets using peak theory in the gaussian limit (dashed black line) and
the exponential approximation (dot-dashed black line; we take $\alpha=0.45$,
and consider non-gaussianities at the level of the leading part of third-order
cumulants). This example is tailored to show that in this case the
gravitational wave signal would drop below the expected stochastic background
due to mergers of astrophysical black holes and neutron stars (the diagonal
gray band). In the same plot we also show the value of $h^{2}\Omega_{\rm GW}$
that we get in the presence of primordial non-gaussianities (label “PNG” this
time computed according to eq. (8)) and in the presence of both primordial
non-gaussianities and non-linearities (label “PNG+NL” computed according to
our discussion in section III.3). As expected, primordial non-gaussianities
(the combination of primordial non-gaussianities and non-linearities) decrease
(increases) the amplitude of the gravitational wave signal compared to the
gaussian result but the net effect in both cases does not exceed the order of
magnitude. As an order-of-magnitude approximation, therefore, the gaussian
limit can be considered trustable (while using the exponential approximation
would hide the signal below the stochastic background).

There is one last point that is worth emphasizing. In the gaussian
approximation, peak theory gives a PBH abundance which is systematically
larger than the one provided by threshold statistics (this is in agreement
with the expectation of ref. Young:2014ana  and discussed in section III.2).
The amplitude of the gravitational wave signal in fig. 6 is indeed smaller
than the one computed with threshold statistics (see ref. Ballesteros:2020qam
) and dangerously closer to the expected stochastic background. For this
reason, we argue that keeping under control the impact of non-gaussianities
could play an important part in the interpretation of possible future
detection.

##  IV Conclusions and outlook

We conclude summarizing the main results and novelties of our work.

  * $\circ$

We have studied the impact of primordial non-gaussianities of local type on
the abundance of PBHs. For this purpose we have focused on the maxima of the
the comoving curvature perturbations. We have shown that those peaks “spiky
enough”, i.e. with a curvature larger than a certain threshold, are good
proxies for the maxima of the density field which undergo gravitational
collapse, and produce PBHs. Exploiting this observation, we have obtained that
the cosmological abundance of PBHs can be computed using eq. (8). Our
calculations extend the gaussian peak theory formalism Bardeen:1985tr  to
include the effect of local non-gaussianities.

We have examined the consequences of our results for models of single-field
inflation with an ultra slow-roll phase. These scenarios have been
investigated for the production of PBHs. In this context, the impact of
primordial non-gaussianities is not dramatic. In fact, the desired PBH
abundance can be obtained with an amplitude of the power spectrum of curvature
perturbations which is only a factor $\lesssim 2$ smaller from the value
inferred with the gaussian approximation.

  * $\circ$

In parallel to the exact result presented in eq. (8), we have developed a
computational strategy that approximates the effect of primordial non-
gaussianities in the form of an expansion in cumulants around the gaussian
result. In order to have full analytical control, we worked out this part in
two (instead of three) spatial dimensions. The presence of local non-
gaussianities affect cumulants at any order, and we have shown what is their
structure. Furthermore, we have shown that the only inclusion of third-order
cumulants is not sufficient to fully capture the effect of local non-
gaussianities.

  * $\circ$

Previous works have studied primordial non-gaussianities in the context of
threshold statistics, finding that non-gaussianities exponentially affect the
PBH abundance as given by eq. (30). We have found that that result is not
correct and largely overestimate the PBH abundance.

  * $\circ$

In addition to non-gaussianities of primordial origin, we have also considered
non-gaussianities that originate from the non-linear relation between
curvature and density perturbations. In this case we have two competing
effects. On the one hand, primordial non-gaussianities alter the statistics of
peaks enhancing the abundance of PBHs as described by eq. (8). On the other
one, primordial non-gaussianities and non-linearities change the condition for
collapse. At the linear level and without accounting for primordial non-
gaussianities, the condition for collapse only involves the laplacian of the
curvature perturbation at the position of the peak of the overdensity field
while in the presence of non-gaussianities it takes the form of eq. (48). This
alteration makes the formation of PBHs harder. When combining the two effects,
we find that non-gaussianities suppress PBH production but the desired PBH
abundance can be obtained with an amplitude of the power spectrum of curvature
perturbations which is only a factor $\lesssim 2$ larger from the value
inferred with the gaussian approximation.

  * $\circ$

En route, we discuss the difference between threshold statistics and peak
theory, and we explain under which conditions (and why) peak theory gives a
PBH abundance which is larger than the one computed by means of threshold
statistics.

Finally, let us mention that our treatment of non-gaussianities can be further
improved. We have computed the threshold for collapse without including the
effect of non-gaussianities (both primordial and non-linear) on the shape of
the collapsing peak (technically speaking, in eq. (295) we used the averaged
density profile—while in principle peak theory also gives informations about
the shape of the peak—and the linear approximation in eq. (294))
191919Formalisms which incorporate the variation of the shape and size of the
peaks have been proposed in refs. Germani:2019zez ; Suyama:2019npc ;
Young:2020xmk  considering gaussian curvature perturbations.. Including these
effects could give an additional modification of $\delta_{c}$. In light of the
motivations discussed before, it would be important to keep under control
these additional effects. We will address this issue in a forthcoming work
focused on the interplay between the computation of the black hole abundance
and the induced gravitational wave signal.

###### Acknowledgements.

The research of A.U. is supported in part by the MIUR under contract 2017
FMJFMW (“New Avenues in Strong Dynamics,” PRIN 2017) and by the INFN grant
“SESAMO – SinergiE di SApore e Materia Oscura.” M.T. acknowledges support from
the INFN grant “LINDARK,” the research grant “The Dark Universe: A Synergic
Multimessenger Approach No. 2017X7X85” funded by MIUR, and the project
“Theoretical Astroparticle Physics (TAsP)” funded by the INFN.

##  Appendix A Scalar perturbations during inflation in the presence of ultra
slow-roll

In this paper we focus on single-field models of inflation. We indicate with
$\phi$ the canonically normalized inflaton field and with $U=U(\phi)$ its
potential. The dynamics of $\phi$ can be obtained solving the equation of
motion

$\displaystyle\frac{d^{2}\phi}{dN^{2}}+3\frac{d\phi}{dN}-\frac{1}{2}\left(\frac{d\phi}{dN}\right)^{3}+\left[3-\frac{1}{2}\left(\frac{d\phi}{dN}\right)^{2}\right]\frac{d\log
U}{d\phi}=0\,,$ (49)

with slow-roll initial conditions; $N$ indicates the number of $e$-folds
defined by $dN=Hdt$, where $H\equiv\dot{a}/a$ is the Hubble rate, $a$ the
scale factor of the Friedmann-Lemaître-Robertson-Walker metric, and $t$ the
cosmic time (with $\dot{}\equiv~{}d/dt$). We will also use in the following
the conformal time $\tau$ defined by means of $dt/d\tau=a$ or, equivalently,
$dN/d\tau=aH$ (notice that, in the limit in which $H$ is constant, we have the
relation $\tau=-1/aH$; the conformal time, therefore, is negative, and late
times towards the end of inflation can be formally identify with the limit
$\tau\to 0^{-}$).

Scalar fluctuations can be efficiently described in terms of the Mukhanov-
Sasaki field variable $u(\tau,\vec{x})$ which is a gauge-invariant combination
of both fluctuations of the inflaton field and scalar fluctuations of the
background Friedmann-Lemaître-Robertson-Walker geometry. The Mukhanov-Sasaki
field variable solves the differential equation (for a review, see ref.
Riotto:2002yw )

$\displaystyle\frac{d^{2}u}{d\tau^{2}}$
$\displaystyle=\left(\triangle+\frac{1}{z}\frac{d^{2}z}{d\tau^{2}}\right)u\,,$
(50) $\displaystyle\frac{1}{z}\frac{d^{2}z}{d\tau^{2}}$
$\displaystyle=a^{2}H^{2}\left[(1+\epsilon-\eta)(2-\eta)+\frac{1}{aH}\left(\frac{d\epsilon}{d\tau}-\frac{d\eta}{d\tau}\right)\right]\,,$
(51)

with $z\equiv(1/H)(d\phi/d\tau)$ and $\triangle$ the laplacian acting on
spatial coordinates. The Hubble parameters are defined by
$\epsilon\equiv-\dot{H}/H^{2}$ and $\eta\equiv-\ddot{H}/2H\dot{H}$.

The field $u$ can be quantized by defining the operator

$\displaystyle\hat{u}(\tau,\vec{x})=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\left[u_{k}(\tau)a_{\vec{k}}e^{+i\vec{k}\cdot\vec{x}}+u_{k}^{*}(\tau)a^{{\dagger}}_{\vec{k}}e^{-i\vec{k}\cdot\vec{x}}\right]\,,$
(52)

with the annihilation and creation operators that satisfy the commutation
relations of bosonic fields

$[a_{\vec{k}},a_{\vec{k}^{\prime}}]=[a_{\vec{k}}^{{\dagger}},a_{\vec{k}^{\prime}}^{{\dagger}}]=0\,,~{}~{}~{}~{}~{}~{}~{}[a_{\vec{k}},a_{\vec{k}^{\prime}}^{{\dagger}}]=(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})\,,~{}~{}~{}~{}~{}~{}~{}a_{\vec{k}}|0\rangle=0\,,$
(53)

where the last condition defines the vacuum. The equation of motion for each
mode $u_{k}(\tau)$ takes the form of a Schrödinger equation

$\displaystyle\frac{d^{2}u_{k}}{d\tau^{2}}+\left(k^{2}-\frac{1}{z}\frac{d^{2}z}{d\tau^{2}}\right)u_{k}=0\,,$
(54)

which can be solved imposing Bunch-Davies initial conditions at some initial
time when $k\gg aH$. This choice defines the initial vacuum state as the
minimum energy eigenstate for an harmonic oscillator with time-independent
frequency in a space-time which is locally flat (because we are at length
scales much smaller than the de Sitter curvature radius). Quantization of
scalar perturbations is easier in terms of the Mukhanov-Sasaki field variable
$u(\tau,\vec{x})$ since the quadratic action corresponding to the equation of
motion in eq. (50) is the action describing a canonically normalized free
field with an effective time-dependent mass. However, the dynamics of the
modes $u_{k}$ in Fourier space is more transparent if we introduce the so-
called comoving curvature perturbation $\mathcal{R}=u/z$; in analogy with eq.
(52), $\mathcal{R}$ can be promoted to a quantum operator

$\displaystyle\hat{\mathcal{R}}(\tau,\vec{x})=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\left[\mathcal{R}_{k}(\tau)a_{\vec{k}}e^{+i\vec{k}\cdot\vec{x}}+\mathcal{R}_{k}^{*}(\tau)a^{{\dagger}}_{\vec{k}}e^{-i\vec{k}\cdot\vec{x}}\right]\,,$
(55)

with $\mathcal{R}_{k}=u_{k}/z$ that satisfies the time-evolution equation

$\displaystyle\frac{d^{2}\mathcal{R}_{k}}{dN^{2}}+\left(3+\epsilon-2\eta\right)\frac{d\mathcal{R}_{k}}{dN}+\frac{k^{2}}{a^{2}H^{2}}\mathcal{R}_{k}=0\,,$
(56)

which is the differential equation of a damped harmonic oscillator.

In the canonical picture of slow-roll inflation (that is for small Hubble
parameters $\epsilon,\eta\ll 1$), eq. (56) can be solved in two complementary
regimes divided by the so-called horizon-crossing condition $k=aH$. At early
times, when $k\gg aH$ (in such case the modes are called sub-horizon), the
last term dominates over the friction one, and the solution oscillates. As
time passes by during inflation, the comoving Hubble radius $1/aH$ shrinks,
the horizon-crossing condition is met, and one eventually enters in the regime
characterized by $k\ll aH$ (in such case the modes are called super-horizon);
the last term in eq. (56) can be neglected and the latter admits the constant
solution $d\mathcal{R}_{k}/dN=0$. In other words, after horizon crossing the
mode $\mathcal{R}_{k}$ freezes to a constant value that it maintains in time
until the end of inflation (more precisely, until the mode re-enters the
conformal Hubble horizon after the end of inflation).

In the canonical picture of slow-roll inflation, after horizon crossing
quantum fluctuations can be regarded as classical, which motivates the
description of cosmological perturbations in terms of classical random fields
Kiefer:1998jk . We will give a more precise definition of random fields in
appendix B. At the conceptual level, the previous statement implies that one
has the schematic relation $\lim_{k/aH\ll
1}\hat{\mathcal{R}}(\tau,\vec{x})=\mathcal{R}(\vec{x})$ meaning that for
super-horizon modes the quantum operator $\hat{\mathcal{R}}$ can be
interpreted as a classical random field $\mathcal{R}(\vec{x})$; notice that
the latter is time-independent because the modes $\mathcal{R}_{k}$ are frozen
in time. The previous relation can be formulated in a more precise form by
saying that vacuum expectation values on the quantum side are interpreted as
statistical averages on the classical side. Formally, in the classical picture
we still have the Fourier decomposition

$\displaystyle\mathcal{R}(\vec{x})=\int\frac{d^{3}\vec{k}}{(2\pi)^{3}}\left(\mathcal{R}_{k}a_{\vec{k}}e^{+i\vec{k}\cdot\vec{x}}+\mathcal{R}_{k}^{*}a^{{\dagger}}_{\vec{k}}e^{-i\vec{k}\cdot\vec{x}}\right)\,,$
(57)

but now creation and annihilation operators are no longer $q$-numbers but
stochastic $c$-numbers which are defined by the statistical averages

$\displaystyle\langle
a_{\vec{k}}a_{\vec{k}^{\prime}}^{{\dagger}}\rangle=\frac{1}{2}(2\pi)^{3}\delta^{(3)}(\vec{k}-\vec{k}^{\prime})=\langle
a_{\vec{k}^{\prime}}^{{\dagger}}a_{\vec{k}}\rangle\,,~{}~{}~{}~{}~{}~{}~{}\langle
a_{\vec{k}}a_{\vec{k}^{\prime}}\rangle=\langle
a^{{\dagger}}_{\vec{k}}a^{{\dagger}}_{\vec{k}^{\prime}}\rangle=0\,.$ (58)

Notice that $\langle
a_{\vec{k}}a_{\vec{k}^{\prime}}^{{\dagger}}\rangle=\langle
a_{\vec{k}^{\prime}}^{{\dagger}}a_{\vec{k}}\rangle$ is only valid for
classical fluctuations, since it implies that the stochastic parameters
commute. Furthermore, the stochastic parameters $a_{\vec{k}}$,
$a_{\vec{k}}^{{\dagger}}$ have zero mean value (consequently,
$\langle\mathcal{R}(\vec{x})\rangle=0$) in order to match the fact that the
quantum operator $\hat{\mathcal{R}}(\tau,\vec{x})$ in eq. (55) has zero vacuum
expectation value.

The random field $\mathcal{R}(\vec{x})$ is fully specified by the entire
hierarchy of its correlation functions. The simplest one is the two-point
correlation function. The latter can be defined by introducing the idea of
power spectrum. The power spectrum $\Delta_{\mathcal{R}}(k)$ is defined by the
Fourier transform of the two-point correlation function

$\displaystyle\langle\mathcal{R}(\vec{x})\mathcal{R}(\vec{x}+\vec{r})\rangle=\int\frac{d^{3}k}{(2\pi)^{3}}e^{i\vec{k}\cdot\vec{r}}\Delta_{\mathcal{R}}(k)=\int_{0}^{\infty}dk\frac{\sin(kr)}{kr}\frac{k^{2}}{2\pi^{2}}\Delta_{\mathcal{R}}(k)\,.$
(59)

The fact that $\Delta_{\mathcal{R}}(k)$ depends only on $k\equiv|\vec{k}|$ and
the explicit angular integrations that we performed in eq. (59) are
consequences of the assumptions of spatial homogeneity and isotropy
(equivalently, spatial homogeneity and isotropy imply that
$\langle\mathcal{R}(\vec{x})\mathcal{R}(\vec{x}+\vec{r})\rangle$ only depends
on the relative distance $r\equiv|\vec{r}|$). The limit $r\to 0$ in eq. (59)
defines the variance $\sigma_{0}^{2}$ of the comoving curvature perturbation

$\displaystyle\sigma_{0}^{2}\equiv\langle\mathcal{R}(\vec{x})\mathcal{R}(\vec{x})\rangle=\lim_{r\to
0}\langle\mathcal{R}(\vec{x})\mathcal{R}(\vec{x}+\vec{r})\rangle=\int_{0}^{\infty}\frac{dk}{k}\,\frac{k^{3}}{2\pi^{2}}\Delta_{\mathcal{R}}(k)\equiv\int_{0}^{\infty}\frac{dk}{k}\,\mathcal{P}_{\mathcal{R}}(k)\,,$
(60)

where we defined the dimensionless power spectrum
$\mathcal{P}_{\mathcal{R}}(k)\equiv(k^{3}/2\pi^{2})\Delta_{\mathcal{R}}(k)$.
This equation gives to the power spectrum an intuitive statistical meaning;
$\mathcal{P}_{\mathcal{R}}(k)$ represents the contribution to the variance of
the field per unit logarithmic bin around the comoving wavenumber $k$. As
stated above, we can also compute eq. (60) by taking the vacuum expectation
value of the quantum operator $\hat{\mathcal{R}}(\tau,\vec{x})$. We have
(using eq. (53))

$\displaystyle\lim_{k/aH\ll
1}\langle\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\rangle=\lim_{k/aH\ll
1}\int_{0}^{\infty}\frac{dk}{k}\,\frac{k^{3}}{2\pi^{2}}|\mathcal{R}_{k}(\tau)|^{2}=\int_{0}^{\infty}\frac{dk}{k}\,\frac{k^{3}}{2\pi^{2}}|\mathcal{R}_{k}|^{2}\,,$
(61)

where the last step means that we are considering a sufficiently late time
(the limit $\lim_{k/aH\ll 1}$, for fixed $k$, represents a time-limit since
$aH$ depends on time) such that the mode with comoving wavenumber $k$ is
frozen to its constant value after horizon crossing (and the time-dependence
drops in the last equality). Eq. (61) gives an operative definition of the
power spectrum. For a given model of inflation, we can solve eq. (56)
(equivalently, eq. (54)) for each $k$ and take a “late-time limit” in the
sense specified above (that is we evaluate $\mathcal{R}_{k}$ at some late time
after it freezes to a constant value). The power spectrum is given by
$\mathcal{P}_{\mathcal{R}}(k)\equiv(k^{3}/2\pi^{2})|\mathcal{R}_{k}|^{2}$ and
fully specifies the two-point correlator of the random field
$\mathcal{R}(\vec{x})$. Equivalently, eq. (61) can be derived from the
computation of $\langle\mathcal{R}(\vec{x})\mathcal{R}(\vec{x})\rangle$ by
means of the decomposition given in eq. (57).

To proceed further, we need to specify the structure of higher-order
correlators. Under the assumption that the random field $\mathcal{R}(\vec{x})$
is gaussian, however, the power spectrum is enough to fully reconstruct
higher-order correlators. This is because the $N$-point correlation function
either vanishes (for odd $N$) or can be expressed in terms of the power
spectrum as a consequence of the Isserlis’ theorem (for even $N$). The
assumption that the statistics of the random field $\mathcal{R}(\vec{x})$ is
gaussian is well-motivated in the context of the canonical picture of slow-
roll inflation. This is because if one tries to compute, on the quantum-side,
the three-point correlator $\lim_{k/aH\ll
1}\langle\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\rangle$
the resulting expression turns out to be suppressed—compared to the two-point
correlator—by additional powers of the Hubble parameters $\epsilon$ and $\eta$
Maldacena:2002vr . This means that non-gaussianities are not relevant during
conventional slow-roll dynamics during which the Hubble parameters take $O(\ll
1)$ values (said differently, this means that any detection of sizable non-
gaussianities at CMB scales will rule out all single field slow-roll models of
inflation).

However, we are interested in a situation which deviates from standard slow-
roll dynamics. We refer to ref. Ballesteros:2020qam  for a detailed (both
analytical and numerical) study, and we summarize here the main points.

Standard slow-roll dynamics takes place at large field values. However, few
$e$-folds before the end of inflation (which ends at the absolute minimum of
the potential located at the origin) the inflaton field crosses an approximate
stationary inflection point. The inflaton field almost stops but it possesses
just enough inertia to overcome the approximate stationary inflection point.
During this part of the dynamics the Hubble parameter $\eta$ transits from
$\eta\simeq 0$ (that is typical of slow-roll) to a large positive value that
is maintained for few $e$-folds until the field crosses the approximate
stationary inflection point. If $\eta\gtrsim 3/2$ (typically one has
$\eta\gtrsim 3$), the friction term in eq. (56) becomes negative. This part of
the dynamics characterized by the presence of negative friction is dubbed
ultra slow-roll. During the negative friction phase, the modes
$\mathcal{R}_{k}$—more precisely, their modulus $|\mathcal{R}_{k}|$—change
exponentially fast, and can be either enhanced or suppressed depending on the
specific value of $k$.

After the end of ultra slow-roll, $\eta$ transits to a phase during which it
takes negative $O(1)$ values. The friction term in eq. (56) turns positive,
and the modes $\mathcal{R}_{k}$—after being enhanced or suppressed by the
negative friction phase—are now free to freeze to their final constant value.
Since
$\mathcal{P}_{\mathcal{R}}(k)\equiv(k^{3}/2\pi^{2})|\mathcal{R}_{k}|^{2}$, the
negative friction phase that modifies exponentially $|\mathcal{R}_{k}|$
produces a distinctive peak in the power spectrum of curvature perturbations.

This peculiar dynamics has important consequences as far as non-gaussianities
are concerned. Because of the presence of the negative friction phase,
classicalization of the modes do not happens after their horizon crossing but
is delayed after the end of ultra slow-roll Ballesteros:2020sre . This means
that the three-point correlator
$\langle\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\rangle$
has to be evaluated after the end of ultra slow-roll. Crucially, after the end
of ultra slow-roll $\eta$ takes sizable negative $O(1)$ values (while we
expect $\epsilon\ll 1$). Let us indicate this value with $\eta_{0}$ (which is
a negative number). This implies that non-gaussianities are no longer
negligible since the expected slow-roll suppression is not valid anymore. The
explicit computation of the three-point correlator in Fourier space gives, for
a triad of comoving wavenumbers $k_{1},k_{2},k_{3}$, the so-called local
bispectrum Atal:2018neu

$\displaystyle
B_{\mathcal{R}}(k_{1},k_{2},k_{3})\simeq-\eta_{0}[\Delta_{\mathcal{R}}(k_{1})\Delta_{\mathcal{R}}(k_{2})+\Delta_{\mathcal{R}}(k_{1})\Delta_{\mathcal{R}}(k_{3})+\Delta_{\mathcal{R}}(k_{2})\Delta_{\mathcal{R}}(k_{3})]\,,$
(62)

where the three-point correlator is given by the Fourier transform

$\displaystyle\lim_{\hskip 42.67912pt{\rm
end\,of\,USR}}\langle\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\hat{\mathcal{R}}(\tau,\vec{x})\rangle=\int\frac{d^{3}\vec{k}_{1}}{(2\pi)^{3}}\frac{d^{3}\vec{k}_{2}}{(2\pi)^{3}}\frac{d^{3}\vec{k}_{3}}{(2\pi)^{3}}(2\pi)^{3}\delta^{(3)}(\vec{k}_{1}+\vec{k}_{2}+\vec{k}_{3})B_{\mathcal{R}}(k_{1},k_{2},k_{3})\,,$
(63)

and where $\lim_{{\rm end\,of\,USR}}$ indicates explicitly that the three-
point correlator has to be evaluated at some late time after the end of ultra
slow-roll, when the modes $\mathcal{R}_{k}$ finally set to their final
constant value (so that the right-hand side of eq. (63) does not depend on
time).

The analysis carried out in refs. Atal:2018neu  shows that non-gaussianities
in the presence of ultra slow-roll are expected to be non-negligible. Eqs.
(62, 63) represent the quantum side of the story, and eq. (62) can be obtained
by means of the so-called “in-in” formalism Maldacena:2002vr . It is important
to understand the implications from the point of view of the random field
$\mathcal{R}(\vec{x})$. Consider the non-gaussian random field
$\mathcal{R}_{\rm NG}(\vec{x})$ defined by

$\displaystyle\mathcal{R}_{\rm NG}(\vec{x})\equiv\mathcal{R}_{\rm
G}(\vec{x})+\frac{(-\eta_{0})}{2}\left[\mathcal{R}_{\rm
G}(\vec{x})^{2}-\langle\mathcal{R}_{\rm G}(\vec{x})^{2}\rangle\right]\,,$ (64)

where $\mathcal{R}_{\rm G}(\vec{x})$ is a gaussian random field—with variance
$\sigma_{0}^{2}$ (see eq. (60))—which admits the decomposition given in eq.
(57); it is a simple exercise to show that the three-point correlator
$\langle\mathcal{R}_{\rm NG}(\vec{x})\mathcal{R}_{\rm
NG}(\vec{x})\mathcal{R}_{\rm NG}(\vec{x})\rangle$ has precisely the form given
in eq. (63). Notice that in eq. (64) the presence of the constant piece
$\langle\mathcal{R}_{\rm G}(\vec{x})^{2}\rangle$ guarantees that the non-
gaussian random field $\mathcal{R}_{\rm NG}(\vec{x})$ has zero mean. This
property is physically motivated by the fact that the background solution is
stable.

The physics-case sketched in this appendix motivates the non-gaussianities
studied in this paper which are of the form given in eq. (64). To avoid
cluttering the notation, in the main part of this work we indicate simply with
$\mathcal{R}$ the gaussian random field $\mathcal{R}_{\rm G}$ in eq. (64) and
with $h$ the non-gaussian one $\mathcal{R}_{\rm NG}$. Furthermore, we set
$\alpha\equiv(-\eta_{0})/2$.

Before proceeding, an important comment is in order. As we have discussed, the
structure of the non-gaussian random field given in eq. (64) is motivated by
the explicit computation of the three-point correlator in eqs. (62, 63) in the
presence of ultra slow-roll. This computation is based on the “in-in”
formalism in which one expands the interaction Hamiltonian up to the cubic
order.

However, computing only the bispectrum is not the end of the story. In
principle, one should compute also the trispectrum (that is the connected part
of the four-point correlator) and check that its contribution does not alter
the form of non-gaussianities derived including only cubic interactions.
Needless to say, the computation of the trispectrum is anything but simple,
and we are not aware of explicit results in the context of ultra slow-roll.
From a more pragmatic phenomenological perspective, one possible way to
proceed—widely used for the analysis of CMB non-Gaussianity, see ref.
Kogo:2006kh —is the following. Instead of eq. (64), one takes the more general
ansatz

$\displaystyle\mathcal{R}_{\rm NG}(\vec{x})\equiv\mathcal{R}_{\rm
G}(\vec{x})+f_{2}\left[\mathcal{R}_{\rm
G}(\vec{x})^{2}-\langle\mathcal{R}_{\rm
G}(\vec{x})^{2}\rangle\right]+f_{3}\mathcal{R}_{\rm G}(\vec{x})^{3}\,,$ (65)

where $f_{2,3}$ are free coefficients that parametrize quadratic and cubic
deviations from the gaussian limit. Notice that eq. (65) generalizes eq. (64)
in the sense that it introduces cubic corrections but it preserves locality
(in the sense that deviations from exact gaussianity at $\vec{x}$ are located
at the same spatial position). As discussed before, in the absence of an
explicit computation there is no guarantee that ultra slow-roll generates
deviations from eq. (64) that have the form given in eq. (65). Nevertheless,
eq. (65) can be considered as a phenomenological parametrization to study
deviations from from eq. (64), as done for instance in ref. Young:2013oia .
For this reason, in our analysis we will try to set the formalism considering
a generic deviation from gaussianity (but always assuming locality) even
though we will present our results for the motivated case of quadratic non-
gaussianities given in eq. (64).

Finally, let us mention that a resummation of local non-gaussianities at all
orders has been presented in ref. Atal:2019cdz  in the context of the $\delta
N$ formalism. The analysis is based on an expansion of the inflaton potential
around a local maximum. It would be interesting to go beyond such
approximation and extend this analysis for the potentials studied here.
Moreover it is also worth investigating the role of higher order contributions
leading to non-local non-gaussianities.

##  Appendix B CliffsNotes on gaussian peak theory

As a warm-up, consider the case of a $n$-dimensional scalar gaussian random
field $\mathcal{R}(\vec{x})$ that we identify with the random field associated
to curvature perturbations. Peak theory in the case of a scalar gaussian
random field is well-known Bardeen:1985tr . However, in this appendix we will
give a detailed discussion. The reason is that in our approach there is an
important conceptual difference compared to the standard results of ref.
Bardeen:1985tr . Ref. Bardeen:1985tr  computes the number density of peaks of
the overdensity field working directly with $\delta(\vec{x},t)$, and without
any reference to curvature perturbations. In this paper, on the contrary, we
aim to compute the same quantity but starting from the distribution of local
maxima of the curvature perturbation. Following this alternative route, we
will be able to find (see appendix D) a generalization that accounts for the
case in which local non-gaussianities in the definition of
$\mathcal{R}(\vec{x})$ are present.

Let us start from basics. A $n$-dimensional scalar random field
$\mathcal{R}(\vec{x})$ is a set of random variables, one for each point
$\vec{x}$ in the $n$-dimensional real space, equipped with a probability
distribution
$p[\mathcal{R}(\vec{x}_{1}),\dots,\mathcal{R}(\vec{x}_{m})]d\mathcal{R}(\vec{x}_{1})\dots
d\mathcal{R}(\vec{x}_{m})$ which measures the probability that the function
$\mathcal{R}$ has values in the range $\mathcal{R}(\vec{x}_{j})$ to
$\mathcal{R}(\vec{x}_{j})+d\mathcal{R}(\vec{x}_{j})$ for each of the
$j=1,\dots,m$, with $m$ an arbitrary integer and
$\vec{x}_{1},\dots,\vec{x}_{m}$ arbitrary points.202020Notice that, as in ref.
Bardeen:1985tr , all spatial separations and length scales are described in
comoving coordinates in the cosmological background. This means that the
number density that we shall compute at the end of this section in eq. (101)
must be understood as a comoving number density.

We are interested in the behavior of the random field for a point in space
that is stationary, and we consider $m=1$ with $\vec{x}_{1}=\vec{x}_{\rm st}$.
We can expand around this point according to

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

from which we have

$\displaystyle\mathcal{R}_{i}(\vec{x})=\sum_{j=1}^{n}\mathcal{R}_{ij}(\vec{x}_{\rm
st})(\vec{x}-\vec{x}_{\rm st})_{j}\,.$ (67)

The goal is to obtain the number density of these stationary points in the
$n$-dimensional space. Not to violate the cosmological principle, we only want
to consider random fields which are statistically homogeneous and isotropic.
Consequently, the specific value of $\vec{x}_{\rm st}$ is irrelevant, and we
can always shift to $\vec{x}_{\rm st}=\vec{0}$. Let us, therefore, drop the
explicit dependence on $\vec{x}_{\rm st}$.

The quantity of central interest for the computation of the number density of
stationary points of $\mathcal{R}$ is the joint probability density
distribution of the field $\mathcal{R}$, its first and second derivatives.
This is intuitively obvious, since identifying maxima (or minima) implies a
set of conditions on field derivatives, and it is thus mandatory to know what
is their probability distribution (derivatives of a random field are also
random variables themselves).

Consider the realistic case with $n=3$. We indicate with
$P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})d\mathcal{R}d^{3}\mathcal{R}_{i}d^{6}\mathcal{R}_{ij}$
the joint probability distribution for the field being in the range
$\mathcal{R}$ to $\mathcal{R}+d\mathcal{R}$, the field gradient being in the
range $\mathcal{R}_{i}$ to $\mathcal{R}_{i}+d\mathcal{R}_{i}$ and the second
derivative matrix elements being in the range $\mathcal{R}_{ij}$ to
$\mathcal{R}_{ij}+d\mathcal{R}_{ij}$, all at the same point in space.
$P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})$ is the joint ten-dimensional
probability density distribution.

If the point is stationary, we can write the joint probability distribution as
$P(\mathcal{R},\mathcal{R}_{i}=0,\mathcal{R}_{ij})\left|{\rm
det}(\mathcal{R}_{ij})\right|d\mathcal{R}d^{3}\vec{x}d^{6}\mathcal{R}_{ij}$.
We set $\mathcal{R}_{i}=0$ since the point is stationary, and we used eq. (66)
to change variables in the gradient volume element. We can, therefore, write
the probability distribution to have a stationary point in a volume
$d^{3}\vec{x}$ with the field being in the range $\mathcal{R}$ to
$\mathcal{R}+d\mathcal{R}$ as

$\displaystyle n_{\rm st}(\mathcal{R})d\mathcal{R}d^{3}\vec{x}\equiv
d\mathcal{R}d^{3}\vec{x}\underbrace{\int
P(\mathcal{R},\mathcal{R}_{i}=0,\mathcal{R}_{ij})\left|{\rm
det}(\mathcal{R}_{ij})\right|d^{6}\mathcal{R}_{ij}}_{\equiv n_{\rm
st}(\mathcal{R})}\,,$ (68)

where the integral is extended to the whole range of variability of the second
derivatives since we are considering generic stationary points. The
probability distribution to have a maximum in a volume $d^{3}\vec{x}$ with the
field being in the range $\mathcal{R}$ to $\mathcal{R}+d\mathcal{R}$ is

$\displaystyle n_{\rm max}(\mathcal{R})d\mathcal{R}d^{3}\vec{x}\equiv
d\mathcal{R}d^{3}\vec{x}\underbrace{\int_{\rm
max}P(\mathcal{R},\mathcal{R}_{i}=0,\mathcal{R}_{ij})\left|{\rm
det}(\mathcal{R}_{ij})\right|d^{6}\mathcal{R}_{ij}}_{\equiv n_{\rm
max}(\mathcal{R})}\,,$ (69)

where now the integral is subject to the conditions on the Hessian matrix that
define a maximum. An equivalent definition holds in the case of a minimum. The
probability density distribution $n_{\rm max}(\mathcal{R})$ represents the
number density of local maxima (where “number density” is defined in a
probabilistic sense) with field value in the range $\mathcal{R}$ to
$\mathcal{R}+d\mathcal{R}$.

In order to extract quantitative informations, we need to compute
$P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})$, set $\mathcal{R}_{i}=0$,
and integrate. This strategy does not depend on the specific statistics of the
random field.

The computation of $P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})$
drastically simplifies in the gaussian case. This is because the joint
probability distribution of a gaussian field and its derivatives is a
multivariate normal distribution. Consider the simplified case with $n=2$ in
which we have more control on analytical formulas;212121The simplest
possibility would be $n=1$. However, the case $n=2$ is the simplest setup in
which a non-trivial discussion about spatial isotropy is possible. we have six
random variables that we collect in the column vector
${R}\equiv(\mathcal{R},\mathcal{R}_{x},\mathcal{R}_{y},\mathcal{R}_{xx},\mathcal{R}_{xy},\mathcal{R}_{yy})^{\rm
T}$. We have

$\displaystyle
P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})=\frac{1}{(2\pi)^{k/2}\sqrt{{\rm
det}\,{C}}}\exp\left[-\frac{1}{2}({R}-\langle{R}\rangle)^{\rm
T}({C}^{-1})({R}-\langle{R}\rangle)\right]\,,$ (70)

where $k=6$ is the dimension of $R$, $\langle{R}\rangle$ is the column vector
of the expectation values of ${R}$ and ${C}$ is the covariance matrix defined
by $C\equiv\langle({R}-\langle{R}\rangle)({R}-\langle{R}\rangle)^{\rm
T}\rangle$ with elements

$\displaystyle
C_{ij}=\langle({R}-\langle{R}\rangle)_{i}({R}-\langle{R}\rangle)_{j}\rangle=\langle
R_{i}R_{j}\rangle-\langle R_{i}\rangle\langle R_{j}\rangle\,.$ (71)

From the computation of the covariance matrix one can fully reconstruct the
joint probability distribution.

We restrict the analysis to zero-mean random fields since this assumption is
physically motivated (see discussion in appendix A). From the explicit
computation of the two-point correlators $\langle R_{i}R_{j}\rangle$, one
finds that the covariance matrix takes the form

$C=\bordermatrix{&{\scriptstyle\mathcal{R}}&{\scriptstyle\mathcal{R}_{x}}&{\scriptstyle\mathcal{R}_{y}}&{\scriptstyle\mathcal{R}_{xx}}&{\scriptstyle\mathcal{R}_{xy}}&{\scriptstyle\mathcal{R}_{yy}}\cr{\scriptstyle\mathcal{R}}&\sigma_{0}^{2}&0&0&-\sigma_{1}^{2}/2&0&-\sigma_{1}^{2}/2\cr{\scriptstyle\mathcal{R}_{x}}&0&\sigma_{1}^{2}/2&0&0&0&0\cr{\scriptstyle\mathcal{R}_{y}}&0&0&\sigma_{1}^{2}/2&0&0&0\cr{\scriptstyle\mathcal{R}_{xx}}&-\sigma_{1}^{2}/2&0&0&3\sigma_{2}^{2}/8&0&\sigma_{2}^{2}/8\cr{\scriptstyle\mathcal{R}_{xy}}&0&0&0&0&\sigma_{2}^{2}/8&0\cr{\scriptstyle\mathcal{R}_{yy}}&-\sigma_{1}^{2}/2&0&0&\sigma_{2}^{2}/8&0&3\sigma_{2}^{2}/8}\,,$
(72)

and we find that $\mathcal{R}_{x}$, $\mathcal{R}_{y}$ and $\mathcal{R}_{xy}$
are completely uncorrelated while $\mathcal{R}$, $\mathcal{R}_{xx}$ and
$\mathcal{R}_{yy}$ are correlated. We introduce the spectral moments

$\sigma_{j}^{2}\equiv\int\frac{dk}{k}\mathcal{P}_{\mathcal{R}}(k)\,k^{2j}\,,$
(73)

where $\mathcal{P}_{\mathcal{R}}(k)$ is the dimensionless power spectrum of
$\mathcal{R}$. Notice that in two spatial dimensions the dimensionless power
spectrum of $\mathcal{R}$ is related to the power spectrum by means of
$\mathcal{P}_{\mathcal{R}}(k)=(k^{2}/2\pi)\Delta_{\mathcal{R}}(k)$ with
$\Delta_{\mathcal{R}}(k)=|\mathcal{R}_{k}|^{2}$.

Consider, as an illustrative example, the computation of
$\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle$. We use the explicit form of
$\mathcal{R}$ given in eq. (57). We find

$\displaystyle\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle$
$\displaystyle=\int\frac{d^{2}\vec{k}_{1}}{(2\pi)^{2}}\frac{d^{2}\vec{k}_{2}}{(2\pi)^{2}}\langle\left(-k_{1,x}^{2}\mathcal{R}_{k_{1}}a_{\vec{k}_{1}}e^{i\vec{k}_{1}\cdot\vec{x}}-k_{1,x}^{2}\mathcal{R}_{k_{1}}^{*}a^{{\dagger}}_{\vec{k}_{1}}e^{-i\vec{k}_{1}\cdot\vec{x}}\right)\left(-k_{2,y}^{2}\mathcal{R}_{k_{2}}a_{\vec{k}_{2}}e^{i\vec{k}_{2}\cdot\vec{x}}-k_{2,y}^{2}\mathcal{R}_{k_{2}}^{*}a^{{\dagger}}_{\vec{k}_{2}}e^{-i\vec{k}_{2}\cdot\vec{x}}\right)\rangle$
$\displaystyle=\int\frac{d^{2}\vec{k}_{1}}{(2\pi)^{2}}\frac{d^{2}\vec{k}_{2}}{(2\pi)^{2}}k_{1,x}^{2}k_{2,y}^{2}\left[e^{i(\vec{k}_{1}-\vec{k}_{2})\cdot\vec{x}}\mathcal{R}_{k_{1}}\mathcal{R}_{k_{2}}^{*}\langle
a_{\vec{k}_{1}}a^{{\dagger}}_{\vec{k}_{2}}\rangle+e^{-i(\vec{k}_{1}-\vec{k}_{2})\cdot\vec{x}}\mathcal{R}_{k_{2}}\mathcal{R}_{k_{1}}^{*}\langle
a^{{\dagger}}_{\vec{k}_{1}}a_{\vec{k}_{2}}\rangle\right]$
$\displaystyle=\int\frac{d^{2}\vec{k}}{(2\pi)^{2}}k_{x}^{2}k_{y}^{2}|\mathcal{R}_{k}|^{2}=\frac{1}{4\pi^{2}}\int
dkd\varphi
k^{5}\cos^{2}\varphi\sin^{2}\varphi\Delta_{\mathcal{R}}(k)=\frac{1}{8}\int\frac{dk}{k}k^{4}\mathcal{P}_{\mathcal{R}}(k)\,,$
(74)

where in the last line we just introduced polar coordinates. All the entries
in eq. (72) can be computed in a similar way.

Interestingly, the pattern of zeros in eq. (72) and the relations among
different non-zero entries—obtained before by means of a direct
computation—can be understood as a consequence of homogeneity and isotropy.

Homogeneity, that is translational invariance, implies that correlators do not
depend on the specific spatial position at which they are computed. For
instance, this means that the spatial derivative of
$\langle\mathcal{R}\mathcal{R}\rangle$ should vanish (this must be true for a
generic correlator evaluated at a given spatial point); from this condition,
one finds
$\partial_{i}(\langle\mathcal{R}\mathcal{R}\rangle)=0\to\langle\mathcal{R}\mathcal{R}_{i}\rangle=0$
so that $\mathcal{R}$, $\mathcal{R}_{x}$ and $\mathcal{R}_{y}$ are
uncorrelated. Similarly, from
$\partial_{x}(\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle)=0$ one gets
$\langle\mathcal{R}_{x}\mathcal{R}_{xx}\rangle=0$ (with similar relations
along other directions) so that first and second derivatives are uncorrelated.
On the contrary, from
$\partial_{x}(\langle\mathcal{R}\mathcal{R}_{x}\rangle)=0$ it follows that
$\langle\mathcal{R}\mathcal{R}_{xx}\rangle=-\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle$
as indeed obtained in eq. (72). Similarly,
$\langle\mathcal{R}\mathcal{R}_{yy}\rangle=-\langle\mathcal{R}_{y}\mathcal{R}_{y}\rangle$.

Isotropy, that is rotational invariance, implies that correlators do not
depend on a particular direction in space. The simplest consequence of
isotropy is that
$\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle=\langle\mathcal{R}_{y}\mathcal{R}_{y}\rangle$
and $\langle\mathcal{R}_{x}\mathcal{R}_{y}\rangle=0$. In order to derive these
two conditions, a very useful trick (that we shall use also in the non-
gaussian computation) is to introduce—instead of the two components $x$ and
$y$ of the two-dimensional vector $\vec{x}$—the complex conjugated variables
$z\equiv x+iy$ and $z^{*}=x-iy$ from which we have
$\partial_{z}=(\partial_{x}-i\partial_{y})/2$ and
$\partial_{z^{*}}=(\partial_{x}+i\partial_{y})/2$. If we now rotate the two-
dimensional vector $\vec{x}$ (without changing its length) the complex number
$z$ changes by a phase factor $e^{i\tau}$, that is $z\to e^{i\tau}z$ (and
$z^{*}\to e^{-i\tau}z^{*}$). Consequently, the derivatives with respect to $z$
and $z^{*}$ change according to $\partial_{z}\to e^{-i\tau}\partial_{z}$ and
$\partial_{z^{*}}\to e^{i\tau}\partial_{z^{*}}$. If we now consider the
correlators $\langle\mathcal{R}_{z}\mathcal{R}_{z}\rangle$ and
$\langle\mathcal{R}_{z^{*}}\mathcal{R}_{z^{*}}\rangle$ they rotate according
to $\langle\mathcal{R}_{z}\mathcal{R}_{z}\rangle\to
e^{-2i\tau}\langle\mathcal{R}_{z}\mathcal{R}_{z}\rangle$ and
$\langle\mathcal{R}_{z^{*}}\mathcal{R}_{z^{*}}\rangle\to
e^{2i\tau}\langle\mathcal{R}_{z^{*}}\mathcal{R}_{z^{*}}\rangle$. Because of
isotropy of the two-dimensional space,
$\langle\mathcal{R}_{z}\mathcal{R}_{z}\rangle$ and
$\langle\mathcal{R}_{z^{*}}\mathcal{R}_{z^{*}}\rangle$ can not depend on
$\tau$ and, therefore, they must vanish. Consequently, the system of equations

$\displaystyle\langle\mathcal{R}_{z}\mathcal{R}_{z}\rangle$
$\displaystyle=\frac{1}{4}\left(\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle-2i\langle\mathcal{R}_{x}\mathcal{R}_{y}\rangle-\langle\mathcal{R}_{y}\mathcal{R}_{y}\rangle\right)=0\,,$
(75) $\displaystyle\langle\mathcal{R}_{z^{*}}\mathcal{R}_{z^{*}}\rangle$
$\displaystyle=\frac{1}{4}\left(\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle+2i\langle\mathcal{R}_{x}\mathcal{R}_{y}\rangle-\langle\mathcal{R}_{y}\mathcal{R}_{y}\rangle\right)=0\,,$
(76)

admits the solution
$\langle\mathcal{R}_{x}\mathcal{R}_{x}\rangle=\langle\mathcal{R}_{y}\mathcal{R}_{y}\rangle$
and $\langle\mathcal{R}_{x}\mathcal{R}_{y}\rangle=0$ which are precisely the
relations we were looking for. More in general, if we consider the rotation
$\partial_{z}\to e^{-i\tau}\partial_{z}$ and $\partial_{z^{*}}\to
e^{i\tau}\partial_{z^{*}}$ a generic correlator will take a phase factor
$e^{i\kappa\tau}$ where $\kappa\equiv(\\#\,z^{*}\,{\rm
derivatives})-(\\#\,z\,{\rm derivatives})$. If $\kappa\neq 0$, then the
correlator must be equal to zero as a consequence of isotropy. For instance,
we have $\langle\mathcal{R}_{zz}\mathcal{R}_{z^{*}}\rangle=0$ (because
$\kappa=-1$) but $\langle\mathcal{R}_{z}\mathcal{R}_{z^{*}}\rangle\neq 0$
(because $\kappa=0$).

Combining homogeneity and isotropy validates the remaining entries in eq.
(72). For instance, from
$\partial_{y}(\langle\mathcal{R}\mathcal{R}_{x}\rangle)=\langle\mathcal{R}_{x}\mathcal{R}_{y}\rangle+\langle\mathcal{R}\mathcal{R}_{xy}\rangle=0$
(homogeneity) we find $\langle\mathcal{R}\mathcal{R}_{xy}\rangle=0$ since
isotropy implies $\langle\mathcal{R}_{x}\mathcal{R}_{y}\rangle=0$. As a final
check, we consider the block of the second derivatives in eq. (72). From the
previous argument, isotropy implies that
$\langle\mathcal{R}_{zz}\mathcal{R}_{zz}\rangle=0$ and
$\langle\mathcal{R}_{zz^{*}}\mathcal{R}_{zz}\rangle=0$ (together with their
complex conjugated
$\langle\mathcal{R}_{z^{*}z^{*}}\mathcal{R}_{z^{*}z^{*}}\rangle=0$ and
$\langle\mathcal{R}_{zz^{*}}\mathcal{R}_{z^{*}z^{*}}\rangle=0$). These two
relations imply
$\langle\mathcal{R}_{xx}\mathcal{R}_{xy}\rangle=\langle\mathcal{R}_{yy}\mathcal{R}_{xy}\rangle=0$
and
$\langle\mathcal{R}_{xx}\mathcal{R}_{xx}\rangle=\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle+2\langle\mathcal{R}_{xy}\mathcal{R}_{xy}\rangle$.
Both these relations are verified by the entries in eq. (72). We can actually
do more since it is possible to show that
$\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle=\langle\mathcal{R}_{xy}\mathcal{R}_{xy}\rangle$.
Let us write

$\displaystyle\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle$
$\displaystyle=\langle\partial_{xx}\mathcal{R}(\vec{x})\partial_{yy}\mathcal{R}(\vec{x})\rangle=\left.\langle\partial_{x_{1}x_{1}}\mathcal{R}(\vec{x}_{1})\partial_{y_{2}y_{2}}\mathcal{R}(\vec{x}_{2})\rangle\right|_{\vec{x}_{1}=\vec{x}_{2}=\vec{x}}$
$\displaystyle=\left.\partial_{x_{1}x_{1}}\partial_{y_{2}y_{2}}\langle\mathcal{R}(\vec{x}_{1})\mathcal{R}(\vec{x}_{2})\rangle\right|_{\vec{x}_{1}=\vec{x}_{2}=\vec{x}}=\left.\partial_{x_{1}}\partial_{y_{2}}\langle\mathcal{R}_{x_{1}}(\vec{x}_{1})\mathcal{R}_{y_{2}}(\vec{x}_{2})\rangle\right|_{\vec{x}_{1}=\vec{x}_{2}=\vec{x}}\,.$
(77)

where the first step is just a more explicit definition of
$\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle$ while in the following ones
we consider two distinct point $\vec{x}_{1}=(x_{1},y_{1})$ and
$\vec{x}_{2}=(x_{2},y_{2})$ that we later set equal again. Now the point is
that because of homogeneity of space the correlator
$\langle\mathcal{R}_{x_{1}}(\vec{x}_{1})\mathcal{R}_{y_{2}}(\vec{x}_{2})\rangle$
depends only on the distance $|\vec{x}_{1}-\vec{x}_{2}|$; we can, therefore,
exchange $1\leftrightarrow 2$ obtaining
$\langle\mathcal{R}_{x_{2}}(\vec{x}_{2})\mathcal{R}_{y_{1}}(\vec{x}_{1})\rangle$
without altering the result. From eq. (77) this means that we have
$\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle=\langle\mathcal{R}_{xy}\mathcal{R}_{xy}\rangle$
as indeed verified in eq. (72). If we combine this result with the previous
relation
$\langle\mathcal{R}_{xx}\mathcal{R}_{xx}\rangle=\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle+2\langle\mathcal{R}_{xy}\mathcal{R}_{xy}\rangle$
we find
$\langle\mathcal{R}_{xx}\mathcal{R}_{xx}\rangle=3\langle\mathcal{R}_{xy}\mathcal{R}_{xy}\rangle$
which is again verified in eq. (72). Finally, we also note that the condition
$\langle\mathcal{R}_{xx}\mathcal{R}_{yy}\rangle=\langle\mathcal{R}_{xy}\mathcal{R}_{xy}\rangle$
implies that
$\langle\mathcal{R}_{zz^{*}}\mathcal{R}_{zz^{*}}\rangle=\langle\mathcal{R}_{zz}\mathcal{R}_{z^{*}z^{*}}\rangle$.
These kind of relations based on homogeneity and isotropy will be useful later
in the non-gaussian case.

Using the properties of the exponential function, eq. (70) takes the form

$\displaystyle
P(\mathcal{R},\mathcal{R}_{i},\mathcal{R}_{ij})=P(\mathcal{R}_{x})P(\mathcal{R}_{y})P(\mathcal{R}_{xy})P(\mathcal{R},\mathcal{R}_{xx},\mathcal{R}_{yy})\,,$
(78)

where

$\displaystyle
P(\mathcal{R}_{x})=\frac{1}{\sqrt{\pi\sigma_{1}^{2}}}\exp\left(-\frac{\mathcal{R}_{x}^{2}}{\sigma_{1}^{2}}\right)\,,\hskip
19.91684ptP(\mathcal{R}_{y})=\frac{1}{\sqrt{\pi\sigma_{1}^{2}}}\exp\left(-\frac{\mathcal{R}_{y}^{2}}{\sigma_{1}^{2}}\right)\,,\hskip
19.91684ptP(\mathcal{R}_{xy})=\frac{2}{\sqrt{\pi\sigma_{2}^{2}}}\exp\left(-\frac{4\mathcal{R}_{xy}^{2}}{\sigma_{2}^{2}}\right)\,,$
(79)

and

$\displaystyle
P(\mathcal{R},\mathcal{R}_{xx},\mathcal{R}_{yy})=\frac{1}{(2\pi)^{3/2}\sqrt{{\rm
det}\tilde{C}}}\exp\left(-\frac{1}{2}\tilde{R}^{\rm
T}\tilde{C}^{-1}\tilde{R}\right)\,,\hskip
28.45274pt\tilde{C}\equiv\bordermatrix{&{\scriptstyle\mathcal{R}}&{\scriptstyle\mathcal{R}_{xx}}&{\scriptstyle\mathcal{R}_{yy}}&\cr{\scriptstyle\mathcal{R}}&\sigma_{0}^{2}&-\sigma_{1}^{2}/2&-\sigma_{1}^{2}/2\cr{\scriptstyle\mathcal{R}_{xx}}&-\sigma_{1}^{2}/2&3\sigma_{2}^{2}/8&\sigma_{2}^{2}/8\cr{\scriptstyle\mathcal{R}_{yy}}&-\sigma_{1}^{2}/2&\sigma_{2}^{2}/8&3\sigma_{2}^{2}/8}\,,$
(80)

with $\tilde{R}\equiv(\mathcal{R},\mathcal{R}_{xx},\mathcal{R}_{yy})^{\rm T}$.
As customary in the gaussian case, from the knowledge of the power spectrum it
is possible to fully reconstruct the statistics of the random field. From eq.
(69) we get the number density of maxima

$\displaystyle n_{\rm max}(\mathcal{R})=\frac{1}{\pi\sigma_{1}^{2}}\int_{\rm
max}d\mathcal{R}_{xx}d\mathcal{R}_{xy}d\mathcal{R}_{yy}\left|\mathcal{R}_{xx}\mathcal{R}_{yy}-\mathcal{R}_{xy}^{2}\right|P(\mathcal{R}_{xy})P(\mathcal{R},\mathcal{R}_{xx},\mathcal{R}_{yy})\,,$
(81)

where we used
$P(\mathcal{R}_{x}=0)=P(\mathcal{R}_{y}=0)=1/\sqrt{\pi\sigma_{1}^{2}}$. The
integration region is defined by the conditions ${\rm
max}=\\{\mathcal{R}_{xx}\mathcal{R}_{yy}-\mathcal{R}_{xy}^{2}>0\land\mathcal{R}_{xx}<0\land\mathcal{R}_{yy}<0\\}$.

The change of variables
$\\{\mathcal{R}_{xx},\mathcal{R}_{yy},\mathcal{R}_{xy}\\}\to\\{r,s,\theta\\}$
defined by

$\displaystyle
r\cos\theta\equiv\frac{1}{2}\left(\mathcal{R}_{xx}-\mathcal{R}_{yy}\right)\,,\hskip
21.33955ptr\sin\theta\equiv\mathcal{R}_{xy}\,,\hskip
21.33955pts\equiv-\frac{1}{2}\left(\mathcal{R}_{xx}+\mathcal{R}_{yy}\right)\,,$
(82)

turns out the be useful. The Jacobian of the transformation is $J=2r$, and we
have $\mathcal{R}_{xx}\mathcal{R}_{yy}-\mathcal{R}_{xy}^{2}=s^{2}-r^{2}$. The
condition $s^{2}-r^{2}>0$ becomes $-s<r<s$ with $s>0$ since
$\mathcal{R}_{xx}<0$ and $\mathcal{R}_{yy}<0$. Furthermore, we restrict to
$r>0$ if $0<\theta<2\pi$. All in all, we have ${\rm max}=\\{0<\theta<2\pi\land
s>0\land 0<r<s\\}$.

The parametrization in terms of $\\{r,s,\theta\\}$ is useful because we have

$\displaystyle
2s=-(\mathcal{R}_{xx}+\mathcal{R}_{yy})=-\triangle\mathcal{R}\,.$ (83)

This implies that any condition that restricts the value of the curvature
$-\triangle\mathcal{R}$ can be implemented through $s$ by imposing $s>s_{\rm