Renormalized Primordial Black Holes

G. Franciolini,1 A. Ianniccari,2 A. Kehagias, 2,3 D. Perrone,2 and A. Riotto,2

1 CERN, Theoretical Physics Department, Esplanade des Particules 1, Geneva
1211, Switzerland

2 Department of Theoretical Physics and Gravitational Wave Science Center
(GWSC),
Université de Genève, CH-1211 Geneva

3 Physics Division, National Technical University of Athens, Athens, 15780,
Greece

Abstract
The formation of primordial black holes in the early universe may happen
through the collapse of large curvature perturbations generated during a non-
attractor phase of inflation or through a curvaton-like dynamics after
inflation. The fact that such small-scale curvature perturbation is typically
non-Gaussian leads to the renormalization of composite operators built up from
the smoothed density contrast and entering in the calculation of the
primordial black abundance. Such renormalization causes the phenomenon of
operator mixing and the appearance of an infinite tower of local, non-local
and higher-derivative operators as well as to a sizable shift in the threshold
for primordial black hole formation. This hints that the calculation of the
primordial black hole abundance is more involved than what generally assumed.

###### Contents

  1. 1 Introduction
  2. 2  A brief summary of the standard PBH abundance calculation
  3. 3 The non-Gaussianity of the curvature perturbation
  4. 4 Renormalization in the PBHs
    1. 4.1 Feynman rules
  5. 5 Renormalization at one-loop: linear operator mixing
  6. 6 Renormalization at one-loop: quadratic operator mixing
    1. 6.1 The impact of renormalization
  7. 7 Comments and conclusions
  8. A Integrals involving the Gegenbauer polynomials
    1. A.1 Integral of one window function
    2. A.2 Integral with two window functions
  9. B Non-local operators

##  1 Introduction

The physics of Primordial Black Holes (PBHs) [1, 2, 3, 4] has recently
attracted much attention thanks to the plethora of gravitational wave
detections from the mergers of BH binaries [5, 6, 7, 8]. Maybe even more
interesting, some of the LIGO/Virgo/KAGRA events might be in fact originating
from the mergers of PBHs [9, 10, 11, 12, 13, 14, 15, 16]. Future gravitational
wave experiments might help in shedding light on the possible existence of
PBHs [17, 18, 19, 20, 21].

In this paper we will start from the rather standard point of view that PBHs
are born in the radiation-dominated era by the collapse of large overdensities
generated on small scales during a non-attractor phase of inflation or through
a curvaton-like dynamics after inflation [1].

Calculating the abundance of PBHs as precisely as possible is a crucial step
in assessing if PBHs may account for a substantial fraction of the dark matter
of the universe or if they satisfy the different observational constraints.
Such computation is rendered difficult as PBHs are rare events and therefore
the formation probability is sensitive to small changes in the various
ingredients, such as the critical threshold of collapse [22, 23], the non-
Gaussian nature of the fluctuations [24, 25], the choice of the window
function to define smoothed observables [26], to mention a few. Furthermore it
has been recently pointed out that non-linear corrections entering in the
calculation of the PBHs abundance from the non-linear radiation transfer
function and the determination of the true physical horizon crossing corrects
the overdensity (on comoving orthogonal slices) by introducing many non-linear
terms, even of non-local type [27]. This makes the calculation of the
formation probability highly nontrivial.

In this paper we point out that the standard way of calculating the PBH
abundance, which we will revise in the next section, overlooks a technical,
but simple point. In the standard calculation one defines a smoothed density
contrast $\delta_{m}(r_{m},\vec{x}_{\rm pk})$, which depends on the position
$\vec{x}_{\rm pk}$ of the peak of the profile of the overdensity and on the
distance $r_{m}$ from it where the compaction function has its maximum. This
field turns out to be the sum of a field $\delta(r_{m},\vec{x}_{\rm pk})$, but
also its square $\delta^{2}(r_{m},\vec{x}_{\rm pk})$. The latter is a
composite operator, i.e., evaluated at the same point of space (and time), and
receives loop contributions from all scales. This happens because
$\delta(r_{m},\vec{x}_{\rm pk})$ is proportional to the gradient of the small-
scale curvature perturbation and the latter is typically non-Gaussian (see
e.g., [28]). In other words the curvature perturbation has cubic and quartic
interactions, leading to the renormalization of the composite operator
$\delta^{2}(r_{m},\vec{x}_{\rm pk})$ and therefore of
$\delta_{m}(r_{m},\vec{x}_{\rm pk})$.

The renormalization procedure leads to the well-known phenomenon of operator
mixing: the renormalized operator $\delta^{2}(r_{m},\vec{x}_{\rm pk})$
receives corrections from all possible operators compatible with the
symmetries of the problem. We will show that indeed one-loop corrections
assume the form of an infinite tower of operators, even non-local in
$\delta(r_{m},\vec{x}_{\rm pk})$ and of the higher-derivative type. It is the
interacting nature of the comoving curvature perturbation which gives rise to
the phenomenon of operator mixing, which leads to a shift in the correlators
and the threshold for PBH formation, as well as a modification of the
formation probability.

The paper is organized as follows. In section 2 we set the stage and describe
briefly the standard way of calculating the PBH abundance. In section 3 we
remind the reader why PBH formation models are characterized by a non-Gaussian
curvature perturbation. In section 4 we introduce the concept of
renormalization of the PBHs and the Feynman rules we will adopt to calculate
the operator mixing at one-loop, which we do in sections 5 and 6. Section 7
contains our conclusions, while we devote the appendices A and B to the
technical details.

##  2  A brief summary of the standard PBH abundance calculation

The goal of this section is to set the stage and briefly recap what is the
standard procedure to calculate the PBH abundance in the literature, see for
instance Refs. [25, 29].

As we already wrote, our focus will be the PBH formation from the collapse of
large overdensities which are generated during inflation and re-enter the
Hubble radius in the radiation-dominated era. The key starting object is the
curvature perturbation $\zeta$ appearing in the metric in the comoving
uniform-energy density gauge

${\rm d}s^{2}=-{\rm d}t^{2}+a^{2}(t)e^{2\zeta(\vec{x})}{\rm d}{\vec{x}}^{2}.$
(2.1)

Here $a(t)$ is the scale factor in terms of the cosmic time. Applying the
gradient expansion on superhoizon scales [30], one can relate the non-linear
density contrast $\delta_{\rm com}(\vec{x},t)$ on comoving orthogonal slicings
and the time independent curvature perturbation $\zeta(\vec{x})$ [31]

$\delta_{\rm
com}(\vec{x},t)=-\frac{8}{9}\frac{1}{a^{2}H^{2}}e^{-5\zeta(\vec{x})/2}\nabla^{2}e^{\zeta(\vec{x})/2},$
(2.2)

where $H$ is the Hubble rate. Cosmological perturbations may gravitationally
collapse to form a PBH depending upon the amplitude measured at the peak of
the compaction function, defined to be the mass excess compared to the
background value in a given radius [31, 22, 32]. On superhorizon scales and
assuming spherical symmetry, it reads

$\mathcal{C}(r)=-\frac{2}{3}r\partial_{r}\zeta(r)\left[2+r\partial_{r}\zeta(r)\right].$
(2.3)

The compaction function has a maximum at the comoving length scale $r_{m}$
such that

$\left.\partial_{r}\zeta(r)+r_{m}\partial_{r}^{2}\zeta(r)\right|_{r=r_{m}}=0,$
(2.4)

and one can define a smoothed perturbation amplitude as the volume average of
the energy density contrast within the scale $r_{m}$ at the Hubble radius
crossing time $t_{H}$ [22]

$\delta_{m}=\frac{3}{\left(r_{m}e^{\zeta(r_{m})}\right)^{3}}\int_{0}^{r_{m}}{\rm
d}r\,\delta_{\rm
com}(r,t_{H})\left(re^{\zeta(r)}\right)^{2}\partial_{r}\left(re^{\zeta(r)}\right),$
(2.5)

where a top-hat window function needs to be used for a proper treatment of the
threshold [26].

The quantity $\delta_{m}$ is crucial in determining the abundance of PBHs. If
computed at the Hubble radius crossing time $t_{H}$, that is

$\epsilon(t_{H})=\frac{r_{H}}{r_{m}e^{\zeta(r_{m})}}=\frac{1}{r_{m}e^{\zeta(r_{m})}\,aH}=1,$
(2.6)

it becomes

$\delta_{m}=\delta_{l}-\frac{3}{8}\delta_{l}^{2},\qquad\delta_{l}=-\frac{4}{3}r_{m}\left.\partial_{r}\zeta(r)\right|_{r=r_{m}}.$
(2.7)

The PBH abundance is subsequently calculated integrating the probability
distribution function of the smoothed density contrast from a threshold value
$\delta_{c}$ on

$\beta=\int_{\delta_{c}}{\rm d}\delta_{m}\left(\frac{M_{\rm
PBH}}{M_{H}}\right)P(\delta_{m}),$ (2.8)

as a function of the PBH mass $M_{\rm PBH}$ and the mass $M_{H}$ enclosed in
the Hubble volume. One can then use the conservation of probability to write
[33, 25, 29, 34, 28, 35]

$P(\delta_{l}){\rm d}\delta_{l}=P(\delta_{m}){\rm d}\delta_{m},$ (2.9)

with threshold given by [25]

$\delta_{l,c}=\frac{4}{3}\left(1-\sqrt{1-\frac{3}{2}\delta_{c}}\right),$
(2.10)

where $\delta_{c}$ is the threshold routinely computed by running numerical
simulations [22].

For Gaussian curvature perturbations, the probability of the linear density
contrast is Gaussian and is exponentially sensitive to the threshold
$\delta_{l,c}$ and the variance $\sigma^{2}_{l}$

$\displaystyle P(\delta_{l})$
$\displaystyle=\frac{1}{\sqrt{2\pi}\sigma_{l}}e^{-\delta_{l}^{2}/2\sigma_{l}^{2}},$
$\displaystyle\sigma_{l}^{2}$ $\displaystyle=\int\frac{{\rm d}k}{k}{\cal
P}_{\delta_{l}}(k),$ (2.11)

where ${\cal P}_{\delta_{l}}$ is the dimensionless power spectrum of the
density contrast $\delta_{l}$.

##  3 The non-Gaussianity of the curvature perturbation

The probability (2) is valid only if the smoothed density contrast
$\delta_{l}$ is linear (hence the l subscript). Given the relation

$\delta_{l}(r_{m})=-\frac{4}{3}r_{m}\left.\partial_{r}\zeta(r)\right|_{r=r_{m}},$
(3.1)

one concludes that $\delta_{l}$ is indeed linear only if the curvature
perturbation is linear. However, in all models in which a large fluctuation in
the curvature perturbation is generated, the curvature perturbation is
typically non-Gaussian (see, for instance, Ref. [28])111We remark that the
non-Gaussianity we are discussing here is induced at small scales and
therefore is not the non-Gaussianity of the curvature perturbation on CMB
scales [36], which are much larger than those interested in the PBH formation.
On CMB scales the non-Gaussianity is severely constrained by CMB data [37]..
Non-Gaussianity among the modes interested in the growth of the curvature
perturbation is generated either by their self-interaction during the ultra
slow-roll phase [38] or after Hubble radius exit when the curvature
perturbation is sourced by a curvaton-like field [39, 40].

Even though in general the exact relation between $\zeta$ and its
corresponding Gaussian component $\zeta_{g}$ can be worked out model by model,
we will adopt a cubic expansion to be model independent,

$\zeta=\zeta_{g}+\frac{3}{5}f_{NL}\zeta_{g}^{2}+\frac{9}{25}g_{NL}\zeta_{g}^{3},$
(3.2)

which is a good approximation if $(3/5)f_{NL}\zeta_{g}$ and
$(9/25)g_{NL}\zeta^{2}_{g}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}1$. For instance, in models with a
sharp transition between the ultra slow-roll phase and the subsequent slow-
roll phase one has $f_{NL}\simeq 5/2$ and $g_{NL}\simeq 25/6$ (see, for
instance, Ref. [38]) and therefore the expansion is motivated as long as
$\zeta_{g}\mathrel{\hbox to0.0pt{\lower 4.0pt\hbox{\hskip
0.5pt$\sim$}\hss}\raise 1.0pt\hbox{$<$}}0.7$. The expansion is also justified
if we perform a perturbative loop expansion and think the non-Gaussianity as
leading to three-, four-point vertices, and so on.

The non-Gaussianity of the curvature perturbation is crucial in changing the
PBH abundance. Indeed, around peaks of the power spectrum of the curvature
perturbation, a positive $f_{NL}$ increases the abundance of the PBHs, while a
negative $f_{NL}$ has the opposite effect, thus helping to interpret the
recent pulsar timing array observation of a stochastic gravitational wave
background [41, 42, 43, 44, 45, 46, 47, 48, 49] by using the scalar-induced
gravitational waves sourced along with PBH formation [50, 51, 52, 53, 54, 55,
56, 57, 58].

The importance of the non-Gaussianity in the curvature perturbation will
become clear in the following section, where we will go back to the smoothed
density contrast $\delta_{m}$ and deal with the fact that it is a composite
operator.

##  4 Renormalization in the PBHs

We start again from the expression (2.7) of the smoothed density contrast

$\delta_{m}(r_{m},\vec{x}_{\rm pk})=\frac{3}{4\pi r_{m}^{3}}\int{\rm
d}^{3}x\,\delta_{\rm com}(\vec{x},t_{H})\,\theta(r_{m}-|\vec{x}-\vec{x}_{\rm
pk}|)=\delta(r_{m},\vec{x}_{\rm pk})-\frac{3}{8}\delta^{2}(r_{m},\vec{x}_{\rm
pk}),$ (4.1)

where this time we removed the subscript l to account for the fact that
$\delta_{l}$ is not a Gaussian field, since the curvature perturbation is not.
Furthermore we have highlighted the fact that $\delta_{m}$ depends on the peak
position $\vec{x}_{\rm pk}$ and from the distance $r_{m}$ from the peak.

Figure 1: Schematic view of operations the leading to the composite operator
renormalization.

The key point is that $\delta_{m}(r_{m},\vec{x}_{\rm pk})$ contains the
composite operator $\delta^{2}(r_{m},\vec{x}_{\rm pk})$ which receives
contributions from all scales, not only from the infrared ones. This is
because two small-scale modes may combine in the loops to form a long mode.
Such contributions may be therefore not small at all, even on large scales.
One should also recall at this stage that, even in the presence of the window
function cutting-off ultraviolet momenta, variances are typically dominated by
integrating momenta corresponding to length scales well beneath the horizon
[27].

One needs to renormalize the composite operator $\delta^{2}(r_{m},\vec{x}_{\rm
pk})$ at the point $\vec{x}_{\rm pk}$ in such a way that all potentially large
contributions can be systematically removed by adding local counterterms
(canceling as well the dependence on an unphysical cut-off scale, if needed)

$\left[\delta^{2}(r_{m},\vec{x}_{\rm
pk})\right]\equiv\delta^{2}(r_{m},\vec{x}_{\rm pk})+\sum_{\cal O}c_{\cal
O}{\cal O}(r_{m},\vec{x}_{\rm pk}),$ (4.2)

where the square brackets indicate the renormalized operator. This means that
the renormalized composite operator will mix with possibly all the other
operators ${\cal O}(r_{m},\vec{x}_{\rm pk})$ which are allowed by the
symmetries of the problem. It will also have an impact on the PBH abundance,
being it so sensitive to minute changes in the threshold and/or in the
variances.

The goal of the subsequent sections is to perform the systematic
renormalization of the smoothed density contrast.

###  4.1 Feynman rules

Our starting point is the expression

$\delta(r_{m},\vec{x}_{\rm
pk})=-\frac{4}{3}r\partial_{r}\zeta(r_{m},\vec{x}_{\rm
pk})=-\frac{4}{9}\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}\;e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}k^{2}r_{m}^{2}W_{3}(kr_{m})\zeta_{\vec{k}},$ (4.3)

whose Fourier transform reads

$\delta_{\vec{k}}(r_{m})=-\frac{4}{9}k^{2}r_{m}^{2}W_{3}(kr_{m})\zeta_{\vec{k}}.$
(4.4)

Here $W_{3}(kr_{m}$) is the Fourier transform of the top-hat window function
in real space with radius $r_{m}$

$W_{3}(x)=3\frac{\sin x-x\cos
x}{x^{3}}=2^{3/2}\Gamma(5/2)\frac{J_{3/2}(x)}{x^{3/2}},$ (4.5)

being $J_{\nu}(x)$ the Bessel function of the first kind. To simplify the
notation we have absorbed the time-dependent radiation transfer function in
the curvature perturbation $\zeta_{\vec{k}}$ (and in its power spectrum). All
power spectra will be intended to be calculated at horizon crossing.

The renormalization of the composite operator $\delta^{2}(r_{m},\vec{x}_{\rm
pk})$ will proceed through the following Feynman rules:

  * •

draw every $N$-point connected graph and conserve momentum in each vertex;

  * •

to each $(n+1)$ point vertex, with a solid line (in-going momentum) and $n$
dashed lines, it corresponds the following vertex

$\displaystyle=(-1)^{n-1}n!\,F_{n}\;\;\frac{(qr_{m})^{2}W_{3}(qr_{m})}{\left(\frac{4}{9}\right)^{n-1}\prod_{i=1}^{n}(q_{i}r_{m})^{2}W_{3}(q_{i}r_{m})},$
(4.6)

where in the numerator there is the factor $(4/9)(qr_{m})^{2}W_{3}(qr_{m})$
for the solid line and at the denominator a product of all the factors
$(4/9)(q_{i}r_{m})^{2}W_{3}(q_{i}r_{m})$ for each dashed lines;

  * •

the propagator, represented with a dashed line, corresponds to

$\parbox[c]{0.0pt}{\includegraphics[scale={0.9}]{Figures/D_propagator_t.pdf}}\
=P_{\delta}(q)=\frac{16}{81}r_{m}^{4}q^{4}\,W^{2}_{3}(qr_{m})P_{\zeta_{g}}(q);$
(4.7)

  * •

integrate for every loop

$\int_{\vec{q}}=\int\frac{d^{3}q}{(2\pi)^{3}};$ (4.8)

  * •

multiply a Dirac delta with all the external momenta, for the global momentum
conservation

$(2\pi)^{3}\delta^{(3)}\left(\sum_{i=1}^{N}\vec{k}_{i}\right);$ (4.9)

Here $F_{n}$ is a function that depends on the number of legs in the vertex,
for example from Eq. (3.2)

$F_{1}=1,\quad F_{2}=\frac{3}{5}f_{NL},\quad F_{3}=\frac{9}{25}g_{NL}.$ (4.10)

Notice that for the $n=1$ vertex we have

$\parbox[c]{0.0pt}{\includegraphics[scale={0.9}]{Figures/D_2-point_t.pdf}}=\frac{(qr_{m})^{2}W_{3}(qr_{m})}{(q_{1}r_{m})^{2}W_{3}(q_{1}r_{m})}=1,$
(4.11)

which is just a consequence of momentum conservation.

##  5 Renormalization at one-loop: linear operator mixing

By indicating the composite operator with the following vertex

$\delta^{2}=\parbox[c]{0.0pt}{\includegraphics[scale={0.8}]{Figures/D_composite_t.pdf}}$
(5.1)

its expectation value is obtained from the Feynman diagram

$\parbox[c]{0.0pt}{\includegraphics[scale={0.8}]{Figures/D_loop_0legs_t.pdf}}\
=(2\pi)^{3}\delta^{(3)}(\vec{k})\int_{\vec{q}}P_{\delta}(q),$ (5.2)

which gives

$\langle\left(\delta^{2}\right)_{\vec{k}}\rangle^{\prime}=\int_{\vec{q}}P_{\delta}(q)=\frac{16}{81}r_{m}^{4}\int_{\vec{q}}q^{4}\,W^{2}_{3}(qr_{m})P_{\zeta_{g}}(q)=\sigma^{2}_{\delta(r_{m})}.$
(5.3)

Here and in the following primes indicate that we remove the factors of
$(2\pi)^{3}$ times the Dirac delta for the momentum conservation. This vacuum
expectation value is removed by adding a constant counterterm

$\left[\delta^{2}(r_{m},\vec{x}_{\rm pk})\right]=\delta^{2}(r_{m},\vec{x}_{\rm
pk})-\sigma^{2}_{\delta(r_{m})}.$ (5.4)

Subtracting this tadpole contribution ensures that the vacuum expectation
value of $\delta_{m}$ vanishes at the loop level. Furthermore, in the
following calculations this constant counterterm will only contribute to
disconnected graphs.

The one-loop contribution to
$\langle(\delta^{2})_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\rangle$ is,
diagramatically222We will ignore from now on diagrams such as  which does not
have power spectra of external momenta and lead to contact terms, with no new
operators.

$\langle(\delta^{2})_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\rangle^{\prime}=2\times\parbox[c]{0.0pt}{\includegraphics[scale={0.8}]{Figures/D_loop_1leg_t.pdf}}\;,$
(5.5)

which gives

$\displaystyle\langle(\delta^{2})_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\rangle^{\prime}$
$\displaystyle=$
$\displaystyle-\frac{12}{5}f_{NL}\int_{\vec{q}}\;\frac{(|\vec{k}_{1}-\vec{q}|r_{m})^{2}W_{3}(|\vec{k}_{1}-\vec{q}|r_{m})}{\frac{4}{9}(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})(qr_{m})^{2}W_{3}(qr_{m})}\;P_{\delta}(k_{2})P_{\delta}(q)$
$\displaystyle=$
$\displaystyle-\frac{12}{5}\cdot\frac{4}{9}\cdot\frac{16}{81}f_{NL}\int_{\vec{q}}r_{m}^{6}q^{2}|\vec{k}_{1}-\vec{q}|^{2}k_{2}^{2}\;W_{3}\left(|\vec{k}_{1}-\vec{q}|r_{m}\right)W_{3}(k_{2}r_{m})W_{3}(qr_{m})P_{\zeta_{g}}(q)P_{\zeta_{g}}(k_{2}).$

We now use Eq. (A.1) to perform the angular integral and obtain

$\displaystyle\langle(\delta^{2})_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\rangle^{\prime}$
$\displaystyle=$ $\displaystyle I^{(1)}_{\rm one-leg}+I^{(2)}_{\rm one-leg},$
$\displaystyle I^{(1)}_{\rm one-leg}$ $\displaystyle=$
$\displaystyle-\frac{256}{1215}f_{NL}\int_{\vec{q}}r_{m}^{6}q^{2}(q^{2}+k_{2}^{2})k_{2}^{2}\;W_{3}^{2}(k_{2}r_{m})W^{2}_{3}(qr_{m})P_{\zeta_{g}}(q)P_{\zeta_{g}}(k_{2})$
$\displaystyle I^{(2)}_{\rm one-leg}$ $\displaystyle=$
$\displaystyle-\frac{256}{1215}f_{NL}\int_{\vec{q}}r_{m}^{6}q^{2}k_{2}^{2}\;W^{2}_{3}(k_{2}r_{m})W^{2}_{3}(qr_{m})\left(\frac{q^{2}}{3}\frac{{\rm
d}\ln W_{3}(k_{2}r_{m})}{{\rm d}\ln
k_{2}r_{m}}+\;\frac{k_{2}^{2}}{3}\frac{{\rm d}\ln W_{3}(qr_{m})}{{\rm d}\ln
qr_{m}}\right)$ $\displaystyle\cdot$ $\displaystyle
P_{\zeta_{g}}(q)P_{\zeta_{g}}(k_{2})\;.$

The first integral gives

$I^{(1)}_{\rm one-
leg}=\frac{12}{5}f_{NL}\left(\sigma^{2}_{\delta(r_{m})}\langle\zeta_{g\vec{k}_{1}}(r_{m})\delta_{\vec{k}_{2}}(r_{m})\rangle^{\prime}+\frac{4}{9}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}\langle\delta_{\vec{k}_{1}}(r_{m})\delta_{\vec{k}_{2}}(r_{m})\rangle^{\prime}\right),$
(5.8)

where

$\zeta_{g}(r_{m},\vec{x}_{\rm pk})=\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}\,e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}\zeta_{g\vec{k}}(r_{m})=\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}\,e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}W_{3}(kr_{m})\zeta_{g\vec{k}}$ (5.9)

and

$\sigma^{2}_{\nabla\zeta_{g}(r_{m})}=\int_{\vec{q}}r_{m}^{2}q^{2}W_{3}^{2}(qr_{m})P_{\zeta_{g}}(q).$
(5.10)

The second integral gives

$\displaystyle I^{(2)}_{\rm one-
leg}=\frac{12}{5}f_{NL}\left[\frac{1}{3}\sigma^{2}_{\delta(r_{m})}\frac{{\rm
d}}{{\rm d}\ln
r_{m}}\langle\zeta_{g\vec{k}_{1}}(r_{m})\delta_{\vec{k}_{2}}(r_{m})\rangle^{\prime}+\frac{4}{9}\cdot\left(\frac{1}{6}\frac{{\rm
d}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}}{{\rm d}\ln
r_{m}}-\frac{1}{3}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}\right)\langle\delta_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\rangle^{\prime}\right].$

Therefore, the renormalized operator at the level of linear mixing becomes
(see Appendix B)

$\displaystyle\left[\delta^{2}(r_{m},\vec{x}_{\rm pk})\right]_{\text{linear
}}$ $\displaystyle=$ $\displaystyle\delta^{2}(r_{m},\vec{x}_{\rm
pk})-\sigma^{2}_{\delta(r_{m})}-\frac{12}{5}f_{NL}\cdot\frac{4}{9}\left(\frac{2}{3}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}+\frac{1}{6}\frac{{\rm
d}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}}{{\rm d}\ln
r_{m}}\right)\delta(r_{m},\vec{x}_{\rm pk})$ $\displaystyle+$
$\displaystyle\frac{12}{5}f_{NL}\cdot\frac{5}{3}\cdot\sigma^{2}_{\delta(r_{m})}\zeta(r_{m},\vec{x}_{\rm
pk}),$ $\displaystyle=$ $\displaystyle\hskip
14.22636pt\parbox[c]{0.0pt}{\includegraphics[scale={0.8}]{Figures/D_loop_1leg_t.pdf}}\hskip
14.22636pt+\hskip
14.22636pt\parbox[c]{0.0pt}{\includegraphics[scale={1}]{Figures/D_counter1_t.pdf}}$
(5.12)

where we have added the corresponding Feynman diagram for the counterterm. We
see that already at this level there appears a non-local operator in the
density smoothed contrast field $\delta(r_{m},\vec{x}_{\rm pk})$ once the
following identity

$\zeta_{g}(r_{m},\vec{x}_{\rm pk})=\frac{9}{4r_{m}^{2}\nabla^{2}_{x_{\rm
pk}}}\delta(r_{m},\vec{x}_{\rm pk})$ (5.13)

is adopted.

##  6 Renormalization at one-loop: quadratic operator mixing

The one-loop contribution to
$\langle[\delta^{2}]_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\delta_{\vec{k}_{3}}\rangle$
is diagramatically

$\langle[\delta^{2}]_{\vec{k}_{1}}\delta_{\vec{k}_{2}}\delta_{\vec{k}_{3}}\rangle^{\prime}\hskip
8.5359pt=\hskip
8.5359pt\parbox[c]{0.0pt}{\includegraphics[scale={0.8}]{Figures/D_loop_2legs_1_t.pdf}}\hskip
8.5359pt+\hskip
8.5359pt2\times\parbox[c]{0.0pt}{\includegraphics[scale={0.8}]{Figures/D_loop_2legs_2_t.pdf}}\hskip
8.5359pt+\hskip
8.5359pt\parbox[c]{0.0pt}{\includegraphics[scale={0.85}]{Figures/D_counter2_t.pdf}}$
(6.1)

Let us start with the first diagram, which gives

$\displaystyle I^{(1)}_{\rm two-legs}$ $\displaystyle=$
$\displaystyle\frac{36}{25}f^{2}_{NL}P_{\delta}(k_{2})P_{\delta}(k_{3})\left(\frac{9}{4}\right)^{2}\int_{\vec{q}}P_{\delta}(q)\frac{(|\vec{k}_{3}+\vec{q}|r_{m})^{2}W_{3}(|\vec{k}_{3}+\vec{q}|r_{m})}{(qr_{m})^{2}W_{3}(qr_{m})(k_{3}r_{m})^{2}W_{3}(k_{3}r_{m})}$
(6.2) $\displaystyle\cdot$
$\displaystyle\frac{(|\vec{k}_{2}-\vec{q}|r_{m})^{2}\;W_{3}(|\vec{k}_{2}-\vec{q}|r_{m})}{(qr_{m})^{2}W_{3}(qr_{m})(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})}.$

Expanding the power spectra we get

$\displaystyle I^{(1)}_{\rm two-legs}$ $\displaystyle=$
$\displaystyle\frac{36}{25}f^{2}_{NL}\left(\frac{16}{81}\right)^{2}(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})P_{\zeta_{g}}(k_{2})(k_{3}r_{m})^{2}W_{3}(k_{3}r_{m})P_{\zeta_{g}}(k_{3})$
(6.3) $\displaystyle\cdot$
$\displaystyle\int_{\vec{q}}P_{\zeta_{g}}(q)\;(|\vec{k}_{3}+\vec{q}|r_{m})^{2}\;W_{3}(|\vec{k}_{3}+\vec{q}|r_{m})(|\vec{k}_{2}-\vec{q}|r_{m})^{2}\;W_{3}(|\vec{k}_{2}-\vec{q}|r_{m}).$

The calculation is rather involved and we do not report it here in full
length, the interested reader can find all the details in the Appendix A.2.
The procedure once again is to disentangle the sum of the vectors in the
window functions using the addition theorem of the Bessel functions. The
integral (6.3) is then written as

$\displaystyle I^{(1)}_{\rm two-legs}$ $\displaystyle=$
$\displaystyle\frac{36}{25}f^{2}_{NL}\left(\frac{4}{9}\right)^{4}(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})P_{\zeta_{g}}(k_{2})(k_{3}r_{m})^{2}W_{3}(k_{3}r_{m})P_{\zeta_{g}}(k_{3})$
$\displaystyle\cdot$ $\displaystyle
r_{m}^{4}\int_{\vec{q}}P_{\zeta_{g}}(q)\;\left[k_{2}^{2}k_{3}^{2}I^{00}(\cos\theta_{23})+2qk_{3}k_{2}^{2}I^{01}(\cos\theta_{23})+2qk_{2}k_{3}^{2}I^{10}(\cos\theta_{23})\right.$
$\displaystyle+$
$\displaystyle\left.q^{2}(k_{2}^{2}+k_{3}^{2})I^{00}(\cos\theta_{23})+4q^{2}k_{2}k_{3}I^{11}(\cos\theta_{23})+2q^{3}k_{3}I^{01}(\cos\theta_{23})\right.$
$\displaystyle+$
$\displaystyle\left.2q^{3}k_{2}I^{10}(\cos\theta_{23})+q^{4}I^{00}(\cos\theta_{23})\right],$

where $\cos\theta_{23}=\vec{k}_{2}\cdot\vec{k}_{3}/k_{2}k_{3}$,

$\displaystyle
I^{ij}(\cos\theta_{23})=\sum_{k=0}^{\infty}\sum_{s=0}^{\infty}N_{ks}I^{ij}_{ks}(\cos\theta_{23}),$
(6.5)

$\displaystyle N_{ks}$
$\displaystyle=\frac{4}{9}\left(\frac{3}{2}+k\right)\left(\frac{3}{2}+s\right)\left[(r_{m}q)^{k}\left(\frac{1}{(r_{m}q)}\frac{{\rm
d}}{{\rm
d}(r_{m}q)}\right)^{k}W_{3}(r_{m}q)\right]\left[(r_{m}q)^{s}\left(\frac{1}{(r_{m}q)}\frac{{\rm
d}}{{\rm d}(r_{m}q)}\right)^{s}W_{3}(r_{m}q)\right]$
$\displaystyle~{}~{}\cdot\left[(r_{m}k_{2})^{k}\left(\frac{1}{(r_{m}k_{2})}\frac{{\rm
d}}{{\rm
d}(r_{m}k_{2})}\right)^{k}W_{3}(r_{m}k_{2})\right]\left[(r_{m}k_{3})^{s}\left(\frac{1}{(r_{m}k_{3})}\frac{{\rm
d}}{{\rm d}(r_{m}k_{3})}\right)^{s}W_{3}(r_{m}k_{3})\right],$ (6.6)

and

$\displaystyle I^{00}_{ks}$ $\displaystyle=$
$\displaystyle\left(\frac{1+(-1)^{k+s}}{2}\right)\dfrac{{\rm d}}{{\rm
d}\cos\theta_{23}}P_{{\rm min}(k,s)+1}(\cos\theta_{23})$ $\displaystyle
I^{10}_{ks}$ $\displaystyle=$
$\displaystyle\left(\frac{1-(-1)^{k+s}}{2}\right)\left[\Theta(k-s)\frac{\rm
d}{{\rm
d}\cos\theta_{23}}P_{k+1}(\cos\theta_{23})+\Theta(s-k)\cos\theta_{23}\frac{\rm
d}{{\rm d}\cos\theta_{23}}P_{s+1}(\cos\theta_{23})\right],$ $\displaystyle
I^{01}_{ks}$ $\displaystyle=$
$\displaystyle\left(\frac{1-(-1)^{k+s}}{2}\right)\left[\Theta(s-k)\frac{\rm
d}{{\rm
d}\cos\theta_{23}}P_{k+1}(\cos\theta_{23})+\Theta(k-s)\cos\theta_{23}\frac{\rm
d}{{\rm d}\cos\theta_{23}}P_{s+1}(\cos\theta_{23})\right],$ $\displaystyle
I^{11}_{ks}$ $\displaystyle=$ $\displaystyle
I^{01}_{k+1,s}-(k+2)\left\\{\left(\frac{1+(-1)^{k+s}}{2}\right)\Theta(s-k-1)P_{k+1}(\cos\theta_{23})+\frac{s+1}{2s+3}\delta_{sk}P_{s+1}(\cos\theta_{23})\right\\}.$

This will contribute to an infinite tower of operators to subtract. The second
diagram in Eq. (6.1) gives

$\displaystyle I^{(2)}_{\rm two-legs}$ $\displaystyle=$ $\displaystyle
g_{NL}\frac{108}{25}\left(\frac{9}{4}\right)^{2}P_{\delta}(k_{3})P_{\delta}(k_{2})\int_{\vec{q}}P_{\delta}(q)\frac{(|\vec{k}_{1}-\vec{q}|r_{m})^{2}\;W_{3}(|\vec{k}_{1}-\vec{q}|r_{m})}{(qr_{m})^{2}W_{3}(qr_{m})(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})(k_{3}r_{m})^{2}W_{3}(k_{3}r_{m})}$
(6.8) $\displaystyle=$ $\displaystyle
g_{NL}\frac{108}{25}\left(\frac{16}{81}\right)^{2}(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})P_{\zeta_{g}}(k_{2})(k_{3}r_{m})^{2}W_{3}(k_{3}r_{m})P_{\zeta_{g}}(k_{3})$
$\displaystyle\cdot$
$\displaystyle\int_{\vec{q}}(qr_{m})^{2}W_{3}(qr_{m})P_{\zeta_{g}}(q)\;(|\vec{k}_{1}-\vec{q}|r_{m})^{2}\;W_{3}(|\vec{k}_{1}-\vec{q}|r_{m}).$

Using again Eq. (A.1), we can perform the angular integration, leading to

$\displaystyle I^{(2)}_{\rm two-legs}$ $\displaystyle=$ $\displaystyle
g_{NL}\frac{108}{25}\left(\frac{16}{81}\right)^{2}(k_{2}r_{m})^{2}W_{3}(k_{2}r_{m})P_{\zeta_{g}}(k_{2})(k_{3}r_{m})^{2}W_{3}(k_{3}r_{m})P_{\zeta_{g}}(k_{3})$
(6.9) $\displaystyle\cdot$
$\displaystyle\int_{\vec{q}}(qr_{m})^{2}W_{3}^{2}(qr_{m})P_{\zeta_{g}}(q)\;W_{3}(k_{1}r_{m})\left[(qr_{m})^{2}+(k_{1}r_{m})^{2}+\frac{(qr_{m})^{2}}{3}\frac{{\rm
d}\ln W_{3}(k_{1}r_{m})}{{\rm d}\ln k_{1}r_{m}}\right.$ $\displaystyle+$
$\displaystyle\left.\frac{(k_{1}r_{m})^{2}}{3}\frac{{\rm d}\ln
W_{3}(qr_{m})}{{\rm d}\ln qr_{m}}\right].$

If we define the composite operator

$\zeta^{2}(r_{m},\vec{x}_{\rm pk})=\int\frac{{\rm
d}^{3}k}{(2\pi)^{3}}\,e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}W_{3}(kr_{m})(\zeta^{2}_{g})_{\vec{k}},$ (6.10)

the operators to be added at this stage are

$-\frac{9}{25}\cdot\frac{16}{81}\cdot
6\;g_{NL}\left(\frac{2}{3}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}+\frac{1}{6}\frac{{\rm
d}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}}{{\rm d}\ln
r_{m}}\right)r_{m}^{2}\nabla^{2}_{x_{\rm pk}}\zeta^{2}(r_{m},\vec{x}_{\rm
pk})-\frac{9}{25}\cdot
6\;g_{NL}\;\sigma^{2}_{\delta(r_{m})}\left(1+\frac{1}{3}\frac{{\rm d}}{{\rm
d}\ln r_{m}}\right)\zeta^{2}(r_{m},\vec{x}_{\rm pk}).$ (6.11)

Finally, from the third diagram in Eq. (6.1) the operators to add are
trivially

$\frac{3}{5}f_{NL}\cdot\frac{4}{9}\left[-\frac{12}{5}f_{NL}\cdot\frac{4}{9}\left(\frac{2}{3}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}+\frac{1}{6}\frac{{\rm
d}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}}{{\rm d}\ln
r_{m}}\right)r_{m}^{2}\nabla^{2}_{x_{\rm pk}}\zeta^{2}(r_{m},\vec{x}_{\rm
pk})+\frac{12}{5}f_{NL}\cdot\frac{5}{3}\cdot\sigma^{2}_{\delta(r_{m})}\zeta^{2}(r_{m},\vec{x}_{\rm
pk})\right].$ (6.12)

Summing up all the contributions from $I^{(2)}_{\rm two-legs}$ and
$I^{(3)}_{\rm two-legs}$, and writing only the first terms for $I^{(1)}_{\rm
two-legs}$, we obtain at the level of quadratic mixing

$\displaystyle\left[\delta^{2}(r_{m},\vec{x}_{\rm pk})\right]_{\text{quadratic
}}$ $\displaystyle=$ $\displaystyle\delta^{2}(r_{m},\vec{x}_{\rm
pk})-\frac{18}{25}f_{NL}^{2}\sigma^{2}_{\zeta_{g}(r_{m})}\delta^{2}(r_{m},\vec{x}_{\rm
pk})$ (6.13) $\displaystyle+$
$\displaystyle\frac{16}{25}f_{NL}^{2}\sigma^{2}_{\nabla\zeta(r_{m})}\zeta(r_{m},\vec{x}_{\rm
pk})\delta(r_{m},\vec{x}_{\rm pk})$ $\displaystyle+$
$\displaystyle\frac{256}{675}f_{NL}^{2}\sigma^{2}_{\nabla\zeta(r_{m})}\nabla_{x_{\rm
pk}^{i}}\zeta(r_{m},\vec{x}_{\rm pk})\nabla^{x_{\rm
pk}^{i}}\zeta(r_{m},\vec{x}_{\rm pk})$ $\displaystyle-$
$\displaystyle\left(\frac{32}{75}f_{NL}^{2}+\frac{64}{225}g_{NL}\right)\left(\frac{2}{3}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}+\frac{1}{6}\frac{{\rm
d}\sigma^{2}_{\nabla\zeta_{g}(r_{m})}}{{\rm d}\ln
r_{m}}\right)r_{m}^{2}\nabla^{2}_{x_{\rm pk}}\zeta^{2}(r_{m},\vec{x}_{\rm
pk})$ $\displaystyle-$
$\displaystyle\left[\frac{54}{25}g_{NL}\;\sigma^{2}_{\delta(r_{m})}\left(1+\frac{1}{3}\frac{{\rm
d}}{{\rm d}\ln
r_{m}}\right)-\frac{26}{75}f^{2}_{NL}\;\sigma^{2}_{\delta(r_{m})}\right]\zeta^{2}(r_{m},\vec{x}_{\rm
pk})$ $\displaystyle+$ $\displaystyle\cdots.$

Eqs. (5.12) and (6.13) are the main results of this paper. In the last
expression, we have indicated by the dots the infinite series of operators
with higher-derivatives applied either to the operators or to the window
functions.

To evaluate the impact of the renormalization and operator mixing, one has to
recall that, calling $k_{\star}$ the typical momentum at which PBHs form, one
has typically (at the linear level) $k_{\star}r_{m}\simeq{\cal O}(3)$ [59,
22]. For instance, for a monochromatic curvature spectrum peaked at
$k_{\star}$, one has $k_{\star}r_{m}\simeq 2.7$. For a broad spectrum
$k_{\star}$ coincides with the maximum momentum scale, as the PBH mass
function peaks at that scale [60], and $k_{\star}r_{m}\simeq 3.5$. One would
therefore expect that the higher-derivative terms are suppressed by powers of
$1/k_{\star}r_{m}$, but having an infinite amount of terms one cannot exclude
a priori that the infinite sum of operators give a sizable contribution. We
will show some numerical insights in the next subsection.

###  6.1 The impact of renormalization

We consider the case of a peaked power spectrum in the curvature perturbation

$\frac{k^{3}}{2\pi^{2}}P_{\zeta_{g}}(k)=A_{\zeta}k_{\star}\delta(k-k_{\star}),$
(6.14)

where the amplitude $A_{\zeta}\simeq 1.8\cdot 10^{-2}$ is fixed by setting a
PBH abundance $\beta\simeq 10^{-10}$, neglecting composite operator
renormalization. With this assumption, the integrand in Eq. (6.3) greatly
simplifies as all momenta are fixed to $k_{\star}$. To study the convergence
of the result as a function of the sum in $k$ and $s$, we have plotted
$N_{ks}$ in Fig. 2 as a function of index $k$ (and analogous scaling is
obtained varying $s$).

The strong hierarchy observed between the different orders in the sum (6.5)
results in a rapidly converging series. We can further investigate this by
computing the angular integration

${\cal I}_{3}=4\int{\rm
d}\theta_{23}\left[I^{00}(\cos\theta_{23})+I^{01}(\cos\theta_{23})+I^{10}(\cos\theta_{23})+I^{11}(\cos\theta_{23})\right],$
(6.15)

obtained from (6.5) in the limit of narrow power spectrum. We show the result
as a function of the maximum index included in the series in Fig. 2. We see
that already at order $k_{\text{\tiny max}}=s_{\text{\tiny max}}={\cal O}(3)$
the series converges towards the asymptotic value within sufficient accuracy.

Figure 2:  Left panel: Coefficients $N_{ks}$ as a function of k, fixing $s=0$.
The red dashed line report a scaling $N_{k0}\propto 1/(k_{*}r_{m})^{3k}$.
Right panel: Values of ${\cal I}_{3}$ as a function of the maximum index
included in the sum (6.5).

Correction to the density variance. The effect of the renormalization of the
composite operator can be quantified as follows. We first consider its impact
on the variance of the field. With a narrow power spectrum, the Gaussian
component of the curvature perturbation is characterised by a mean peak
profile $\zeta_{g}(r)=\zeta_{\rm pk}\sin(k_{*}r)/(k_{*}r)$. We fix the
amplitude of the curvature perturbation in such a way that it corresponds to a
realisation of density contrast with the expectation value
$\delta_{l}\approx\sigma_{\delta(r_{m})}$. We find

$\left[\delta^{2}(r_{m},\vec{x}_{\rm
pk})\right]=0.022+0.015\left(\frac{f_{NL}}{5/2}\right)+0.019\left(\frac{f_{NL}}{5/2}\right)^{2}-0.0016\left(\frac{g_{NL}}{25/6}\right),$
(6.16)

which should be compared with the bare quantity
$\delta^{2}(r_{m})\approx\sigma_{\delta(r_{m})}^{2}=0.022$. This propagates on
the smoothed density contrast, which can be estimated as

$\delta_{m}(r_{m},\vec{x}_{\rm
pk})=0.14+0.00091\left(\frac{f_{NL}}{5/2}\right)-0.0070\left(\frac{f_{NL}}{5/2}\right)^{2}+0.00075\left(\frac{g_{NL}}{25/6}\right).$
(6.17)

Correction to the threshold. In order to evaluate the correction to the PBH
threshold for collapse, one should estimate the renormalization of the
composite operators assuming the fields amplitude reach the threshold
$\delta_{c}$. We fix the amplitude of the Gaussian curvature perturbation
$\zeta_{g}$ assuming it leads to a collapsing overdensity with threshold
amplitude $\delta_{m}=\delta_{c}$, where for a narrow spectrum we take
$\delta_{c}=0.59$ [23]. This corresponds to $\delta_{l}=0.88$ and $\zeta_{\rm
pk}=0.68$. In this case, we find

$\left[\delta^{2}(r_{m},\vec{x}_{\rm
pk})\right]=0.78+0.091\left(\frac{f_{NL}}{5/2}\right)+0.66\left(\frac{f_{NL}}{5/2}\right)^{2}-0.031\left(\frac{g_{NL}}{25/6}\right),$
(6.18)

while its bare value is $\delta^{2}(r_{m})\approx 0.78$. In terms of smoothed
density contrast, the renormalisation of the quadratic operators gives

$\delta_{m}(r_{m},\vec{x}_{\rm
pk})=0.59+0.20\left(\frac{f_{NL}}{5/2}\right)-0.25\left(\frac{f_{NL}}{5/2}\right)^{2}+0.020\left(\frac{g_{NL}}{25/6}\right).$
(6.19)

##  7 Comments and conclusions

We have shown that in the renormalization of the quadratic power of the
smoothed density contrast, a composite operator entering in the calculation of
the PBH abundance, leads to the well-known phenomenon of operator mixing due
to the non-Gaussian nature of the curvature perturbation. The mixing gives to
an infinite tower of operators due to the necessary operation of smoothing out
with the top-hat window function in real space. There are some comments we can
offer at this stage:

  1. 1.

The calculation of the PBH abundance is not as straightforward as standardly
assumed. Not only non-linear corrections entering in the calculation of the
PBH from the non-linear radiation transfer function and the determination of
the true physical horizon crossing are important and lead to large
uncertainties in the final result [27], but also the renormalization procedure
leads to an infinite tower of operators, making the calculation of the
formation probability of the PBH a difficult task. Our findings indicate a
large impact on the threshold for PBH formation, depending on the sign of the
non-Gaussianity parameters. Since in ultra-slow inflation the sign of the
parameter $f_{NL}$ is always positive [61], our results indicate that the
threshold might in fact decrease due to quadratic corrections ${\cal
O}(f_{NL}^{2})$, thus increasing the PBH abundance and reversing the ${\cal
O}(f_{NL})$ impact.

  2. 2.

The one-loop operator mixing leads to an infinite tower of operators, made of
quadratic powers of the density contrast $\delta$ and higher-derivatives
thereof. Furthermore, the mixing occurs as well with non-local operators if
considered from the point of view of the smoothed density contrast (or local
if expressed in terms of the original curvature perturbation). Maybe one can
intuitively understand why an infinite tower of operators is needed by the
following argument. Around the peak of the curvature perturbation, one can
perform a rotation of the coordinate axes to be aligned with the principal
axes of length $\lambda_{i}$ ($i=1,2,3$) of the constant-curvature
perturbation ellipsoid and Taylor expanding up to second-order gives [62]

$\zeta(r)\simeq\zeta_{\rm
pk}-\frac{1}{2}\sum_{i=1}^{3}\lambda_{i}\left(x_{i}-x_{i}^{\rm
pk}\right)^{2},$ (7.1)

and

$\delta(r_{m},\vec{x}_{\rm
pk})=-\frac{4}{3}r_{m}\zeta^{\prime}(r_{m})\simeq\frac{8}{3}\left[\zeta_{\rm
pk}-\zeta(r_{m})\right].$ (7.2)

This shows that already at the linear level, the statistics of the smoothed
density contrast calculated in a volume of radius $r_{m}$ demands knowing all
correlations of the curvature perturbation in two different spatial points
[34] and therefore all its gradients.

  3. 3.

In an idealized scenario, one might envisage to calculate the PBH abundance by
involving only superhorizon physics to avoid the complications arising at
horizon crossing, and thus relying on the fact that the comoving number of
peaks which eventually will collapse into PBHs at horizon crossing remains
constant on superhorizon scales. The renormalization procedure however must be
applied as well on superhorizon scales as composite operators probe short
scales and the latter may combine in the loops to generate long modes.
Therefore, even on superhorizon scales calculating the PBH formation
probability might be more involved than naively thought.

  4. 4.

The probability of forming a PBH should not depend on any smoothing procedure,
and therefore should be independent of any window function we decide to force
into the calculation of the PBH abundance. The renormalization procedure is
exactly performed to leave no sign of such arbitrary choice. Being the non-
linear density contrast $\delta_{\rm com}(\vec{x},t)$ on comoving orthogonal
slicings (2.2) a composite operator in terms of the time independent curvature
perturbation $\zeta(\vec{x})$, one concludes that its renormalization should
be performed in order to lead to a result independent of the choice of the
smoothing procedure.

We intend to elaborate about these points in future work.

Acknowledgments.

We thank V. De Luca for useful discussions. A.I. and A.R. acknowledge support
from the Swiss National Science Foundation (project number CRSII5_213497). D.
P. and A.R. are supported by the Boninchi Foundation for the project “PBHs in
the Era of GW Astronomy”.

##  Appendix A Integrals involving the Gegenbauer polynomials

###  A.1 Integral of one window function

The window function $W_{3}(|\vec{k}-\vec{q}|r_{m})$, defined as

$W_{3}(x)=2^{3/2}\Gamma(5/2)\frac{J_{3/2}(x)}{x^{3/2}},$ (A.1)

can be treated in the following way. First we use the property

$\frac{J_{\nu}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{\nu}}=2^{\nu}\Gamma(\nu)\sum_{k=0}^{\infty}(\nu+k)\frac{J_{\nu+k}(l_{1})}{l_{1}^{\nu}}\frac{J_{\nu+k}(l_{2})}{l_{2}^{\nu}}C_{k}^{\nu}(\cos\theta),$
(A.2)

where $C_{k}^{\nu}(\cos\theta)$ are the Gegenbauer polynomials and $\theta$ is
the angle between the vectors $\vec{l}_{1}$ and $\vec{l}_{2}$. In general

$\int_{0}^{\pi}{\rm
d}\theta\,C^{\nu}_{k}(\cos\theta)\sin^{2\nu}\theta=0,\quad\forall\;k\neq 0.$
(A.3)

Therefore,

$\displaystyle\int\frac{{\rm
d}\Omega}{4\pi}\frac{J_{1/2}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{1/2}}$
$\displaystyle=$ $\displaystyle
2^{1/2}\Gamma(1/2)\frac{1}{2}\frac{J_{1/2}(l_{1})}{l_{1}^{1/2}}\frac{J_{1/2}(l_{2})}{l_{2}^{1/2}}\int\frac{{\rm
d}\Omega}{4\pi}C_{0}^{1/2}(\cos\theta)\;+$ (A.4) $\displaystyle+$
$\displaystyle
2^{1/2}\Gamma(1/2)\sum_{k=1}^{\infty}(1/2+k)\frac{J_{1/2+k}(l_{1})}{l_{1}^{1/2}}\frac{J_{1/2+k}(l_{2})}{l_{2}^{1/2}}\int\frac{{\rm
d}\Omega}{4\pi}C_{k}^{1/2}(\cos\theta).$

Only the first term in the sum is not zero, leading to the simple relation

$\int\frac{{\rm
d}\Omega}{4\pi}\frac{J_{1/2}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{1/2}}=2^{1/2}\Gamma(1/2)\frac{1}{2}\frac{J_{1/2}(l_{1})}{l_{1}^{1/2}}\frac{J_{1/2}(l_{2})}{l_{2}^{1/2}}.$
(A.5)

Now we can use the general relations

$l\frac{{\rm d}}{{\rm
d}l}\frac{J_{D/2-1}(l)}{l^{D/2-1}}=-\frac{J_{D/2}(l)}{l^{D/2-2}},$ (A.6)
$\frac{J_{D/2-1}(l)}{l^{D/2-1}}=l\frac{{\rm d}}{{\rm
d}l}\frac{J_{D/2}(l)}{l^{D/2}}+D\frac{J_{D/2}(l)}{l^{D/2}},$ (A.7)

to transform the equation above from $J_{1/2}$ to $J_{3/2}$, using the
relation

$l_{1}\partial_{l_{1}}\frac{J_{1/2}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{1/2}}=-(l_{1}^{2}+\vec{l}_{1}\cdot\vec{l}_{2})\frac{J_{3/2}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{3/2}}.$
(A.8)

We take the derivative with respect to $l_{1}$ of Eq. (A.5), leading to

$\int\frac{{\rm
d}\Omega}{4\pi}W_{3}(|\vec{l}_{1}-\vec{l}_{2}|)\left(1-\frac{\vec{l}_{1}\cdot\vec{l}_{2}}{l_{1}^{2}}\right)=W_{3}(l_{1})W_{3}(l_{2})+W_{3}(l_{1})\;\frac{l_{2}}{3}\frac{{\rm
d}}{{\rm d}l_{2}}W_{3}(l_{2}).$ (A.9)

Multiplying both sides with $l_{1}^{2}$ and summing the same relation
exchanging $l_{1}$ with $l_{2}$, we get

$\displaystyle\int\frac{{\rm
d}\Omega}{4\pi}W_{3}(|\vec{l}_{1}-\vec{l}_{2}|)\left(l_{1}^{2}+l_{2}^{2}-2\vec{l}_{1}\cdot\vec{l}_{2}\right)=\int\frac{{\rm
d}\Omega}{4\pi}W_{3}(|\vec{l}_{1}-\vec{l}_{2}|)|\vec{l}_{1}-\vec{l}_{2}|^{2}=$
$\displaystyle=W_{3}(l_{1})W_{3}(l_{2})(l_{1}^{2}+l_{2}^{2})+W_{3}(l_{1})\;\frac{l_{1}^{2}l_{2}}{3}\frac{{\rm
d}}{{\rm d}l_{2}}W_{3}(l_{2})+W_{3}(l_{2})\;\frac{l_{2}^{2}l_{1}}{3}\frac{{\rm
d}}{{\rm d}l_{1}}W_{3}(l_{1})=$
$\displaystyle=W_{3}(l_{1})W_{3}(l_{2})\left(l_{1}^{2}+l_{2}^{2}+\frac{l_{1}^{2}l_{2}}{3W_{3}(l_{2})}\frac{{\rm
d}}{{\rm d}l_{2}}W_{3}(l_{2})+\;\frac{l_{2}^{2}l_{1}}{3W_{3}(l_{1})}\frac{{\rm
d}}{{\rm d}l_{1}}W_{3}(l_{1})\right).$ (A.10)

###  A.2 Integral with two window functions

For the composite operator with two external legs, we have to evaluate the
following momentum integral

$\displaystyle I$
$\displaystyle=\int_{\vec{l}_{1}}|\vec{l}_{1}+\vec{l}_{2}|^{2}W_{3}(|\vec{l}_{1}+\vec{l}_{2}|)|\vec{l}_{1}+\vec{l}_{3}|^{2}W_{3}(|\vec{l}_{1}+\vec{l}_{3}|)$
$\displaystyle=\int_{\vec{l}_{1}}(l_{1}^{2}+l_{2}^{2}+2l_{1}l_{2}\cos\theta_{12})W_{3}(|\vec{l}_{1}+\vec{l}_{2}|)(l_{1}^{2}+l_{3}^{2}+2l_{1}l_{3}\cos\theta_{13})W_{3}(|\vec{l}_{1}+\vec{l}_{3}|)$
$\displaystyle=\int_{\vec{l}_{1}}\Big{(}A_{00}+A_{10}\cos\theta_{12}+A_{01}\cos\theta_{13}+A_{11}\cos\theta_{12}\cos\theta_{13}\Big{)}W_{3}(|\vec{l}_{1}+\vec{l}_{2}|)W_{3}(|\vec{l}_{1}+\vec{l}_{3}|)$
$\displaystyle=2^{3}\Gamma^{2}(5/2)\int_{\vec{l}_{1}}\Big{(}A_{00}+A_{10}\cos\theta_{12}+A_{01}\cos\theta_{13}+A_{11}\cos\theta_{12}\cos\theta_{13}\Big{)}\frac{J_{3/2}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{3/2}}\frac{J_{3/2}(|\vec{l}_{1}+\vec{l}_{3}|)}{|\vec{l}_{1}+\vec{l}_{3}|^{3/2}},$
(A.11)

where

$\displaystyle
A_{00}=l_{1}^{4}+l_{1}^{2}(l_{2}^{2}+l_{3}^{2})+l_{2}^{2}l_{3}^{2},~{}~{}~{}~{}~{}A_{10}=2l_{1}l_{2}(l_{1}^{2}+l_{3}^{2}),$
$\displaystyle
A_{01}=2l_{1}l_{3}(l_{1}^{2}+l_{2}^{2}),~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}A_{11}=4l_{1}^{2}l_{2}l_{3}.$
(A.12)

We have parametrized the vectors $\vec{l_{i}}$ as

$\displaystyle\vec{l}_{2}=l_{2}\;(0,\;0,\;1),$
$\displaystyle\vec{l}_{3}=l_{3}\;(0,\;\sqrt{1-y^{2}},\;y),$
$\displaystyle\vec{l}_{1}=l_{1}\;(\cos\beta\sqrt{1-x^{2}},\;\sin\beta\sqrt{1-x^{2}},\;x),$
(A.13)

so that

$\cos\theta_{12}=x,~{}~{}~{}~{}~{}\cos\theta_{13}=xy+\sqrt{1-y^{2}}\sqrt{1-x^{2}}\sin\beta.$
(A.14)

Using the addition theorem

$\displaystyle\frac{J_{\nu}(|\vec{l}_{1}+\vec{l}_{2}|)}{|\vec{l}_{1}+\vec{l}_{2}|^{\nu}}\frac{J_{\nu}(|\vec{l}_{1}+\vec{l}_{3}|)}{|\vec{l}_{1}+\vec{l}_{3}|^{\nu}}$
$\displaystyle=(2^{\nu}\Gamma(\nu))^{2}\sum_{k=0}^{\infty}\sum_{s=0}^{\infty}\frac{J_{\nu+k}(l_{1})}{l_{1}^{\nu}}\frac{J_{\nu+k}(l_{2})}{l_{2}^{\nu}}\frac{J_{\nu+s}(l_{1})}{l_{1}^{\nu}}\frac{J_{\nu+s}(l_{3})}{l_{3}^{\nu}}$
$\displaystyle\cdot(\nu+k)(\nu+s)C_{k}^{\nu}(\cos\theta_{12})C_{s}^{\nu}(\cos\theta_{13}),$
(A.15)

we can express the integral (A.11) as

$\displaystyle I$
$\displaystyle=\int_{\vec{l}_{1}}\Big{(}A_{00}+A_{10}\cos\theta_{12}+A_{01}\cos\theta_{13}+A_{11}\cos\theta_{12}\cos\theta_{13}\Big{)}\sum_{k=0}^{\infty}\sum_{s=0}^{\infty}9\pi^{2}\left(\frac{3}{2}+k\right)$
$\displaystyle\cdot\left(\frac{3}{2}+s\right)\frac{J_{3/2+k}(l_{1})}{l_{1}^{3/2}}\frac{J_{3/2+k}(l_{2})}{l_{2}^{3/2}}\frac{J_{3/2+s}(l_{1})}{l_{1}^{3/2}}\frac{J_{3/2+s}(l_{3})}{l_{3}^{3/2}}C_{k}^{3/2}(\cos\theta_{12})C_{s}^{3/2}(\cos\theta_{13}).$
(A.16)

Therefore, the integral (A.11) can be written as

$\displaystyle I=\sum_{i,j=0}^{1}\int\frac{l^{2}{\rm
d}l}{2\pi^{2}}A_{ij}I^{ij},$ (A.17)

where

$\displaystyle
I^{ij}=\sum_{k=0}^{\infty}\sum_{s=0}^{\infty}N_{ks}I^{ij}_{ks}(y),$ (A.18)

$\displaystyle N_{ks}$
$\displaystyle=9\pi^{2}\left(\frac{3}{2}+k\right)\left(\frac{3}{2}+s\right)\frac{J_{3/2+k}(l_{1})}{l_{1}^{3/2}}\frac{J_{3/2+k}(l_{2})}{l_{2}^{3/2}}\frac{J_{3/2+s}(l_{1})}{l_{1}^{3/2}}\frac{J_{3/2+s}(l_{3})}{l_{3}^{3/2}}$
$\displaystyle=\frac{4}{9}\left(\frac{3}{2}+k\right)\left(\frac{3}{2}+s\right)\left[l_{1}^{k}\left(\frac{1}{l_{1}}\frac{{\rm
d}}{{\rm
d}l_{1}}\right)^{k}W_{3}(l_{1})\right]\left[l_{1}^{s}\left(\frac{1}{l_{1}}\frac{{\rm
d}}{{\rm d}l_{1}}\right)^{s}W_{3}(l_{1})\right]$
$\displaystyle~{}~{}\times\left[l_{2}^{k}\left(\frac{1}{l_{2}}\frac{{\rm
d}}{{\rm
d}l_{2}}\right)^{k}W_{3}(l_{2})\right]\left[l_{3}^{s}\left(\frac{1}{l_{3}}\frac{{\rm
d}}{{\rm d}l_{3}}\right)^{s}W_{3}(l_{3})\right],$ (A.19)

and

$\displaystyle I^{ij}_{ks}(y=\cos\theta_{23})=\int\frac{{\rm
d}\Omega}{4\pi}(\cos\theta_{12})^{i}(\cos\theta_{13})^{j}C_{k}^{3/2}(\cos\theta_{12})C_{s}^{3/2}(\cos\theta_{13}).$
(A.20)

Therefore, in order to calculate Eq. (A.11), we need to calculate the four
integrals

$\displaystyle
I^{00}_{ks}(y),~{}I^{10}_{ks}(y),~{}I^{01}_{ks}(y),~{}~{}~{}\mbox{and}~{}~{}~{}I^{11}_{ks}(y),$
(A.21)

which we will do in the following. Inserting the angles (A.14) and using the
fact that

$\displaystyle\int_{0}^{2\pi}\frac{{\rm
d}\beta}{2\pi}\,C_{k}^{3/2}\left(xy+\sqrt{1-y^{2}}\sqrt{1-x^{2}}\sin\beta\right)$
$\displaystyle=\left\\{\begin{array}[]{lll}{\displaystyle{\sum_{n={\rm
even}}^{k}}}(2n+1)P_{n}(x)P_{n}(y)&~{}{\rm if}&~{}k=\mbox{even},\\\ &&\\\
{\displaystyle{\sum_{n={\rm odd}}^{k}}}(2n+1)P_{n}(x)P_{n}(y)&~{}{\rm
if}&~{}k=\mbox{odd},\\\ \end{array}\right.$ (A.25)

and

$\displaystyle C_{k}^{3/2}(x)$ $\displaystyle=\frac{{\rm d}}{{\rm
d}x}P_{k+1}(x)$
$\displaystyle=(2k+1)P_{k}(x)+\Big{(}2(k-2)+1\Big{)}P_{k-2}(x)+\Big{(}2(k-4)+1\Big{)}P_{k-4}(x)+\cdots,$
(A.26)

the integration over the angles gives

$\displaystyle I^{00}_{ks}$ $\displaystyle=\int\frac{{\rm
d}\Omega}{4\pi}C_{k}^{3/2}(x)C_{s}^{3/2}\left(xy+\sqrt{1-y^{2}}\sqrt{1-x^{2}}\sin\beta\right)=$
$\displaystyle\hskip 113.81102pt=\left\\{\begin{array}[]{lll}\dfrac{{\rm
d}}{{\rm d}y}P_{{\rm min}(k,s)+1}(y)&{\rm if}&(k,s)\in 2{\mathbb{N}},\\\ &&\\\
\dfrac{{\rm d}}{{\rm d}y}P_{{\rm min}(k,s)+1}(y)&{\rm if}&(k,s)\in
2{\mathbb{N}}+1,\\\ &&\\\ 0&{\rm otherwise},&\\\ \end{array}\right.,$ (A.32)

or

$\displaystyle I^{00}_{ks}=\left(\frac{1+(-1)^{k+s}}{2}\right)\dfrac{{\rm
d}}{{\rm d}y}P_{{\rm min}(k,s)+1}(y).$ (A.33)

The next integral to evaluate is

$I_{ks}^{10}(y)=\int\frac{{\rm
d}\Omega}{4\pi}\,x\,C_{k}^{3/2}(x)C_{s}^{3/2}(\cos\theta_{13}).$ (A.34)

To proceed, we need the recursive relation of Legendre polynomials

$\displaystyle x\frac{{\rm d}P_{k+1}(x)}{{\rm d}x}=\frac{{\rm
d}P_{k+2}(x)}{{\rm d}x}-(k+2)P_{k+1}(x),$ (A.35)

which can be written as

$\displaystyle x\,C_{k}^{3/2}(x)=C^{3/2}_{k+1}-(k+2)P_{k+1}(x).$ (A.36)

Then, using Eq. (A.25), we can easily see that

$\displaystyle\int_{0}^{2\pi}\frac{{\rm
d}\Omega}{4\pi}\,P_{k+1}(x)C_{s}^{3/2}(\cos\theta_{13})=\left\\{\begin{array}[]{lll}\Theta(s-k-1)\,P_{k+1}(y)&{\rm
if}&(k,s)\in 2{\mathbb{N}},\\\ &&\\\ \Theta(s-k-1)\,P_{k+1}(y)&{\rm
if}&(k,s)\in 2{\mathbb{N}}+1,\\\ &&\\\ 0&{\rm otherwise},&\\\
\end{array}\right.$ (A.42)

where $\Theta(n)=0$ for $n<0$ and $\Theta(n)=1$ for $n\geq 0$. Therefore,
using Eqs. (A.36) and (A.42), we find that

$\displaystyle
I^{10}_{ks}=\left(\frac{1-(-1)^{k+s}}{2}\right)\left[\Theta(k-s)\frac{\rm
d}{{\rm d}y}P_{k+1}(y)+\Theta(s-k)y\frac{\rm d}{{\rm d}y}P_{s+1}(y)\right]$
(A.43)

To calculate $I^{01}$ we need the relation

$\displaystyle\int_{0}^{2\pi}\frac{{\rm d}\beta}{2\pi}\,$
$\displaystyle\cos\theta_{13}C_{k}^{3/2}\left(\cos\theta_{13}\right)=$
$\displaystyle=\left\\{\begin{array}[]{lll}2(k+1)P_{k+1}(x)P_{k+1}(y)+{\displaystyle{\sum_{n={\rm
odd}}^{k}}}2(2n+1)P_{n}(x)P_{n}(y)&~{}{\rm if}&~{}k=\mbox{even},\\\ &&\\\
2(k+1)P_{k+1}(x)P_{k+1}(y)+{\displaystyle{\sum_{n={\rm
even}}^{k}}}2(2n+1)P_{n}(x)P_{n}(y)&~{}{\rm if}&~{}k=\mbox{odd},\\\
\end{array}\right.$ (A.47)

Then, by using Eq. (A.26), we find that the integral

$\displaystyle I^{01}_{ks}$ $\displaystyle=$ $\displaystyle\int\frac{{\rm
d}\Omega}{4\pi}\,C^{3/2}_{k}(x)\cos\theta_{13}C_{s}^{3/2}(\cos\theta_{13})$
(A.48)

turns out to be

$I^{01}_{ks}=\left(\frac{1-(-1)^{k+s}}{2}\right)\left[\Theta(s-k)\frac{\rm
d}{{\rm d}y}P_{k+1}(y)+\Theta(k-s)y\frac{\rm d}{{\rm d}y}P_{s+1}(y)\right].$
(A.49)

The final integral $I^{11}_{ks}$ can be calculated using Eq. (A.36). The
result is

$\displaystyle
I^{11}_{ks}=I^{01}_{k+1,s}-(k+2)\left\\{\left(\frac{1+(-1)^{k+s}}{2}\right)\Theta(s-k-1)P_{k+1}(y)+\frac{s+1}{2s+3}\delta_{sk}P_{s+1}(y)\right\\}.$
(A.50)

##  Appendix B Non-local operators

To derive the expression (5.12), we start Eq. (4.3) which can be rewritten as

$\frac{\delta(r_{m},\vec{x}_{\rm
pk})}{r_{m}^{2}}=\frac{4}{9}\nabla^{2}_{x_{\rm
pk}}\int_{\vec{k}}\;e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}W_{3}(kr_{m})\zeta_{\vec{k}}.$ (B.1)

We now take the derivative $r_{m}\partial_{r_{m}}$ of both sides to get

$\frac{1}{r_{m}}\partial_{r}\delta(r,\vec{x}_{\rm
pk})\Big{|}_{r=r_{m}}=-2\frac{\delta(r_{m},\vec{x}_{\rm
pk})}{r_{m}^{2}}+\frac{1}{r_{m}}\partial_{r}\delta(r,\vec{x}_{\rm
pk})\Big{|}_{r=r_{m}}=\frac{4}{9}\nabla^{2}_{x_{\rm
pk}}\int_{\vec{k}}\;e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}\zeta_{\vec{k}}\,r_{m}\partial_{r}W_{3}(kr)\Big{|}_{r=r_{m}}.$ (B.2)

By using Eq. (2.4) we get

$\partial_{r}\delta(r,\vec{x}_{\rm pk})\Big{|}_{r=r_{m}}=0,$ (B.3)

and consequently

$-2\frac{\delta(r_{m},\vec{x}_{\rm
pk})}{r_{m}^{2}}=\frac{4}{9}\nabla^{2}_{x_{\rm
pk}}\int_{\vec{k}}\;e^{i\vec{k}\cdot\vec{x}_{\rm
pk}}\zeta_{\vec{k}}\,r_{m}\partial_{r}W_{3}(kr)\Big{|}_{r=r_{m}}$ (B.4)

or

$r_{m}\partial_{r}W_{3}(kr)\Big{|}_{r=r_{m}}\zeta_{\vec{k}}=-\frac{9}{2r_{m}^{2}}\int_{\vec{x}_{\rm
pk}}\;e^{-i\vec{k}\cdot\vec{x}_{\rm pk}}\;\nabla^{-2}_{x_{\rm
pk}}\delta(r_{m},\vec{x}_{\rm pk})$ (B.5)

which is a non-local operator.

## References

  * [1] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, _Primordial black holes—perspectives in gravitational wave astronomy_ , _Class. Quant. Grav._ 35 (2018) 063001, [1801.05235].
  * [2] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on primordial black holes_ , _Rept. Prog. Phys._ 84 (2021) 116902, [2002.12778].
  * [3] A. M. Green and B. J. Kavanagh, _Primordial Black Holes as a dark matter candidate_ , _J. Phys. G_ 48 (2021) 043001, [2007.10722].
  * [4] LISA Cosmology Working Group collaboration, E. Bagui et al., _Primordial black holes and their gravitational-wave signatures_ , 2310.19857.
  * [5] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102, [1602.03837].
  * [6] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., _GWTC-1: A Gravitational-Wave Transient Catalog of Compact Binary Mergers Observed by LIGO and Virgo during the First and Second Observing Runs_ , _Phys. Rev. X_ 9 (2019) 031040, [1811.12907].
  * [7] LIGO Scientific, Virgo collaboration, R. Abbott et al., _GWTC-2: Compact Binary Coalescences Observed by LIGO and Virgo During the First Half of the Third Observing Run_ , _Phys. Rev. X_ 11 (2021) 021053, [2010.14527].
  * [8] LIGO Scientific, VIRGO, KAGRA collaboration, R. Abbott et al., _GWTC-3: Compact Binary Coalescences Observed by LIGO and Virgo During the Second Part of the Third Observing Run_ , 2111.03606.
  * [9] S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz et al., _Did LIGO detect dark matter?_ , _Phys. Rev. Lett._ 116 (2016) 201301, [1603.00464].
  * [10] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, _Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914_ , _Phys. Rev. Lett._ 117 (2016) 061101, [1603.08338].
  * [11] S. Clesse and J. García-Bellido, _The clustering of massive Primordial Black Holes as Dark Matter: measuring their mass distribution with Advanced LIGO_ , _Phys. Dark Univ._ 15 (2017) 142–147, [1603.05234].
  * [12] Y. Ali-Haïmoud, E. D. Kovetz and M. Kamionkowski, _Merger rate of primordial black-hole binaries_ , _Phys. Rev. D_ 96 (2017) 123523, [1709.06576].
  * [13] G. Hütsi, M. Raidal, V. Vaskonen and H. Veermäe, _Two populations of LIGO-Virgo black holes_ , _JCAP_ 03 (2021) 068, [2012.02786].
  * [14] V. De Luca, G. Franciolini, P. Pani and A. Riotto, _Bayesian Evidence for Both Astrophysical and Primordial Black Holes: Mapping the GWTC-2 Catalog to Third-Generation Detectors_ , _JCAP_ 05 (2021) 003, [2102.03809].
  * [15] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti et al., _Searching for a subpopulation of primordial black holes in LIGO-Virgo gravitational-wave data_ , _Phys. Rev. D_ 105 (2022) 083526, [2105.03349].
  * [16] G. Franciolini, I. Musco, P. Pani and A. Urbano, _From inflation to black hole mergers and back again: Gravitational-wave data-driven constraints on inflationary scenarios with a first-principle model of primordial black holes across the QCD epoch_ , _Phys. Rev. D_ 106 (2022) 123526, [2209.05959].
  * [17] Z.-C. Chen and Q.-G. Huang, _Distinguishing Primordial Black Holes from Astrophysical Black Holes by Einstein Telescope and Cosmic Explorer_ , _JCAP_ 08 (2020) 039, [1904.02396].
  * [18] O. Pujolas, V. Vaskonen and H. Veermäe, _Prospects for probing gravitational waves from primordial black hole binaries_ , _Phys. Rev. D_ 104 (2021) 083521, [2107.03379].
  * [19] V. De Luca, G. Franciolini, P. Pani and A. Riotto, _The minimum testable abundance of primordial black holes at future gravitational-wave detectors_ , _JCAP_ 11 (2021) 039, [2106.13769].
  * [20] S. Barsanti, V. De Luca, A. Maselli and P. Pani, _Detecting Subsolar-Mass Primordial Black Holes in Extreme Mass-Ratio Inspirals with LISA and Einstein Telescope_ , _Phys. Rev. Lett._ 128 (2022) 111104, [2109.02170].
  * [21] S. S. Bavera, G. Franciolini, G. Cusin, A. Riotto, M. Zevin and T. Fragos, _Stochastic gravitational-wave background as a tool for investigating multi-channel astrophysical and primordial black-hole mergers_ , _Astron. Astrophys._ 660 (2022) A26, [2109.05836].
  * [22] I. Musco, _Threshold for primordial black holes: Dependence on the shape of the cosmological perturbations_ , _Phys. Rev. D_ 100 (2019) 123524, [1809.02127].
  * [23] I. Musco, V. De Luca, G. Franciolini and A. Riotto, _Threshold for primordial black holes. II. A simple analytic prescription_ , _Phys. Rev. D_ 103 (2021) 063538, [2011.03014].
  * [24] V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Ünal, _The Ineludible non-Gaussianity of the Primordial Black Hole Abundance_ , _JCAP_ 07 (2019) 048, [1904.00970].
  * [25] S. Young, I. Musco and C. T. Byrnes, _Primordial black hole formation and abundance: contribution from the non-linear relation between the density and curvature perturbation_ , _JCAP_ 11 (2019) 012, [1904.00984].
  * [26] S. Young, _The primordial black hole formation criterion re-examined: Parametrisation, timing and the choice of window function_ , _Int. J. Mod. Phys. D_ 29 (2019) 2030002, [1905.01230].
  * [27] V. De Luca, A. Kehagias and A. Riotto, _How well do we know the primordial black hole abundance: The crucial role of nonlinearities when approaching the horizon_ , _Phys. Rev. D_ 108 (2023) 063531, [2307.13633].
  * [28] G. Ferrante, G. Franciolini, A. Iovino, Junior. and A. Urbano, _Primordial non-Gaussianity up to all orders: Theoretical aspects and implications for primordial black hole models_ , _Phys. Rev. D_ 107 (2023) 043520, [2211.01728].
  * [29] M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, _The formation probability of primordial black holes_ , _Phys. Lett. B_ 820 (2021) 136602, [2105.07810].
  * [30] M. Shibata and M. Sasaki, _Black hole formation in the Friedmann universe: Formulation and computation in numerical relativity_ , _Phys. Rev. D_ 60 (1999) 084002, [gr-qc/9905064].
  * [31] T. Harada, C.-M. Yoo, T. Nakama and Y. Koga, _Cosmological long-wavelength solutions and primordial black hole formation_ , _Phys. Rev. D_ 91 (2015) 084057, [1503.03934].
  * [32] A. Escrivà, C. Germani and R. K. Sheth, _Universal threshold for primordial black hole formation_ , _Phys. Rev. D_ 101 (2020) 044022, [1907.13311].
  * [33] C. Germani and R. K. Sheth, _Nonlinear statistics of primordial black holes from Gaussian curvature perturbations_ , _Phys. Rev. D_ 101 (2020) 063520, [1912.07072].
  * [34] V. De Luca and A. Riotto, _A note on the abundance of primordial black holes: Use and misuse of the metric curvature perturbation_ , _Phys. Lett. B_ 828 (2022) 137035, [2201.09008].
  * [35] A. D. Gow, H. Assadullahi, J. H. P. Jackson, K. Koyama, V. Vennin and D. Wands, _Non-perturbative non-Gaussianity and primordial black holes_ , _EPL_ 142 (2023) 49001, [2211.08348].
  * [36] N. Bartolo, E. Komatsu, S. Matarrese and A. Riotto, _Non-Gaussianity from inflation: Theory and observations_ , _Phys. Rept._ 402 (2004) 103–266, [astro-ph/0406398].
  * [37] Planck collaboration, Y. Akrami et al., _Planck 2018 results. IX. Constraints on primordial non-Gaussianity_ , _Astron. Astrophys._ 641 (2020) A9, [1905.05697].
  * [38] Y.-F. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D.-G. Wang and Z. Wang, _Revisiting non-Gaussianity from non-attractor inflation models_ , _JCAP_ 05 (2018) 012, [1712.09998].
  * [39] N. Bartolo, S. Matarrese and A. Riotto, _On nonGaussianity in the curvaton scenario_ , _Phys. Rev. D_ 69 (2004) 043503, [hep-ph/0309033].
  * [40] M. Sasaki, J. Valiviita and D. Wands, _Non-Gaussianity of the primordial perturbation in the curvaton model_ , _Phys. Rev. D_ 74 (2006) 103003, [astro-ph/0607627].
  * [41] NANOGrav collaboration, G. Agazie et al., _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8, [2306.16213].
  * [42] NANOGrav collaboration, G. Agazie et al., _The NANOGrav 15 yr Data Set: Observations and Timing of 68 Millisecond Pulsars_ , _Astrophys. J. Lett._ 951 (2023) L9, [2306.16217].
  * [43] EPTA collaboration, J. Antoniadis et al., _The second data release from the european pulsar timing array iii. search for gravitational wave signals_ , 2306.16214.
  * [44] EPTA collaboration, J. Antoniadis et al., _The second data release from the european pulsar timing array i. the dataset and timing analysis_ , 2306.16224.
  * [45] EPTA collaboration, J. Antoniadis et al., _The second data release from the european pulsar timing array: V. implications for massive black holes, dark matter and the early universe_ , 2306.16227.
  * [46] D. J. Reardon et al., _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L6, [2306.16215].
  * [47] PPTA collaboration, A. Zic et al., _The parkes pulsar timing array third data release_ , 2306.16230.
  * [48] D. J. Reardon et al., _The Gravitational-wave Background Null Hypothesis: Characterizing Noise in Millisecond Pulsar Arrival Times with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L7, [2306.16229].
  * [49] H. Xu et al., _Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I_ , _Res. Astron. Astrophys._ 23 (2023) 075024, [2306.16216].
  * [50] G. Franciolini, A. Iovino, Junior., V. Vaskonen and H. Veermae, _The recent gravitational wave observation by pulsar timing arrays and primordial black holes: the importance of non-gaussianities_ , 2306.17149.
  * [51] L. Liu, Z.-C. Chen and Q.-G. Huang, _Implications for the non-Gaussianity of curvature perturbation from pulsar timing arrays_ , 2307.01102.
  * [52] S. Wang, Z.-C. Zhao, J.-P. Li and Q.-H. Zhu, _Exploring the Implications of 2023 Pulsar Timing Array Datasets for Scalar-Induced Gravitational Waves and Primordial Black Holes_ , 2307.00572.
  * [53] Y.-F. Cai, X.-C. He, X. Ma, S.-F. Yan and G.-W. Yuan, _Limits on scalar-induced gravitational waves from the stochastic background by pulsar timing array observations_ , 2306.17822.
  * [54] K. Inomata, K. Kohri and T. Terada, _The Detected Stochastic Gravitational Waves and Sub-Solar Primordial Black Holes_ , 2306.17834.
  * [55] D. G. Figueroa, M. Pieroni, A. Ricciardone and P. Simakachorn, _Cosmological Background Interpretation of Pulsar Timing Array Data_ , 2307.02399.
  * [56] Z. Yi, Q. Gao, Y. Gong, Y. Wang and F. Zhang, _The waveform of the scalar induced gravitational waves in light of Pulsar Timing Array data_ , 2307.02467.
  * [57] Q.-H. Zhu, Z.-C. Zhao and S. Wang, _Joint implications of BBN, CMB, and PTA Datasets for Scalar-Induced Gravitational Waves of Second and Third orders_ , 2307.03095.
  * [58] J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki et al., _What is the source of the PTA GW signal?_ , 2308.08546.
  * [59] C. Germani and I. Musco, _Abundance of Primordial Black Holes Depends on the Shape of the Inflationary Power Spectrum_ , _Phys. Rev. Lett._ 122 (2019) 141302, [1805.04087].
  * [60] V. De Luca, G. Franciolini and A. Riotto, _On the Primordial Black Hole Mass Function for Broad Spectra_ , _Phys. Lett. B_ 807 (2020) 135550, [2001.04371].
  * [61] H. Firouzjahi and A. Riotto, _The Sign of non-Gaussianity and the Primordial Black Holes Abundance_ , 2309.10536.
  * [62] J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, _The Statistics of Peaks of Gaussian Random Fields_ , _Astrophys. J._ 304 (1986) 15–61.


# Caption Enriched Samples for Improving Hateful Memes Detection

Efrat Blaier
Tel-Aviv University
&Itzik Malkiel∗
Tel-Aviv University
&Lior Wolf
Tel-Aviv University  Denotes equal contribution.

###### Abstract

The recently introduced hateful meme challenge demonstrates the difficulty of
determining whether a meme is hateful or not. Specifically, both unimodal
language models and multimodal vision-language models cannot reach the human
level of performance. Motivated by the need to model the contrast between the
image content and the overlayed text, we suggest applying an off-the-shelf
image captioning tool in order to capture the first. We demonstrate that the
incorporation of such automatic captions during fine-tuning improves the
results for various unimodal and multimodal models. Moreover, in the unimodal
case, continuing the pre-training of language models on augmented and original
caption pairs, is highly beneficial to the classification accuracy. Our code
is publicly available 111https://github.com/efrat-safanov/caption-enriched-
samples-research.

##  1 Introduction

Multimodal transformers, including VisualBERT Li et al. (2019), VilBERT Lu et
al. (2019), and UNITER Chen et al. (2020) are currently the state of the art
methods in tasks that involve both text and images, such as visual question
answering Antol et al. (2015) and image captioning Chen et al. (2015). One
particular example is that of the “hateful memes” challenge Kiela et al.
(2020), in which the task is to classify whether a given meme is hateful or
not. By construction (and similar to real-world memes), the classifier has to
apply complex reasoning that involves a multimodal text-image analysis from
both image and text. The prevalence of sarcasm adds another layer of
difficulty to both machines and humans. Kiela et al. (2020) results indicate
that the human accuracy is only about 85%. Recent state-of-the-art models
perform considerably less accurately, achieving up to 64.73% accuracy.

In this work, we adopt an off-the-shelf caption generator, and show that
combining such captions as a third input (in addition to the image and the
overlayed text) can effectively improve the performance of both multimodal and
unimodal models. Two techniques are used: continued pre-training of the
baseline model with the additional augmented captions, and using the added
information as part of the fine-tuned model used for the hateful memes task.

Figure 1: An example of memes in which the classifiction depends on both the
text and the image. Specifically, in this Benign Confounder example, the same
text appears and the classification changes depending on the image.

##  2 Related Work

Transformer-based Language Models (LMs) have shown significant performance
gains in various natural language understanding tasks, such as text-similarity
and language inference Devlin et al. (2019); Liu et al. (2019). These models
are usually pre-trained on the Masked Language Model (MLM) objective followed
by a task-specific fine-tuning process Wang et al. (2018).

Recently, vision and language tasks have gained a lot of traction, where
transformer-based multimodality models were introduced, showing great promise
in solving multimodality problems, such as visual reasoning Suhr et al.
(2019), image caption generation Chen et al. (2015) and more.

In Kiela et al. (2020), the authors introduced the Hateful Memes challenge - a
dataset of 10,000 memes, each associated with an image and text. Each meme
sample image-text pair is binary labeled as hateful or not hateful. The
dataset is split into 8,500 train memes, 500 and 1,000 validation and test
samples, respectively. The test labels are not public and are used for
leaderboard scoring of models.

The dataset was artificially constructed based on real-world memes. To ensure
that the dataset requires multimodal reasoning, image and text confounders
were crafted to enrich the dataset with contrastive examples. Confounders
samples are memes that appear in two variations: same image with different
texts that have different labels or vice versa (see Fig. 1).

In Kiela et al. (2020), the authors show that both unimodal models (such as
Bert Devlin et al. (2019)) and multimodal models (e.g. VisualBERT Li et al.
(2019)) seem to produce relatively poor performance, that is significantly
inferior to the human accuracy (reported to be 84.7%). This can be attributed
to the task difficulty, which requires an advanced reasoning process, where
some samples are found to also confuse humans.

The task of caption generation has been shown to yield promising results in
various works, such as the show and tell Vinyals et al. (2015), show tell and
attend Xu et al. (2016), and more. Among the different caption generator
techniques, the Neural Image Caption (NIC) Vinyals et al. (2017) has been
shown to yield promising results on the various caption tasks. The NIC model
is an encoder-decoder network, trained to generate captions from images. We
use this generator to augment the hateful memes dataset, by associating each
meme image-text pair with an additional auto-generated caption. In our work,
we adopt the implementation of the IBM Image Caption
Generator222github.com/IBM/MAX-Image-Caption-Generator.

Gururangan et al. (2020) show that continued pre-training of a pre-trained
language model on a given dataset at hand improves the performance of the
fine-tuning procedure. In a similar spirit, we demonstrate that for the
hateful meme task, it can be beneficial to continue pre-training on text-pairs
of meme text and augmented captions.

Figure 2: Two representative samples of unimodal and multimodal predictions,
with and without the CES approach.

##  3 Method

The hateful memes challenge is a binary classification problem in which each
meme sample $x=(m,c)$ consists of an image in some image domain $M$ and an
associated text in a language domain $C$. In both cases, we employ an off-the-
shelf image-caption generator $G:M\rightarrow C$, which receives an image from
domain $M$ and returns a caption from $C$. Every sample thus becomes
$x^{*}=(m,c,c^{*})$, where $c^{*}=G(m)$.

We next present the training strategy, termed Caption Enriched Samples (CES),
in either the unimodal or the multimodal settings. In the former, the network
classifier considers only text input, and in the latter image-text inputs.

| Validation | Test
---|---|---
| Baseline | CES | Baseline | CES
Human | - |  | 82.65 | -
BERT | 64.65 | 69.32 | 65.08 | 70.70
RoBERTa | 65.77 | 72.95 | 65.77 | 72.37
VisualBERT | 73.97 | 75.55 | 71.41 | 72.93
ViLBERT | 70.07 | 75.98 | 70.03 | 72.59

Table 1: AUROC for various models reported on validation and test sets of the
hateful memes dataset.

###  3.1 The unimodal case

In unimodal settings, we employ a dual-phase training approach, which utilizes
caption-pairs composed of the original caption and the generated caption
$(c,c^{*})$. In the first training phase, we consider the caption-pairs of
each sample, tokenize each element, concatenate them using the special [SEP]
token, and apply a standard masking procedure. Similar to BERT, the training
objective is to optimize a standard masked language model (MLM). Notably, in
contrast to Gururangan et al. (2020), our method applies this pre-training on
an augmented text, where half of the data was generated by an image-to-text
generator.

In the second phase, we fine-tune the same model on the task at hand, by
feeding the model with the pairs $(c,c^{*})$. During this phase, we employ a
standard fine-tuning approach, where we initialize a classification head on
top of the pre-trained model and minimize a binary cross-entropy loss.

###  3.2 The multimodal case

In multimodal settings, in order to reduce training costs, we do not continue
the pre-training phase of the multimodal backbones on the triplets
$(m,c,c^{*})$, as done in the unimodal case, and employ only a fine-tuning
strategy.

We, therefore, adopt a multimodal backbone network that was pre-trained on a
general vision-language task, and initialize a classification head composed of
two fully connected layers on top. We fine-tune the models, by propagating the
triplets $(m,c,c^{*})$ and optimizing a standard binary cross entropy loss for
the prediction of the meme labels (hateful/non hateful).

The classification head is employed over the pooled representation of the
embeddings. The pooled embeddings are the ones proposed by each of the
multimodal models at hand, as detailed in Sec.4. For example, in the ViLBERT
model Lu et al. (2019), the learned classification head is employed on top of
the element-wise product of the image and text representations, which are the
embeddings of the first tokens obtained from each of the image sequence and
the textual sequence. Note that when employing CES on the ViLBERT, the textual
sequence is a pair of sentences separated by the special [SEP].

Model | Validation | Test
---|---|---
BERT | 64.65 | 65.08
i CES-BERT | 69.32 | 70.7
ii CES-BERT | 65.19 | *
iii CES-BERT | 69.25 | *
CES-BERT | 69.76 | *
RoBERTa | 65.44 | 65.77
i CES-RoBERTa | 70.47 | 70.91
ii CES-RoBERTa | 64.8 | *
iii CES-RoBERTa | 71.11 | 71.08
CES-RoBERTa | 72.95 | 72.37
UNITER | 77.01 | 78.68
UNITER+RoBERTa | 78.04 | *
UNITER+i.CES-RoBERTa | 77.53 | *
UNITER+ii.CES-RoBERTa | 78.31 | *
UNITER+iii.CES-RoBERTa | 78.57 | *
UNITER+CES-RoBERTa | 78.29 | 78.90

Table 2: Ablation study results for the unimodal models. Reported are the
AUROC scores on the validation and test sets. (*) For some ablations, we do
not have the test results, since the Hateful Memes challenge server went
offline on May 1st, 2021.

##  4 Results

In this section, we evaluate and report the performance of CES applied to the
following models:

#### BERT and RoBERTa

are the BERT${}_{\text{BASE}}$ Devlin et al. (2019) and
RoBERTa${}_{\text{BASE}}$Liu et al. (2019) models, fine-tuned on the meme text
using a standard fine-tuning approach.

#### ViLBERT

is a multimodal Transformer with a two-stream architecture that embeds image
regions and language separately, employed with co-attentional transformer
layers that allow both streams to interact with each other Lu et al. (2019).

#### VisualBERT

is another multimodal Transformer that was pre-trained on image-text inputs,
by aligning elements of the text and regions in the image through cross-
attention and self-attention operations Li et al. (2019). In contrast to
ViLBERT, this model employs a single-stream architecture, embedding image
regions and language tokens through the same transformer blocks.

For the image-caption generator $G$ we employ the IBM MaxCaption model Vinyals
et al. (2017) for generating captions for all memes in the dataset. The length
of the generated captions is at most 15 words (see Fig.2).

In Tab.1, we report the performance of the above models applied with and
without our CES strategy. As can be seen, CES is shown to improve the
performance of all models by a sizeable margin. Specifically, for the BERT and
RoBERTa unimodal models, we observe that CES improves the test performance by
8.6% and 10%, respectively. This can be attributed to the fact that in CES,
the unimodal receives a textual representation of the image (the generated
caption $c^{*}$), which is a proxy to the image itself. In the multimodal
settings, we observe that CES improves the performance of ViLBERT and
VisualBERT by 2.1% and 3.6%, respectively. A similar table reporting the
accuracy scores can be found in the appendix.

Fig. 2 depicts two representative samples demonstrating the effectiveness of
CES applied to the VisualBERT model. As can be seen in the figure, the CES-
VisualBERT model was able to correctly predict the right label for both
images, while the baseline VisualBERT model was confused.

###  4.1 Ablation Study

We perform an ablation study for the CES technique, by applying multiple
variants of it to the BERT and RoBERTa models. We consider three variants: (i)
skipping the first phase of pre-training on the pairs $(c,c^{*})$ of meme text
and generated captions. (ii) applying the continued pre-training followed by
fine-tuning, without utilizing the generated captions in both phases (i.e. we
solely use the original meme text). (iii) applying continued pre-training
solely on the original memes text (without the augmented captions), followed
by fine-tuning that utilizes the pairs $(c,c^{*})$.

The results, shown in Tab. 2, indicate that it is important to continue the
pre-training on the text-pairs of original meme text and augmented captions,
and that the use of the augmented captions during fine-tuning improves the
accuracy of the models.

To further demonstrate the importance of each of the CES components, we
evaluate the performance of each RoBERTa variant in an ensemble that includes
the state-of-the-art UNITER model. The UNITER model is multi-modal
Transformer-based network employing a single-stream architecture Chen et al.
(2020) for both vision and textual input. The UNITER pre-training utilizes a
combination of four tasks: Masked Language Modeling, Masked Region Modeling,
Image-Text Matching (ITM), and Word-Region Alignment (WRA).

We enhance the UNITER model by either of the RoBERTa ablation variants, by
utilizing the late fusion technique form Kiela et al. (2020). Specifically, we
concatenate the CLS embeddings of the UNITER and each RoBERTa model, and train
a two-layer MLP on top. For the UNITER model, we use the Vilio Muennighoff
(2020) implementation. As can be seen in Tab. 2, all ensemble variants improve
the UNITER performance, while the ones that include a continued pre-trained
RoBERTa (UNITER+ii.CES-RoBERTa, UNITER+iii.CES-RoBERTa, UNITER+CES-RoBERTa)
yield the largest gains. Following a t-test on the validation set, we observe
that those three ensembles yield significant improvement compared to the
UNITER and UNITER+RoBERTa baselines (p-value<0.02), while the differences
between them are not statistically significant.

Figure 3: Limitations of augmented captions. (A) A representative sample of an
inaccurate augmented caption (wrong animal), for which CES was able to recover
and improve over the baseline. (B) A representative sample of a wrong caption
(detecting a dog playing instead of a shoe), for which CES failed to classify
the meme as hateful.

###  4.2 Limitations and Discussion

CES relies on augmented captions, and, therefore, it also depends on the
performance of the caption generator. In Fig. 3(A), we observe a failure case
of the caption generator, where the pig in the image was interpreted as a cow.
In this case, we observe that CES could recover from the non-accurate caption,
perhaps since the cow could be negatively interpreted in this context. Fig.
3(B) demonstrates a case where a wrong augmented caption can confuse CES. In
the figure, we observe that the caption generator describes the image as a dog
with a frisbee, instead of a shoe getting out of the ocean. As can be seen,
CES misclassified the image as not hateful.

##  5 Summary

Detecting hateful memes is a challenging task that requires complex image-text
reasoning. In this work, we introduce CES, a training approach for both
multimodal and unimodal models that leverages an off-the-shelf image-caption
generator to enhance the accuracy of hateful meme classifiers. Our results
indicate that augmented captions are highly beneficial for transformer-based
models, which seem to be able to effectively attend between the augmented
caption, meme text, and image.

## Acknowledgment

This project has received funding from the European Research Council (ERC)
under the European Union Horizon 2020 research and innovation programme (grant
ERC CoG 725974). The contribution of the second author is part of a PhD thesis
research conducted at Tel Aviv University.

## References

  * Antol et al. (2015) Stanislaw Antol, Aishwarya Agrawal, Jiasen Lu, Margaret Mitchell, Dhruv Batra, C Lawrence Zitnick, and Devi Parikh. 2015.  Vqa: Visual question answering.  In _Proceedings of the IEEE international conference on computer vision_ , pages 2425–2433.
  * Chen et al. (2015) Xinlei Chen, Hao Fang, Tsung-Yi Lin, Ramakrishna Vedantam, Saurabh Gupta, Piotr Dollar, and C. Lawrence Zitnick. 2015.  Microsoft coco captions: Data collection and evaluation server.
  * Chen et al. (2020) Yen-Chun Chen, Linjie Li, Licheng Yu, Ahmed El Kholy, Faisal Ahmed, Zhe Gan, Yu Cheng, and Jingjing Liu. 2020.  Uniter: Universal image-text representation learning.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  Bert: Pre-training of deep bidirectional transformers for language understanding.
  * Gururangan et al. (2020) Suchin Gururangan, Ana Marasović, Swabha Swayamdipta, Kyle Lo, Iz Beltagy, Doug Downey, and Noah A. Smith. 2020.  Don’t stop pretraining: Adapt language models to domains and tasks.
  * Kiela et al. (2020) Douwe Kiela, Hamed Firooz, Aravind Mohan, Vedanuj Goswami, Amanpreet Singh, Pratik Ringshia, and Davide Testuggine. 2020.  The hateful memes challenge: Detecting hate speech in multimodal memes.  _ArXiv_ , abs/2005.04790.
  * Li et al. (2019) Liunian Harold Li, Mark Yatskar, Da Yin, Cho-Jui Hsieh, and Kai-Wei Chang. 2019\.  Visualbert: A simple and performant baseline for vision and language.
  * Liu et al. (2019) Yinhan Liu, Myle Ott, Naman Goyal, Jingfei Du, Mandar Joshi, Danqi Chen, Omer Levy, Mike Lewis, Luke Zettlemoyer, and Veselin Stoyanov. 2019.  Roberta: A robustly optimized bert pretraining approach.
  * Lu et al. (2019) Jiasen Lu, Dhruv Batra, Devi Parikh, and Stefan Lee. 2019.  Vilbert: Pretraining task-agnostic visiolinguistic representations for vision-and-language tasks.  In _Advances in Neural Information Processing Systems_ , volume 32. Curran Associates, Inc.
  * Muennighoff (2020) Niklas Muennighoff. 2020.  Vilio: State-of-the-art visio-linguistic models applied to hateful memes.
  * Suhr et al. (2019) Alane Suhr, Stephanie Zhou, Ally Zhang, Iris Zhang, Huajun Bai, and Yoav Artzi. 2019\.  A corpus for reasoning about natural language grounded in photographs.
  * Vinyals et al. (2015) Oriol Vinyals, Alexander Toshev, Samy Bengio, and Dumitru Erhan. 2015.  Show and tell: A neural image caption generator.
  * Vinyals et al. (2017) Oriol Vinyals, Alexander Toshev, Samy Bengio, and Dumitru Erhan. 2017.  Show and tell: Lessons learned from the 2015 mscoco image captioning challenge.  _IEEE Transactions on Pattern Analysis and Machine Intelligence_ , 39(4):652–663.
  * Wang et al. (2018) Alex Wang, Amanpreet Singh, Julian Michael, Felix Hill, Omer Levy, and Samuel Bowman. 2018.  GLUE: A multi-task benchmark and analysis platform for natural language understanding.  In _Proceedings of the 2018 EMNLP Workshop BlackboxNLP: Analyzing and Interpreting Neural Networks for NLP_ , pages 353–355, Brussels, Belgium. Association for Computational Linguistics.
  * Xu et al. (2016) Kelvin Xu, Jimmy Ba, Ryan Kiros, Kyunghyun Cho, Aaron Courville, Ruslan Salakhutdinov, Richard Zemel, and Yoshua Bengio. 2016.  Show, attend and tell: Neural image caption generation with visual attention.

Supplementary Appendices333Put here for the reader’s convenience.

Baseline Source | Model | Validation | Test
---|---|---|---
|  | Baseline | Captions | Baseline | Captions
|  | Acc | Acc | Acc | Acc
Hateful Memes | human |  |  | 84.7 |
| Bert | 58.8 | 63.4 | 59.2 | 63.2
| Roberta | 59.8 | 62.2 | 60.8 | 62.7
| VisualBERT COCO | 65.06 | 67.8 | 61.7 | 64.44
| ViLBERT CC | 61.4 | 65.6 | 61.1 | 64.4

Table 3: Accuracy of captions usage for various models

##  Appendix A Additional Results

Tab. 3 depicts the accuracy of all models reported in Table 1 from the main
text.

Fig.4 presents an additional sample with better classification in most CES
models, although the generated caption inaccurate.

Figure 4: An additional sample of inaccurate caption (cellphone vs
microphone), for which CES is able to produce better classification in three
out of four models Figure 5: RoBERTa results.

Additional qualitative results for RoBERTa, BERT, ViLBERT and VisualBert are
presented in Fig.5, Fig.6, Fig.7 and Fig.8, respectively.

##  Appendix B Datasets

We use the publicly available hateful memes dataset along with Vilio’s
features dataset for UNITER implementation. Locations:

  * •

Hateful Memes Dataset can be downloaded from here: competition page. We use
Phase 1 data.

  * •

Features for UNITER can be found on Kaggle: hmtsvfeats and hmfeatureszipfin

We augment the dataset with generated captions. The augmentation code can be
found in the software supplementary.

##  Appendix C Code explanation

The code contains the following:

  * •

mmf - changes done to mmf codebase to utilize captions for BERT, RoBERTa,
VisualBert, ViLBERT.

  * •

vilio_with_captions - changes done to Vilio codebase to support UNITER with
captions.

  * •

MAX-Image-Caption-Generator - provided as is for caption generations

  * •

my_hateful_memes - code that uses MAX-Image-Caption-Generator to create
captions for the dataset and the generated captions in 2 formats - CSVs for
the dataset and datasets for BERT and RoBERTa’s pretraining (using huggingface
transformers code, not supplied)

  * •

Kaggle-Ensemble-Guide - changes to the Kaggle-Ensemble-Guide codebase to
generate ensembles for the results.

##  Appendix D Runtime and implementation details

Fine-tuning models using MMF codebase: we re-used best settings as for the
original Hateful memes challenge: 22000 updates for number of updates. We use
weighted Adam with cosine learning rate schedule and fixed 2000 warmup steps
for optimization without gradient clipping. Batch size: 32 for all models.
Learning rate: 5e-5 for all models.

Finetuning after pretraining for the BERT and RoBERTA models - we use 5e-6 vs
5e-5 learning rate.

All models are trained on a single GPU GeForce GTX TITAN X with 12GB of RAM
(same GPU in all runs below). Different models have different training times
but mostly it took 10 to 13 hours for each model.

To continue the pre-training of the language models (BERT and RoBERTA) we used
the transformers library, each pretraining was with batch size: 32 and, spread
across 4 GPUs (effective batch size: 128). Pretraining runs in about 2-3 hours
(100 epochs).

For the UNITER baseline, we use batch size:8 on a single GPU (learning rate:
1e-5). For ensembling it with a pre-trained RoBERTA model, we use batch size:6
with gradient accumulation: 2 (same learning rate). Each run takes roughly 1
hour. We used 36 features for the images.

Figure 6: BERT results. Figure 7: ViLBERT results. Figure 8: VisualBert
results.