aainstitutetext: Laboratoire de Physique Théorique et Hautes Energies (LPTHE),
UMR 7589 CNRS & Sorbonne Université, 4 Place Jussieu, F-75252, Paris,
Francebbinstitutetext: Dipartimento di Fisica, Università di Torino & INFN,
Sezione di Torino, via P. Giuria 1, I-10125 Torino, Italyccinstitutetext:
Theoretical Astrophysics Department, Fermi National Accelerator Laboratory,
Batavia, Illinois, 60510, USAddinstitutetext: Nikhef, Science Park 105, 1098
XG Amsterdam, The Netherlandseeinstitutetext: Institut d’Astrophysique de
Paris, UMR 7095 CNRS & Sorbonne Université, 98 bis boulevard Arago, F-75014
Paris, Franceffinstitutetext: Department of Physics and Astronomy, The Johns
Hopkins University, 3400 N. Charles Street, Baltimore, MD 21218,
U.S.A.gginstitutetext: Beecroft Institute for Particle Astrophysics and
Cosmology, University of Oxford, Keble Road, Oxford OX1 3RH, U.K.

# Baryogenesis from ultra-slow-roll inflation

Yi-Peng Wu a, b, c Elena Pinetti a, d Kalliopi Petraki e, f, g and Joseph Silk
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

The ultra-slow-roll (USR) inflation represents a class of single-field models
with sharp deceleration of the rolling dynamics on small scales, leading to a
significantly enhanced power spectrum of the curvature perturbations and
primordial black hole (PBH) formation. Such a sharp transition of the
inflationary background can trigger the coherent motion of scalar condensates
with effective potentials governed by the rolling rate of the inflaton field.
We show that a scalar condensate carrying (a combination of) baryon or lepton
number can achieve successful baryogenesis through the Affleck-Dine mechanism
from unconventional initial conditions excited by the USR transition. Viable
parameter space for creating the correct baryon asymmetry of the Universe
naturally incorporates the specific limit for PBHs to contribute significantly
to dark matter, shedding light on the cosmic coincidence problem between the
baryon and dark matter densities today.

##  1 Introduction

Cosmic inflation naturally gives rise to coherent production of scalar
condensates with extremely large vacuum expectation values (VEVs). These VEVs
are characterized by the Hubble scale during inflation, $H_{\ast}$, which
could be around the order of $10^{13}$ GeV Tristram:2020wbi ; Akrami:2018odb .
If some of these scalar condensates carry baryon or lepton numbers from the
Standard Model or dark sectors of the particle theories, their relaxation and
decay after inflation ends can create a non-negligible baryon (lepton)
asymmetry to the Universe. This is usually referred to as the Affleck-Dine
(AD) mechanism of baryogenesis Affleck:1984fy ; Linde:1985gh ; Dolgov:1991fr ;
Dine:2003ax .

In order to obtain a quantitative prediction about the amount of baryon
asymmetry produced by the AD mechanism from the energy scale of inflation,
there are at least two essential questions requiring clarification at the
beginning of the process, as pointed out in Dine:1995kz : (1) What are the
origins of the baryon/lepton number violating interactions in the scalar
potential, and (2) what are the initial conditions for the relaxation of
scalar condensates?

So far, the mainstream research on AD baryogenesis has been devoted to the
presence of “flat directions” in the early Universe arising from, typically,
the supersymmetric extension of the Standard Model field theories Dine:1995kz
. The necessary supersymmetry breaking at the energy scale of inflation
provides answers to both questions (1) and (2). Firstly, self-couplings of
flat directions in the Kahler potential (or the superpotential) generically
foster baryon/lepton number (B/L) violating terms that are non-renormalizable.
These B/L violating terms can have complex coefficients to serve as the source
of CP violation in the AD mechanism, yet this is not a necessary condition
since CP violation can emerge from the stochastic pick-up of scalar VEVs
during inflation Dine:1995kz ; Wu:2020pej ; Hook:2015foa . Secondly, flat
directions are in general lifted by soft supersymmetry violating terms and
higher-order non-renormalizable terms in the superpotential with parameters of
order $H_{\ast}$. These parameters determine the global minimum of the
effective potential (and thus the VEVs of flat directions). Thanks to the
tachyonic instability at the potential origin caused by soft terms, initial
VEVs for the AD mechanism may be larger than $H_{\ast}$ by $10^{2}$ $-$
$10^{3}$ times, ensuring successful baryogenesis over a wide parameter space.

In this work, we explore a novel class of initial conditions for AD
baryogenesis in the unusual limit with VEVs of the charged scalar (the AD
field) provided that they are much smaller than $H_{\ast}$. Our scenario does
not rely on a tachyonic mass term, so that the average value of the AD field
decays exponentially in time. This is similar to that of the flat direction
scenario with a minimal Kahler potential in which case the AD mechanism is
expected to fail Dine:1995kz . The key difference with the mainstream scenario
is that the AD field is dynamically excited around a later phase of inflation
so that the scalar condensate is not in an equilibrium state
Starobinsky:1994bd . With a slow but non-zero coherent motion of the AD field
by the end of inflation, we show that initial VEVs around $10^{-2}$ \-
$10^{-1}$ in units of $H_{\ast}$ are sufficient to produce the baryon
asymmetry observed today.

As a concrete demonstration of the basic idea, we consider that the background
spacetime of inflation is created by single-field models that exhibit multi-
stage transitions of the rolling dynamics of the inflaton field
Yokoyama:1998pt ; Saito:2008em ; Garcia-Bellido:2017mdw ; Germani:2017bcs ;
Motohashi:2017kbs ; Cicoli:2018asa ; Ozsoy:2018flq ; Byrnes:2018txb ;
Liu:2020oqe ; Cheng:2018qof ; Ozsoy:2019lyy ; Ballesteros:2020sre ; Ng:2021hll
; Carrilho:2019oqg ; Kannike:2017bxn . This class of models is proposed for
primordial black hole (PBH) formation on small scales with significantly
enhanced curvature perturbation induced by a sharp deceleration of inflaton
(from the primary slow-roll phase on large scales). The so-called ultra-slow-
roll (USR) inflation Kinney:2005vj ; Tsamis:2003px ; Martin:2012pe ;
Motohashi:2014ppa ; Anguelova:2017djf  is a special case of interest, which in
particular provides a convenient limit for illustration.

We show that the presence of derivatives coupling between the inflaton and the
AD field in general introduces soft effective mass terms governed by the
rolling, and more importantly, the rate of rolling of inflaton. Due to the
sharp transition of the inflaton rolling rate in the USR scenario for PBH
formation, coherent motion of the AD field condensate can be triggered during
inflation, leading to successful baryogenesis from initial conditions in the
parameter space that were usually omitted in previous investigations.

There is a crucial asset in considering baryogenesis triggered by inflation
scenarios of the USR class, since PBHs are currently under close examination
of their feasibility for contributing all (or a significant fraction) of dark
matter. Suppose that PBHs occupy more than 10 percent of the dark matter
density, their energy density today is of the same order as that of the
baryons. This could be considered as a special version of the “cosmic
coincidence problem.”

Baryogenesis initiated by USR inflation is a mechanism in parallel with that
of PBH formation, and can be realized even without the presence of any PBHs
today. Thus it is not identical to the idea where baryon asymmetry of the
Universe is created by existing PBHs Turner:1979bt ; Barrow:1990he ;
Majumdar:1995yr ; Baumann:2007yr ; Bambi:2008hp ; Hamada:2016jnq ;
Hooper:2020otu ; Ambrosone:2021lsx ; JyotiDas:2021shi , providing a different
insight to the cosmic coincidence between dark matter and baryon densities
sourced by PBHs Fujita:2014hha ; Garcia-Bellido:2019vlf  or originated from
the asymmetric dark sector in particle physics Petraki:2013wwa ;
vonHarling:2012yn ; Bell:2011tn ; Flores:2020drq . Based on the assumption of
PBH dark matter, the scenario investigated in this work allows a more
systematic study on the viable parameter space for the cosmic coincidence
problem as model parameters of the USR inflation simultaneously control the
dark matter density and the baryon asymmetry of the Universe.

This work is organized as follows: in Section 2, we outline the structure of
the inflaton dynamics in the (quasi) USR limit, the template of the
corresponding power spectrum of the curvature perturbation, and the coherent
motion of a massive scalar field under such an inflationary background. In
Section 3, we provide a simple realization of the AD field with effective mass
terms controlled by parameters of the USR inflation. We then compute the AD
field evolution during inflation and estimate the amount of baryon asymmetry
created at the end of inflation. These results are initial conditions for
Section 4 to obtain the final baryogenesis in the radiation-dominated era of
the Universe. Finally, in Section 5 we take into account the associated PBH
abundance in terms of the USR parameters based on several statistical
approaches that address different uncertainties about the PBH formation from
inflation. We explore the viable parameter space for PBHs to be a significant
contributor to dark matter with the correct amount of baryon asymmetry.
Conclusions for baryogenesis and implications for the cosmic coincidence
problem are given in Section 6.

##  2 Inflation and primordial black holes

Models of inflation investigated in this work are motivated by primordial
black hole (PBH) production via an enhanced power spectrum $P_{\zeta}$ of the
curvature perturbation $\zeta$ on small scales. In single-field inflationary
framework, the enhancement of $P_{\zeta}$ can be realized by a sharp
deceleration of the inflaton field, $\phi$, in its potential $V(\phi)$ that
causes the decrease of the first slow-roll parameter
$\epsilon_{H}\equiv-\dot{H}/H^{2}\approx\epsilon_{V}$ by many orders of
magnitude, where $\epsilon_{V}\equiv(V_{,\phi}/V)^{2}M_{P}^{2}$ denotes the
first slow-roll parameter for the inflaton potential. Throughout this work, we
will keep in mind the key parameter that controls such a deceleration, defined
as

$\delta=\frac{\ddot{\phi}}{H\dot{\phi}},$ (1)

where $\delta\approx\epsilon_{V}-\eta_{V}$ and $\eta_{V}\equiv
V_{,\phi\phi}/VM_{P}^{2}$ denotes the second slow-roll parameter for the
potential. $\delta$ can be $\mathcal{O}(1)$ or larger due to the significant
violation of the slow-roll conditions from $\eta_{V}$. We restrict ourselves
to models that satisfy the condition $\epsilon_{H}\approx\epsilon_{V}\ll 1$
where the background is well described by the de Sitter space.

###  2.1 Ultra-slow-roll (USR) templates

Figure 1: Analytic template for the power spectrum of the curvature
perturbation for a 3-stage constant-rate inflation with Phase 1 $\delta_{1}=0$
of the primary slow-roll and Phase 2 $\delta_{2}=-3.05$ in the (quasi-)ultra-
slow-roll limit.

To illustrate the basic idea for baryogenesis, we will assume that inflation
is composed of three essential phases with a nearly constant rate-of-rolling
($\delta\approx$ const.) during each period. The analytic structure of
$P_{\zeta}$ with the constant-rate condition $\dot{\delta}=0$ has been well
studied in multi-stage inflation Byrnes:2018txb ; Liu:2020oqe ; Cheng:2018qof
; Ozsoy:2019lyy ; Ballesteros:2020sre ; Ng:2021hll ; Carrilho:2019oqg . The
first stage of inflation is the primary slow-roll phase with a nearly scale-
invariant power spectrum probed by cosmic microwave background (CMB)
observations and we take $\delta_{1}=\delta(t<t_{0})\rightarrow 0$ for
simplicity. Superhorizon enhancement of $P_{\zeta}$ occurs when the second
stage ($t_{0}\leq t<t_{\ast}$) has a constant rate $\delta_{2}<-3/2$
Liu:2020oqe . Except for further notification, in the following discussion we
focus on a second stage of inflation close to (but not exactly in) the ultra-
slow-roll (USR) limit Kinney:2005vj ; Tsamis:2003px ; Martin:2012pe ;
Motohashi:2014ppa ; Anguelova:2017djf , where

$\delta_{2}=\frac{\ddot{\phi}}{H\dot{\phi}}\lesssim-3,\quad t_{0}<t<t_{\ast}.$
(2)

To suppress the non-Gaussian tail in the high-sigma limit of $\zeta$ led by
quantum diffusion in the exact USR limit ($\delta_{2}\rightarrow-3$)
Biagetti:2018pjj ; Ezquiaga:2019ftu ; Pattison:2021oen ; Figueroa:2020jkf ;
Biagetti:2021eep , we put an upper bound to $\delta_{2}$ as discussed in
Appendix A. Note that for $\delta\leq-3$, the curvature perturbation $\zeta$
is led by the late-time scaling dimension (or conformal weight) $\Delta\equiv
3/2-\sqrt{9/4+3\delta+\delta^{2}}=3/2-|3/2+\delta|$, and thus the
dimensionless power spectrum $P_{\zeta}\sim k^{2\Delta}\rightarrow k^{0}$ as
$\delta\rightarrow-3$. The USR limit is therefore a specific case in the slow-
roll violation regime that can give rise to a scale-invariant power spectrum
as that from the standard slow-roll inflation ($\delta\rightarrow 0$).

The momentum scaling of the power spectrum $P_{\zeta}(k)$ from a 3-stage
constant rate inflation can be summarized as Liu:2020oqe ; Ng:2021hll :

$\displaystyle P_{\zeta}=\left\\{\begin{array}[]{ll}A_{\rm CMB}&k<k_{\rm
min},\\\ A_{\rm PBH}(k/k_{0})^{4},&k_{\rm min}<k<k_{0},\\\ A_{\rm
PBH}(k/k_{0})^{6+2\delta_{2}},&k>k_{0},\end{array}\right.$ (6)

where $A_{\rm CMB}\approx H_{\ast}^{2}/(\epsilon_{\rm CMB}M_{P}^{2})$ measured
on CMB scales has negligible contribution for PBH formations, and the enhanced
spectral amplitude

$A_{\rm PBH}\approx\frac{H_{\ast}^{2}}{\epsilon_{\rm
CMB}M_{P}^{2}}\left(\frac{k_{0}}{k_{\ast}}\right)^{6+4\delta_{2}}.$ (7)

Here $M_{P}$ is the reduced Planck mass and $H_{\ast}$ is the Hubble parameter
during inflation and $0\leq\epsilon_{\rm CMB}<0.0063$ Akrami:2018odb . In
Figure 1 we plot the case in the USR limit with $\delta_{2}=-3$. We choose
$N=\ln a=0$ at the time $t=t_{0}$ where $k_{0}=a(t_{0})H_{\ast}$ crosses the
horizon so that $A_{\rm PBH}/A_{\rm CMB}=e^{-N_{\ast}(6+4\delta_{2})}$. Models
of inflation that realize the power spectrum given in Figure 1 are sometimes
called “punctuated inflation.” Ragavendra:2020sop ; Jain:2008dw ;
Allahverdi:2006we ; Jain:2009pm  See also Byrnes:2018txb ; Biagetti:2018pjj ;
Motohashi:2019rhu ; Ragavendra:2020sop  for PBHs from inflationary scenarios
that exhibit a secondary slow-roll phase.

The transition from the primary slow-roll phase (Phase 1) to the USR phase
(Phase 2) shows an apparent violation of the continuity of the time derivative
for the leading $\zeta$ mode with $\Delta_{1}=0$ on superhorizon scales. This
is due to the efficient entropy production driven by the sharp deceleration of
inflaton that violates the adiabaticity of the curvature perturbation
Ragavendra:2020sop ; Ng:2021hll ; Leach:2000yw . Subleading entropy modes with
the next-to-lowest scaling dimension come to dominate the power spectrum in
the range of $k_{\rm min}\leq k<k_{0}$ with the special power $P_{\zeta}\sim
k^{4}$, where $k_{\rm min}\approx k_{0}(k_{0}/k_{\ast})^{3/2}$ for
$\delta_{2}=-3$.111The leading scaling dimension for the curvature
perturbation $\zeta$ in slow-roll is the same as in USR ($\Delta_{\rm
SR}=\Delta_{\rm USR}=0$). However, the leading scaling dimension for the time
derivative $\dot{\zeta}$ in USR is $\Delta_{\rm USR}-1$, yet it is
$\Delta_{\rm SR}+1$ in slow-roll. Thus the steepest growth $P_{\zeta}\sim
k^{4}$ occurs only for the slow-roll $\rightarrow$ USR transition
Byrnes:2018txb  and it does not occur for the inverse USR $\rightarrow$ slow-
roll transition Cai:2017bxr . For a more detailed discussion, see Ng:2021hll .
The $k^{4}$ scaling of $P_{\zeta}$ is firstly identified in USR inflation by
Leach-Sasaki-Wands-Liddle Leach:2001zf , and this mechanism is also called as
the steepest growth of the power spectrum Byrnes:2018txb ; Carrilho:2019oqg .

The transition from the USR phase to the final phase (Phase 3) with a non-
negative rate $\delta_{3}\geq 0$ is necessary for stopping the growth of the
power spectrum. For inflation models that can realized a stable negative rate,
$\delta_{2}\leq-3$ with $\dot{\delta}_{2}\approx 0$, in Phase 2, the value of
$\delta_{3}$ is constrained by the continuity of the scaling dimension
$\Delta$, where the condition $\Delta_{2}=\Delta_{3}=-\delta_{3}$ indicates
that $\delta_{3}=-\delta_{2}-3$. This is due to the non-violation of the
adiabatic condition in the acceleration phase of inflaton, and the effective
mass of inflaton retains continuous Ng:2021hll . Thus if Phase 2 goes to the
USR limit with $\delta_{2}\rightarrow-3$, one finds that
$\delta_{3}\rightarrow 0$ where Phase 3 must be a secondary slow-roll phase.

The termination of the secondary slow-roll phase in this class of models
relies an additional assumptions Ragavendra:2020sop ; Jain:2008dw ;
Allahverdi:2006we ; Jain:2009pm . We will consider that inflation is
terminated by (1) the sudden decay of inflaton into radiation at $N_{\rm end}$
and drives instantaneous reheating (Section 3) or (2) the coherent
oscillations of the inflaton as it rolls into a deep valley and obtains a
large mass (Section 4).

###  2.2 Massive scalars under USR transition

Before going to the concrete scenario for baryogenesis, we derive in this
section some useful formulae that describe the coherent motion of massive
scalar fields during the multi-stage constant-rate inflation illustrated
above, focusing on the limit where Phase 2 is USR. We are interested in a
scenario where a spectator field, $\sigma$, with negligible energy density
during inflation receives different effective masses in each phase, where
$m_{\sigma}=m_{i}$ for $i=1,2,3$ in Phase $i$. The vacuum expectation value
(VEV) in Phase 1 is developed naturally by the stochastic effect in
equilibrium Starobinsky:1994bd . This gives the initial conditions for Phase 2
as

$\displaystyle\sigma_{1}(t_{0})=\sqrt{\frac{3}{8\pi^{2}}}\frac{H_{\ast}^{2}}{m_{1}},\quad\dot{\sigma}_{1}(t_{0})=0,$
(8)

where the VEV is the result of a massive non-interacting scalar.

With the transition from the primary slow-roll phase to USR, $\sigma$ gains a
different mass $m_{2}\neq m_{1}$ which suddenly drives the coherent condensate
out of equilibrium in the potential. This leads to the evolution of the VEV
with respect to a simplified equation of motion as

$\displaystyle\ddot{\sigma}_{2}+3H_{\ast}\dot{\sigma}_{2}+m_{2}^{2}\,\sigma_{2}=0,\qquad
t_{0}\leq t<t_{\ast}.$ (9)

According to the initial conditions given by (8), the solution of $\sigma_{2}$
reads Antoniadis:2011ib :

$\displaystyle\sigma_{2}(N)=\sqrt{\frac{3}{8\pi^{2}}}\frac{H_{\ast}^{2}}{m_{1}}\frac{1}{2\nu_{2}}\left(\Delta_{2}^{+}e^{-\Delta_{2}^{-}N}-\Delta_{2}^{-}e^{-\Delta_{2}^{+}N}\right),$
(10)

where we use $N=H_{\ast}t$ with $N_{0}=H_{\ast}t_{0}=0$, and the scaling
dimensions (conformal weights) for a massive scalar in the de Sitter space
Strominger:2001pn ; Ng:2012xp ; Jafferis:2013qia ; Antoniadis:2011ib  are
given by

$\displaystyle\Delta_{2}^{\pm}\equiv\frac{3}{2}\pm\nu_{2},\qquad\nu_{2}=\sqrt{\frac{9}{4}-\frac{m_{2}^{2}}{H_{\ast}^{2}}}.$
(11)

Note that this is the leading scaling dimension for the coherent motion of
$\sigma$ (to be distinguished from those for the curvature perturbation
$\zeta$ in Section 2.1). With
$\Delta_{2}^{-}\Delta_{2}^{+}=9/4-\nu_{2}^{2}=m_{2}^{2}/H_{\ast}^{2}$, the
time-derivative of $\sigma_{2}$ reads

$\displaystyle\frac{d\sigma_{2}}{dN}=\frac{\dot{\sigma}_{2}}{H_{\ast}}=\sqrt{\frac{3}{8\pi^{2}}}\frac{m_{2}^{2}}{m_{1}}\frac{1}{2\nu_{2}}\left(e^{-\Delta_{2}^{+}N}-e^{-\Delta_{2}^{-}N}\right),$
(12)

which gives non-zero initial conditions for Phase 3 at $N=N_{\ast}$.

After the scalar condensate is excited by the slow-roll to USR transition, the
equation of motion for $\sigma$ from USR to the secondary slow-roll phase
(Phase 3) is unchanged but the non-zero time-derivative in Phase 2 asks the
solution to take the general form as

$\displaystyle\sigma_{3}(N)=Ce^{-\Delta_{3}^{+}N}+De^{-\Delta_{3}^{-}N},$ (13)

where similarly, the leading scaling dimension is fixed by the mass in Phase 3
as

$\displaystyle\Delta_{3}^{\pm}\equiv\frac{3}{2}\pm\nu_{3},\qquad\nu_{3}=\sqrt{\frac{9}{4}-\frac{m_{3}^{2}}{H_{\ast}^{2}}}.$
(14)

Matching boundary conditions at the time slice $N=N_{\ast}$ at the end of USR,
one finds

$\displaystyle C=\sqrt{\frac{3}{8\pi^{2}}}\frac{H_{\ast}^{2}}{m_{1}}$
$\displaystyle\frac{\Delta_{2}^{-}}{\Delta_{2}^{-}-\Delta_{2}^{+}}\frac{e^{\Delta_{3}^{+}N_{\ast}}}{\Delta_{3}^{-}-\Delta_{3}^{+}}$
$\displaystyle\times\left[\left(\Delta_{3}^{-}-\Delta_{2}^{+}\right)e^{-\Delta_{2}^{+}N_{\ast}}+\Delta_{2}^{+}\left(1-\frac{\Delta_{3}^{-}}{\Delta_{2}^{-}}\right)e^{-\Delta_{2}^{-}N_{\ast}}\right],$
(15) $\displaystyle D=\sqrt{\frac{3}{8\pi^{2}}}\frac{H_{\ast}^{2}}{m_{1}}$
$\displaystyle\frac{\Delta_{2}^{-}}{\Delta_{2}^{-}-\Delta_{2}^{+}}\frac{e^{\Delta_{3}^{-}N_{\ast}}}{\Delta_{3}^{+}-\Delta_{3}^{-}}$
$\displaystyle\times\left[\left(\Delta_{3}^{+}-\Delta_{2}^{+}\right)e^{-\Delta_{2}^{+}N_{\ast}}+\Delta_{2}^{+}\left(1-\frac{\Delta_{3}^{+}}{\Delta_{2}^{-}}\right)e^{-\Delta_{2}^{-}N_{\ast}}\right].$
(16)

For a constant mass $m_{i}/H_{\ast}\leq 3/2$, $\nu_{i}$ is real so that the
positive branch $\Delta^{+}$ decays much faster than the negative branch
$\Delta^{-}$. For $m_{i}/H_{\ast}>3/2$, $\nu_{i}$ or $\Delta^{\pm}_{i}$
contain imaginary parts so that both $\pm$ branches contribute to the late-
time dynamics of the mass scalar. In this work, we focus on the mass range
$m_{i}/H_{\ast}\leq 3/2$ and thus only the negative branch survives at the end
of inflation as

$\displaystyle\sigma_{3}(N_{\rm end})\rightarrow De^{-\Delta_{3}^{-}N_{\rm
end}},\qquad\left.\frac{d\sigma_{3}}{dN}\right|_{N_{\rm
end}}\rightarrow-\Delta_{3}^{-}De^{-\Delta_{3}^{-}N_{\rm end}}.$ (17)

The results in (17) are the two most relevant quantities for the computation
of baryon production from a charged scalar field at or after the end of
inflation.

##  3 Baryogenesis during inflation

In most of the cases, baryon or lepton number production during inflation is
considered negligible due to the exponential dilution of the particle number
density in time. However, for baryon production triggered by the phase
transition of USR at $N=N_{0}$ on scales much smaller than those constrained
by CMB observations, the duration $\Delta N=N_{\rm end}-N_{0}$ is less
constrained as long as the total duration $N_{\rm end}-N_{\rm
CMB}\gtrsim\mathcal{O}(50)$.

Baryogenesis investigated in this section is driven by a charged scalar
$\sigma$ that possesses a global $U(1)$ particle number and develops non-zero
expectation values during inflation. The baryon or lepton number excitation is
triggered by the phase transition of the USR dynamics and it is essentially
described by the Affleck-Dine (AD) mechanism Affleck:1984fy , except for the
fact that we evaluate the process in the (quasi-)de Sitter space.

###  3.1 Asymmetric scalar masses

To illustrate how the transition of inflaton rolling rate $\delta$ given by
(1) can lead to a transition of effective masses for a charged scalar $\sigma$
(the AD field), let us consider

$\displaystyle\mathcal{L}=\mathcal{L}_{\phi}+|\partial\sigma|^{2}+m_{\sigma}^{2}|\sigma|^{2}+\Delta\mathcal{L},$
(18)

where $\mathcal{L}_{\phi}$ describes the multi-stage constant-rate (including
the USR) model for PBH formation. $\Delta\mathcal{L}$ includes self-
interactions of $\sigma$ and derivative couplings with $\phi$ of the form

$\displaystyle\Delta\mathcal{L}\supset\frac{c_{1}}{\Lambda}\left|\sigma^{2}\right|\square\phi+\frac{c_{2}}{\Lambda}\partial_{\mu}\phi\left[\sigma\partial^{\mu}\sigma+\sigma^{\ast}\partial^{\mu}\sigma^{\ast}\right]+\frac{c_{3}}{\Lambda^{2}}\left(\partial\phi\right)^{2}|\sigma|^{2},$
(19)

where we have imposed C and CP to the system so that $c_{1}$, $c_{2}$ and
$c_{3}$ are $\mathcal{O}(1)$ real constants.222In general, $c_{1}$, $c_{2}$
and $c_{3}$ can have imaginary components given by non-zero phases. However,
the combination of $c_{1}$ and $c_{3}$ terms with their complex conjugates
only gives a redefinition of $c_{1}$ and $c_{3}$, while the phase of $c_{2}$
can be absorbed into $\sigma$ via a constant shift of its phase term. The
derivative couplings used in (19) are similar to those of the inflaton-induced
chemical potential for an enhanced charged scalar production Wang:2019gbi ;
Wang:2020ioa ; Bodas:2020yho . However, CP-violating couplings, such as
$\partial_{\mu}\phi(\sigma^{\ast}\partial^{\mu}\sigma-\sigma\partial^{\mu}\sigma^{\ast})$
Bodas:2020yho  is forbidden for simplicity. The necessary C/CP violation for
baryogenesis is spontaneously realized by the initial condensation of the
charged scalar at the VEVs during inflation.

One can rewrite the Lagrangian for $\sigma$ with respect to the mass
eigenstates $\sigma_{\pm}$, based on the decomposition
$\sigma\equiv(\sigma_{-}+i\sigma_{+})/\sqrt{2}$, to obtain a system for two
real scalars as

$\displaystyle\mathcal{L}_{\sigma}=$
$\displaystyle\frac{1}{2}\left(\partial\sigma_{-}\right)^{2}+\frac{1}{2}\left[m_{\sigma}^{2}+\frac{c_{1}-c_{2}}{\Lambda}\square\phi+\frac{c_{3}}{\Lambda^{2}}\left(\partial\phi\right)^{2}\right]\sigma_{-}^{2}$
$\displaystyle+\frac{1}{2}\left(\partial\sigma_{+}\right)^{2}+\frac{1}{2}\left[m_{\sigma}^{2}+\frac{c_{1}+c_{2}}{\Lambda}\square\phi+\frac{c_{3}}{\Lambda^{2}}\left(\partial\phi\right)^{2}\right]\sigma_{+}^{2}.$
(20)

The mass eigenstates provide a convenient basis for the study of baryon number
production, and, unless with a fine tuning of the coefficients, they in
general have non-degenerated effective masses

$\displaystyle m_{\pm}^{2}=m_{\sigma}^{2}+\frac{c_{1}\pm
c_{2}}{\Lambda}\square\phi+\frac{c_{3}}{\Lambda^{2}}\left(\partial\phi\right)^{2}.$
(21)

$\Lambda$ is the cutoff scale for the effective field theory that satisfies
$H_{\ast}\ll\Lambda\leq M_{P}$.

We now apply the multi-stage constant-rate inflation described in Section 2.1
to specify the spatially homogeneous dynamics of $\phi$. With a constant
$\delta$, we have
$\square\phi=-\ddot{\phi}-3H\dot{\phi}\approx-(\delta+3)\sqrt{2\epsilon_{H}}M_{P}H_{\ast}^{2}$
and $(\partial\phi)^{2}=\dot{\phi}^{2}\approx
2\epsilon_{H}M_{P}^{2}H_{\ast}^{2}$. If $\delta\neq 0$, the first slow-roll
parameter $\epsilon_{H}$ is evolving in each phase where
$\epsilon_{2}=\epsilon_{\rm CMB}e^{2\delta_{2}N}$ in Phase 2 and
$\epsilon_{3}=\epsilon_{\ast}e^{2\delta_{3}(N-N_{\ast})}$ in Phase 3 with
$\epsilon_{\ast}=\epsilon_{\rm CMB}e^{2\delta_{2}N_{\ast}}$. As a result, the
effective masses for $\sigma_{\pm}$ change with the rolling rate of inflaton
in different phases:

$\displaystyle m_{i\pm}^{2}=m_{\sigma}^{2}+\frac{c_{1}\pm
c_{2}}{\Lambda}[-(\delta_{i}+3)]\sqrt{2\epsilon_{i}}M_{P}H_{\ast}^{2}+\frac{c_{3}}{\Lambda^{2}}2\epsilon_{i}M_{P}^{2}H_{\ast}^{2}.$
(22)

For the USR scenario of our interest ($\delta_{1}=0$, $\delta_{2}=-3$ and
$\delta_{3}=-\delta_{2}-3=0$), the masses $m_{1\pm}$ and $m_{3\pm}$ in Phase 1
& 3 are always constants. The time-dependence of $m_{2\pm}$ replies on the
choices of parameters, and there are two common assumptions for the mass
spectrum of the AD field:

Heavy AD field. Charged scalars arising as the superpartners of the Standard
Model fermions are natural candidates for the AD field Dine:1995kz . The
necessary supersymmetry breaking generically induces soft masses for those
charged scalars of order the Hubble constant $H_{\ast}$ due to inflation.
Baryon-violating terms residing in the superpotential, a non-minimal coupling
with gravity or the stochastic effect of quantum fluctuations during inflation
can also introduce effective masses related to $H_{\ast}$ Linde:1985gh ;
Dolgov:1991fr ; Dine:2003ax . If the scalar mass is associated with the Hubble
scale of inflation, namely $m_{\sigma}\sim\mathcal{O}(H_{\ast})$, we can find
nearly constant masses in each phase where

$\displaystyle m_{1\pm}^{2}$ $\displaystyle=m_{\sigma}^{2}-3\frac{c_{1}\pm
c_{2}}{\Lambda}\sqrt{2\epsilon_{\rm
CMB}}M_{P}H_{\ast}^{2}+\frac{c_{3}}{\Lambda^{2}}2\epsilon_{\rm
CMB}M_{P}^{2}H_{\ast}^{2},$ $\displaystyle m_{2\pm}^{2}$
$\displaystyle=m_{\sigma}^{2}+\frac{c_{3}}{\Lambda^{2}}2\epsilon_{2}M_{P}^{2}H_{\ast}^{2}\approx
m_{\sigma}^{2},$ (23) $\displaystyle m_{3\pm}^{2}$
$\displaystyle=m_{\sigma}^{2}+\delta_{2}\frac{c_{1}\pm
c_{2}}{\Lambda}\sqrt{2\epsilon_{\ast}}M_{P}H_{\ast}^{2}+\frac{c_{3}}{\Lambda^{2}}2\epsilon_{\ast}M_{P}^{2}H_{\ast}^{2},$

where $\delta_{3}=-\delta_{2}-3$ is used for $m_{3\pm}$. The approximation
$m_{2\pm}\approx m_{\sigma}^{2}$ is valid as long as
$m_{\sigma}^{2}/H_{\ast}^{2}\gg\epsilon_{\rm CMB}M_{P}^{2}/\Lambda^{2}$, which
implies a lower bound for the cutoff $\Lambda\gg\epsilon_{\rm CMB}M_{P}$. In
this case the time evolution of $\sigma_{\pm}$ is directly given by the
results found in Section 2.2. Note that $m_{\sigma}\sim H_{\ast}$ also ensures
that the density fluctuation of $\sigma$ during inflation merely contributes
to negligible isocurvature perturbations in the CMB background radiations.

Massless AD field. Another usual assumption for the bare mass parameter
$m_{\sigma}$ is of order the electroweak scale ($m_{\sigma}=m_{3/2}$), which
has negligible contribution at the energy scale of inflation. The special
massless limit ($m_{\sigma}\rightarrow 0$) can harvest a non-negligible baryon
asymmetry at the end of inflation which is investigated in Appendix B.

A possible scenario towards the UV completion of (19) would be to consider a
vector field $A^{\mu}$ that couples to the AD field and a scalar $\varphi$
through the interactions

$\displaystyle\Delta\mathcal{L}_{A}\supset c_{A1}\varphi
A^{\mu}\partial_{\mu}\varphi+c_{A2}\sigma
A^{\mu}\partial_{\mu}\sigma+c_{A3}\sigma^{\ast}A^{\mu}\partial_{\mu}\sigma+c.c.,$
(24)

where $\varphi$ should be identified with the inflaton via the redefinition
$\Lambda_{A}\phi\equiv\varphi^{2}/2$. Since inflaton is real scalar, (24)
cannot be a gauge theory. $A^{\mu}$ shall be heavy with a mass
$m_{A}\sim\Lambda$ so that it can be integrated out at the energy scale of
inflation. This results in the effective theory $c_{A1}\varphi
A^{\mu}\partial_{\mu}\varphi\equiv
c_{A1}\Lambda_{A}A^{\mu}\partial_{\mu}\phi\supset-
c_{A1}c_{A2}\frac{\Lambda_{A}}{m_{A}^{2}}\sigma\partial^{\mu}\sigma\partial_{\mu}\phi-
c_{A1}c_{A3}\frac{\Lambda_{A}}{m_{A}^{2}}\sigma^{\ast}\partial^{\mu}\sigma\partial_{\mu}\phi+\cdots$.
The only choice for the effective theory to reproduce (19) is to make sure
$\Lambda_{A}=\Lambda$ where one can recognize $c_{1}\sim c_{A1}c_{A3}$,
$c_{2}\sim c_{A1}c_{A2}$ and the $c_{3}$ term becomes higher-order
corrections. Although there are analogous terms to
$\frac{1}{\Lambda^{2}}\sigma^{2}(\partial\sigma)^{2}$,
$\frac{1}{\Lambda^{2}}|\sigma|^{2}|\partial\sigma|^{2}$, $\cdots$ arise in
this scenario that introduce non-canonical couplings to the mass eigenstates,
one can check that these terms of $\Lambda^{-2}$ are indeed subdominant
corrections to (19) with the proper choice of the cutoff
$H_{\ast}\ll\Lambda\leq M_{P}$.

There is a simple constraint $\mathcal{L}_{\phi}\approx
3M_{P}^{2}H_{\ast}^{2}\gg\Delta\mathcal{L}$ to ensure the de Sitter background
of single-field inflation. The background assumption
$\dot{\phi}\sim\sqrt{\epsilon_{H}}M_{P}H_{\ast}$ and the initial condition
$\sigma_{\pm}(t_{0})\sim H_{\ast}^{2}/m_{\pm}$ indicate that the $c_{1}$ or
$c_{2}$ term has
$\Delta\mathcal{L}\sim\sqrt{\epsilon_{H}}\frac{M_{P}}{\Lambda}\frac{H_{\ast}^{6}}{m_{\pm}^{2}}$,
and the $c_{3}$ term has
$\Delta\mathcal{L}\sim\epsilon_{H}\frac{M_{P}^{2}}{\Lambda^{2}}\frac{H_{\ast}^{6}}{m_{\pm}^{2}}$.
Thus $\Delta\mathcal{L}\ll H_{\ast}^{4}\ll\mathcal{L}_{\phi}$ is naturally
consistent with the constant mass approximation
$m_{\sigma}^{2}/H_{\ast}^{2}\gg\epsilon_{\rm CMB}M_{P}^{2}/\Lambda^{2}$ used
in (3.1).

###  3.2 Baryon number production

The sudden transition of scalar masses induced by the change of inflaton
rolling rate $\delta$ from slow-roll to USR triggers the coherent motion of
the $\sigma$ VEV. As $\sigma$ starts in motion, the conserved current
$j^{\mu}=i(\sigma^{\ast}\partial_{\mu}\sigma-\sigma\partial_{\mu}\sigma^{\ast})$
associated with a global $U(1)$ quantum number in the limit of
$\Delta\mathcal{L}\rightarrow 0$ is produced from zero. We will assign such a
$U(1)$ to the baryon number for a typical example, despite that in general the
$U(1)$ can be a combination of baryon, lepton or even quantum numbers in the
dark sector Petraki:2013wwa ; Bell:2011tn ; vonHarling:2012yn . In terms of
mass eigenstates, the baryon number is

$\displaystyle
n_{B}=j^{0}=\sigma_{+}\dot{\sigma}_{-}-\sigma_{-}\dot{\sigma}_{+}=H_{\ast}\left(\sigma_{+}\frac{d\sigma_{-}}{dN}-\sigma_{-}\frac{d\sigma_{+}}{dN}\right).$
(25)

Note that it is only the $c_{2}$ term in (19) that violates the baryon number.

To estimate the baryon asymmetry at the end of inflation, in this section we
assume instantaneous reheating soon after $N_{\rm end}$, where the radiation
energy density $\rho_{r}(t_{\rm end})\approx 3M_{P}^{2}H_{\ast}^{2}$. The
temperature $T(t)=[\frac{30}{\pi^{2}g_{\ast}}\rho_{r}(t)]^{1/4}$ at the
beginning of radiation domination is therefore given by $T_{\ast}\equiv
T(H_{\ast})$ and the corresponding entropy production
$s_{\ast}=2\pi^{2}g_{\ast}T_{\ast}^{3}/45$, where the number of relativistic
degrees of freedom above $300$ GeV is $g_{\ast}=106.75$.

Based on the results given by (17) for each mass eigenstate, we can compute
$n_{B}$ to find out the baryon asymmetry as

$\displaystyle Y_{B}=\frac{n_{B}(N_{\rm
end})}{s_{\ast}}=\frac{H_{\ast}}{s_{\ast}}D_{+}D_{-}\left(\Delta_{3+}^{-}-\Delta_{3-}^{-}\right)e^{-(\Delta_{3+}^{-}+\Delta_{3-}^{-})N_{\rm
end}},$ (26)

where $D_{\pm}$ and $\Delta_{3\pm}^{-}$ are the coefficient (2.2) and the
Phase 3 scaling dimension (14) with respect to the mass eigenstates
$m_{3\pm}$. For a large duration number $N_{\rm end}\gg
1/(\Delta_{3+}^{-}+\Delta_{3-}^{-})$, the exponential suppression of $Y_{B}$
illustrates the conventional expectation for a negligible baryon asymmetry
from long-term inflation.

Figure 2: The parameter scan for the baryon asymmetry $|Y_{B}|$ at the end of
inflation from a heavy AD field ($m_{\sigma}=H_{\ast}/2$) in the USR limit
$\delta_{2}=-3$ with $\\{c_{1},c_{2},c_{3}\\}=\\{-2,1,1\\}$, $N_{\ast}=0.5$
and $H_{\ast}=2.37\times 10^{13}$ GeV. Dashed lines are contours of $|Y_{B}|$.
Colored regions are excluded by validity of the late-time approximation for
$m_{3\pm}/H_{\ast}<3/2$ and the constant-mass formalism (3.1).

Let us find out the typical $|Y_{B}|$ at the end of inflation from a heavy AD
field. For $m_{\sigma}\sim H_{\ast}$, we apply constant-mass approximation
(3.1) to compute the evolution of each mass eigenstate, where the solutions
for $\sigma_{\pm}$ at the end of inflation are readily given in Section 2.2.
With $H_{\ast}\sim 10^{13}$ GeV, the instantaneous reheating gives
$H_{\ast}^{3}/s_{\ast}\sim 10^{-9}$ and with $\Lambda\sim M_{P}$,
$D_{+}D_{-}/H_{\ast}^{2}\sim 10^{-2}$ for arbitrary choices of $N_{\ast}$ and
$N_{\rm end}$. $|\Delta_{3+}^{-}-\Delta_{3-}^{-}|\gtrsim 10^{-1}$ if
$N_{\ast}\lesssim 0.5$ and $|\Delta_{3+}^{-}-\Delta_{3-}^{-}|\lesssim 10^{-3}$
if $N_{\ast}\gtrsim 2.5$. (Note that $N_{\ast}>2.5$ is usually required for
the byproduct PBHs from the USR inflation to be important dark matter, see the
discussion Section 5.)

As a result, for heavy AD mass $0.2<m_{\sigma}/H_{\ast}<1.5$, our numerical
results show that the maximal baryon asymmetry given by (26) is realized at
the order of

$\displaystyle\left|Y_{B}\right|(m_{\sigma},N_{\ast},N_{\rm end})\lesssim
10^{-11},$ (27)

for short term USR inflation with $N_{\ast}<0.5$ and $N_{\rm end}<5$ that
gives $e^{-(\Delta_{3+}^{-}+\Delta_{3-}^{-})N_{\rm end}}\sim\mathcal{O}(1)$. A
parameter scan for $|Y_{B}|$ at the end of inflation from a heavy AD field
with $m_{\sigma}=H_{\ast}/2$ is shown in Figure 2. We conclude that for a
heavy charged scalar with a bare mass $m_{\sigma}\sim H_{\ast}$, the observed
baryon asymmetry cannot be purely explained by its dynamical excitation
triggered by the USR phase transition during inflation. However, when taking
into account the coherent oscillation of inflaton after the end of inflation
in the usual (non-instantaneous) reheating scenario, a heavy scalar can still
generate sufficient baryon asymmetry through the standard AD mechanism. We
will explore baryogenesis after inflation via reheating in Section 4.

The correlation length. Here we address the correlation length (or the
variance) of the local baryon asymmetry associated with the initial
conditions. First of all, it is remarkable that the choice of (8) for the mass
eigenstates in Phase 1 correspond to VEVs with maximal CP violation. The
averaged baryon asymmetry among the probability distribution of all VEVs is
expected to be dominated by those values with maximal CP violation Wu:2020pej
. To see what is the meaning of maximal CP violation, let us apply the polar
representation $\sigma=Re^{i\theta}/\sqrt{2}$ to the AD field, where the
radial VEV is given by $\langle
R^{2}\rangle=\langle\sigma_{-}^{2}\rangle+\langle\sigma_{+}^{2}\rangle$ and
the angular VEV can be computed via the relation
$\tan\theta=\sigma_{+}/\sigma_{-}$. For heavy AD field considered here,
$m_{1\pm}$ given by (3.1) is dominated by $m_{\sigma}$ so that
$|\sigma_{+}/\sigma_{-}|=m_{1-}/m_{1+}\approx 1$, which implies that the
initial VEV of $\theta$ is approximately at $\pm\pi/4$ or $\pm 3\pi/4$.

In terms of the polar representation, the equation of motion reads

$\displaystyle\ddot{R}+3H\dot{R}+\left[m_{\sigma}^{2}+\frac{c_{1}}{\Lambda}\square\phi+\frac{c_{3}}{\Lambda^{2}}(\partial\phi)^{2}-\dot{\theta}^{2}\right]R=\frac{c_{2}}{\Lambda}\square\phi
R\cos(2\theta),$ (28)
$\displaystyle\ddot{\theta}+\left(3H+2\frac{\dot{R}}{R}\right)\dot{\theta}=-\frac{c_{2}}{\Lambda}\square\phi\sin(2\theta).$
(29)

One can see that initial conditions with $\theta=\theta_{\rm
max}\equiv\pm(2n+1)\pi/4$ are the values that provide the maximal source for
the $\theta$-motion on the right-hand-side of (29). Given that
$n_{B}=R\dot{\theta}$ and that the radial VEV is always positive, a local
Universe filled with matter or antimatter is in fact determined by the initial
VEV of $\theta$, which controls the sign of $\dot{\theta}$.

We can derive the effective potential $V(R,\theta)=m_{R}(\theta)^{2}R^{2}/2$
for the AD field from (28), where
$m_{R}(\theta)^{2}=m_{\sigma}^{2}+(c_{1}-c_{2}\cos(2\theta))\square\phi/\Lambda+c_{3}(\partial\phi)^{2}/\Lambda^{2}$.
For a heavy AD field with $m_{R}\sim\mathcal{O}(H_{\ast})$, the perturbation
of the radial mode on top of the coherent value $R_{0}\equiv\sqrt{\langle
R^{2}\rangle}$ is negligible. In contrast to the presence of “flat directions”
in higher-order non-renormalizable potential Dine:1995kz , $\theta$ may not
have a well-defined coherent value in this scenario, yet one can easily check
that the effective mass of the angular mode at the maximal CP violation is
given by
$m_{\theta}^{2}\sim\partial_{\theta}\partial_{\theta}V(R,\theta)|_{\theta=\theta_{\rm
max}}\rightarrow 0$. The correlation length $x_{c}$ for a massive scalar in de
Sitter, defined from $G(x_{c})=G(0)/2$, is $x_{c}/x_{\rm
ref}=2^{3H_{\ast}^{2}/(2m^{2})}$ Wu:2020ilx ; Starobinsky:1994bd , where $m$
is the scalar mass, $x_{\rm ref}$ is a reference length scale and $G(x)\equiv
G(|\vec{x}_{1}-\vec{x}_{2}|)$ is the two-point spatial correlation function.
Since we are interested in long wavelength modes that have exited the horizon
by the time of USR transition (namely $x_{\rm ref}\sim 1/k_{\rm CMB}$ can be a
good choice), the condition $m_{\theta}\rightarrow 0$ implies that baryon
asymmetry from $\theta$ picked up by the choice (8) for the mass eigenstates
has a correlation length much larger than the Hubble scale, which ensures that
a local patch of the Universe is left with a pure (anti)matter (see also the
discussion in Section 4.1 of Ref. Dine:1995kz ).

##  4 Baryogenesis after inflation

The presence of a coherent oscillation of the inflaton following the end of
the (secondary) slow-roll phase is a conventional assumption to terminate
inflation. Once the inflaton $\phi$ rolls into a deep potential valley, it
acquires a large mass $m_{I}\gg H_{\ast}$ and starts to oscillate around the
potential minimum, leading to an effective matter-domination epoch. Reheating
starts with the decay of inflation into radiative degrees of freedom and
finally the energy density of radiation dominates the Universe. In this
section, we investigate the subsequent evolution of the AD field from inflaton
coherent oscillations to completed reheating (or namely the beginning of
radiation domination).

Figure 3:  The inflaton energy density
$\Omega_{I}\equiv\rho_{I}/(3M_{p}^{2}H^{2})$ as a function of time based on
the smoothed and analytic approximations with an inflaton mass
$m_{I}=20H_{\ast}$. The dotted-dashed line depicts the energy density of
radiation $\Omega_{r}\equiv\rho_{r}/(3M_{p}^{2}H^{2})$. The vertical dashed
line is $t=1/\Gamma_{I}$.

We assume the existence of a deep valley
$\partial^{2}V(\phi)/\partial\phi^{2}\approx m_{I}^{2}\gg H_{\ast}^{2}$
without specifying the shape of the potential $V(\phi)$. The decay of $\phi$
is governed by a perturbative channel, featured by the decay width
$\Gamma_{I}$ of the inflaton, where the energy density of inflaton and
radiation ($\rho_{I}$ and $\rho_{r}$) evolve as

$\displaystyle\dot{\rho}_{I}+3H\rho_{I}=-\Gamma_{I}\rho_{I},\quad\dot{\rho}_{r}+4H\rho_{r}=\Gamma_{I}\rho_{I}.$
(30)

The evolution of each density species satisfies the constraint
$3M_{P}^{2}H^{2}=\rho_{I}+\rho_{r}$. It is convenient to defined the energy
scale of inflation $\Lambda_{I}^{4}=3M_{P}^{2}H_{\ast}^{2}=\rho_{I0}$ as the
initial condition for the post-inflationary dynamics, where
$\rho_{I0}=\rho_{I}(t_{\rm end})$. The energy density of radiation is created
from zero: $\rho_{r0}=0$.

To mimic the background expansion in the effective matter domination Universe
driven by the rapid oscillation of $\phi$, we impose a simplified analytic
representation Wu:2019ohx  as

$\displaystyle\phi(t)=\phi_{\rm max}a^{-3/2}\cos\left[m_{I}(t-t_{\rm
end})\right]e^{-\Gamma_{I}(t-t_{\rm end})/2},$ (31)

where $\phi_{\rm max}\equiv\Lambda_{I}^{2}/m_{I}$ is the maximal oscillatory
amplitude for $\phi$ at the beginning of oscillation. This analytic
representation is a good approximation as long as $m_{I}\gg H_{\ast}$ and
$H_{\ast}(t-t_{\rm end})\gg 1$. For clearness and simplicity we take $t_{\rm
end}=0$ from now on. In the limit of $m_{I}\rightarrow\infty$, one can
reproduce the smoothed energy density of the inflaton as $\rho_{I}=\langle
m_{I}^{2}\phi^{2}\rangle\rightarrow\Lambda_{I}^{4}a^{-3}e^{-\Gamma_{I}t}$. We
plot in Figure 3 an analytic example with $m_{I}=20H_{\ast}$. One can see that
$t_{r}\equiv 1/\Gamma_{I}$ is the approximated time scale at the beginning of
radiation domination for arbitrary choice of $m_{I}$.

Based on the analytic representation (31) we can express the mass terms (21)
of the mass eigenstates $\sigma_{\pm}$ as explicit functions of the background
dynamics through $\square\phi=-\ddot{\phi}-3H\dot{\phi}=m_{I}^{2}\phi$ and
$(\partial\phi)^{2}=\dot{\phi}^{2}$. Therefore we are now ready to study the
time evolution of $\sigma_{\pm}$ in the post-inflationary epochs with initial
conditions at the end of inflation investigated in the previous
sections.333The AD field has a maximal energy density at the end of inflation
where the $c_{1}$ or $c_{2}$ term has
$\Delta\mathcal{L}\lesssim\frac{1}{\Lambda}\frac{H_{\ast}^{4}}{m_{\sigma}^{2}}m_{I}^{2}\phi_{\rm
max}\sim\frac{m_{I}\Lambda_{I}^{2}}{\Lambda m_{\sigma}^{2}}H_{\ast}^{4}$ and
the $c_{3}$ term has
$\Delta\mathcal{L}\lesssim\frac{\Lambda_{I}^{4}}{\Lambda^{2}}\frac{H_{\ast}^{4}}{m_{\sigma}^{2}}$
as $\sigma_{\pm}\sim H_{\ast}^{2}/m_{\sigma}$ is slightly decay during
inflation. For a heavy AD field with $m_{\sigma}\sim H_{\ast}$,
$\Delta\mathcal{L}\ll\mathcal{L}_{\phi}$ holds in an apparent way so that
reheating driven by the coherent oscillation of $\phi$ is not interrupted by
the AD field back reactions.

Dynamical initial conditions. As seen from (17), for the ending time of
inflation $N_{\rm end}\lesssim 1/\Delta_{3\pm}^{-}$, VEVs of the AD field
excited by the USR phase transition remains in motion, providing dynamical
initial conditions for their relaxation during reheating.

Figure 4: The time evolution of mass eigenstates $\sigma_{\pm}$ from the end
of inflation to radiation domination with $m_{I}=20H_{\ast}$, $\Lambda=M_{P}$,
$H_{\ast}=2.37\times 10^{13}$ GeV and $\Gamma_{I}=10^{13}$ GeV. In both panels
$N_{\ast}=1$ and $N_{\rm end}=5$ are used. The vertical dotted lines are
$t=\pi/m_{I}$ and the vertical dashed lines are $t=1/\Gamma_{I}$. Left Panel:
The case of tachyonic initial masses with
$\\{c_{1},c_{2},c_{3}\\}=\\{-2,1,1\\}$ and $m_{\sigma}=H_{\ast}/2$. Right
Panel: The case of positive masses with $\\{c_{1},c_{2},c_{3}\\}=\\{2,-1,1\\}$
and $m_{\sigma}=H_{\ast}$.

Let us focus on a heavy AD field with $m_{\sigma}\sim H_{\ast}$ and
$\Lambda\sim M_{P}$ which guarantees initial conditions of $\sigma_{\pm}$ that
are well computed by the constant-mass approximation (3.1). For
$m_{\sigma}\lesssim H_{\ast}$ and $m_{I}\gg H_{\ast}$, the initial masses at
$t_{\rm end}$ are dominated by $c_{1}$ and $c_{2}$ terms where
$m_{\pm}^{2}\approx(c_{1}\pm c_{2})m_{I}\Lambda_{I}^{2}/\Lambda\sim(c_{1}\pm
c_{2})m_{I}H_{\ast}\gg H_{\ast}^{2}$. This is due to the transition of
potential energy to kinetic terms for the inflaton at the end of inflation,
which gives the opportunity for $m_{\pm}$ to become tachyonic at the beginning
of reheating, depending on the choices of $c_{1}$ and $c_{2}$. In Figure 4, we
show examples for the evolution of $\sigma_{\pm}$ from tachyonic initial
masses $m_{\pm}^{2}(t_{\rm end})<0$ and positive initial masses
$m_{\pm}^{2}(t_{\rm end})>0$ with different choices of the coefficients
$c_{i}$. For both choices of $c_{1}$ and $c_{2}$ in Figure 4, the effective
masses $m_{i\pm}^{2}$ during inflation given by (3.1) are all positively
defined. It is remarkable that these initial VEVs are in the regime of
$\sigma_{\pm}\ll H_{\ast}$, while the conventional AD baryogenesis from flat
directions Dine:1995kz  relies on initial VEVs much greater than the scale of
$H_{\ast}$.

Figure 5: The baryon asymmetry $|Y_{B}|$ generated by a heavy AD field
$m_{\sigma}=H_{\ast}/2$ with tachyonic initial masses led by the choice of
$\\{c_{1},c_{2},c_{3}\\}=\\{-2,1,1\\}$. Upper Panel: Initial conditions of the
mass eigenstates (left) and their time derivatives (right) with respect to
$N_{\rm end}$. Lower Left Panel: The time evolution of $|Y_{B}|$ from the end
of inflation to reheating complete with respect to various choices of
$N_{\ast}$ and $N_{\rm end}=5$. The vertical dotted line is $t=\pi/m_{I}$ and
the vertical dashed line is $t=1/\Gamma_{I}$. Lower Right Panel: The final
$|Y_{B}|$ in radiation domination with respect to various choices of $N_{\rm
end}$. In all panels, $\Lambda=M_{P}$, $H_{\ast}=2.37\times 10^{13}$ GeV and
$\Gamma_{I}=10^{13}$ GeV are used.

With the numerical solutions of the mass eigenstates $\sigma_{\pm}$, we can
obtain the baryon asymmetry $Y_{B}(t)=n_{B}(t)/s(t)$ during reheating where
$n_{B}=\sigma_{+}\dot{\sigma}_{-}-\sigma_{-}\dot{\sigma}_{+}$. For arbitrary
choices of the initial conditions, the ending time of inflation, $N_{\rm
end}$, is the most important parameter that determines the final $Y_{B}$ in
radiation domination. Some examples based on the tachyonic initial masses of
$m_{\pm}$ at the end of inflation are given in Figure 5. The decay of
$|Y_{B}|$ due to the increase of $N_{\rm end}$ in the initial conditions is
easy to understand since the initial VEVs of $\sigma_{\pm}$ and their time
derivative are exponentially diluted by $N_{\rm end}$.

On the other hand, $Y_{B}$ is not sensitive to the duration of the USR phase,
$N_{\ast}$, from arbitrary initial conditions. This is due to the fact that
with sufficiently large $N_{\ast}$, $\epsilon_{\ast}=\epsilon_{\rm
CMB}e^{2\delta_{2}N_{\ast}}$ becomes too small so that the effective masses
$m_{2\pm}=m_{3\pm}\approx m_{\sigma}$, as seen from (3.1). As a result,
initial conditions of the mass eigenstates at the end of inflation are nearly
the same for $N_{\ast}\gtrsim 1$ when $\delta_{2}=-3.05$. Such an asymptotic
constant behavior for the final baryon asymmetry towards the large $N_{\ast}$
limit has an important implication to the PBH formation from the USR inflation
models. As will be discussed in Section 5, $N_{\ast}>2$ is generally required
for the by product PBHs from USR inflation to act as important dark matter
today.

The decay of AD condensate due to scattering with thermalized particles during
reheating could be an important issue, especially for baryogenesis following
high-scale inflation considered in this work. It is usually assumed that the
finite temperature effect induced a thermal mass $m_{T}$ to the AD field led
by the smallest Yukawa coupling $y_{s}$ as $m_{T}\sim y_{s}T$ Anisimov:2000wx
. For the scale of inflation $\Lambda_{I}\sim 10^{16}$ GeV used in our
examples, the maximal temperature during reheating is $T_{\rm
max}\sim(M_{P}\Lambda_{I}^{2}\Gamma_{I})^{1/4}\sim 10^{15}$ GeV. However, one
can check that if $y_{s}$ is of order $10^{-2}-10^{-4}$, $m_{T}$ has
negligible contribution to the relaxation process since
$m_{\pm}\lesssim(m_{I}\Lambda_{I}^{2}/\Lambda)^{1/2}$ is dominated by inflaton
kinetic energy at the beginning of reheating. Therefore, unless $y_{s}$ could
be of order $10^{-1}$ or larger, we find that the finite temperature effect
does not play an important role in the AD mechanism.

Long-term inflation. For models of inflation that experience a long duration
until the end of the secondary slow-roll phase with $N_{\rm end}\gg
1/\Delta_{3\pm}^{-}$, the coherent motion of the AD field triggered by the
Phase 1 to 2 transition becomes negligibly small. In this case the
accumulation of long wavelength perturbations that exit the horizon in Phase 3
can develop new condensates for the mass eigenstates $\sigma_{\pm}$. As the
subsequent scalar condensate reaches equilibrium under the stochastic effect
Starobinsky:1994bd , the VEV of $\sigma_{\pm}$ can be estimated by

$\displaystyle\sigma_{\rm
end\pm}=\sqrt{\frac{3}{8\pi}}\frac{H_{\ast}^{2}}{m_{\rm
end\pm}},\qquad\dot{\sigma}_{\rm end\pm}=0,$ (32)

where $\epsilon_{\rm end}\sim\mathcal{O}(1)$, and for convenience we take
$\epsilon_{\rm end}=1/2$ so that (21) at $N_{\rm end}$ gives

$\displaystyle m_{\rm end\pm}^{2}=m_{3\pm}^{2}(t_{\rm
end})=m_{\sigma}^{2}+\frac{c_{1}\pm
c_{2}}{\Lambda}\delta_{2}M_{P}H_{\ast}^{2}+\frac{c_{3}}{\Lambda^{2}}M_{P}^{2}H_{\ast}^{2}.$
(33)

Note that $m_{\rm end\pm}$ is a constant in constant-rate inflation with
$\delta_{2}\leq-3$.

Baryogenesis from the static VEVs given by (32) does not rely on the constant-
mass approximation (3.1) nor the late-time expression (17), so that the
applicable parameter space for $m_{\sigma}$ and the cutoff scale $\Lambda$ is
less constrained. This allows us to explore cases away from the USR limit with
$\delta_{2}<-3$. In Figure 6 we compare the resulting $Y_{B}$ with different
choices of $\delta_{2}$ away from the USR limit. The final $Y_{B}$ approaches
to a constant value in the massless limit where $m_{\sigma}\ll m_{\rm
end\pm}$. The enhancement of $Y_{B}$ in the large AD mass limit $m_{\sigma}\gg
H_{\ast}$ is due to the resonance amplification led by the choice of inflaton
mass $m_{I}/H_{\ast}\sim\mathcal{O}(10)$. Such a resonance between inflaton
and the AD field could happen in reality but is not a necessary condition for
the scenario to reach enough baryon asymmetry. We remark that the purpose of
this section is to show a complete analysis of the possible initial conditions
for baryogenesis triggered by constant-rate inflation. However, information
related to parameters for PBH formation, such as $N_{\ast}$ or $N_{\rm end}$,
will be washed away by the long-term inflation. For the discussion of the PBH-
baryon correlation in Section 5, we shall focus on dynamic initial conditions.

Figure 6: The baryon asymmetry $|Y_{B}|$ from static initial conditions (32)
after long-term inflation with respect to the mass parameter $m_{\sigma}$ of
the AD field and various choices of $\delta_{2}$. In this plot, we use
$\\{c_{1},c_{2},c_{3}\\}=\\{-2,1,1\\}$, $m_{I}=20H_{\ast}$,
$\Lambda=0.3M_{P}$, where $H_{\ast}=2.37\times 10^{13}$ GeV and
$\Gamma_{I}=10^{13}$ GeV.

##  5 Baryon-PBH correlation

So far we have investigated baryon production during or after multi-stage USR
inflation without taking into account PBH formation is generated by the
enhanced power spectrum $P_{\zeta}$ on scales much smaller than those of the
CMB. As the inflaton decays during reheating, the enhanced small-scale
curvature perturbation $\zeta$ is inherited by the density perturbations of
the radiation, and high peaks of these superhorizon density perturbations
collapse to form PBHs when they reenter the horizon. For those byproduct PBHs
from USR inflation to occupy a significant fraction of the dark matter density
$\Omega_{\rm DM}$ at matter-radiation equality ($t=t_{\rm eq}$), the power
spectrum $P_{\zeta}$ has to satisfy several constraints set by current
observations. In this section, we investigate the viable parameter space for
the multi-stage USR inflation to produce significant PBH dark matter and a
large enough baryon asymmetry.

###  5.1 PBH dark matter

We focus on the viable window for PBHs to be all dark matter in the asteroid-
mass range $M_{\rm PBH}/M_{\odot}\simeq 10^{-16}-10^{-12}$. This indicates
that the pivot scale $k_{0}/k_{\rm eq}\approx$ $(M_{\rm eq}/M_{\rm
PBH})^{1/2}$ for the initiation of the USR phase is subject to the range of
$k_{0}\simeq 10^{12}-10^{14}$ Mpc-1, where $k_{\rm eq}=0.01$ Mpc-1 and $M_{\rm
eq}\simeq 2.94\times 10^{17}M_{\odot}$ are the horizon wavenumber and the
horizon mass at matter-radiation equality, respectively. The comoving horizon
length $R_{H}\equiv 1/(aH)$ as a function of the horizon mass
$M_{H}=4\pi\rho(H)H^{-3}/3$ is given by

$\displaystyle R_{H}(M_{H})=\frac{1}{k_{\rm eq}}\left(\frac{M_{H}}{M_{\rm
eq}}\right)^{1/2}\left(\frac{g_{\ast}}{g_{\rm eq}}\right)^{1/6},$ (34)

where $g_{\rm eq}=3$ accounts for the relativistic degrees of freedom at
$t_{\rm eq}$. Note that $R_{\rm eq}=R_{H}(M_{\rm eq})\approx 1.57\times
10^{40}$ GeV-1.

We define $\Omega_{\rm PBH}(R_{H})=\Omega_{\rm PBH}(M_{H})$ as the density
parameter of PBHs at a given scale $R_{H}$ (or namely $M_{H}$) in radiation
domination. The mass fraction, $\beta(M_{\rm PBH},M_{H})\equiv d\Omega_{\rm
PBH}/d\ln M_{\rm PBH}$, describes the distribution of PBH masses $M_{\rm PBH}$
at the given scale $M_{H}$. The statistics of the PBH density at the comoving
scale $R_{H}$ can be summarized as

$\displaystyle\Omega_{\rm
PBH}(R_{H})=\int\cdots\int_{F_{c}}^{\infty}\frac{M_{\rm
PBH}}{M_{H}}f_{c}(y_{i})P(F,y_{i},\sigma_{Fi})dFdy_{1}\cdots dy_{i},$ (35)

where a Gaussian random field is an ideal choice for $F$ and $y_{i}$ denotes
spatial or temporal derivatives (constraints) of $F$. $f_{c}(y_{i})$ describes
spacetime constraints to ensure that the selected $F$ meets the criteria of
PBH formation Bardeen:1985tr ; Young:2014ana ; Green:2004wb ; Suyama:2019npc ;
Wu:2020ilx . $\sigma_{Fi}$ describes the $i$-th spectral moment of $F$ which
can be computed by the power spectrum of $F$ as

$\sigma_{Fi}^{2}(R_{H})=\int_{0}^{\infty}k^{2i}W^{2}(kR_{H})P_{F}(k)d\ln k,$
(36)

where $W(kR_{H})$ is a window function smoothing over $R_{H}$ to prevent the
divergence in the limit of $k\rightarrow\infty$. The probability distribution
function $P(F,y_{i},\sigma_{Fi})$ is in general a multi-variable structure due
to the non-trivial conditional distribution of $F$. The density contrast of
radiation is the most widely adopted choice of $F$ for PBH formation
Carr:1975qj . The curvature perturbation $\zeta$ is another natural option for
the fundamental random field in the statistics of PBH abundance from models of
inflation Yoo:2018esr ; DeLuca:2019qsy ; Kalaja:2019uju . The compaction
function is a good choice of $F$ to explore the $P(F,y_{i},\sigma_{Fi})$
correlation with the threshold value $F_{c}$ Germani:2019zez ; Young:2020xmk ,
and the inflaton perturbation, $\delta\phi$, is a convenient candidate of $F$
to study the effect of quantum diffusion in ultra-slow-roll scenarios
Biagetti:2021eep .

####  5.1.1 The fiducial statistics

The simplest one-variable Gaussian statistics, referred to as the Press-
Schechter method Carr:1975qj , is the standard approach to obtain a reference
PBH abundance from models of inflation. As the fiducial method to forecast the
PBH abundance from the constant-rate model (6), we adopt the Press-Schechter
formula for the mass fraction of the form:

$\displaystyle\beta_{\rm
PS}(M_{H})=2\int_{\nu_{c}}\frac{e^{-\nu^{2}/2}}{\sqrt{2\pi}}d\nu=\textrm{erfc}\left(\nu_{c}/\sqrt{2}\right),$
(37)

where $F$ is identified as the density contrast of radiation, $\nu\equiv
F/\sigma_{F0}$ is the peak value for a Gaussian random field $F$ over its
variance $\sigma_{F0}$ (or the zeroth spectral moment), and $\nu_{c}\equiv
F_{c}/\sigma_{F0}$ is the critical peak value defined at the threshold $F_{c}$
of gravitational collapse. The analytic expression of (37) is given in
Appendix C.1.

Given that PBHs behave as dust-like matter, the mass fraction grows as
$\beta(M_{H})=\beta(M_{\rm eq})(a/a_{\rm eq})$ in the radiation-dominated
Universe. The final PBH density at matter-radiation equality is therefore
given by

$\Omega_{\rm PBH,eq}=\int\beta(M_{H})\left(\frac{M_{\rm
eq}}{M_{H}}\right)^{1/2}d\ln M_{H}\equiv\Omega_{\rm DM,eq}\int f(M_{H})d\ln
M_{H}$ (38)

where the scaling relation $M_{H}/M_{\rm eq}\sim a/a_{\rm eq}$ is used. The
mass function $f(M_{H})=\beta(M_{H})$ $(M_{\rm eq}/M_{H})^{1/2}/\Omega_{\rm
DM,eq}$ shows the PBH ratio in dark matter density at each horizon scale
$M_{H}$ and $f_{\rm PBH}=\Omega_{\rm PBH,eq}/\Omega_{\rm DM,eq}$ measures the
contribution from all PBHs in the dark matter density.

Let us explore the viable parameter space for realizing $f_{\rm PBH}>0.1$,
where PBHs occupy a significant fraction of the dark matter density. A perfect
USR inflation with $\delta_{2}=-3$ will launch a scale-invariant power
spectrum from $k_{0}$ all the way to $k_{\rm end}$. This means that the
dominant contribution to the mass function $f(M_{H})$ at matter-radiation
equality comes from the mass range around $M_{H}\gtrsim 10^{15}$ gram, as
lower mass PBHs have been fully evaporated via Hawking radiation. The perfect
USR case with $f_{\rm PBH}>0.1$ is thus immediately ruled out by the severe
constraints on PBH evaporation in the CMB, and the extra and inner galactic
background radiation.

Figure 7:  The PBH mass function $f(M_{\rm PBH})$ (blue) based on the fiducial
(Press-Schechter) approach (C.1, solid line) and extended mass functions given
by the effect of critical collapse (C.2, dashed), the finite-size effect (C.3,
dash-dotted) and non-linear density contrast (C.4, dotted). The pivot scale
$k_{0}=7.3\times 10^{12}$ Mpc-1 for the USR transition is used, and $N_{\ast}$
in each function is fixed by $f_{\rm PBH}=1$ with $\delta_{2}=-3.05$ and
$N_{\rm end}=20$. The existing observational bounds are shown from: Galactic
(red) and extragalactic (green) probes, as well as CMB and Edges (black).

Fortunately, a slightly red-tilted USR inflation with $\delta_{2}\lesssim-3$
can lead to sufficient decay of $f(M_{\rm PBH})$ towards the small mass limit.
In Figure 7 we provide an example of the mass function computed by the
fiducial method (37) (solid blue line) with $\delta_{2}=-3.05$ and $f_{\rm
PBH}=1$. For a comparison, we also provide extended mass functions based on
several statistical uncertainties discussed in Appendix C. In all cases, one
can see that $f(M_{\rm PBH})<10^{-6}$ for $M_{\rm PBH}/M_{\odot}\leq 10^{-18}$
so that the mass functions given in Figure 7 are compatible with current
constraints444The plot of the observational constraints was generated using
the publicly available code PlotBounds.. Notably, the green lines denote the
Galactic bounds from: the Voyager1 data Boudaud:2018hqb , the MeV diffuse flux
observed by the Integral/Spi detector Laha:2020ivk , the 511 keV line observed
in our Galaxy Laha:2019ssq ; DeRocco:2019fjq  and the Comptel measurements
Coogan:2020tuf . The red lines indicate the extragalactic bounds from Super-
Kamiokande Dasgupta:2019cae  and from the extragalactic background radiation
Carr:2009jm . Finally, the CMB and Edges constraints are displayed in black.

There is a lower bound of $\delta_{2}$ that comes from the constant-mass
approximation (3.1) for the AD field, given that time-dependence appears in
the Phase 3 mass eigenstates $m_{3\pm}$ when $\delta_{2}$ is away from the USR
limit at $-3$. For example, with $m_{\sigma}/H_{\ast}=1/2$, $N_{\ast}=2$ and
$N_{\rm end}=20$, the constant-mass approximation $m_{3\pm}\approx m_{\sigma}$
holds for $\Lambda>0.2M_{P}$ at which the lower bound reads
$\delta_{2}>-3.12$. For $\Lambda=M_{P}$, the constant-mass holds for a larger
parameter space with $\delta_{2}>-3.3$. We remark that, however, the violation
of constant-mass approximation does not imply a failure of the baryogenesis
but just invokes an extension of the formalism in Section 2.2.

####  5.1.2 Uncertainties

The PBH abundance from inflationary models of the USR class involves several
uncertainties in the statistical approaches. The purpose of this section is to
clarify the possible deviations from the fiducial parameter space led by
various uncertain effects for the PBH formation. Here we summarize these
uncertainties in terms of the difference in $N_{\ast}$ with respect to the
fiducial value reported by $\beta_{\rm PS}$ defined in (37), where $N_{\ast}$
is obtained by $f_{\rm PBH}(\delta_{2},N_{\ast},N_{\rm end})=1$ with fixed
$\delta_{2}$ and $N_{\rm end}$.

  * •

Quantum diffusion. The leading uncertainties for the PBH abundance in this
scenario. In the exact USR limit ($\delta_{2}=-3$), investigations in
Biagetti:2021eep ; Figueroa:2020jkf  suggest that the real abundance enhanced
by quantum diffusion might lead to a deviation $\beta/\beta_{\rm PS}\sim
10^{5}-10^{12}$, which can be translated into an
$\mathcal{O}(10^{-2}-10^{-1})$ deviation in $N_{\ast}$. However,
generalization of the method to the quasi-USR case with $\delta_{2}<-3$ is
required for the current scenario, where we expect a significant suppression
of the effect of quantum diffusion due to the large effective inflaton mass of
order $H_{\ast}$ (see Appendix A). As discussed in Section 5.2, the viable
parameter space for PBH as important dark matter with large enough baryon
asymmetry allows a deviation $\beta/\beta_{\rm PS}\gg 10^{100}$ or namely an
$\mathcal{O}(1)$ shift in $N_{\ast}$.

  * •

Non-linear density. The leading uncertainty for PBH formation arises from
models of inflation. The deviation could be up to $\mathcal{O}(10^{-1})$
difference in $N_{\ast}$ (Appendix C.4). The non-linear effect moves
$N_{\ast}$ towards a larger value since the required spectral amplitude
$A_{\rm PBH}$ becomes larger Young:2019yug . The shape of the resulting mass
function is similar to that of the critical collapse, see Figure 7.

  * •

Finite-size effect. PBHs are not zero-size objects in real space and they only
form at the local maximum (peaks) of the density contrast Bardeen:1985tr .
These spatial constraints introduce a $\nu_{c}^{3}$ factor in addition to the
Press-Schechter method Wu:2020ilx ; Young:2014ana . As shown in Appendix C.3,
the effect gives a $\mathcal{O}(10^{-2})$ correction towards a smaller
$N_{\ast}$ for a fixed PBH abundance.

  * •

Critical collapse. The effect of critical collapse Niemeyer:1997mt ;
Yokoyama:1998xd ; Musco:2012au ; Musco:2008hv  is an intrinsic uncertainty for
PBH formation from the reenter of large density perturbation into horizon. The
extended mass function under critical collapse (Appendix C.2) shows an
evidently increased low-mass tail and a lowered peak value, yet it has a
negligible correction of $\mathcal{O}(10^{-3})$ to the parameter $N_{\ast}$
for a fixed PBH abundance.

Note that corrections from the transfer function that describes the time
evolution of the density contrast after horizon reentry have been neglected in
our results. There may be more uncertainties that arise from the choices of
window functions, the profile dependence of the curvature perturbation, and
the exact threshold value of gravitational collapse. We expect the combination
of all uncertainties can only lead to a difference $\beta/\beta_{\rm PS}\ll
10^{100}$, which means at most an $\mathcal{O}(10^{-1})$ uncertainty for
$N_{\ast}$.

###  5.2 The cosmic coincidence

The cold dark matter density today $\Omega_{\rm CDM0}=0.265$ in the standard
$\Lambda$CDM Universe Aghanim:2018eyx  indicates that at matter-radiation
equality $\Omega_{\rm CDMeq}=0.42$ and the baryon density $\Omega_{\rm
Beq}=m_{B}n_{\rm Beq}=0.08$, where $m_{B}=0.938$ GeV is the averaged nucleon
mass and $n_{\rm Beq}=|Y_{B}|s(t_{\rm eq})$ is the baryon number density. The
specific ratio between the two density parameters:

$\displaystyle\Omega_{\rm CDMeq}/\Omega_{\rm Beq}\approx 5,$ (39)

is referred to as the “cosmic coincidence problem.” This cosmic coincidence
also implies that $\Omega_{\rm PBHeq}/\Omega_{\rm Beq}\sim\mathcal{O}(1)$ if
PBHs occupy $10$ $-$ $100$ percent of the dark matter density.

Taking $H_{0}=67.36$ km s-1 Mpc-1, $\Omega_{\Lambda 0}=1-\Omega_{m0}=0.6847$,
and the redshift $z_{\rm eq}=3402$ from Aghanim:2018eyx , we find that $H_{\rm
eq}=H_{0}\sqrt{\Omega_{\Lambda 0}+2\Omega_{m0}(a_{0}/a_{\rm eq})^{3}}\approx
2.264\times 10^{-37}$ GeV and thus the baryon asymmetry at equality is
$|Y_{B}|=6.25\times 10^{-11}$. Based on this estimate, we seek a viable
parameter space for baryogenesis to satisfy $|Y_{B}|\gtrsim 10^{-10}$.

Figure 8:  The parameter scan of PBH-baryon production from USR inflation with
$\delta_{2}=-3.05$ in the $\\{N_{\ast},N_{\rm end}\\}$ plane, where $N=0$ is
chosen at the onset of USR transition in Figure 1. The PBH-to-dark matter
ratio $f_{\rm PBH}(N_{\ast},N_{\rm end})=0.1$, $0.5$, $1$ is estimated by the
monochromatic relation $M_{H}=M_{\rm PBH}$. Initial conditions for
$\sigma_{\pm}(N_{\ast},N_{\rm end})$ are computed by (17) based on the
constant-mass approximation with $\\{c_{1},c_{2},c_{3}\\}=\\{-2,1,1\\}$,
$m_{\sigma}=H_{\ast}/2$, and $\Lambda=0.3M_{P}$. The contours are the final
baryon asymmetry $|Y_{B}|$ resulting from the given initial conditions
$\\{N_{\ast},N_{\rm end}\\}$, where $m_{I}=20H_{\ast}$ and
$\Gamma_{I}=10^{13}$ GeV are taken in the reheating scenario. Colored regions
in the left panel are excluded by the minimal $e$-fold number of single-field
inflation and the condition $m_{\pm}/H_{\ast}<3/2$ for the late-time
approximation of (17). Right Panel: A zoom-in for the ratio $f_{\rm
PBH}(N_{\ast},N_{\rm end})$ on the $N_{\ast}$ axis.

In Figure 8, we scan the parameter space $\\{N_{\ast},N_{\rm end}\\}$ for the
USR inflation with a slightly tilted spectrum $\delta_{2}=-3.05$ that can be
realized $f_{\rm PBH}>0.1$ and keep the constant-mass formalism (3.1) as a
good approximation. The results in Figure 8 show that baryogenesis triggered
by the USR inflation with $|Y_{B}|\gtrsim 10^{-10}$ is compatible with PBH
dark matter in the range of $19<N_{\rm end}<24$ and this viable parameter
space is not sensitive to a small variation in $N_{\ast}$. On the other hand,
the PBH abundance is very sensitive to a perturbation in $N_{\ast}$ so that a
specific value in the range of $0.1\leq f_{\rm PBH}\leq 1$ will fix a precise
value of $N_{\ast}$ in the region of $N_{\ast}>2$. In other words, if PBHs
compose an important dark matter fraction with $0.1\leq f_{\rm PBH}\leq 1$,
the inflationary parameter $N_{\ast}$ is fixed to the space that gives
$|Y_{B}|\gtrsim 10^{-10}$ so that $\Omega_{\rm PBHeq}/\Omega_{\rm
Beq}\sim\mathcal{O}(1)$ becomes a natural outcome in the scenario. This
provides a potential explanation to the cosmic coincidence problem for PBH
dark matter.

The above conclusion is unchanged even if one takes into account the effects
of critical collapse (C.2) or the non-linear density corrections (C.4) in the
PBH mass function. As one can see from Figure 7, the effect of critical
collapse slightly decreases the peak amplitude of $f(M_{\rm PBH})$ and lifts
the population of PBHs towards the small mass limit. As a result, the values
of $N_{\ast}$ required by $0.1\leq f_{\rm PBH}\leq 1$ are slightly larger than
those without critical collapse. Similarly, the non-linearity of the density
contrast asks a larger $N_{\ast}$ by $\mathcal{O}(0.1)$, yet a constant shift
of $N_{\ast}$ towards the region of $N_{\ast}>2$ does not affect the resulting
$|Y_{B}|$ in this scenario.

Note that in the paradigm of slow-roll single-field inflation, the exclusion
of a perfect scale-invariant $P_{\zeta}$ in general requires a finite $e$-fold
number $50<N_{\rm end}-N_{0.002}<60$, where $N_{0.002}$ is the number at the
pivot scale of CMB measurements $k_{\rm CMB}=0.002$ Mpc-1 Akrami:2018odb . The
non-observation of spatial curvature requires an additional $4.9$ $e$-folds to
produce the pivot scale Akrami:2018odb . These two constraints imply a minimal
$e$-fold number from CMB scales to the smallest scales at the end of
inflation: $\ln(k_{\rm end}/k_{\rm CMB})\approx\ln(a_{\rm end}/a_{\rm
CMB})>54.9$, which translates into a constraint $N_{\rm end}=\ln(k_{\rm
end}/k_{0})>54.9-\ln(k_{0}/k_{\rm CMB})$.

Before closing this section, we remark that the final baryon asymmetry
triggered by the USR transition during inflation always approaches a constant
value in the regime of $N_{\ast}>2$ for a heavy AD field with
$m_{\sigma}\sim\mathcal{O}(H_{\ast})$, which we have numerically confirmed
with various choices of the $\mathcal{O}(1)$ parameters, $c_{1}$, $c_{2}$,
$c_{3}$ (but avoided the fine cancellation between $c_{1}$ and $c_{2}$), and
$m_{I}>10H_{\ast}$. This conclusion relies on the constant-mass approximation
(3.1) which sets a lower bound to the cutoff as $\Lambda>\sqrt{\epsilon_{\rm
CMB}}M_{P}H_{\ast}/m_{\sigma}$. However, a $\Lambda$ below this bound does not
necessarily mean a failure of the baryogenesis, but only requires a
modification of the solutions in Section 2.2 to that of a time-dependent mass
term. In general, a larger $Y_{B}$ is obtained with a lower $\Lambda$ since it
gives a larger source of the baryon number violation in the equation of motion
(namely the $c_{2}$ term).

The “asymptotic constant” behavior of $|Y_{B}|$ implies that the viable
parameter space for explaining the cosmic coincidence problem holds as long as
$A_{\rm PBH}/A_{\rm CMB}>5\times 10^{2}$ for $\delta_{2}\lesssim-3$, which can
be translated into the viable range $10^{-6}<A_{\rm PBH}<10^{-1}$ with the
upper bound set by the perturbativity of the inflation model.555For a
reference in terms of $f_{\rm PBH}$ based on the fiducial method, the value
$N_{\ast}=2.3$ with $\delta_{2}=-3.05$ gives $f_{\rm PBH}\sim 10^{-130}$ and
$A_{\rm PBH}\sim 3\times 10^{-3}$, and $N_{\ast}=2.3$ with $\delta_{2}=-3.1$
gives $f_{\rm PBH}\sim 10^{-83}$ and $A_{\rm PBH}\sim 5\times 10^{-3}$.  As a
result, even if there could be some $\mathcal{O}(1-10)$ uncertainty for the
necessary amplitude $A_{\rm PBH}$ to realize $0.1\leq f_{\rm PBH}\leq 1$ due
to any uncovered effect involved in the PBH formation from USR inflation, the
prediction of the ratio $\Omega_{\rm PBHeq}/\Omega_{\rm Beq}$ remains
unaffected.

##  6 Conclusion

The sharp deceleration of the slow-roll dynamics on small scales into an
“ultra-slow-roll (USR)” phase is a generic mechanism to enhance the primordial
power spectrum for PBH formation in single-field inflation. If PBHs indeed
play an important role in dark matter, the cosmic coincidence problem among
the energy densities of the Universe today might hint a correlated origin for
baryon and dark matter. In this work, we have explored and confirmed the
viability of baryogenesis, based on the AD mechanism, with modified initial
conditions and baryon-number-violating operators triggered by the generic USR
transition for PBH formation.

Baryogenesis driven by the USR transition during inflation significantly
extends the viable parameter space for the initial conditions from preceding
conclusions. It does not rely on the presence of “flat directions” nor a
tachyonic soft mass term induced by non-minimal Kahler couplings. Scalar
fields with initial VEVs much smaller than the Hubble scale of inflation can
create sufficient baryon asymmetry as long as they were dynamically excited
during inflation. In other words, the AD mechanism of baryogenesis can be
realized with a minimal Kahlar potential or without the assumption of
supersymmetry. Albeit the formalism presented in this work has been focused on
the constant-mass assumption for the AD field near the USR limit
($\delta_{2}\simeq-3$) with a bound $0<m_{\pm}/H_{\ast}<3/2$ on the effective
masses, we expect successful baryogenesis beyond those constraints. The
coherent production of scalar motion away from the USR limit and the modified
late-time behavior for $m_{\pm}/H_{\ast}>3/2$ are interesting topics for
further investigations.

The asymptotic constant behavior of the final baryon asymmetry towards the
large USR duration limit ($N_{\ast}\gg 1$) is a generic feature of the
presented scenario and this fact has important implications to the cosmic
coincidence problem. For PBHs to occupy a significant fraction of the dark
matter density $\Omega_{\rm PBH}$, the value of $N_{\ast}$ must be precisely
fixed, whereas the exact value depends on various uncertainties involved in
the statistics of the PBH abundance. However, all statistical methods tested
in this work suggest that the allowed space of $N_{\ast}$ for PBH dark matter
lies deep inside the constant plateau of the correct baryon asymmetry, where
the coincidence between $\Omega_{\rm PBHeq}$ and $\Omega_{\rm Beq}$ can
incorporate at least $10^{100}$ orders of uncertainties in the ratio of
$f_{\rm PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM}$. Nevertheless, despite the
fact that the USR transition provides a quantitative connection between baryon
and PBH densities with a viable answer to the cosmic coincidence problem
independent of the statistical uncertainties, the existing fine-tuning problem
for the inflaton potential to realize PBH dark matter inevitably enters the
story of cosmology. On the other hand, the asteroid-mass window for PBH dark
matter, one of the fundamental assumption in this work, could be tested by
near future astrophysical or gravitational-wave experiments Nemiroff:1995ak ;
Katz:2018zrn ; Jung:2019fcs ; Ballesteros:2019exr ; Coogan:2020tuf ;
Capela:2013yf ; Graham:2015apa ; Ray:2021mxu ; Dutta:2020lqc ; Carr:2020gox ;
Green:2020jor ; Ali-Haimoud:2019khd ; Bai:2018bej .

###### Acknowledgements.

We thank Julien Lesgourgues and Vincent Vennin for helpful discussions. Y.-P.
Wu and K. Petraki were supported by the Agence Nationale de la Recherche (ANR)
Accueil de Chercheurs de Haut Niveau (ACHN) 2015 grant (“TheIntricateDark”
project). K. Petraki was also supported by the NWO Vidi grant “Self-
interacting asymmetric dark matter.” E. Pinetti is supported by: the Fermi
Research Alliance, LLC under Contract No. DE-AC02-07CH11359 with the U.S.
Department of Energy, Office of High Energy Physics; Department of Excellent
grant 2018-2022, awarded by the Italian Ministry of Education, University and
Research (MIUR); Research grant of the Italo-French University, under Bando
Vinci 2020.

##  Appendix A The inflaton mass in quasi-USR inflation

The exact USR limit with $\delta_{2}=-3$ corresponds to a region of completely
flat potential in which the inflaton field $\phi$ is effectively massless. In
this region stochastic effect driven by modes well inside the horizon can play
an important role in the VEV of $\phi$, leading to non-linear translation of
the curvature perturbation $\zeta$ from $\phi$ beyond the standard Gaussian
approximation Biagetti:2018pjj ; Ezquiaga:2019ftu ; Pattison:2021oen ;
Figueroa:2020jkf ; Biagetti:2021eep . In this section, we drive an upper bound
for $\delta_{2}$ to prevent the breakdown of the Gaussian assumption for
$\zeta$.

The constant-rate condition $\dot{\delta}_{2}=0$ used for the USR template in
Section 2.1 is in fact a constraint on the inflaton mass. To see this, one
uses $\delta=\ddot{\phi}/(H\dot{\phi})$ to replace the classical motion
$\ddot{\phi}+3H\dot{\phi}+V_{\phi}=0$ to obtain
$\delta+3=-V_{\phi}/(H\dot{\phi})$. Taking the time derivative of $\delta$ one
obtains the relation with the inflaton mass $V_{\phi\phi}$ Ng:2021hll , where
the constant-rate condition $\dot{\delta}=0$ gives

$\displaystyle\frac{V_{\phi\phi}}{H^{2}}=\left(\epsilon_{H}-\delta\right)\left(\delta+3\right),$
(40)

where $\epsilon_{H}\rightarrow 0$ manifests a de Sitter background. One can
see that $\delta=0$ (slow-roll) and $\delta=-3$ (USR) are two special cases
for the massless inflaton. The condition $\delta_{3}=-\delta_{2}-3$ with
$\delta_{2}\leq-3$ implies a continuous inflaton mass from Phase 2 to 3 so
that the scaling power of $P_{\zeta}$ is also continuous Ng:2021hll .

To avoid a large quantum diffusion led by the stochastic fluctuation of
$\phi$, we ask the effective mass $m_{\phi}^{2}\equiv|V_{\phi\phi}|$ to
satisfy the condition, $m_{\phi}^{2}/H>H/(2\pi)$, so that the classical
deviation is larger than the size of the quantum fluctuation. With
$V_{\phi\phi}=-\delta(\delta+3)H^{2}$, one finds that
$m_{\phi}/H>1/\sqrt{2\pi}$ gives the upper bound $\delta_{2}<-3.05$.

The statistics of PBH abundance in the exact USR inflation with
$\delta_{2}=-3$ is very sensitive to the effect of quantum diffusion
Ezquiaga:2019ftu ; Pattison:2021oen . For more physical cases with
$\delta_{2}<-3$, the inflaton mass (40) led by the constant-rate condition
$\dot{\delta}_{2}=0$ in fact has fixed solutions of the inflaton mode
functions (with the de Sitter condition $\epsilon_{H}\rightarrow 0$), see
Ng:2021hll . In other words, power spectrum of the Gaussian inflaton
perturbation $\delta\phi$ has been uniquely determined by the background
evolution of $\phi$, where in the USR limit with $\delta_{2}\rightarrow-3$ the
spectrum $P_{\delta\phi}$ is exactly scale invariant. There is no degree of
freedom to impose an additional delta-like spectrum for $\delta\phi$ as in
Biagetti:2021eep  and the effect of quantum diffusion in the probability
distribution function of the density contrast away from the USR limit is left
for future effort.

##  Appendix B Baryogenesis from massless AD field

Let us explore in this section the baryogenesis during inflation triggered by
a charged scalar with a negligible bare mass. In the limit of
$m_{\sigma}\rightarrow 0$, $m_{1\pm}$ and $m_{3\pm}$ given in (3.1) are still
constants yet the scalar mass in Phase 2 only given by the $c_{3}$ term
becomes decaying with time, where

$\displaystyle m_{1\pm}^{2}=$ $\displaystyle-3\frac{c_{1}\pm
c_{2}}{\Lambda}\sqrt{2\epsilon_{\rm CMB}}M_{P}H_{\ast}^{2}$ $\displaystyle
m_{2\pm}^{2}=$ $\displaystyle\frac{c_{3}}{\Lambda^{2}}2\epsilon_{\rm
CMB}M_{P}^{2}H_{\ast}^{2}e^{2\delta_{2}N}=\bar{m}^{2}e^{-6N},$ (41)
$\displaystyle m_{3\pm}^{2}=$ $\displaystyle-3\frac{c_{1}\pm
c_{2}}{\Lambda}\sqrt{2\epsilon_{\ast}}M_{P}H_{\ast}^{2}+\frac{c_{3}}{\Lambda^{2}}2\epsilon_{\ast}M_{P}^{2}H_{\ast}^{2},$
(42)

and we denote $\bar{m}^{2}=2c_{3}\epsilon_{\rm
CMB}M_{P}^{2}H_{\ast}^{2}/\Lambda^{2}$ for convenience. To study the
baryogenesis in this limit, the coefficients $c_{i}$ are chosen such that
$m_{i\pm}$ in each phase is always positive defined. We ask
$m_{1\pm}\sim\mathcal{O}(H_{\ast})$ for the AD field to be compatible with the
isocurvature constraint in Phase 1, which gives an upper bound to the cutoff
as $\Lambda<3\sqrt{2\epsilon_{\rm CMB}}M_{P}$.

For $m_{\sigma}\ll H_{\ast}$, we adopt the massless approximation (B) to the
mass eigenstates, where $m_{2\pm}$ are degenerated but they have explicit
time-dependence. The Phase 2 equation of motion for $\sigma_{\pm}$ becomes

$\displaystyle\ddot{\sigma}_{2\pm}+3H_{\ast}\dot{\sigma}_{2\pm}+\bar{m}^{2}e^{2\delta_{2}H_{\ast}t}\,\sigma_{2\pm}=0,\qquad
t_{0}\leq t<t_{\ast},$ (43)

where $\bar{m}^{2}=2c_{3}\epsilon_{\rm CMB}M_{P}^{2}H_{\ast}^{2}/\Lambda^{2}$,
and in the USR limit $\delta_{2}=-3$ the solution takes a simple form of

$\displaystyle\sigma_{2\pm}(N)=\sqrt{\frac{3}{8\pi^{2}}}\frac{H_{\ast}^{2}}{m_{1\pm}}\cos\left[\frac{\bar{m}}{3H_{\ast}}\left(1-e^{-3N}\right)\right].$
(44)

Note that initial conditions given by (8) at $N=0$ have been used. The
decaying-mass solutions in Phase 2 modify the coefficients of the constant-
mass Phase 3 solutions
$\sigma_{3\pm}=C_{\pm}e^{-\Delta_{3\pm}^{+}N}+D_{\pm}e^{-\Delta_{3\pm}^{-}N}$,
by matching boundary conditions at $N=N_{\ast}$, as

$\displaystyle
C_{\pm}=\frac{e^{\Delta_{3\pm}^{+}N_{\ast}}}{\Delta_{3\pm}^{-}-\Delta_{3\pm}^{+}}\left(\sigma_{\ast\pm}\Delta_{3\pm}^{-}-A_{\pm}\right),\quad
D_{\pm}=\frac{e^{\Delta_{3\pm}^{-}N_{\ast}}}{\Delta_{3\pm}^{-}-\Delta_{3\pm}^{+}}\left(A_{\pm}-\sigma_{\ast\pm}\Delta_{3\pm}^{+}\right),$
(45)

where $\sigma_{\ast\pm}=\sigma_{2\pm}(N_{\ast})$,
$A_{\pm}=\sigma_{0\pm}(\bar{m}/H_{\ast})e^{-3N_{\ast}}\sin[\bar{m}(1-e^{-3N_{\ast}})/(3H_{\ast})]$
and
$\sigma_{0\pm}\equiv\sigma_{\pm}(t_{0})=\sqrt{3/(8\pi^{2})}H_{\ast}^{2}/m_{1\pm}$.

Figure 9: The parameter scan for the baryon asymmetry $|Y_{B}|$ at the end of
inflation from massless AD field ($m_{\sigma}=0$) in the USR limit
$\delta_{2}=-3$ with $\\{c_{1},c_{2},c_{3}\\}=\\{-2,1,1\\}$, $N_{\ast}=0.5$
and $H_{\ast}=2.37\times 10^{13}$ GeV. Dashed lines are contours of $|Y_{B}|$.
Colored regions are excluded by validity of the late-time approximation for
$m_{3\pm}/H_{\ast}<3/2$ and the existence of the USR phase ($N_{\rm
end}>N_{\ast}$).

The resulting baryon asymmetry estimated at the end of inflation from a
massless AD field ($m_{\sigma}\rightarrow 0$) can marginally reach the
observational value $|Y_{B}|\lesssim 10^{-10}$ in a very limited parameter
space. An example using the late-time approximation (26) for
$m_{\pm}/H_{\ast}<3/2$ with $N_{\ast}=0.5$ and $H_{\ast}=2.37\times 10^{13}$
GeV is given in Figure 9. By definition, the ending time of inflation $N_{\rm
end}>N_{\ast}$ is necessary for the existence of the USR phase (Phase 2).

##  Appendix C The statistics of PBH abundance

In this section we provide a more detailed formulation of the uncertain
effects in the statistics of PBH abundance considered in Section 5.1.2. In
Table 1, a summary is given about the condition $f_{\rm PBH}=1$ for PBH as all
dark matter from different statistical approaches in terms of the duration of
the USR phase $N_{\ast}$ and the maximal spectral amplitude $A_{\rm
PBH}=A_{\rm CMB}e^{-N_{\ast}(6+4\delta_{2})}$. The corresponding mass function
in each approach given in Figure 7 is compatible with existing observational
constraints from galactic Boudaud:2018hqb ; DeRocco:2019fjq ; Laha:2019ssq ;
Laha:2020ivk , extra-galactic Carr:2009jm ; Dasgupta:2019cae , and CMB or 21cm
experiments Poulin:2016anj ; Clark:2016nst ; Clark:2018ghm ; Hektor:2018qqw ;
Mittal:2021egv .

|  $\beta_{\rm PS}$ |  $\beta_{c}$ |  $\beta_{\rm 3d}$ |  $\beta_{\rm NL}$
---|---|---|---|---
$\delta_{2}=-3.05$ |  2.6645 |  2.6687 |  2.6423 |  2.7473
$A_{\rm PBH}$ |  0.0329 |  0.0286 |  0.0337 |  0.0549
$\delta_{2}=-3.1$ |  2.5837 |  2.5877 |  2.5628 |  2.6638
$A_{\rm PBH}$ |  0.0334 |  0.0292 |  0.0342 |  0.0558
Table 1: The duration of the USR phase, $N_{\ast}$, and the corresponding
spectral amplitude, $A_{\rm PBH}$, for realizing PBH as all dark matter
$f_{\rm PBH}=1$ from various statistical approaches with $N_{\rm end}=20$.

###  C.1 The fiducial approach

In this work, we forecast the PBH abundance by assuming the monochromatic
relation $M_{\rm PBH}=M_{H}$ (no critical collapse)666Due to the effect of
critical collapse Niemeyer:1997mt ; Yokoyama:1998xd ; Musco:2012au ;
Musco:2008hv , $M_{\rm PBH}$ in general can be much smaller than the horizon
mass $M_{H}$ at the epoch of formation. at each comoving epoch so that the
ratio $M_{\rm PBH}/M_{H}=1$ in (35). Assuming that $F$ is the only variable
relevant to the PBH formation without further constraints, we have
$f_{c}(y_{i})=1$ and
$P(F,y_{i},\sigma_{Fi})=e^{F^{2}/(2\sigma_{F0}^{2})}/\sqrt{2\pi}$ so that (35)
simply gives the Press-Schechter formula (37) (up to a factor of 2).

Conventionally, the Gaussian field $F$ is identified as the density
perturbation of radiation converted from the decay of inflaton during
reheating. At linear order, the power spectrum of $F$ is related to
$P_{\zeta}$ according to

$\displaystyle
P_{F}=\frac{4(1+w)^{2}}{(5+3w)^{2}}\left(\frac{k}{aH}\right)^{4}P_{\zeta},$
(46)

where $w=1/3$ is the equation of state of the Universe in radiation
domination. The linear simplification (46) of the density spectrum inevitably
loses non-linear and non-Gaussian contributions from $P_{\zeta}$
DeLuca:2019qsy ; Kalaja:2019uju ; Kawasaki:2019mbl ; Young:2019yug ;
Saito:2008em ; Byrnes:2012yx ; Tada:2015noa ; Franciolini:2018vbk ;
Atal:2018neu ; Taoso:2021uvl ; Germani:2019zez ; Young:2015cyn . However,
given that the PBH abundance is exponentially sensitive to the peak value
$\nu$, the statistical results based on non-linear or non-Gaussian approaches
can usually be mimicked by the linear Gaussian statistics with a mild change
in the input parameters (or, to be more specifically, the spectral amplitude
$A_{\rm PBH}$), see also DeLuca:2019qsy ; Taoso:2021uvl .

The gain from using the linear relation (46) is that one can obtain explicit
analytic expressions for the spectral moments (36) given by the USR template
(6). Importing the Gaussian window function $W(kR)=\exp[-k^{2}R^{2}/2]$, one
can decompose the template (6) into one blue-tilted trapezoidal spectrum (from
$k_{\rm min}$ to $k_{0}$) plus one step/trapezoidal spectrum (from $k_{0}$ to
$k_{\rm end}$) to obtain the analytic result of $\sigma_{Fi}$ based on the
method in Wu:2020ilx , and it reads

$\displaystyle\sigma_{Fi}^{2}(R_{H})$
$\displaystyle=\frac{4(1+w)^{2}}{(5+3w)^{2}}\frac{A_{\rm PBH}}{2R_{H}^{2i}}$
(47) $\displaystyle\times$
$\displaystyle\left[\frac{\Gamma(5+i+\delta_{2},r_{0}^{2})}{r_{0}^{6+2\delta}}-\frac{\Gamma(5+i+\delta_{2},r_{\rm
end}^{2})}{r_{0}^{6+2\delta}}+\frac{\Gamma(4+i,r_{\rm
min}^{2})}{r_{0}^{4}}-\frac{\Gamma(4+i,r_{0}^{2})}{r_{0}^{4}}\right],$

where $r_{0}=k_{0}R_{H}$, $r_{\rm min}=k_{\rm min}R_{H}$ and $r_{\rm
end}=k_{\rm end}R_{H}$. Note that $A_{\rm PBH}=A_{\rm
CMB}e^{-N_{\ast}(6+4\delta_{2})}$, $k_{\rm min}=(A_{\rm CMB}/A_{\rm
PBH})^{1/4}k_{0}$, and $k_{\rm end}=k_{0}e^{N_{\rm end}}$ are all expressed by
free parameters $\\{\delta_{2},N_{\ast},N_{\rm end}\\}$ of the USR inflation.
$k_{0}$ is fixed by the targeting mass range of PBHs.

###  C.2 The effect of critical collapse

In this section, we derive the PBH mass function with extended $M_{\rm PBH}$
to $M_{H}$ relation led by the effect of critical collapse Niemeyer:1997mt ;
Yokoyama:1998xd ; Musco:2012au ; Musco:2008hv . This effect is conventionally
described by the scaling relation of the form:

$\displaystyle M_{\rm PBH}=KM_{H}(F-F_{c})^{\gamma_{c}},$ (48)

where $K=3.3$ and $\gamma_{c}=0.35$ are two numerical constants. $F_{c}$ is
the threshold density perturbation for gravitational collapse. The scaling
formula generalizes the peak value $\nu=F/\sigma_{F0}$ to be a function of
both $M_{\rm PBH}$ and $M_{H}$ so that Wu:2020ilx ; Byrnes:2018clq :

$\displaystyle\beta_{c}\left(M_{\rm
PBH},M_{H}\right)=\frac{K}{\sqrt{2\pi}\gamma_{c}\sigma_{F0}(M_{H})}\left(\frac{M_{\rm
PBH}}{KM_{H}}\right)^{1+1/\gamma_{c}}\exp\left[-\frac{F^{2}}{2\sigma_{F0}^{2}(M_{H})}\right],$
(49)

where $F=[M_{\rm PBH}/(KM_{H})]^{1/\gamma_{c}}+F_{c}$.

The PBH mass function based on the extended mass relation at matter-radiation
equality is therefore a sum of the contribution at each epoch of the horizon
mass $M_{H}$, where

$\displaystyle f_{c}(M_{\rm PBH})=\frac{1}{\Omega_{\rm DM,eq}}\int_{\ln M_{\rm
min}}^{\ln M_{\rm eq}}\left(\frac{M_{\rm
eq}}{M_{H}}\right)^{1/2}\beta_{c}\left(M_{\rm PBH},M_{H}\right)d\ln M_{H}.$
(50)

The lower limit $M_{\rm min}$ comes from the upper bound for $F$ as over-large
density perturbation would have already collapsed to form PBHs in earlier
epoch. We adopt a conservative upper bound $F_{\rm max}\leq 2F_{c}$ so that
one can derive the lower limit $M_{\rm min}\geq M_{\rm
PBH}/(KF_{c}^{\gamma_{c}})$ from the scaling relation (48).

Finally, the PBH ratio in dark matter is the sum of all masses at equality:
$f_{\rm PBH}\equiv\int f_{c}(M_{\rm PBH})d\ln M_{\rm PBH}$. As a comparison,
for $f_{\rm PBH}=1$ with $k_{0}=7.3\times 10^{12}$ Mpc-1, we find that
$N_{\ast}=2.6687$ ($2.6645$) with (or without) the effect of critical
collapse.

###  C.3 The finite-size effect

The fact that PBHs are not point particles in real space and they only form at
the local maximum of the density contrast introduces spatial constraints of
the peak value $\nu=F/\sigma_{F0}$ with respect to its first and second
spatial derivatives. These constraints in 3-dimensional space lead to non-
trivial functions $f_{c}(y_{i})$ among the spectral moments in the PBH density
(35) and invoke a system of 10-variable Gaussian statistics Bardeen:1985tr .
In the high-peak limit $\nu_{c}=F_{c}/\sigma_{F0}\gg 1$, the mass fraction
with 3-dimensional spatial constraints is found as Wu:2020ilx ; Young:2014ana
:

$\displaystyle\beta_{\rm
3d}(M_{H})\simeq\frac{1}{2}Q^{3/2}\nu_{c}^{3}\beta_{\rm PS}(M_{H}),$ (51)

where $Q\equiv R^{2}\sigma_{F1}^{2}/(3\sigma_{F0}^{2})$ is an $\mathcal{O}(1)$
factor near the pivot scale $k_{0}$ of the USR transition, and $\nu_{c}^{3}$
is nothing but the finite-size effect for PBHs in 3-dimensional real space.
The detailed derivation of $\beta_{\rm 3d}$ and the extended version with the
effective of critical collapse can be found in Wu:2020ilx .

The finite-size effect usually becomes negligible if $P_{F}$ is a broad
spectrum spanning a wide mass range. By fixing $f_{\rm
PBH}(\delta_{2},N_{\ast},N_{\rm end})=1$ at $N_{\rm end}=20$ for a broad
spectrum with $\delta_{2}=-3.05$, we find that the result from $\beta_{\rm
3d}$ ($\beta_{\rm PS}$) gives $N_{\ast}=2.6423$ ($2.6645$). For a comparison,
a narrower spectrum with $\delta_{2}=-3.1$ at $N_{\rm end}=20$, we find that
$N_{\ast}=2.563$ from $\beta_{\rm 3d}$ and $N_{\ast}=2.584$ from $\beta_{\rm
PS}$.

###  C.4 Non-linear density contrast

The general solution of the Einstein equation for the relation between the
curvature perturbation and the density contrast of radiation is non-linear.
Recent studies Musco:2018rwt ; Young:2019yug ; Harada:2015yda  suggested that
the density perturbation of radiation, $\delta\rho_{r}/\rho_{r}$, in terms of
the compaction function provides a well-defined criterion for PBH formation,
which is related to the profile of the curvature perturbation, $\zeta(r)$,
through the non-linear equation as

$\displaystyle\frac{\delta\rho_{r}}{\rho_{r}}\equiv
F_{NL}=F-\frac{3}{8}F^{2},$ (52)

where $F=-4r_{m}\partial_{r}\zeta(r_{m})$ is the linear density contrast
imposed in the fiducial approach and $r_{m}$ is the local maximum of the
compaction function. Such a non-linear relation promotes the extended mass
relation of critical collapse (48) to be

$\displaystyle M_{\rm PBH}=KM_{H}(F_{\rm
NL}-F_{c})^{\gamma_{c}}=KM_{H}\left(F-\frac{3}{8}F^{2}-F_{c}\right)^{\gamma_{c}}.$
(53)

Completing the square of the linear perturbation (with $2/3>F>F_{c}$
Young:2019yug ), one can derive the modified relation

$\displaystyle
dF=\frac{4}{3\gamma_{c}}\left[\frac{16}{9}-\frac{8}{3}F_{c}-\frac{8}{3}\left(\frac{M_{\rm
PBH}}{KM_{H}}\right)^{1/\gamma_{c}}\right]^{-1/2}\left(\frac{M_{\rm
PBH}}{KM_{H}}\right)^{1/\gamma_{c}}d\ln M_{\rm PBH}.$ (54)

Given that our goal is only to clarify the effect of non-linearity for the
mass function based on the fiducial approach, we compute the PBH density by
using the one-variable formula with Gaussian probability distribution function
as

$\displaystyle\Omega_{\rm
PBH}(M_{H})=\frac{1}{\sqrt{2\pi}\sigma_{F0}}\int_{F_{c}}\frac{M_{\rm
PBH}}{M_{H}}e^{-F^{2}/(2\sigma_{F0}^{2})}dF.$ (55)

Taking the modified relations from (53) and (54), one can derive the mass
fraction from non-linear density contrast of the form

$\displaystyle\beta_{\rm NL}(M_{\rm PBH},M_{H})=$
$\displaystyle\frac{4K}{3\sqrt{2\pi}\gamma_{c}\sigma_{F0}}\left(\frac{M_{\rm
PBH}}{KM_{H}}\right)^{1+1/\gamma_{c}}$ (56)
$\displaystyle\times\left[\frac{16}{9}-\frac{8}{3}F_{c}-\frac{8}{3}\left(\frac{M_{\rm
PBH}}{KM_{H}}\right)^{1/\gamma_{c}}\right]^{-1/2}\exp\left[-\frac{F^{2}}{2\sigma_{F0}^{2}}\right].$

The unified result of non-linear density contrast with the finite-size effect
of density peaks is provided in Young:2019yug .

The extended mass function from non-linear density perturbation can be
computed similarly as (50) by using $\beta_{\rm NL}$. Comparing with the
resulting $N_{\ast}$ for $f_{\rm PBH}(\delta_{2},N_{\ast},N_{\rm end})=1$ with
fixed $\delta_{2}$ and $N_{\rm end}$ based on the fiducial method, the mass
function from non-linear density contrast asks a larger $N_{\ast}$ by
$\mathcal{O}(10^{-1})$.

## References

  * (1) Y. Akrami et al. [Planck], [arXiv:1807.06211 [astro-ph.CO]].
  * (2) M. Tristram, A. J. Banday, K. M. Górski, R. Keskitalo, C. R. Lawrence, K. J. Andersen, R. B. Barreiro, J. Borrill, H. K. Eriksen and R. Fernandez-Cobos, et al. [arXiv:2010.01139 [astro-ph.CO]].
  * (3) I. Affleck and M. Dine, Nucl. Phys. B 249, 361 (1985).
  * (4) A. D. Linde, Phys. Lett. 160B, 243 (1985).
  * (5) A. D. Dolgov, Phys. Rept. 222, 309 (1992).
  * (6) M. Dine and A. Kusenko, Rev. Mod. Phys. 76, 1 (2003) [arXiv:hep-ph/0303065 [hep-ph]].
  * (7) M. Dine, L. Randall and S. D. Thomas, Nucl. Phys. B 458, 291-326 (1996) [arXiv:hep-ph/9507453 [hep-ph]].
  * (8) A. Hook, JHEP 1511, 143 (2015) [arXiv:1508.05094 [hep-ph]].
  * (9) Y. P. Wu and K. Petraki, [arXiv:2008.08549 [hep-ph]].
  * (10) J. Yokoyama, Phys. Rev. D 58, 083510 (1998) [astro-ph/9802357].
  * (11) R. Saito, J. Yokoyama and R. Nagata, JCAP 0806, 024 (2008) [arXiv:0804.3470 [astro-ph]].
  * (12) J. Garcia-Bellido and E. Ruiz Morales, Phys. Dark Univ. 18, 47 (2017) [arXiv:1702.03901 [astro-ph.CO]].
  * (13) K. Kannike, L. Marzola, M. Raidal and H. Veermäe, JCAP 09, 020 (2017) [arXiv:1705.06225 [astro-ph.CO]].
  * (14) C. Germani and T. Prokopec, Phys. Dark Univ. 18, 6 (2017) [arXiv:1706.04226 [astro-ph.CO]].
  * (15) H. Motohashi and W. Hu, Phys. Rev. D 96, no.6, 063503 (2017) [arXiv:1706.06784 [astro-ph.CO]].
  * (16) M. Cicoli, V. A. Diaz and F. G. Pedro, JCAP 1806, no. 06, 034 (2018) [arXiv:1803.02837 [hep-th]].
  * (17) O. Özsoy, S. Parameswaran, G. Tasinato and I. Zavala, JCAP 07, 005 (2018) [arXiv:1803.07626 [hep-th]].
  * (18) C. T. Byrnes, P. S. Cole and S. P. Patil, JCAP 06, 028 (2019) [arXiv:1811.11158 [astro-ph.CO]].
  * (19) S. L. Cheng, W. Lee and K. W. Ng, Phys. Rev. D 99, no.6, 063524 (2019) [arXiv:1811.10108 [astro-ph.CO]].
  * (20) P. Carrilho, K. A. Malik and D. J. Mulryne, Phys. Rev. D 100, no.10, 103529 (2019) [arXiv:1907.05237 [astro-ph.CO]].
  * (21) O. Özsoy and G. Tasinato, JCAP 04, 048 (2020) [arXiv:1912.01061 [astro-ph.CO]].
  * (22) J. Liu, Z. K. Guo and R. G. Cai, Phys. Rev. D 101, no.8, 083535 (2020) [arXiv:2003.02075 [astro-ph.CO]].
  * (23) G. Ballesteros, J. Rey, M. Taoso and A. Urbano, JCAP 08, 043 (2020) [arXiv:2006.14597 [astro-ph.CO]].
  * (24) K. W. Ng and Y. P. Wu, [arXiv:2102.05620 [astro-ph.CO]].
  * (25) S. M. Leach and A. R. Liddle, Phys. Rev. D 63, 043508 (2001) [arXiv:astro-ph/0010082 [astro-ph]].
  * (26) S. M. Leach, M. Sasaki, D. Wands and A. R. Liddle, Phys. Rev. D 64, 023512 (2001) [arXiv:astro-ph/0101406 [astro-ph]].
  * (27) N. C. Tsamis and R. P. Woodard, Phys. Rev. D 69, 084005 (2004) [arXiv:astro-ph/0307463 [astro-ph]].
  * (28) W. H. Kinney, Phys. Rev. D 72, 023515 (2005) [arXiv:gr-qc/0503017 [gr-qc]].
  * (29) J. Martin, H. Motohashi and T. Suyama, Phys. Rev. D 87, no.2, 023514 (2013) [arXiv:1211.0083 [astro-ph.CO]].
  * (30) H. Motohashi, A. A. Starobinsky and J. Yokoyama, JCAP 09, 018 (2015) [arXiv:1411.5021 [astro-ph.CO]].
  * (31) L. Anguelova, P. Suranyi and L. C. R. Wijewardhana, JCAP 02, 004 (2018) [arXiv:1710.06989 [hep-th]].
  * (32) M. S. Turner, Phys. Lett. B 89, 155-159 (1979)
  * (33) J. D. Barrow, E. J. Copeland, E. W. Kolb and A. R. Liddle, Phys. Rev. D 43, 984-994 (1991) doi:10.1103/PhysRevD.43.984
  * (34) A. S. Majumdar, P. Das Gupta and R. P. Saxena, Int. J. Mod. Phys. D 4, 517-529 (1995)
  * (35) D. Baumann, P. J. Steinhardt and N. Turok, [arXiv:hep-th/0703250 [hep-th]].
  * (36) C. Bambi, A. D. Dolgov and A. A. Petrov, JCAP 09, 013 (2009) [arXiv:0806.3440 [astro-ph]].
  * (37) Y. Hamada and S. Iso, PTEP 2017, no.3, 033B02 (2017) [arXiv:1610.02586 [hep-ph]].
  * (38) D. Hooper and G. Krnjaic, Phys. Rev. D 103, no.4, 043504 (2021) [arXiv:2010.01134 [hep-ph]].
  * (39) S. Jyoti Das, D. Mahanta and D. Borah, [arXiv:2104.14496 [hep-ph]].
  * (40) A. Ambrosone, R. Calabrese, D. F. G. Fiorillo, G. Miele and S. Morisi, [arXiv:2106.11980 [hep-ph]].
  * (41) T. Fujita, M. Kawasaki, K. Harigaya and R. Matsuda, Phys. Rev. D 89, no.10, 103501 (2014) [arXiv:1401.1909 [astro-ph.CO]].
  * (42) J. García-Bellido, B. Carr and S. Clesse, [arXiv:1904.11482 [astro-ph.CO]].
  * (43) K. Petraki and R. R. Volkas, Int. J. Mod. Phys. A 28, 1330028 (2013) [arXiv:1305.4939 [hep-ph]].
  * (44) N. F. Bell, K. Petraki, I. M. Shoemaker and R. R. Volkas, Phys. Rev. D 84, 123505 (2011) [arXiv:1105.3730 [hep-ph]].
  * (45) B. von Harling, K. Petraki and R. R. Volkas, JCAP 05, 021 (2012) [arXiv:1201.2200 [hep-ph]].
  * (46) M. M. Flores and A. Kusenko, Phys. Rev. Lett. 126, no.4, 041101 (2021) [arXiv:2008.12456 [astro-ph.CO]].
  * (47) R. K. Jain, P. Chingangbam, J. O. Gong, L. Sriramkumar and T. Souradeep, JCAP 01, 009 (2009) [arXiv:0809.3915 [astro-ph]].
  * (48) R. K. Jain, P. Chingangbam, L. Sriramkumar and T. Souradeep, Phys. Rev. D 82, 023509 (2010) [arXiv:0904.2518 [astro-ph.CO]].
  * (49) R. Allahverdi, K. Enqvist, J. Garcia-Bellido, A. Jokinen and A. Mazumdar, JCAP 06, 019 (2007) [arXiv:hep-ph/0610134 [hep-ph]].
  * (50) H. V. Ragavendra, P. Saha, L. Sriramkumar and J. Silk, [arXiv:2008.12202 [astro-ph.CO]].
  * (51) M. Biagetti, G. Franciolini, A. Kehagias and A. Riotto, JCAP 07, 032 (2018) [arXiv:1804.07124 [astro-ph.CO]].
  * (52) H. Motohashi, S. Mukohyama and M. Oliosi, JCAP 03, 002 (2020) [arXiv:1910.13235 [gr-qc]].
  * (53) Y. F. Cai, X. Chen, M. H. Namjoo, M. Sasaki, D. G. Wang and Z. Wang, JCAP 05, 012 (2018) [arXiv:1712.09998 [astro-ph.CO]].
  * (54) A. A. Starobinsky and J. Yokoyama, Phys. Rev. D 50, 6357 (1994) [astro-ph/9407016].
  * (55) I. Antoniadis, P. O. Mazur and E. Mottola, JCAP 09, 024 (2012) [arXiv:1103.4164 [gr-qc]].
  * (56) A. Strominger, JHEP 10, 034 (2001) [arXiv:hep-th/0106113 [hep-th]].
  * (57) G. S. Ng and A. Strominger, Class. Quant. Grav. 30, 104002 (2013) [arXiv:1204.1057 [hep-th]].
  * (58) D. L. Jafferis, A. Lupsasca, V. Lysov, G. S. Ng and A. Strominger, JHEP 01, 004 (2015) [arXiv:1305.5523 [hep-th]].
  * (59) L. T. Wang and Z. Z. Xianyu, JHEP 02, 044 (2020) [arXiv:1910.12876 [hep-ph]].
  * (60) L. T. Wang and Z. Z. Xianyu, JHEP 11, 082 (2020) [arXiv:2004.02887 [hep-ph]].
  * (61) A. Bodas, S. Kumar and R. Sundrum, JHEP 02, 079 (2021) [arXiv:2010.04727 [hep-ph]].
  * (62) Y. P. Wu, L. Yang and A. Kusenko, JHEP 12, 088 (2019) [arXiv:1905.10537 [hep-ph]].
  * (63) A. Anisimov and M. Dine, Nucl. Phys. B 619, 729-740 (2001) [arXiv:hep-ph/0008058 [hep-ph]].
  * (64) B. J. Carr, Astrophys. J. 201, 1 (1975).
  * (65) J. M. Bardeen, J. R. Bond, N. Kaiser and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * (66) A. M. Green, A. R. Liddle, K. A. Malik and M. Sasaki, Phys. Rev. D 70, 041502 (2004) [astro-ph/0403181].
  * (67) S. Young, C. T. Byrnes and M. Sasaki, JCAP 07, 045 (2014) [arXiv:1405.7023 [gr-qc]].
  * (68) T. Suyama and S. Yokoyama, arXiv:1912.04687 [astro-ph.CO].
  * (69) Y. P. Wu, Phys. Dark Univ. 30, 100654 (2020) [arXiv:2005.00441 [astro-ph.CO]].
  * (70) S. Young and M. Musso, arXiv:2001.06469 [astro-ph.CO].
  * (71) C. M. Yoo, T. Harada, J. Garriga and K. Kohri, PTEP 2018, no. 12, 123E01 (2018) [arXiv:1805.03946 [astro-ph.CO]].
  * (72) V. De Luca, G. Franciolini, A. Kehagias, M. Peloso, A. Riotto and C. Ünal, JCAP 1907, 048 (2019) [arXiv:1904.00970 [astro-ph.CO]].
  * (73) A. Kalaja, N. Bellomo, N. Bartolo, D. Bertacca, S. Matarrese, I. Musco, A. Raccanelli and L. Verde, arXiv:1908.03596 [astro-ph.CO].
  * (74) M. Kawasaki and H. Nakatsuka, Phys. Rev. D 99, no. 12, 123501 (2019) [arXiv:1903.02994 [astro-ph.CO]].
  * (75) S. Young, I. Musco and C. T. Byrnes, arXiv:1904.00984 [astro-ph.CO].
  * (76) C. T. Byrnes, E. J. Copeland and A. M. Green, Phys. Rev. D 86, 043512 (2012) [arXiv:1206.4188 [astro-ph.CO]].
  * (77) Y. Tada and S. Yokoyama, Phys. Rev. D 91, no. 12, 123534 (2015) [arXiv:1502.01124 [astro-ph.CO]].
  * (78) S. Young, D. Regan and C. T. Byrnes, JCAP 1602, no. 02, 029 (2016) [arXiv:1512.07224 [astro-ph.CO]].
  * (79) G. Franciolini, A. Kehagias, S. Matarrese and A. Riotto, JCAP 1803, no. 03, 016 (2018) [arXiv:1801.09415 [astro-ph.CO]].
  * (80) V. Atal and C. Germani, Phys. Dark Univ. 24, 100275 (2019) [arXiv:1811.07857 [astro-ph.CO]].
  * (81) M. Taoso and A. Urbano, [arXiv:2102.03610 [astro-ph.CO]].
  * (82) C. Germani and R. K. Sheth, Phys. Rev. D 101, no.6, 063520 (2020) [arXiv:1912.07072 [astro-ph.CO]].
  * (83) C. T. Byrnes, M. Hindmarsh, S. Young and M. R. S. Hawkins, JCAP 1808, no. 08, 041 (2018) [arXiv:1801.06138 [astro-ph.CO]].
  * (84) J. C. Niemeyer and K. Jedamzik, Phys. Rev. Lett. 80, 5481 (1998) [astro-ph/9709072].
  * (85) J. Yokoyama, Phys. Rev. D 58, 107502 (1998) [gr-qc/9804041].
  * (86) I. Musco, J. C. Miller and A. G. Polnarev, Class. Quant. Grav. 26, 235001 (2009) [arXiv:0811.1452 [gr-qc]].
  * (87) I. Musco and J. C. Miller, Class. Quant. Grav. 30, 145009 (2013) [arXiv:1201.2379 [gr-qc]].
  * (88) N. Aghanim et al. [Planck], Astron. Astrophys. 641, A6 (2020) [arXiv:1807.06209 [astro-ph.CO]].
  * (89) J. M. Ezquiaga, J. García-Bellido and V. Vennin, JCAP 03, 029 (2020) doi:10.1088/1475-7516/2020/03/029 [arXiv:1912.05399 [astro-ph.CO]].
  * (90) D. G. Figueroa, S. Raatikainen, S. Rasanen and E. Tomberg, [arXiv:2012.06551 [astro-ph.CO]].
  * (91) C. Pattison, V. Vennin, D. Wands and H. Assadullahi, [arXiv:2101.05741 [astro-ph.CO]].
  * (92) M. Biagetti, V. De Luca, G. Franciolini, A. Kehagias and A. Riotto, [arXiv:2105.07810 [astro-ph.CO]].
  * (93) T. Harada, C. M. Yoo, T. Nakama and Y. Koga, Phys. Rev. D 91, no.8, 084057 (2015) [arXiv:1503.03934 [gr-qc]].
  * (94) I. Musco, Phys. Rev. D 100, no.12, 123524 (2019) [arXiv:1809.02127 [gr-qc]].
  * (95) M. Boudaud and M. Cirelli, Phys. Rev. Lett. 122, no.4, 041104 (2019) [arXiv:1807.03075 [astro-ph.HE]].
  * (96) W. DeRocco and P. W. Graham, Phys. Rev. Lett. 123, no.25, 251102 (2019) [arXiv:1906.07740 [astro-ph.CO]].
  * (97) R. Laha, Phys. Rev. Lett. 123, no.25, 251101 (2019) [arXiv:1906.09994 [astro-ph.HE]].
  * (98) R. Laha, J. B. Muñoz and T. R. Slatyer, Phys. Rev. D 101, no.12, 123514 (2020) [arXiv:2004.00627 [astro-ph.CO]].
  * (99) B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, Phys. Rev. D 81, 104019 (2010) [arXiv:0912.5297 [astro-ph.CO]].
  * (100) B. Dasgupta, R. Laha and A. Ray, Phys. Rev. Lett. 125, no.10, 101101 (2020) [arXiv:1912.01014 [hep-ph]].
  * (101) V. Poulin, J. Lesgourgues and P. D. Serpico, JCAP 03, 043 (2017) [arXiv:1610.10051 [astro-ph.CO]].
  * (102) S. Clark, B. Dutta, Y. Gao, L. E. Strigari and S. Watson, Phys. Rev. D 95, no.8, 083006 (2017) [arXiv:1612.07738 [astro-ph.CO]].
  * (103) S. Clark, B. Dutta, Y. Gao, Y. Z. Ma and L. E. Strigari, Phys. Rev. D 98, no.4, 043006 (2018) [arXiv:1803.09390 [astro-ph.HE]].
  * (104) A. Hektor, G. Hütsi, L. Marzola, M. Raidal, V. Vaskonen and H. Veermäe, Phys. Rev. D 98, no.2, 023503 (2018) [arXiv:1803.09697 [astro-ph.CO]].
  * (105) S. Mittal, A. Ray, G. Kulkarni and B. Dasgupta, [arXiv:2107.02190 [astro-ph.CO]].
  * (106) B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, [arXiv:2002.12778 [astro-ph.CO]].
  * (107) A. Katz, J. Kopp, S. Sibiryakov and W. Xue, JCAP 12, 005 (2018) [arXiv:1807.11495 [astro-ph.CO]].
  * (108) A. Ray, R. Laha, J. B. Muñoz and R. Caputo, [arXiv:2102.06714 [astro-ph.CO]].
  * (109) F. Capela, M. Pshirkov and P. Tinyakov, Phys. Rev. D 87, no.12, 123524 (2013) [arXiv:1301.4984 [astro-ph.CO]].
  * (110) P. W. Graham, S. Rajendran and J. Varela, Phys. Rev. D 92, no.6, 063007 (2015) [arXiv:1505.04444 [hep-ph]].
  * (111) A. M. Green and B. J. Kavanagh, [arXiv:2007.10722 [astro-ph.CO]].
  * (112) R. J. Nemiroff and A. Gould, Astrophys. J. Lett. 452, L111 (1995) [arXiv:astro-ph/9505019 [astro-ph]].
  * (113) S. Jung and T. Kim, Phys. Rev. Res. 2, no.1, 013113 (2020) [arXiv:1908.00078 [astro-ph.CO]].
  * (114) Y. Bai and N. Orlofsky, Phys. Rev. D 99, no.12, 123019 (2019) [arXiv:1812.01427 [astro-ph.HE]].
  * (115) B. Dutta, A. Kar and L. E. Strigari, JCAP 03, 011 (2021) [arXiv:2010.05977 [astro-ph.HE]].
  * (116) A. Coogan, L. Morrison and S. Profumo, Phys. Rev. Lett. 126, no.17, 171101 (2021) [arXiv:2010.04797 [astro-ph.CO]].
  * (117) G. Ballesteros, J. Coronado-Blázquez and D. Gaggero, Phys. Lett. B 808, 135624 (2020) [arXiv:1906.10113 [astro-ph.CO]].
  * (118) A. Kashlinsky, Y. Ali-Haimoud, S. Clesse, J. Garcia-Bellido, L. Wyrzykowski, A. Achucarro, L. Amendola, J. Annis, A. Arbey and R. G. Arendt, et al. [arXiv:1903.04424 [astro-ph.CO]].