# Baryogenesis from ultralight primordial black holes and strong gravitational
waves from cosmic strings

Satyabrata Datta<EMAIL_ADDRESS>Saha Institute of Nuclear
Physics, HBNI, 1/AF Bidhannagar, Kolkata 700064, India  Ambar Ghosal
<EMAIL_ADDRESS>Saha Institute of Nuclear Physics, HBNI, 1/AF
Bidhannagar, Kolkata 700064, India  Rome Samanta<EMAIL_ADDRESS>CEICO,
Institute of Physics of the Czech Academy of Sciences, Na Slovance 1999/2, 182
21 Prague 8, Czech Republic

###### Abstract

Ultralight primordial black holes (PBHs)($\lesssim 10^{9}$g) completely
evaporate via Hawking radiation (HR) and produce all the particles in a given
theory regardless of their other interactions. If the right handed (RH)
neutrinos are produced from PBH evaporation, successful baryogenesis via
leptogenesis predicts mass scale of RH neutrinos as well as black holes. We
show that, given the lepton number violation (generation of RH neutrino
masses) in the theory is a consequence of a gauged $U(1)$ breaking which is
then followed by the formation of PBHs, a network of cosmic strings naturally
gives rise to strong stochastic gravitational wave (GW) signal at the
sensitivity level of pulsar timing arrays (PTA) and LIGO5. Besides, due to a
transient period of black hole domination in the early universe, for which
baryon asymmetry is independent of initial PBH density, a break in the GW
spectra occurs around MHz frequency. Therefore, to observe the break along
with the usual GW signal by the emission of gravitons via HR, GW detectors at
higher frequencies are called for. The recent finding by the NANOGrav PTA of a
stochastic common spectrum process (interpreted as GWs) across many pulsars is
in tension with PBH baryogenesis for large cosmic string loops ($\alpha\simeq
0.1$).

##  I Introduction

Physics of Primordial Black Holes (PBH)[1, 2, 3] that entails rich
cosmological consequences, has gained a renewed interest after the discovery
of Gravitational Waves (GW) from black hole mergers by LIGO and Virgo
collaboration [4, 5, 6, 7, 8, 9, 10]. Starting from the mass $M_{BH}\simeq
0.1$g up to few hundreds of solar mass, PBHs constitute seemingly intriguing
physical aspects. Long-lived PBHs with mass ($M_{BH}\gtrsim 10^{15}$g) that
survive until today may contribute to a significant portion of Dark Matter
(DM)[11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21], whereas PBHs with
$M_{BH}\lesssim 10^{9}$g completely evaporate before Big Bang Nucleosynthesis
(BBN) ($T\sim 4$ MeV) via Hawking Radiation (HR)[22]. On the other hand, the
intermediate mass black holes $10^{9}{\rm g}\lesssim M_{BH}\lesssim 10^{15}$g
are subjected to several constraints, e.g., from BBN [23, 24, 25]. Though the
ultralight PBHs ($0.1{\rm g}\lesssim M_{BH}\lesssim 10^{9}$g) completely
evaporate, the particle emission being democratic in nature and poor
constraints on initial PBH density at the formation111can be constrained,
e.g., with the production of gravitinos in supersymmetric models[26, 27].,
they leave plenty of room that can accommodate new physics. Perhaps this is
why parallel to the usual study of PBH as Dark matter (DM) ($\gtrsim
10^{15}$g), physics of ultralight PBHs have gained a lot of interest in recent
time. This includes, for example, production of DM particles[28, 29, 30, 31,
32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42], generation of baryon
asymmetry[43, 44, 28, 45, 46, 47, 48, 49], change in GW spectrum induced by
scalar perturbations[50], GWs induced by ultralight primordial black holes[51,
52] and study of vacuum stability of the Standard Model (SM) Higgs[53, 54]. In
this work, we confine ourselves with the detailed study of baryogenesis via
non-thermal leptogenesis by ultralight PBH evaporation and its possible test
with GWs.

Amongst several mechanisms of baryogenesis, the simplest one is baryogenesis
via leptogenesis[55, 56, 57, 58, 59, 60, 61, 62, 63] from the decays of heavy
Right Handed (RH) neutrinos which are introduced in the SM to obtain light
neutrino masses via Type-I seesaw mechanism[64, 65, 66, 67, 68]. However, to
test standard leptogenesis within Type-I seesaw, one has to either lower the
RH mass scales for collider searches[69, 70, 71, 72] or invoke restrictions on
the relevant parameters, e.g., considering discrete symmetries[73, 74, 75, 76]
and theories like SO(10) grand unification (GUT)[77, 78, 79, 80, 81].
Recently, a pathway to probe leptogenesis with GWs has been introduced in
Ref.[82]. The idea is based on formation of cosmic string network[83, 84, 85]
upon breaking of an Abelian symmetry and consequent emission of GWs. Though
emission of GWs from cosmic string loops still remains under debate, numerical
simulations based on the Nambu–Goto action[86, 87] indicate that cosmic string
loops loose energy dominantly via GW radiation, if the underlying broken
symmetry corresponds to a local gauge symmetry, e.g., $U(1)$, which in the
context of seesaw, can be readily promoted to a $U(1)_{B-L}$[88, 89, 90]. Once
$U(1)_{B-L}$ gets spontaneously broken, RH neutrinos get mass and cosmic
strings are generated. The scale of leptogenesis (RH masses) and amplitude of
the GWs, both are related to the symmetry breaking scale ($\Lambda_{CS}$).
Therefore, one may probe the scale of leptogenesis through the detection of
GWs. With this idea, two other variations of leptogenesis have also been
studied[91, 92]. In fact, due to the most distinguishable feature of strong GW
signal from cosmic strings across a wide range of frequencies, there is a
rapidly growing interest in this field[93, 94, 95, 96, 97, 98, 99, 100, 101,
102, 103, 104, 105, 106]. Moreover, the recent finding of stochastic common
spectrum processes across many pulsars by the NANOGrav collaboration[107] has
strengthened the possibility to consider GW emission from cosmic strings,
since the new data can be explained better with cosmic string models[92, 108,
109, 110, 111, 112, 113] than the single value power spectral density as
suggested by the models of supermassive black holes (SMBHs).

In this work, we follow a similar idea as in Ref.[82, 91, 92] to study ‘GW
friendly leptogenesis (dynamical generation of RH mass scale)’ from PBH
evaporation wherein massive RH neutrinos are emitted by the PBHs and then
decay CP asymmetrically to create lepton asymmetry[58]. However, as we explain
below, unlike Ref.[82], here we give more emphasis on the PBH mass rather than
RH neutrino masses. We discuss a scenario where after the end of inflation,
there is a $U(1)_{B-L}$ phase-transition (at $\Lambda_{CS}$) so that RH
neutrinos become massive[82] which is then followed by the formation of PBHs
(at $T_{Bf}$) that emit massive RH neutrinos thereafter, i.e.,
$\Lambda_{CS}\gtrsim T_{Bf}$. In this scenario, i.e., in the non-thermal
leptogenesis from PBH evaporation, the requirement to obtain correct baryon
asymmetry, strongly constrains RH neutrino masses as well as the mass of the
PBHs. The most interesting outcome is, one obtains an upper bound on $M_{BH}$
which translates into a lower bound on the formation temperature
($T_{Bf}\gtrsim 10^{15}$ GeV) of PBHs and hence on $\Lambda_{CS}$. Since with
the decrease of $M_{BH}$, the formation temperature increases[16], one
naturally obtains GWs with large amplitude throughout the parameter space of a
successful leptogenesis, given the fact that the GW amplitude increases with
the breaking scale $\Lambda_{CS}$[94] (for example, $\Lambda_{CS}\gtrsim
10^{15}$ GeV would correspond to $\Omega_{GW}h^{2}(f)\gtrsim 10^{-9}$).
Therefore, the mechanism is fully testable in most of the planned and present
GW detectors over a wide range of frequencies. We focus more on the case in
which after the formation, PBHs dominate the energy density of the universe
until they completely evaporate, which is also motivated by the fact that in
this case, the baryon asymmetry is independent of PBH density at
formation[43]. With such a non-standard cosmological evolution (a matter
domination before the usual radiation domination) one would expect a break in
the GW spectrum[94, 91]. We show that, for a successful leptogenesis scenario,
this break happens at much higher frequencies compare to the sensitivity, e.g,
of LIGO ($\sim 25$ Hz). Finally, we discuss the relevance of PBH baryogenesis
in the light of recent 12.5 yrs NANOGrav PTA data and conclude that the
predicted GW amplitude for large cosmic string loops $\alpha^{max}\simeq
0.1$[114, 115] is in tension with the NANOGrav data.

The rest of the paper is organised as follows: In Sec.II, we discuss the
dynamics of PBH energy density and its impact on thermal leptogenesis. Sec.III
is devoted to the discussion on baryogenesis via non-thermal leptogenesis in
radiation domination and in PBH domination. In Sec.IV, we review the GW
emission from cosmic string and present the numerical results which also
include the fitting to the 12.5 yrs NANOGrav data. We conclude in Sec.V.

##  II PBH dynamics: Radiation domination Vs. Black hole domination

In this section, we first derive the condition on PBH mass and its density at
the formation to understand the scenarios like PBH domination and radiation
domination after the time of black hole formation (Ref.[16]: a review that
discusses PBH formation mechanisms). For the radiation energy density to
dominate the universe until black hole evaporation, the ratio of the energy
density of black holes and the radiation at the time of black hole evaporation
would satisfy

$\displaystyle r(t_{ev})=\frac{\rho_{BH}(t_{ev})}{\rho_{\rm R}(t_{ev})}<1,$
(II.1)

where $t_{ev}$ is the time at evaporation, $\rho_{BH}$ and $\rho_{R}$ are the
black hole and radiation densities respectively. Thus the ratio of the $r$
parameters at the time of black hole formation and evaporation can be
calculated as

$\displaystyle\frac{r(t_{ev})}{r(t_{Bf})}=\frac{a(t_{ev})}{a(t_{Bf})}=\left(\frac{t_{ev}}{t_{Bf}}\right)^{1/2}$
(II.2)

with $t_{Bf}$ being the time at black hole formation. Recalling the Friedmann
equation

$\displaystyle H(t)^{2}=\frac{8\pi}{3M_{Pl}^{2}}\rho_{R}(t),~{}~{}{\rm
with}~{}~{}M_{Pl}=1.22\times 10^{19}~{}{\rm GeV}$ (II.3)

and using the Hubble parameter in radiation domination and the radiation
density as

$\displaystyle
H=\frac{1}{2t},~{}~{}\rho_{R}(T)=\frac{\pi^{2}g_{*}(T)T^{4}}{30},$ (II.4)

where $g_{*}(T)$ ($\simeq$ 106.75 in SM) is the effective degrees of freedom
that contribute to the radiation, Eq.II.2 can be re-written as

$\displaystyle\frac{r(t_{ev})}{r(t_{Bf})}=\left(\frac{g_{*}(T_{Bf})}{g_{*}(T_{ev})}\right)^{1/4}\frac{T_{Bf}}{T_{ev}}.$
(II.5)

Therefore Eq.II.1 translates to the condition

$\displaystyle
r(t_{Bf})\equiv\beta<\left(\frac{g_{*}(T_{ev})}{g_{*}(T_{Bf})}\right)^{1/4}\frac{T_{ev}}{T_{Bf}}.$
(II.6)

We now need to calculate the black hole formation temperature ($T_{Bf}$) and
the evaporation temperature $(T_{ev})$. The formation temperature can be
calculated once we know the formation time $t_{Bf}$. Assuming radiation
domination at black hole formation, the mass of the black hole from
gravitational collapse is typically close to the value enclosed by the post-
inflation particle horizon and is given by[16]

$\displaystyle M_{BH}=\gamma\frac{4}{3}\pi(H_{Bf}^{-1})^{3}\rho_{Bf}~{}~{}{\rm
with}~{}~{}\rho_{Bf}=\frac{3H_{Bf}^{2}M_{Pl}^{2}}{8\pi},~{}~{}H_{Bf}=\frac{1}{2t_{Bf}}.$
(II.7)

The quantity $\gamma\simeq 0.2$ is a numerical factor which depends on the
details of the gravitational collapse mechanism. Thus the time $t_{Bf}$ can be
calculated as

$\displaystyle t_{Bf}=\frac{M_{BH}}{M_{Pl}^{2}\gamma}.$ (II.8)

Now using the Friedmann equation (Eq.II.3) and radiation energy density
(Eq.II.4) one obtains the black hole formation temperature as

$\displaystyle
T_{Bf}=\left(\frac{45\gamma^{2}}{16\pi^{3}g_{*}(T_{Bf})}\right)^{1/4}\left(\frac{M_{Pl}}{M_{BH}}\right)^{1/2}M_{Pl}.$
(II.9)

We now calculate the black hole evaporation temperature $T_{ev}$. Assuming
radiation dominated universe throughout the evolution of the black holes, we
can write the Hubble parameter as

$\displaystyle H(t_{ev})^{2}=\frac{1}{4t_{ev}^{2}}\simeq\frac{1}{4\tau^{2}},$
(II.10)

where the $\tau$ is the lifetime of the PBH. Therefore, from the Friedmann
equation in Eq.II.3 and the expression for the radiation energy density in
Eq.II.4, we can calculate the temperature $T_{ev}$ as

$\displaystyle
T_{ev}=\left(\frac{45M_{Pl}^{2}}{16\pi^{3}g_{*}(T_{ev})\tau^{2}}\right)^{1/4}.$
(II.11)

The lifetime of the PBH can be calculated from the dynamics of the mass loss
of a PBH via Hawking radiation[22]. The rate at which PBH looses its mass is
given by

$\displaystyle-\frac{dM_{BH}}{dt}=f_{ev}(4\pi r_{BH}^{2})\frac{dE}{dt},$
(II.12)

where the rate of energy loss can be calculated as

$\displaystyle\frac{dE}{dt}=2g_{BH}(T_{BH})\pi^{2}\int_{0}^{\infty}d\nu\frac{\nu^{3}}{{\rm
exp}(2\pi\nu/T)-1}=\frac{\pi^{2}}{120}g_{BH}(T_{BH})T_{BH}^{4}.$ (II.13)

The quantity $f_{ev}$ is the efficiency of PBH evaporation and $r_{BH}$ is the
Schwarzschild radius of the PBH given by

$\displaystyle r_{BH}=2GM_{BH}.$ (II.14)

The quantity $g_{BH}(T_{BH})$ counts the bosonic and fermionic degrees of
freedom below $T_{BH}$.

Figure 1: Left: $\beta$ vs. $M_{BH}$. The solid black line separates the
regions of PBH domination and radiation domination. With the red stars, we
indicate three benchmark points for which we show the evolution of the energy
densities by solving the Boltzmann equations. Right: PBH formation temperature
and the Hawking temperature as a function of PBH mass. In each of the plots,
the light blue region is excluded by the constraint from CMB on Hubble scale
of inflation $H_{\rm inf}\lesssim 3\times 10^{14}$ GeV at 95$\%$ CL[116].

Recalling the PBH temperature ($T_{BH}=M_{Pl}^{2}/8\pi M_{BH}$[22]) and using
the expression for $r_{BH}$, Eq.II.12 can be recast as

$\displaystyle\frac{dM_{BH}}{dt}=-\frac{\mathcal{G}g_{*B}(T_{BH})}{30720\pi}\frac{M_{Pl}^{4}}{M^{2}_{BH}},$
(II.15)

where we have used $f_{ev}g_{BH}(T_{BH})=\mathcal{G}g_{*B}(T_{BH})$ with
$\mathcal{G}\simeq 3.8$ being the graybody factor and $g_{*B}\simeq 100$
(counts all the particle species below $T_{BH}$) in the SM with three left
handed neutrinos[117]. It is now trivial to find the lifetime $\tau$ as

$\displaystyle\tau=\int_{t_{Bf}}^{t_{ev}}dt=-\int_{M_{BH}}^{0}dM_{BH}\frac{30720\pi
M_{BH}^{2}}{\mathcal{G}g_{*B}(T_{BH})M_{Pl}^{4}}=\frac{10240\pi
M_{BH}^{3}}{\mathcal{G}g_{*B}(T_{BH})M_{Pl}^{4}}.$ (II.16)

Therefore, combining Eq.II.16, Eq.II.11 and Eq.II.9, we can recast Eq.II.6 as

$\displaystyle\beta<\gamma^{-1/2}\left(\frac{\mathcal{G}g_{*B}(T_{BH})}{10240\pi}\right)^{1/2}\frac{M_{Pl}}{M_{BH}}.$
(II.17)

On the left-hand side of figure 1, the relative PBH density $\beta$ as a
function of PBH mass has been shown. As one sees, for very light PBHs to
dominate the energy density, the relative PBH density at formation should be
much higher. Thus in the context of gravitino production[26, 27] very light
PBH domination may not be possible. On the right-hand side, we plot the PBH
formation temperature and the PBH temperature as a function of PBH mass. From
this figure one infers that due to the chosen breaking scale of $U(1)_{B-L}$,
i.e., $\Lambda_{CS}\gtrsim T_{Bf}$, the RH mass scale could be higher as well
as lower than $T_{BH}$, depending upon the perturbative Yukawa coupling
($y_{N_{i}}N_{iR}N_{iR}\phi_{B-L}$) of the $U(1)_{B-L}$ scalar and the RH
neutrino fields. Therefore, in the computation of leptogenesis, we need to
consider RH neutrino production from PBH in both the scenarios, i.e.,
$T_{BH}>M_{R}$ as well as $T_{BH}<M_{R}$. To understand the dynamics of black
hole vs. radiation energy density throughout the evolution of the universe,
one needs to solve the following Boltzmann equations[118]

$\displaystyle\frac{d\rho_{R}}{dt}+4H\rho_{R}=-\frac{\dot{M}_{BH}}{M_{BH}}\rho_{BH},$
(II.18)
$\displaystyle\frac{d\rho_{BH}}{dt}+3H\rho_{BH}=+\frac{\dot{M}_{BH}}{M_{BH}}\rho_{BH},$
(II.19)
$\displaystyle\frac{ds}{dt}+3Hs=-\frac{\dot{M}_{BH}}{M_{BH}}\frac{\rho_{BH}}{T},$
(II.20)

where Eq.II.20 represents the non-conservation of entropy due to black hole
evaporation and $\dot{M}_{BH}\equiv dM_{BH}/dt$. It is convenient to recast
the above set of equations in terms of temperature. From Eq.II.20, the
relation between temperature and time can be derived as

$\displaystyle\frac{1}{T}\frac{dT}{dt}=-\left(H+\frac{1}{3g_{*s}(T)}\frac{dg_{*s}(T)}{dt}+\frac{\dot{M}_{BH}}{M_{BH}}\frac{\rho_{BH}}{4\rho_{R}}\right)=-\tilde{H}.$
(II.21)

where we have used

$\displaystyle s(t)=\frac{2\pi^{2}g_{*s}T^{3}}{45}$ (II.22)

and in the right-hand side of Eq.II.21, we have introduced the quantity
$\rho_{R}$ after considering $g_{*}(T)\equiv g_{*s}(T)$. Using Eq.II.21,
Eq.II.18 and Eq.II.19 can be re-written as

$\displaystyle\frac{d\rho_{R}}{dz}+\frac{4}{z}\rho_{R}$ $\displaystyle=$
$\displaystyle 0$ (II.23)
$\displaystyle\frac{d\rho_{BH}}{dz}+\frac{3}{z}\frac{H}{\tilde{H}}\rho_{BH}$
$\displaystyle-$
$\displaystyle\frac{\dot{M}_{BH}}{M_{BH}}\frac{1}{z\tilde{H}}\rho_{BH}=0,$
(II.24)

where $z=M_{0}/T$ with $M_{0}$ being an arbitrary mass scale which we consider
as $M_{0}\equiv T_{Bf}$. On the left-hand side of figure 2, we show the
normalised energy densities of radiation and PBHs for different values of
$\beta$, as indicated in figure 1. We take a benchmark value $M_{BH}=10$ g for
computation. It is clear that with the increase of $\beta$, one gets a larger
epoch of PBH domination which has a great impact on the thermally produced
lepton asymmetry. As we mentioned in the introduction, we shall consider only
non-thermal production of RH neutrinos by PBH evaporation. This requires a
justification that the non-thermal contribution from PBH dominates over the
thermal contribution from the plasma.

Figure 2: Left: Evolution of radiation and black hole energy densities
(normalised to the total energy density). Right: A representative figure to
show the washout (via entropy production by PBH evaporation) of the thermally
produced lepton asymmetry (solid line). The solid black line representing
$N_{B-L}$ indicates the thermally produced asymmetry with no PBHs. We choose
$\Delta=10^{-3}$, $\varepsilon\sim 10^{-6}$ and $K_{1}=50.$

As we explain below, this is indeed the case for a significant period of PBH
domination (large $\beta$). In presence of PBHs, using the time-temperature
relation in Eq.II.21, the standard Boltzmann equations (BEs) for thermal
leptogenesis[61] ($N_{1}$-dominated scenario and without flavour effects[119,
120, 121, 122, 123, 124]) can be generalised to

$\displaystyle\frac{dN_{N_{1}}}{dz}$ $\displaystyle=$
$\displaystyle\frac{3}{\tilde{H}z}(\tilde{H}-H)N_{N_{1}}-\frac{\Gamma_{D_{N_{1}}}}{\tilde{H}z}(N_{N_{1}}-N_{N_{1}}^{\rm
eq}),$ (II.25) $\displaystyle\frac{dN_{B-L}}{dz}$ $\displaystyle=$
$\displaystyle-\varepsilon_{1}\frac{\Gamma_{D_{N_{1}}}}{\tilde{H}z}(N_{N_{1}}-N_{N_{1}}^{\rm
eq})-\left(\frac{\Gamma_{ID_{N_{1}}}}{\tilde{H}z}-\frac{3}{\tilde{H}z}(\tilde{H}-H)\right)N_{B-L},$
(II.26)

where

$\displaystyle\Gamma_{D_{N_{1}}}=K_{1}\left(\frac{m^{*}}{8\pi
v^{2}}\right)\Delta^{2}T_{Bf}^{2}\frac{\mathcal{K}_{1}(z_{1})}{\mathcal{K}_{2}(z_{1})},~{}~{}\Gamma_{ID_{N_{1}}}=\frac{1}{4}K_{1}\left(\frac{m^{*}}{8\pi
v^{2}}\right)\Delta^{2}T_{Bf}^{2}z_{1}^{2}\mathcal{K}_{1}(z_{1})$ (II.27)

In Eq.II.27, $N_{N_{i}}$s are number densities normalised to the ultra-
relativistic number density of $N_{1}$ ($\propto T^{3}$), $\mathcal{K}_{i}$s
are the modified Bessel functions, $N_{N_{1}}^{\rm
eq}=\frac{1}{2}z_{1}^{2}\mathcal{K}_{2}(z_{1})$ and $\Delta=M_{1}/T_{Bf}$
which implies $z_{1}(=M_{1}/T)=z\Delta$. The quantity $K_{1}$ depends on the
Yukawa couplings of SM Higgs, lepton doublets and $N_{1}$ which in the strong
washout scenario has a value ranging typically from $K_{1}\sim 3-50$[61]. The
equilibrium neutrino mass is given by $m^{*}\simeq 10^{-3}$ eV. The quantity
$\varepsilon$ is the CP violation responsible for leptogenesis and has an
upper bound in the hierarchical limit of the RH neutrinos as (see Ref.[125]
for a derivation)

$\displaystyle\varepsilon_{1}\lesssim\frac{3M_{1}m_{3}}{8\pi v^{2}},$ (II.28)

where $v=174$ GeV is the vacuum expectation value of the SM Higgs and we
assume normal mass ordering of the light neutrinos with $m_{3}$ being the
heaviest light neutrino. Note in Eq.II.25 and Eq.II.26, that in the limit
$\tilde{H}\rightarrow H$ we recover the standard BEs for thermal
leptogenesis[61] (no black holes), i.e.,

$\displaystyle\frac{dN_{N_{1}}}{dz}$ $\displaystyle=$ $\displaystyle-
D_{N_{1}}(N_{N_{1}}-N_{N_{1}}^{\rm eq}),$ (II.29)
$\displaystyle\frac{dN_{B-L}}{dz}$ $\displaystyle=$
$\displaystyle-\varepsilon_{1}D_{N_{1}}(N_{N_{1}}-N_{N_{1}}^{\rm
eq})-W^{ID}_{N_{1}}N_{B-L},$ (II.30)

where

$\displaystyle\frac{\Gamma_{D_{N_{1}}}}{Hz}\equiv
D_{N_{N_{1}}}=K_{1}\Delta^{2}z\frac{\mathcal{K}_{1}(z_{1})}{\mathcal{K}_{2}(z_{1})},~{}~{}~{}\frac{\Gamma_{ID_{N_{1}}}}{Hz}\equiv
W^{ID}_{N_{1}}=\frac{1}{4}K_{1}\Delta^{4}z^{3}\mathcal{K}_{1}(z_{1}).$ (II.31)

The main effect that is caused by introducing PBHs, that in Eq.II.26, in
addition to the standard inverse decays ($\Gamma_{ID_{N_{1}}}$), one has a new
‘washout term’ which at the time of PBH evaporation could be strong enough to
erase any asymmetry produced by the thermal plasma. We solve the four coupled
equations presented in Eq.II.23, Eq.II.24, Eq.II.25 and Eq.II.26 to understand
the dynamics of the thermally produced lepton asymmetry in presence of PBHs.
On the right hand side of figure 2, it can be seen that as one goes from a
radiation dominated phase to a PBH dominated phase (i.e., when $\beta$
increases), the contribution from the thermal plasma becomes negligible. For
$\beta=10^{-7}$, the red dotted line and the black solid line (solution of
Eq.II.29 and Eq.II.30) practically coincide. This means for lower values of
$\beta$, the effect of PBHs is negligible, as expected. The produced entropy
can be calculated by solving

$\displaystyle\frac{da}{dz}=\left(1-\frac{\dot{M}_{BH}}{M_{BH}}\frac{\rho_{BH}}{4\rho_{R}\tilde{H}}\right)\frac{a}{z},$
(II.32)

where we have neglected the time variation of $g_{*s}$. On the right-hand side
of figure 2, the quantity $\tilde{S}\propto a^{3}/z^{3}$ shows how for a
significant period of PBH domination, a large amount of entropy is dumped into
plasma and the thermally produced asymmetry gets diluted. One can analytically
quantify this dilution as well. Assuming the washout by the inverse decays and
PBHs do not coincide in time, the magnitude of the final asymmetry produced by
the thermal plasma in presence of PBHs is then given by

$\displaystyle N_{B-L}^{\rm Final}\simeq N_{B-L}^{\rm Thermal}{\rm
Exp}\left[-\int_{\infty}^{0}dz\frac{3}{\tilde{H}z}\frac{\dot{M}_{BH}}{M_{BH}}\frac{\rho_{BH}(z)}{4\rho_{R}(z)}\right].$
(II.33)

Thus overall, with the discussion above, we now understand why a long period
of PBH domination is preferred to study a purely non-thermal leptogenesis and
as we will see in the next section, in PBH domination the non-thermally
produced asymmetry is independent of $\beta$. Let’s now calculate the number
of a particular species of particles (say $`X$’, in our case these are RH
neutrinos) emitted by a black hole which is necessary to compute the final
baryon asymmetry. The differential number of particles emitted by a black hole
can be obtained as[44]

$\displaystyle
dN=dE/3T_{BH}=\frac{M_{Pl}^{2}}{24\pi}\frac{1}{T_{BH}^{3}}dT_{BH},$ (II.34)

where we have used

$\displaystyle
dE\equiv-d(M_{BH})=\frac{M_{Pl}^{2}}{8\pi}\frac{dT_{BH}}{T_{BH}^{2}}$ (II.35)

and consider mean energy of the radiated particles as $\bar{E}=3T$. Assuming
the black holes completely evaporate, the total number of $X$ particle emitted
can be calculated as

$\displaystyle
N_{X}=\frac{g_{X}}{g_{*B}}\int_{T_{BH}}^{\infty}dN=\frac{4\pi}{3}\frac{g_{X}}{g_{*B}}\left(\frac{M_{BH}}{M_{Pl}}\right)^{2}~{}~{}{\rm
for}~{}~{}T_{BH}>M_{X},$ (II.36) $\displaystyle
N_{X}=\frac{g_{X}}{g_{*B}}\int_{M_{X}}^{\infty}dN=\frac{1}{48\pi}\frac{g_{X}}{g_{*B}}\left(\frac{M_{Pl}}{M_{X}}\right)^{2}~{}~{}{\rm
for}~{}~{}T_{BH}<M_{X}.$ (II.37)

Having set up all the necessary pre-requisites, we now proceed to the detailed
calculation of leptogenesis from PBH evaporation.

##  III Leptogenesis from black hole evaporation in radiation and in black
hole domination

In this section, we first discuss non-thermal leptogenesis in radiation
domination. However, as discussed in the previous section, while in radiation
domination, one has to be bit more careful, since the thermal contribution to
the lepton asymmetry pops up. Therefore, either one has to distinguish the
parameter space of thermal and non-thermal leptogenesis or considering lepton
asymmetry produced by PBHs as an additional contribution, re-analyse scenarios
of thermal leptogenesis. That could be interesting to some extent, given the
fact that in that case, the predictions of thermal leptogenesis may change
significantly.

Figure 3: Upper Panel: Allowed parameter space for successful leptogenesis in
radiation domination (white region) for $\beta=\beta_{\max}\times 10^{-2}$.
The same (gray region) for $\beta=\beta_{\max}\times 10^{-1}$. All other
regions are excluded by the constraints discussed in the text. Lower panel:
Allowed parameter space for successful leptogenesis in PBH domination (white
region).

Assuming no further entropy production after PBH evaporation, the final
asymmetry produced non-thermally by PBH evaporation can be computed as

$\displaystyle Y_{B}=a_{\rm sph}N_{\nu}\varepsilon\frac{n_{BH}(T_{ev})}{s},$
(III.1)

where $a_{\rm sph}\simeq 1/3$ is the sphaleron conversion coefficient. Now
given the upper bound on the CP asymmetry parameter in Eq.II.28, one gets the
following constraints on the parameter space

$\displaystyle M_{R}<\left(\frac{m_{3}M_{Pl}^{2}}{384Y_{B}^{\rm
Obs}\pi^{2}v^{2}}\right)\left(\frac{g_{R}}{g_{*B}}\right)\frac{n_{BH}(T_{ev})}{s}~{}~{}{\rm
for}~{}~{}M_{R}>T_{BH},$ (III.2)

and

$\displaystyle 1/M_{R}<\left(\frac{m_{3}M_{BH}^{2}}{6Y_{B}^{\rm
Obs}M_{Pl}^{2}v^{2}}\right)\left(\frac{g_{R}}{g_{*B}}\right)\frac{n_{BH}(T_{ev})}{s}~{}~{}{\rm
for}~{}~{}M_{R}<T_{BH},$ (III.3)

where $g_{R}$ is the number of degrees of freedom associated with the RH
neutrino and $Y_{B}^{\rm Obs}\simeq 8.75\times 10^{-11}$[116]. Potentially,
there could be another constraint if the RH scale is below the PBH evaporation
temperature. In that case, the produced lepton asymmetry is washed out by the
inverse decays. This constraint can be worked out after properly computing the
inverse decay out-of-equilibrium temperature as[61, 126]

$\displaystyle M_{R}/T^{\rm Inv}_{\rm out}\equiv
M_{R}/T_{ev}>\Phi(K)~{}~{}{\rm
with}~{}~{}\Phi(K_{1})=2+4K_{1}^{0.13}e^{\frac{-2.5}{K_{1}}}$ (III.4)

which translates into the constraint on $M_{R}$ as

$\displaystyle
M_{R}>\left(\frac{45M_{Pl}^{2}}{16\pi^{3}g_{*}(T_{ev})\tau^{2}}\right)^{1/4}\Phi(K_{1}).$
(III.5)

With no fine-tuning in the seesaw formula, $K_{1}\simeq 10^{3}m_{1}$, with
$m_{1}$ being bounded from above from cosmology as $m_{1}<51$ meV[116]. We use
this value in the numerical computation. The factor $n_{BH}/s$ factor in
radiation domination can be calculated as

$\displaystyle
Y_{BH}^{R}(T_{ev})\equiv\frac{n_{BH}}{s}\left(T_{ev}\right)=\frac{1}{M_{BH}}\frac{\rho_{R}(T_{ev})}{s(T_{ev})}r(T_{ev}).$
(III.6)

Now using Eq.II.4, Eq.II.5, Eq.II.9 and Eq.II.22, $Y_{BH}^{\rm R}(T_{ev})$ can
be calculated as

$\displaystyle Y_{BH}^{\rm
R}(T_{ev})=\frac{3\beta}{4}\left(\frac{g_{*}(T_{Bf})}{g_{*}(T_{ev})}\right)^{1/4}\frac{g_{*}(T_{ev})}{g_{*s}(T_{ev})}\left(\frac{45\gamma^{2}}{16\pi^{3}g_{*}(T_{Bf})}\right)^{1/4}\left(\frac{M_{Pl}}{M_{BH}}\right)^{3/2}.$
(III.7)

On the other hand, for a significant period of PBH domination, the number
density $n_{BH}$ can be calculated as

$\displaystyle n_{BH}(T_{ev})=\frac{M_{Pl}^{2}}{6\pi M_{BH}\tau^{2}},$ (III.8)

where we have assumed $H(T_{ev})\simeq 2/3\tau$. Now deriving $T_{ev}$ for PBH
domination (cf. Eq.II.11 for radiation domination) and using Eq.II.16 as well
as Eq.II.22, one arrives at

$\displaystyle Y_{BH}^{\rm
B}(T_{ev})=\frac{15}{4\pi^{3}g_{*s}(T_{ev})}\left(\frac{\pi^{3}g_{*}(T_{ev})}{5}\right)^{3/4}\left(\frac{\mathcal{G}g_{*B}(T_{BH})}{10240\pi}\right)^{1/2}\left(\frac{M_{Pl}}{M_{BH}}\right)^{5/2}.$
(III.9)

In figure 3, we show the allowed parameter space for successful leptogenesis
for radiation (upper panel) and black hole domination (lower panel). In the
plots, the green, blue and the red regions are excluded due to the constraints
derived in Eq.III.2, Eq.III.3 and Eq.III.5. For radiation domination, the
white region is allowed when we consider $\beta=10^{-1}\beta_{\rm max}$ (cf.
Eq.II.17) while the gray region is allowed for $\beta=10^{-2}\beta_{\rm max}$.
We note that the allowed parameter space vanishes for $\beta=10^{-3}\beta_{\rm
max}$. As mentioned earlier, since we are mostly interested in PBH domination,
we discuss the parameter space in a bit more detail. The solid black line and
the dashed line represent the PBH formation temperature and the Hawking
temperature. Therefore, the $U(1)_{B-L}$ breaking scale always lies above the
solid black line. On the other hand, comparing Eq.III.2 and Eq.III.3, one
obtains the upper bound on $M_{BH}$ as

$\displaystyle M_{BH}^{\rm max}=\left[2304\left(\frac{Y_{B}^{\rm Obs}\pi
v^{2}}{m_{3}}\right)^{2}\left(\frac{g_{*B}}{g_{R}}\right)^{2}\frac{1}{(\bar{Y}^{B}_{BH}(T_{ev}))^{2}}\right]^{-\frac{1}{3}}\simeq
13{\rm g}~{}(m_{3}=0.05{\rm eV}),$ (III.10)

where $\bar{Y}^{B}_{BH}(T_{ev})=Y^{B}_{BH}(T_{ev})M_{BH}^{5/2}$. Consequently,
$\Lambda_{CS}^{\rm min}$ is calculated as

$\displaystyle\Lambda_{CS}^{\rm
min}=\left(\frac{45\gamma^{2}}{16\pi^{3}g_{*}(T_{Bf})}\right)^{1/4}\left(\frac{M_{Pl}}{M_{BH}^{\rm
max}}\right)^{1/2}M_{Pl}\simeq 1.18\times 10^{15}{\rm GeV}$ (III.11)

which can further be lowered by considering the full range of neutrino
$3\sigma$ oscillation data. For smaller $M_{BH}$, one gets $\Lambda_{CS}^{\rm
min}$ until the maximum reheat temperature $T_{RH}\sim 10^{16}$ GeV is
reached. We take the value of $\Lambda_{CS}^{\rm min}$ in Eq.III.11 and
compute the GW spectra in the next section. We conclude this section with the
following remarks: In the computation of GW spectra from cosmic strings, we do
not consider possible formation of black hole-string networks which for PBHs
of solar mass level with $\Lambda_{CS}\gtrsim T_{Bf}$, may leave observable
consequences on the GW spectra[127]. Though we are primarily dealing with the
non-thermal production, one would like also to understand the overall status
of thermal leptogenesis in presence of PBHs. Firstly, as we see from the
figure 3, that the non-thermal contribution becomes prominent for very light
black holes and with $M_{1}\gtrsim 10^{11}$ GeV, whereas for the RH masses,
say $M\in[10^{9},10^{10}]$, the non-thermal contribution is sub-dominant.
However, for $M\in[10^{9},10^{10}]$, correct baryon asymmetry can be generated
by the standard thermal leptogenesis, since in that case, the light PBHs
already evaporate before the asymmetry production from the thermal bath stops.
On the left side of figure 4, we show the dynamics of the asymmetry production
(standard hierarchical $N_{1}$-dominated scenario) for $M_{BH}=1$g with
$M_{1}=4\times 10^{12}$ GeV and $M_{1}=4\times 10^{9}$ GeV. It is clear that
for the latter case, the effect of PBH on the final asymmetry is negligible
while for the former case, the thermally produced asymmetry is washed out. On
the other hand, heavy PBHs evaporate much later in time and therefore they
erase the thermally produced asymmetry unless the typical RH mass scale is
below the evaporation temperature. On the right-hand side, we show the
dynamics for $M_{BH}=100$g with $M_{1}=4\times 10^{11}$ GeV and $M_{1}=4\times
10^{8}$ GeV. This shows that for heavier PBHs, one has to go beyond the
hierarchical scenario for a successful thermal leptogenesis, given the fact
that for hierarchical leptogenesis, $M_{1}\gtrsim 10^{9}$ GeV (without flavour
effect)[128] and $M_{1}\gtrsim 10^{8}$ GeV (with flavour effect)[122]. A more
precise numerical relation between $M_{BH}$ and the RH scale can be found in
Ref.[48] (which we found in the middle stage of this work) which solves
density matrix equations for the lepton asymmetry matrix (treating flavour
effects rigorously).

##  IV Gravitational waves from cosmic string and phenomenological discussion

Figure 4: Left: Evolution of the $B-L$ asymmetry for $M_{BH}=1$g,
$\Delta=10^{-3}(M_{1}\simeq 4\times 10^{12}{\rm GeV})$ and
$\Delta=10^{-6}(M_{1}\simeq 4\times 10^{9}{\rm GeV})$. Right: Evolution of the
$B-L$ asymmetry for $M_{BH}=100$g, $\Delta=10^{-3}(M_{1}\simeq 4\times
10^{11}{\rm GeV})$ and $\Delta=10^{-6}(M_{1}\simeq 4\times 10^{8}{\rm GeV})$.

Now we proceed to the computation of GW spectrum originating from cosmic
string network. Cosmic strings which are nearly one-dimensional objects arise
generically with topologically-stable field configurations in BSM theories
featuring spontaneously broken U(1) symmetry [84, 129]. They are considered to
be one-dimensional objects characterised macroscopically with their string
tension $\mu$ which is typically taken to be of the order of the square of the
symmetry breaking scale $\Lambda_{CS}$. Normalised string tension $G\mu$, with
$G$ being the Newton’s constant, is directly constrained by CMB observation as
$G\mu\lesssim 1.1\times 10^{-7}$[130]. After formation, cosmic strings are
expected to reach a scaling regime in which their net energy density tracks
the total energy density of the universe with a relative fraction $G\mu$. This
regime consists of many closed loops as well as Hubble-length long strings
which intersect to form new loops as the universe expands. All these loops
oscillate and emit radiation, including GWs. We consider stochastic
gravitational background (SGWB) from cosmic string network by considering
Nambu-Goto strings which radiate energy dominantly in the form of GW
radiation[86, 87]. We follow Ref.[94] to calculate SGWB from cosmic string.
Once the loops are formed, they radiate energy in the form of gravitational
radiation at a constant rate as

$\displaystyle\frac{dE}{dt}=-\Gamma G\mu^{2},$ (IV.1)

where $\Gamma=50$[131, 132]. Therefore, the initial length $l_{i}=\alpha
t_{i}$ of the loop decreases as

$\displaystyle l(t)=\alpha t_{i}-\Gamma G\mu(t-t_{i})$ (IV.2)

until it disappears completely. The loop size $\alpha$ has a distribution and
for the largest loop one typically has $\alpha=0.1$[133, 134]. We consider
$\alpha=0.1$ in the numerical calculation and also discuss the effects of
smaller loops in the GW spectrum. The total energy loss from a loop is
decomposed into a set of normal-mode oscillations at frequencies
$\tilde{f}=2k/l$, where $k=1,2,3..$. The GW density parameter is given by

$\displaystyle\Omega_{GW}=\frac{f}{\rho_{c}}\frac{d\rho_{GW}}{df},$ (IV.3)

where $f$ is the red-shifted frequency and $\rho_{c}=3H_{0}^{2}/8\pi G$. The
density parameter $\Omega_{GW}$ can be written as a sum over all relic
densities corresponding to a mode $k$ as

$\displaystyle\Omega_{GW}(f)=\sum_{k}\Omega_{GW}^{(k)}(f),$ (IV.4)

where 222Please see Appendix:A for a derivation using the loop number density
obtained from the Velocity-dependent One-Scale (VOS) model[135, 136, 95].

$\displaystyle\Omega_{GW}^{(k)}(f)=\frac{1}{\rho_{c}}\frac{2k}{f}\frac{\mathcal{F}_{\alpha}\Gamma^{(k)}G\mu^{2}}{\alpha(\alpha+\Gamma
G\mu)}\int_{t_{F}}^{t_{0}}d\tilde{t}\frac{C_{\rm
eff}(t_{i}^{(k)})}{t_{i}^{(k)4}}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{5}\left[\frac{a(t_{i}^{k})}{a(\tilde{t})}\right]^{3}\Theta(t_{i}^{(k)}-t_{F})$
(IV.5)

and the integration runs over the emission time with $t_{F}$ as the time when
the network reach scaling regime shortly after formation. The numerical values
of $C_{\rm eff}$ are found to be 5.4 and 0.5 at radiation and matter
domination and $\mathcal{F}_{\alpha}$ has a value $\sim 0.1$[133, 134]. The
quantity $t_{i}^{(k)}$ is the formation time of the loops contributing to the
mode $k$ and is given by

$\displaystyle t_{i}^{(k)}(\tilde{t},f)=\frac{1}{\alpha+\Gamma
G\mu}\left[\frac{2k}{f}\frac{a(\tilde{t})}{a(t_{0})}+\Gamma
G\mu\tilde{t}\right].$ (IV.6)

Assuming GW emission is dominated by the cusps like structures, the relative
emission rate per mode is given by

$\displaystyle\Gamma^{(k)}=\frac{\Gamma
k^{-4/3}}{\sum_{m=1}^{\infty}m^{-4/3}}$ (IV.7)

with $\zeta(4/3)=\sum_{m=1}^{\infty}m^{-4/3}\simeq 3.6$ and
$\sum_{k}\Gamma^{(k)}=\Gamma$. On the left side of figure 5, we show the GW
spectrum for $\Lambda_{CS}^{\rm min}$ with loop size $\alpha=0.1$ and
$k_{max}\simeq 10^{3}$. As one sees, that even for the $\Lambda_{CS}^{\rm
min}$, the mechanism predicts GW with a strong amplitude with a flat plateau
which will be fully tested in the space-based interferometers such as Taiji
[138], TianQin[139], LISA[137], BBO[140], DECIGO[141], ground-based
interferometers like Einstein Telescope (ET)[142] and Cosmic Explorer
(CE)[143], and atomic interferometers MAGIS[144], AEDGE[145].

Figure 5: Top left: GWs from cosmic strings for $M_{BH}^{\rm max}$ and
$\Lambda_{CS}^{\rm min}$ for $\alpha=0.1$. Top right: GWs from cosmic strings
for $\Lambda_{CS}^{\rm min}$ and $\alpha=0.1$ (red dashed) and graviton
emission for $M_{BH}^{\rm max}$ (blue dashed). Bottom left: Effect of
including larger number modes on the GW spectrum beyond the turning point
frequency. Bottom right: GW spectra for smaller initial loop sizes
($\alpha=10^{-3},10^{-4}$ and $10^{-5}$) considering the analytic expression
obtained from the VOS model.

Note that the amplitude-frequency relation in Eq.IV.5 is explicitly dependent
on the cosmological history of the universe. Therefore, any non-standard
cosmological scenario (we have a PBH dominated early universe) would leave its
imprints on the GW spectra. Any observable break in the GW spectra due to this
non-standard expansion history could be a strong probe of PBH baryogenesis,
given that the duration of the non-standard expansion is fixed by the
requirement to obtain correct baryon asymmetry. Given $T_{\Delta}$ is the
temperature at which the most recent radiation era begins (before BBN one
should have a radiation domination), the frequency (the turning point
frequency, see Appendix:B) at which the flat plateau (see top panel of figure
5) breaks can be estimated as[94]

$\displaystyle f_{\Delta}=\sqrt{\frac{8}{z_{\rm eq}\alpha\Gamma
G\mu}}\left[\frac{g_{*}(T_{\Delta})}{g_{*}(T_{0})}\right]^{1/4}\left(\frac{T_{\Delta}}{T_{0}}\right)t_{0}^{-1},$
(IV.8)

where the red-shift at the standard matter radiation equality is given by
$z_{\rm eq}\simeq 3387$ and $T_{0}=2.7$K is the present temperature. For an
early matter era, the typical amplitude of the GW for the fundamental mode
($k=1$) goes with frequency as

$\displaystyle\Omega_{\rm GW}(f)\simeq\Omega_{\rm GW}^{\rm
plt}\left(\frac{f}{f_{\Delta}}\right)^{\kappa},$ (IV.9)

where $\kappa\simeq 0$ for $f\lesssim f_{\Delta}$ and $\kappa\simeq-1$ for
$f\gtrsim f_{\Delta}$. In the case of PBH baryogenesis $T_{\Delta}\simeq
T_{ev}$. Now using Eq.II.11, Eq.III.10 and Eq.IV.8, the frequency $f_{\Delta}$
can be calculated as $f_{\Delta}\simeq 3\times 10^{6}$ Hz, for $\alpha=0.1$.
Therefore, to observe this break one needs GW detectors at MHz frequencies. In
addition, let us mention that beyond the turning point frequency, the fall of
the spectrum as $\Omega_{GW}\propto f^{-1}$ is true only for the fundamental
or first few modes[146, 147, 91]. Once larger number of modes are included,
the spectrum falls as $f^{-1/3}$ (see Appendix:C for a derivation with a
cartoon diagram)–a feature that can be observed in the bottom-left panel of
Fig.5. There could be another gravitational wave background via graviton
emission by the PBHs[148, 149, 150]. However, for light PBHs as in the present
scenario, GW background via graviton emission typically peaks at frequencies
more than $10^{12}$ Hz. An analytical relation between energy density and
frequency is given by[148]

$\displaystyle\Omega_{GW}^{GR}(f)=\frac{128\pi^{4}f^{4}}{\rho_{c}M_{Pl}^{4}}\int_{t_{Bf}}^{t_{ev}}\frac{M(t)^{2}n_{PBH}(t)}{\exp\left[16\pi^{2}fM(t)a(t_{0})/a(t)M_{Pl}^{2}\right]-1}dt.$
(IV.10)

In figure 5 (top), we have shown the full spectrum of gravitational waves
relevant to probe PBH baryogenesis. In the left panel we have shown the GW
spectrum from CS only, for $\alpha=0.1$, whereas in the right panel, we show
the GW spectrum from CS for $\alpha=0.1$ (red dashed) as well as due to the
graviton emission (blue dashed) for $M_{BH}=13$g. The thick solid red line on
the right hand side represents the combined signal that would manifest a
successful baryogenesis. There have been proposals to detect GWs at higher
frequencies[151, 152, 153] which may shed light on the behaviour of GW
spectrum at higher frequencies compatible with successful PBH baryogenesis.
Since we discuss GWs with strong amplitude, probably it is worthwhile to have
some comments on the prospect of PBH baryogenesis in the light of recent
finding by the NANOGrav collaboration[107]. The NANOGrav collaboration with
their recently released 12.5 yrs data set, reports strong evidence for a
stochastic common-spectrum process over independent-red noises across 45
pulsars [107]. However, the detection has not been claimed as GWs, since the
time residuals do not show characteristic spatial relation which is described
by the Hellings–Downs (HD) curve[154]. On top of that, other systematics such
as pulsar spin noise[155] and solar system effects[156] may affect the signal
and therefore, the analysis requires proper handling of these two
effects[107]. Nonetheless, if the signal is interpreted as GWs, it would
unequivocally open up a new direction to probe Early Universe cosmology. Note
that, the new NANOGrav 12.5 yr data[107] is nearly consistent with previous
EPTA data[157], but they are in tension with previous limits from PPTA [158]
and a previous NANOGrav analysis of their 11 yrs (older) data[159]. According
to Ref.[107], this tension would be reduced using improved prior to the
intrinsic pulsar red noise to the older data. The 12.5 yrs NANOGrav data are
expressed in terms of power-law signal with characteristic strain given by

$\displaystyle h_{c}(f)=A\left(\frac{f}{f_{yr}}\right)^{(3-\gamma)/2}$ (IV.11)

with $f_{yr}=1yr^{-1}$ and $A$ being the characteristic strain amplitude. The
abundance of GWs has the standard form and can be recast as:

$\displaystyle\Omega_{GW}(f)=\Omega_{yr}\left(\frac{f}{f_{yr}}\right)^{5-\gamma},~{}~{}~{}~{}~{}{\rm
with}~{}~{}~{}~{}~{}\Omega_{yr}=\frac{2\pi^{2}}{3H_{0}^{2}}A^{2}f^{2}_{yr}.$
(IV.12)

Given our predicted amplitude for $\alpha=0.1$ within the NANOGrav frequency
range, we find a fitting point (see Ref.[108, 92] for a procedure) on the
spectral index ($\gamma$)-amplitude ($A$) plane well outside the $2\sigma$
contour which also evident from the Fig.5. It can be seen on the bottom panel
in Fig.5, that the GW amplitude could be consistent with the NANOGrav data for
smaller loop sizes (see also Fig.6 for a fit). As mentioned earlier, recent
numerical simulations of cosmic string networks report a population of large
loops that effectively contribute to the GWs, with initial loop size parameter
peaked around $\alpha\simeq 0.1$[133, 134]. Around 10$\%$ of the energy
released by the long string network goes to the large loops, whereas the
remaining $90\%$ goes to the kinetic energy of highly boosted smaller loops.
The kinetic energy simply redshifts away and is not transferred to GWs. This
fact is taken into account in VOS model by the normalization factor
$\mathcal{F}_{\alpha}=0.1$ which for large loops ($\sim\alpha\simeq 0.1$),
provides a good description of the results obtained from numerical simulation
(see Appendix A). To compute GW spectrum for smaller loops, we use the loop
number density obtained analytically from VOS model without including the
normalisation factor $\mathcal{F}_{\alpha}=0.1$ (see e.g., Ref.[95]).

Figure 6: A fit to the NANOGrav data for small loops ($\alpha<0.1$ spectrum)
and $\mathcal{F}_{\alpha}=1$. Three benchmark values of $\alpha$ are shown in
red ($\alpha=10^{-3}$), purple ($\alpha=10^{-4}$) and green
($\alpha=10^{-5}$).

##  V Summary

In this work, we have studied baryogenesis via non-thermal leptogenesis from
ultralight primordial black hole evaporation and its test with cosmic string
induced gravitational waves. We show that if the right handed neutrino masses
are generated dynamically which is then followed by black hole formation,
successful baryogenesis gives rise to GWs with strong amplitude which will be
tested within a wide range of frequencies by the future GW detectors. We focus
more on the scenario where the black holes dominate the energy density of the
universe before they evaporate. In this case, the final baryon asymmetry is
independent of black hole density at the formation and there is a break in the
GW spectra which occurs around MHz frequency. Therefore, scenarios like PBH
baryogenesis could be a good motivation to think of GW detectors beyond the
frequency sensitivity of LIGO. We also consider the recent finding by the
NANOGrav PTA of a stochastic common spectrum process across many pulsars which
could be a GW signal. The predicted GW amplitude by successful PBH
baryogenesis is in tension with the reported cross-power spectral density and
characteristic strain amplitude contours for large string loops as predicted
by numerical simulations of cosmic string network. However, an agnostic
analytical consideration of GWs produced by smaller loops can explain NANOGrav
data.

## acknowledgement

R. Samanta is supported by the project MSCA-IF IV FZU -
CZ.02.2.69/0.0/0.0/$20\\_079$/0017754 and acknowledges European Structural and
Investment Fund and the Czech Ministry of Education, Youth and Sports.

##  Appendix A Derivation of the Gravitational waves spectrum

The normalised energy density of gravitational waves at present time is
expressed as

$\displaystyle\Omega_{GW}(t_{0},f)=\frac{f}{\rho_{c}}\frac{d\rho_{GW}}{df}=\sum_{k}\Omega_{GW}^{(k)}(t_{0},f).$
(A.1)

The frequency derivative is given by[134]

$\displaystyle\frac{d\rho_{GW}^{(k)}}{df}=\int_{t_{F}}^{t_{0}}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{4}P_{GW}(\tilde{t},f_{k})\frac{dF}{df}d\tilde{t},$
(A.2)

where the factor $\frac{dF}{df}=f\left[\frac{a(t_{0})}{a(\tilde{t})}\right]$
is introduced to take into account the red-shifting of the frequency and the
quantity $P_{GW}(\tilde{t},f_{k})$ represents the power emitted by the loops
and is given by

$\displaystyle P_{GW}(\tilde{t},f_{k})=G\mu^{2}\Gamma_{k}\int
n(l,\tilde{t})\delta\left(f_{k}-\frac{2k}{l}\right)dl.$ (A.3)

Integrating Eq.A.3 over the loop lengths gives

$\displaystyle
P_{GW}(\tilde{t},f_{k})=\frac{2kG\mu^{2}\Gamma_{k}}{f_{k}^{2}}n(\tilde{t},f_{k})=\frac{2kG\mu^{2}\Gamma_{k}}{f^{2}\left[\frac{a(t_{0})}{a(\tilde{t})}\right]^{2}}n\left(\tilde{t},\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)$
(A.4)

From Eq.A.4 and Eq.A.2 one gets

$\displaystyle\frac{d\rho_{GW}^{(k)}}{df}=\frac{2kG\mu^{2}\Gamma_{k}}{f^{2}}\int_{t_{F}}^{t_{0}}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{5}n\left(\tilde{t},\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)d\tilde{t}$
(A.5)

Therefore the energy density corresponding to the mode ‘$k$’ is given by

$\displaystyle\Omega_{GW}^{(k)}(t_{0},f)=\frac{2kG\mu^{2}\Gamma_{k}}{f\rho_{c}}\int_{t_{F}}^{t_{0}}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{5}n\left(\tilde{t},\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)d\tilde{t}.$
(A.6)

The most important aspect to obtain gravitational wave spectrum from cosmic
string loops is to compute the number density
$n\left(\tilde{t},l_{k}(\tilde{t})\equiv\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)$.
In this regard, two different approaches can be considered. The first approach
could be obtaining the number density directly from numerical solution[133]
which reads (considering the loops created during radiation domination,
similar argument also holds true for matter domination)

$\displaystyle
n(\tilde{t},l_{k}(\tilde{t}))=\frac{0.18}{\left[l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}\right]^{5/2}\tilde{t}^{3/2}}.$ (A.7)

On the other hand one may opt for an analytic approach by considering the
Velocity dependent One Scale (VOS) model which gives

$\displaystyle
n(\tilde{t},l_{k}(\tilde{t}))=\frac{A_{r}}{\alpha}\frac{(\alpha+\Gamma
G\mu)^{3/2}}{\left[l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}\right]^{5/2}\tilde{t}^{3/2}}\equiv\frac{A_{r}N_{\alpha}}{\left[l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}\right]^{5/2}\tilde{t}^{3/2}},$ (A.8)

where $A_{r}=5.4$[95]. Note that the VOS model assumes all the loops are of
same length at creation. However, at the moment of production, the loops may
follow a distribution depending on $\alpha$. In that case the above formula
should be modified as

$\displaystyle n(\tilde{t},l_{k}(\tilde{t}))=\frac{A_{r}\int
w(\alpha)N_{\alpha}d\alpha}{\left[l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}\right]^{5/2}\tilde{t}^{3/2}}.$ (A.9)

Therefore, to use the VOS formula in Eq.A.8 and at the same time to be
consistent with the numerical result, one has to normalise Eq.A.8, i.e.,

$\displaystyle
n(\tilde{t},l_{k}(\tilde{t}))=\frac{\mathcal{F}_{\alpha}A_{r}N_{\alpha}}{\left[l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}\right]^{5/2}\tilde{t}^{3/2}},~{}~{}{\rm
with}~{}~{}\mathcal{F}_{\alpha}=N_{\alpha}^{-1}\int
w(\alpha)N_{\alpha}d\alpha$ (A.10)

As one can see that for $\alpha=0.1$, Eq.A.7 and Eq.A.8 is consistent for
$\mathcal{F}_{\alpha}\sim 0.18/(A_{r}\sqrt{\alpha})\sim 0.1$ which we consider
in the numerical computation. The most general formula for the number density
in an expanding background that scales as $a\sim t^{\beta}$ is given by

$\displaystyle
n(\tilde{t},l_{k}(\tilde{t}))=\frac{A_{\beta}}{\alpha}\frac{(\alpha+\Gamma
G\mu)^{3(1-\beta)}}{\left[l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}\right]^{4-3\beta}\tilde{t}^{3\beta}}.$ (A.11)

Now having initial time $t_{i}^{(k)}$ as

$\displaystyle t_{i}^{(k)}=\frac{l_{k}(\tilde{t})+\Gamma
G\mu\tilde{t}}{\alpha+\Gamma G\mu},$ (A.12)

Eq.A.11 can be re-expressed as

$\displaystyle
n(\tilde{t},l_{k}(\tilde{t}))=\frac{A_{\beta}(t_{i}^{(k)})}{\alpha(\alpha+\Gamma
G\mu)t_{i}^{(k)4}}\left[\frac{a(t_{i}^{k})}{a(\tilde{t})}\right]^{3}.$ (A.13)

Putting the value of $n(\tilde{t},l_{k}(\tilde{t}))$ from Eq.A.13 into Eq.A.6,
one gets

$\displaystyle\Omega_{GW}^{(k)}(t_{0},f)=\frac{2kG\mu^{2}\Gamma_{k}}{f\rho_{c}\alpha(\alpha+\Gamma
G\mu)}\int_{t_{F}}^{t_{0}}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{5}\frac{C_{\rm
eff}(t_{i}^{(k)})}{t_{i}^{(k)4}}\left[\frac{a(t_{i}^{k})}{a(\tilde{t})}\right]^{3}\Theta(t_{i}^{(k)}-t_{F})d\tilde{t},$
(A.14)

where we have renamed $A_{\beta}$ as $C_{\rm eff}$ and the quantity
$\Gamma^{(k)}$ is given by

$\displaystyle\Gamma^{(k)}=\frac{\Gamma
k^{-\delta}}{\sum_{m=1}^{\infty}m^{-\delta}},~{}~{}\delta=\begin{cases}4/3~{}~{}{\rm
cusps}\\\ 5/3~{}~{}{\rm kinks}.\end{cases}$ (A.15)

##  Appendix B Derivation of the flat plateau and the turning point frequency

The expression for the fundamental mode in Eq.A.6 can be re-written as

$\displaystyle\Omega_{GW}^{(k)}(t_{0},f)=\frac{2G\mu^{2}\Gamma}{f\zeta(\delta)\rho_{c}}\int_{t_{F}}^{t_{0}}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{5}n\left(\tilde{t},\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)d\tilde{t}.$
(B.1)

It is convenient to do the integration with respect to the scale factor.
Eq.B.1 then becomes

$\displaystyle\Omega_{GW}^{(k)}(t_{0},f)=\frac{16\pi}{3\zeta(\delta)}\left(\frac{G\mu}{H_{0}}\right)^{2}\frac{\Gamma}{fa(t_{0})}\int_{a_{F}}^{a_{eq}}H(a)^{-1}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]^{4}n\left(\tilde{t},\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)da,$
(B.2)

where

$\displaystyle
H=H_{0}\Omega_{r}^{1/2}\left(\frac{a(\tilde{t})}{a(t_{0})}\right)^{-2}~{}~{}{\rm
with}~{}~{}\Omega_{r}\simeq 9\times 10^{-5}.$ (B.3)

The number density
$n\left(\tilde{t},l_{k}(\tilde{t})\equiv\frac{2k}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]\right)$
in Eq.A.8 can also be expressed in terms of the scale factor and is given by

$\displaystyle
n(\tilde{t},l_{k}(\tilde{t}))=\frac{A_{r}}{\alpha}\frac{(\alpha+\Gamma
G\mu)^{3/2}}{\left[\frac{2}{f}\left[\frac{a(\tilde{t})}{a(t_{0})}\right]+\Gamma
G\mu/2H\right]^{5/2}(2H)^{-3/2}.}$ (B.4)

Putting Eq.B.4 in Eq.B.2 and after performing the integration one gets

$\displaystyle\Omega_{GW}^{(1)}(f)=\frac{128\pi
G\mu}{9\zeta(\delta)}\frac{A_{r}}{\epsilon_{r}}\Omega_{r}(1+\epsilon_{r})^{3/2}\left[\left(\frac{f}{f+\epsilon_{r}f_{\rm
min}\left(\frac{t_{i}}{t_{\rm
eq}}\right)^{1/2}}\right)^{3/2}-\left(\frac{f}{f+\epsilon_{r}f_{\rm
min}}\right)^{3/2}\right],$ (B.5)

where we define $\epsilon_{r}=\alpha/\Gamma G\mu$ the $f_{min}=\frac{2}{\alpha
t_{i}}\frac{a_{i}}{a_{0}}=\frac{4H_{0}\Omega_{r}^{1/2}}{\alpha}\frac{a_{0}}{a_{i}}$
is the minimum frequency emitted by a given loop. Now given the scaling
solution of the loop production rate, which decreases with the fourth power in
time, $f\simeq f_{\rm min}$ is a reasonable assumption. In that case, with
$t_{i}\ll t_{\rm eq}$ one has

$\displaystyle\Omega_{GW}^{(1)}(f)=\frac{128\pi
G\mu}{9\zeta(\delta)}\frac{A_{r}}{\epsilon_{r}}\Omega_{r}\left[(1+\epsilon_{r})^{3/2}-1\right].$
(B.6)

Note that the expression for the flat plateau matches with Ref.[95] apart from
the factor $\zeta(\delta)$ in the denominator. This is simply due to the fact
that our definition of $\Omega_{GW}^{(k)}(f)$ in Eq.A.14 is inclusive of
$\Gamma^{k}$. In the above, we have assumed that the dominant emission comes
from the very earliest epoch of loop creation. In practice, a precise value of
the time can be computed by maximizing the integral in Eq.B.1 with respect to
$\tilde{t}$. This gives

$\displaystyle\tilde{t}_{M}\simeq\frac{2}{f\Gamma
G\mu}\frac{a_{M}}{a_{0}}\equiv\frac{1}{2\Gamma
G\mu}\frac{4a_{M}}{fa_{0}}\equiv\frac{l_{i}}{2\Gamma G\mu}$ (B.7)

where $f$ is the frequency observed today that was emitted at time $t_{M}$
when the a given loop $l_{i}=\alpha t_{i}$ reached to the half of its size
$l_{i}/2$, i.e., $t_{M}$ is eventually the half life of a given loop. Given
the time $t_{\Delta}$ at which the most recent radiation domination begins,
the approximate frequency up to which the spectrum shows a flat plateau is
given by

$\displaystyle f_{\Delta}=\sqrt{\frac{8}{\Gamma
G\mu}}t_{\Delta}^{-1/2}t_{0}^{-2/3}t_{\rm eq}^{1/6}\simeq\sqrt{\frac{8z_{\rm
eq}}{\Gamma G\mu}}\left(\frac{t_{\rm eq}}{t_{\Delta}}\right)^{1/2}t_{0}^{-1}.$
(B.8)

Eq.B.8 when expressed in terms of temperature gives Eq.IV.8.

##  Appendix C Effect of large-k modes on the spectral shape beyond the
turning point

Figure 7: A cartoon diagram showing how the contributions from large-k modes
are accounted for to analyse spectral shape beyond the turning point.

As pointed out in sec.IV, when the number of modes increases in the sum, the
spectral behaviour beyond the turning point deviates from that of the
fundamental mode[146, 147, 91]. The reason being the following:

First of all, from Eq.A.14 it is evident that

$\displaystyle\Omega_{GW}(f)=\sum_{k}\Omega_{GW}^{(k)}(f)=\sum_{k}k^{-\delta}\Omega^{(1)}(f/k).$
(C.1)

Now to perform the sum or in other words to compute the GW amplitude at a
given frequency (called test point in the cartoon diagram (Fig.7)) let us
expand the RHS of Eq.C.1 for some first few benchmark modes, i.e.,

$\displaystyle\Omega_{GW}(f)$ $\displaystyle=$
$\displaystyle\sum_{k}k^{-\delta}\Omega^{(1)}(f/k)$ (C.2) $\displaystyle=$
$\displaystyle
1^{-\delta}\Omega^{(1)}(f/1)+m^{-\delta}\Omega^{(1)}(f/m)+n^{-\delta}\Omega^{(1)}(f/n)+r^{-\delta}\Omega^{(1)}(f/r)+...,$

where the integers obey $1<m<n<r$. This suggests, if we keep on increasing the
mode numbers, there will be critical value $k\equiv k_{\Delta}$ for which the
amplitude $\Omega_{GW}^{(1)}(f_{\Delta}=f/k_{\Delta})$ contributes to the test
point. Therefore we can split the sum into two parts. The first one is from
$k=1$ up to $k_{\Delta}$ for which the test point receives contributions from
the non-flat part of the fundamental spectrum and the second one is from
$k_{\Delta}$ to $k_{max}(k_{\Delta}<\bar{r}<\bar{n}<\bar{m}<k_{max})$ for
which the test point receives contribution from the flat part of the
fundamental spectrum, i.e.,

$\displaystyle\Omega_{GW}(f)$ $\displaystyle=$
$\displaystyle\sum_{k=1}^{k=k_{\Delta}}k^{-\delta}\Omega_{GW}^{(1)}(f/k>f_{\Delta})+\sum_{k=k_{\Delta}}^{k=k_{max}}k^{-\delta}\Omega_{GW}^{(1)}(f/k<f_{\Delta})$
(C.3) $\displaystyle=$
$\displaystyle\sum_{k=1}^{k=k_{\Delta}}k^{-\delta}\Omega_{GW}^{\rm
plt}\left(\frac{f_{\Delta}}{f/k}\right)+\sum_{k=k_{\Delta}}^{k=k_{max}}k^{-\delta}\Omega_{GW}^{\rm
plt}.$ (C.4)

This sum can easily be evaluated in the large $k_{\Delta}$ limit using
asymptotic expansion of Euler-Maclaurin series for the first term and Hurwitz
zeta function for the second term. The final result is (up to some pre-factor
depending on $\delta$)

$\displaystyle\Omega_{GW}(f)\simeq\Omega_{GW}^{\rm
plt}\left(\frac{f_{\Delta}}{f}\right)^{\delta-1}~{}~{}{\rm
i.e.,}~{}~{}\Omega_{GW}(f)\propto\begin{cases}f^{-1/3}~{}~{}{\rm cusps}\\\
f^{-2/3}~{}~{}{\rm kinks}.\end{cases}$ (C.5)

As an aside, let us mention that in the scenarios like kination[160, 161]
where the equation of state is stiffer ($w=1$), beyond the turning point, the
energy density goes as $f^{1}$[94]. In that case, the first term (swapping
numerator and denominator) in Eq.C.4 gives the dominant contribution. In the
large $k_{\Delta}$ limit, the sum is therefore

$\displaystyle\Omega_{GW}(f)\simeq\Omega_{GW}^{\rm
plt}\left(\frac{f}{f_{\Delta}}\right)\sum_{k=1}^{k=k_{\Delta}}k^{-(\delta+1)}\simeq\Omega_{GW}^{\rm
plt}\left(\frac{f}{f_{\Delta}}\right)\zeta(\delta+1).$ (C.6)

Thus for stiffer equation of states like kination, the spectral shape is quite
similar to the fundamental mode even after taking in account the contribution
from the higher k-modes.

## References

  * [1] S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * [2] B. J. Carr and S. W. Hawking, Mon. Not. Roy. Astron. Soc. 168, 399 (1974).
  * [3] B. J. Carr, Astrophys. J. 201, 1 (1975). doi:10.1086/153853
  * [4] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 6, 061102 (2016) doi:10.1103/PhysRevLett.116.061102 [arXiv:1602.03837 [gr-qc]].
  * [5] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 24, 241103 (2016) doi:10.1103/PhysRevLett.116.241103 [arXiv:1606.04855 [gr-qc]].
  * [6] B. P. Abbott et al. [LIGO Scientific and VIRGO Collaborations], Phys. Rev. Lett. 118, no. 22, 221101 (2017) Erratum: [Phys. Rev. Lett. 121, no. 12, 129901 (2018)] doi:10.1103/PhysRevLett.118.221101, 10.1103/PhysRevLett.121.129901 [arXiv:1706.01812 [gr-qc]].
  * [7] B. . P. .Abbott et al. [LIGO Scientific and Virgo Collaborations], Astrophys. J. Lett. 851, L35 (2017) doi:10.3847/2041-8213/aa9f0c [arXiv:1711.05578 [astro-ph.HE]].
  * [8] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 119, no. 14, 141101 (2017) doi:10.1103/PhysRevLett.119.141101 [arXiv:1709.09660 [gr-qc]].
  * [9] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 119 (2017) no.16, 161101 doi:10.1103/PhysRevLett.119.161101 [arXiv:1710.05832 [gr-qc]].
  * [10] R. Abbott et al. [LIGO Scientific and Virgo Collaborations], Astrophys. J. Lett. 896, no. 2, L44 (2020) doi:10.3847/2041-8213/ab960f [arXiv:2006.12611 [astro-ph.HE]].
  * [11] B. Carr, F. Kuhnel and M. Sandstad, Phys. Rev. D 94, no. 8, 083504 (2016) doi:10.1103/PhysRevD.94.083504 [arXiv:1607.06077 [astro-ph.CO]].
  * [12] S. Clesse and J. García-Bellido, Phys. Dark Univ. 15, 142 (2017) doi:10.1016/j.dark.2016.10.002 [arXiv:1603.05234 [astro-ph.CO]].
  * [13] Y. Ali-Haïmoud and M. Kamionkowski, Phys. Rev. D 95, no. 4, 043534 (2017) doi:10.1103/PhysRevD.95.043534 [arXiv:1612.05644 [astro-ph.CO]].
  * [14] H. Niikura et al., Nature Astron. 3, no. 6, 524 (2019) doi:10.1038/s41550-019-0723-1 [arXiv:1701.02151 [astro-ph.CO]].
  * [15] B. Carr, M. Raidal, T. Tenkanen, V. Vaskonen and H. Veermäe, Phys. Rev. D 96, no. 2, 023514 (2017) doi:10.1103/PhysRevD.96.023514 [arXiv:1705.05567 [astro-ph.CO]].
  * [16] M. Sasaki, T. Suyama, T. Tanaka and S. Yokoyama, Class. Quant. Grav. 35, no. 6, 063001 (2018) doi:10.1088/1361-6382/aaa7b4 [arXiv:1801.05235 [astro-ph.CO]].
  * [17] L. Barack et al., Class. Quant. Grav. 36, no. 14, 143001 (2019) doi:10.1088/1361-6382/ab0587 [arXiv:1806.05195 [gr-qc]].
  * [18] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, arXiv:2002.12778 [astro-ph.CO].
  * [19] A. M. Green and B. J. Kavanagh, arXiv:2007.10722 [astro-ph.CO].
  * [20] A. L. Miller, S. Clesse, F. De Lillo, G. Bruno, A. Depasse and A. Tanasijczuk, arXiv:2012.12983 [astro-ph.HE].
  * [21] R. g. Cai, S. Pi and M. Sasaki, Phys. Rev. Lett. 122, no.20, 201101 (2019) doi:10.1103/PhysRevLett.122.201101 [arXiv:1810.11000 [astro-ph.CO]].
  * [22] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975) Erratum: [Commun. Math. Phys. 46, 206 (1976)]. doi:10.1007/BF02345020, 10.1007/BF01608497
  * [23] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, arXiv:2002.12778 [astro-ph.CO].
  * [24] S. K. Acharya and R. Khatri, JCAP 2006, 018 (2020) doi:10.1088/1475-7516/2020/06/018 [arXiv:2002.00898 [astro-ph.CO]].
  * [25] C. Keith, D. Hooper, N. Blinov and S. D. McDermott, Phys. Rev. D 102, no. 10, 103512 (2020) doi:10.1103/PhysRevD.102.103512 [arXiv:2006.03608 [astro-ph.CO]].
  * [26] A. G. Polnarev and M. Y. Khlopov, Sov. Phys. Usp. 28, 213 (1985) [Usp. Fiz. Nauk 145, 369 (1985)]. doi:10.1070/PU1985v028n03ABEH003858
  * [27] M. Y. Khlopov, A. Barrau and J. Grain, Class. Quant. Grav. 23, 1875 (2006) doi:10.1088/0264-9381/23/6/004 [astro-ph/0406621].
  * [28] T. Fujita, M. Kawasaki, K. Harigaya and R. Matsuda, Phys. Rev. D 89, no. 10, 103501 (2014) doi:10.1103/PhysRevD.89.103501 [arXiv:1401.1909 [astro-ph.CO]].
  * [29] R. Allahverdi, J. Dent and J. Osinski, Phys. Rev. D 97, no. 5, 055013 (2018) doi:10.1103/PhysRevD.97.055013 [arXiv:1711.10511 [astro-ph.CO]].
  * [30] O. Lennon, J. March-Russell, R. Petrossian-Byrne and H. Tillim, JCAP 1804, 009 (2018) doi:10.1088/1475-7516/2018/04/009 [arXiv:1712.07664 [hep-ph]].
  * [31] L. Morrison, S. Profumo and Y. Yu, JCAP 1905, 005 (2019) doi:10.1088/1475-7516/2019/05/005 [arXiv:1812.10606 [astro-ph.CO]].
  * [32] D. Hooper, G. Krnjaic and S. D. McDermott, JHEP 1908, 001 (2019) doi:10.1007/JHEP08(2019)001 [arXiv:1905.01301 [hep-ph]].
  * [33] C. Lunardini and Y. F. Perez-Gonzalez, JCAP 2008, 014 (2020) doi:10.1088/1475-7516/2020/08/014 [arXiv:1910.07864 [hep-ph]].
  * [34] A. Chaudhuri and A. Dolgov, arXiv:2001.11219 [astro-ph.CO].
  * [35] I. Masina, Eur. Phys. J. Plus 135, no. 7, 552 (2020) doi:10.1140/epjp/s13360-020-00564-9 [arXiv:2004.04740 [hep-ph]].
  * [36] I. Baldes, Q. Decant, D. C. Hooper and L. Lopez-Honorez, JCAP 2008, 045 (2020) doi:10.1088/1475-7516/2020/08/045 [arXiv:2004.14773 [astro-ph.CO]].
  * [37] P. Gondolo, P. Sandick and B. Shams Es Haghi, Phys. Rev. D 102, no. 9, 095018 (2020) doi:10.1103/PhysRevD.102.095018 [arXiv:2009.02424 [hep-ph]].
  * [38] N. Bernal and Ó. Zapata, arXiv:2010.09725 [hep-ph].
  * [39] N. Bernal and Ó. Zapata, arXiv:2011.02510 [hep-ph].
  * [40] N. Bernal and Ó. Zapata, arXiv:2011.12306 [astro-ph.CO].
  * [41] J. Auffinger, I. Masina and G. Orlando, arXiv:2012.09867 [hep-ph].
  * [42] T. Kitabayashi, [arXiv:2101.01921 [hep-ph]].
  * [43] J. D. Barrow, E. J. Copeland, E. W. Kolb and A. R. Liddle, Phys. Rev. D 43, 984 (1991). doi:10.1103/PhysRevD.43.984
  * [44] D. Baumann, P. J. Steinhardt and N. Turok, hep-th/0703250 [HEP-TH].
  * [45] A. Hook, Phys. Rev. D 90, no. 8, 083535 (2014) doi:10.1103/PhysRevD.90.083535 [arXiv:1404.0113 [hep-ph]].
  * [46] Y. Hamada and S. Iso, PTEP 2017, no. 3, 033B02 (2017) doi:10.1093/ptep/ptx011 [arXiv:1610.02586 [hep-ph]].
  * [47] D. Hooper and G. Krnjaic, arXiv:2010.01134 [hep-ph].
  * [48] Y. F. Perez-Gonzalez and J. Turner, arXiv:2010.03565 [hep-ph].
  * [49] G. Aliferis and V. Zarikas, Phys. Rev. D 103, no.2, 023509 (2021) doi:10.1103/PhysRevD.103.023509 [arXiv:2006.13621 [gr-qc]].
  * [50] K. Inomata, M. Kawasaki, K. Mukaida, T. Terada and T. T. Yanagida, Phys. Rev. D 101, no. 12, 123533 (2020) doi:10.1103/PhysRevD.101.123533 [arXiv:2003.10455 [astro-ph.CO]].
  * [51] T. Papanikolaou, V. Vennin and D. Langlois, arXiv:2010.11573 [astro-ph.CO].
  * [52] G. Domènech, C. Lin and M. Sasaki, arXiv:2012.08151 [gr-qc].
  * [53] P. Burda, R. Gregory and I. Moss, Phys. Rev. Lett. 115, 071303 (2015) doi:10.1103/PhysRevLett.115.071303 [arXiv:1501.04937 [hep-th]].
  * [54] P. Burda, R. Gregory and I. Moss, JHEP 1606, 025 (2016) doi:10.1007/JHEP06(2016)025 [arXiv:1601.02152 [hep-th]].
  * [55] A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5, 32 (1967) [JETP Lett. 5, 24 (1967)] [Sov. Phys. Usp. 34, no. 5, 392 (1991)] [Usp. Fiz. Nauk 161, no. 5, 61 (1991)]. doi:10.1070/PU1991v034n05ABEH002497
  * [56] V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. 155B, 36 (1985). doi:10.1016/0370-2693(85)91028-7
  * [57] E. W. Kolb and M. S. Turner, Front. Phys. 69, 1 (1990).
  * [58] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45 (1986). doi:10.1016/0370-2693(86)91126-3
  * [59] A. Riotto and M. Trodden, Ann. Rev. Nucl. Part. Sci. 49, 35 (1999) doi:10.1146/annurev.nucl.49.1.35 [hep-ph/9901362].
  * [60] A. Pilaftsis and T. E. J. Underwood, Nucl. Phys. B 692, 303 (2004) doi:10.1016/j.nuclphysb.2004.05.029 [hep-ph/0309342].
  * [61] W. Buchmuller, P. Di Bari and M. Plumacher, Annals Phys. 315, 305 (2005) doi:10.1016/j.aop.2004.02.003 [hep-ph/0401240].
  * [62] S. Davidson, E. Nardi and Y. Nir, Phys. Rept. 466, 105 (2008) doi:10.1016/j.physrep.2008.06.002 [arXiv:0802.2962 [hep-ph]].
  * [63] D. Bodeker and W. Buchmuller, arXiv:2009.07294 [hep-ph].
  * [64] P. Minkowski, Phys. Lett. 67B, 421 (1977). doi:10.1016/0370-2693(77)90435-X
  * [65] M. Gell-Mann, P. Ramond and R. Slansky, Conf. Proc. C 790927, 315 (1979) [arXiv:1306.4669 [hep-th]].
  * [66] T. Yanagida, Prog. Theor. Phys. 64, 1103 (1980). doi:10.1143/PTP.64.1103
  * [67] R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44, 912 (1980). doi:10.1103/PhysRevLett.44.912
  * [68] A variant of Type-1 seesaw: R. N. Mohapatra, Phys. Rev. Lett. 56, 561 (1986). doi:10.1103/PhysRevLett.56.561
  * [69] E. K. Akhmedov, V. A. Rubakov and A. Y. Smirnov, Phys. Rev. Lett. 81, 1359 (1998) doi:10.1103/PhysRevLett.81.1359 [hep-ph/9803255].
  * [70] T. Hambye and D. Teresi, Phys. Rev. Lett. 117, no. 9, 091801 (2016) doi:10.1103/PhysRevLett.117.091801 [arXiv:1606.00017 [hep-ph]].
  * [71] P. S. Bhupal Dev, P. Millington, A. Pilaftsis and D. Teresi, Nucl. Phys. B 886 (2014) 569 doi:10.1016/j.nuclphysb.2014.06.020 [arXiv:1404.1003 [hep-ph]].
  * [72] P. S. Bhupal Dev, P. Millington, A. Pilaftsis and D. Teresi, Nucl. Phys. B 891, 128 (2015) doi:10.1016/j.nuclphysb.2014.12.003 [arXiv:1410.6434 [hep-ph]].
  * [73] G. Altarelli and F. Feruglio, Rev. Mod. Phys. 82, 2701 (2010) doi:10.1103/RevModPhys.82.2701 [arXiv:1002.0211 [hep-ph]].
  * [74] H. Ishimori, T. Kobayashi, H. Ohki, Y. Shimizu, H. Okada and M. Tanimoto, Prog. Theor. Phys. Suppl. 183, 1 (2010) doi:10.1143/PTPS.183.1 [arXiv:1003.3552 [hep-th]].
  * [75] S. F. King, Prog. Part. Nucl. Phys. 94, 217 (2017) doi:10.1016/j.ppnp.2017.01.003 [arXiv:1701.04413 [hep-ph]].
  * [76] S. T. Petcov, Eur. Phys. J. C 78, no. 9, 709 (2018) doi:10.1140/epjc/s10052-018-6158-5 [arXiv:1711.10806 [hep-ph]].
  * [77] P. Di Bari and A. Riotto, Phys. Lett. B 671, 462 (2009) doi:10.1016/j.physletb.2008.12.054 [arXiv:0809.2285 [hep-ph]].
  * [78] P. Di Bari and A. Riotto, JCAP 1104, 037 (2011) doi:10.1088/1475-7516/2011/04/037 [arXiv:1012.2343 [hep-ph]].
  * [79] C. S. Fong, D. Meloni, A. Meroni and E. Nardi, JHEP 1501, 111 (2015) doi:10.1007/JHEP01(2015)111 [arXiv:1412.4776 [hep-ph]].
  * [80] P. Di Bari, L. Marzola and M. Re Fiorentin, Nucl. Phys. B 893, 122 (2015) doi:10.1016/j.nuclphysb.2015.02.005 [arXiv:1411.5478 [hep-ph]].
  * [81] P. Di Bari and R. Samanta, arXiv:2005.03057 [hep-ph].
  * [82] J. A. Dror, T. Hiramatsu, K. Kohri, H. Murayama and G. White, Phys. Rev. Lett. 124, no. 4, 041804 (2020) doi:10.1103/PhysRevLett.124.041804 [arXiv:1908.03227 [hep-ph]].
  * [83] T. W. B. Kibble, J. Phys. A 9, 1387 (1976). doi:10.1088/0305-4470/9/8/029
  * [84] R. Jeannerot, J. Rocher and M. Sakellariadou, Phys. Rev. D 68, 103514 (2003) doi:10.1103/PhysRevD.68.103514 [hep-ph/0308134].
  * [85] H. B. Nielsen and P. Olesen, Nucl. Phys. B 61, 45 (1973). doi:10.1016/0550-3213(73)90350-7
  * [86] C. Ringeval, M. Sakellariadou and F. Bouchet, JCAP 0702, 023 (2007) doi:10.1088/1475-7516/2007/02/023 [astro-ph/0511646].
  * [87] J. J. Blanco-Pillado, K. D. Olum and B. Shlaer, Phys. Rev. D 83, 083514 (2011) doi:10.1103/PhysRevD.83.083514 [arXiv:1101.5173 [astro-ph.CO]].
  * [88] A. Davidson, Phys. Rev. D 20, 776 (1979). doi:10.1103/PhysRevD.20.776
  * [89] R. E. Marshak and R. N. Mohapatra, Phys. Lett. 91B, 222 (1980). doi:10.1016/0370-2693(80)90436-0
  * [90] R. N. Mohapatra and R. E. Marshak, Phys. Rev. Lett. 44, 1316 (1980) Erratum: [Phys. Rev. Lett. 44, 1643 (1980)]. doi:10.1103/PhysRevLett.44.1644.2, 10.1103/PhysRevLett.44.1316
  * [91] S. Blasi, V. Brdar and K. Schmitz, arXiv:2004.02889 [hep-ph].
  * [92] R. Samanta and S. Datta, arXiv:2009.13452 [hep-ph].
  * [93] Y. Cui, M. Lewicki, D. E. Morrissey and J. D. Wells, Phys. Rev. D 97, no. 12, 123505 (2018) doi:10.1103/PhysRevD.97.123505 [arXiv:1711.03104 [hep-ph]].
  * [94] Y. Cui, M. Lewicki, D. E. Morrissey and J. D. Wells, JHEP 1901, 081 (2019) doi:10.1007/JHEP01(2019)081 [arXiv:1808.08968 [hep-ph]].
  * [95] P. Auclair et al., JCAP 2004, 034 (2020) doi:10.1088/1475-7516/2020/04/034 [arXiv:1909.00819 [astro-ph.CO]].
  * [96] S. Jung and T. Kim, JCAP 2007, 068 (2020) doi:10.1088/1475-7516/2020/07/068 [arXiv:1810.04172 [astro-ph.CO]].
  * [97] D. Bettoni, G. Domènech and J. Rubio, JCAP 1902, 034 (2019) doi:10.1088/1475-7516/2019/02/034 [arXiv:1810.11117 [astro-ph.CO]].
  * [98] C. F. Chang and Y. Cui, Phys. Dark Univ. 29, 100604 (2020) doi:10.1016/j.dark.2020.100604 [arXiv:1910.04781 [hep-ph]].
  * [99] Y. Gouttenoire, G. Servant and P. Simakachorn, JCAP 2007, 032 (2020) doi:10.1088/1475-7516/2020/07/032 [arXiv:1912.02569 [hep-ph]].
  * [100] Y. Gouttenoire, G. Servant and P. Simakachorn, JCAP 2007, 016 (2020) doi:10.1088/1475-7516/2020/07/016 [arXiv:1912.03245 [hep-ph]].
  * [101] L. Sousa, P. P. Avelino and G. S. F. Guedes, Phys. Rev. D 101, no. 10, 103508 (2020) doi:10.1103/PhysRevD.101.103508 [arXiv:2002.01079 [astro-ph.CO]].
  * [102] M. A. Fernandez, S. Bird and Y. Cui, Phys. Rev. D 102, no. 4, 043509 (2020) doi:10.1103/PhysRevD.102.043509 [arXiv:2004.13752 [astro-ph.CO]].
  * [103] R. Zhou and L. Bian, arXiv:2006.13872 [hep-ph].
  * [104] W. Buchmuller, V. Domcke, H. Murayama and K. Schmitz, Phys. Lett. B 809, 135764 (2020) doi:10.1016/j.physletb.2020.135764 [arXiv:1912.03695 [hep-ph]].
  * [105] S. F. King, S. Pascoli, J. Turner and Y. L. Zhou, arXiv:2005.13549 [hep-ph].
  * [106] B. Fornal and B. Shams Es Haghi, Phys. Rev. D 102, no. 11, 115037 (2020) doi:10.1103/PhysRevD.102.115037 [arXiv:2008.05111 [hep-ph]].
  * [107] Z. Arzoumanian et al. [NANOGrav Collaboration], arXiv:2009.04496 [astro-ph.HE].
  * [108] J. Ellis and M. Lewicki, arXiv:2009.06555 [astro-ph.CO].
  * [109] S. Blasi, V. Brdar and K. Schmitz, arXiv:2009.06607 [astro-ph.CO].
  * [110] W. Buchmuller, V. Domcke and K. Schmitz, arXiv:2009.10649 [astro-ph.CO].
  * [111] B. C. Barish, S. Bird and Y. Cui, arXiv:2012.07874 [gr-qc].
  * [112] L. Bian, J. Liu and R. Zhou, arXiv:2009.13893 [astro-ph.CO].
  * [113] S. Chigusa, Y. Nakai and J. Zheng, arXiv:2011.04090 [hep-ph].
  * [114] J. J. Blanco-Pillado and K. D. Olum, Phys. Rev. D 96, no. 10, 104046 (2017) doi:10.1103/PhysRevD.96.104046 [arXiv:1709.02693 [astro-ph.CO]].
  * [115] J. J. Blanco-Pillado, K. D. Olum and B. Shlaer, Phys. Rev. D 89, no. 2, 023512 (2014) doi:10.1103/PhysRevD.89.023512 [arXiv:1309.6637 [astro-ph.CO]].
  * [116] Y. Akrami et al. [Planck Collaboration], Astron. Astrophys. 641, A10 (2020) doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].
  * [117] J. H. MacGibbon, Phys. Rev. D 44, 376 (1991). doi:10.1103/PhysRevD.44.376
  * [118] G. F. Giudice, E. W. Kolb and A. Riotto, Phys. Rev. D 64, 023508 (2001) doi:10.1103/PhysRevD.64.023508 [hep-ph/0005123].
  * [119] A. Abada, S. Davidson, A. Ibarra, F.-X. Josse-Michaux, M. Losada and A. Riotto, JHEP 0609, 010 (2006) doi:10.1088/1126-6708/2006/09/010 [hep-ph/0605281].
  * [120] E. Nardi, Y. Nir, E. Roulet and J. Racker, JHEP 0601, 164 (2006) doi:10.1088/1126-6708/2006/01/164 [hep-ph/0601084].
  * [121] R. Samanta, R. Sinha and A. Ghosal, JHEP 1910, 057 (2019) doi:10.1007/JHEP10(2019)057 [arXiv:1805.10031 [hep-ph]].
  * [122] S. Blanchet and P. Di Bari, JCAP 0703, 018 (2007) doi:10.1088/1475-7516/2007/03/018 [hep-ph/0607330].
  * [123] P. S. B. Dev, P. Di Bari, B. Garbrecht, S. Lavignac, P. Millington and D. Teresi, Int. J. Mod. Phys. A 33, 1842001 (2018) doi:10.1142/S0217751X18420010 [arXiv:1711.02861 [hep-ph]].
  * [124] R. Sinha, R. Samanta and A. Ghosal, JHEP 1712, 030 (2017) doi:10.1007/JHEP12(2017)030 [arXiv:1706.00946 [hep-ph]].
  * [125] R. Samanta, A. Biswas and S. Bhattacharya, arXiv:2006.02960 [hep-ph].
  * [126] R. Samanta and M. Sen, JHEP 2001, 193 (2020) doi:10.1007/JHEP01(2020)193 [arXiv:1908.08126 [hep-ph]].
  * [127] A. Vilenkin, Y. Levin and A. Gruzinov, JCAP 1811, 008 (2018) doi:10.1088/1475-7516/2018/11/008 [arXiv:1808.00670 [astro-ph.CO]].
  * [128] S. Davidson and A. Ibarra, Phys. Lett. B 535, 25 (2002) doi:10.1016/S0370-2693(02)01735-5 [hep-ph/0202239].
  * [129] H. B. Nielsen and P. Olesen, Nucl. Phys. B 61, 45 (1973). doi:10.1016/0550-3213(73)90350-7
  * [130] T. Charnock, A. Avgoustidis, E. J. Copeland and A. Moss, Phys. Rev. D 93, no. 12, 123503 (2016) doi:10.1103/PhysRevD.93.123503 [arXiv:1603.01275 [astro-ph.CO]].
  * [131] T. Vachaspati and A. Vilenkin, Phys. Rev. D 31, 3052 (1985). doi:10.1103/PhysRevD.31.3052
  * [132] A. Vilenkin, Phys. Lett. 107B, 47 (1981). doi:10.1016/0370-2693(81)91144-8
  * [133] J. J. Blanco-Pillado and K. D. Olum, Phys. Rev. D 96, no. 10, 104046 (2017) doi:10.1103/PhysRevD.96.104046 [arXiv:1709.02693 [astro-ph.CO]].
  * [134] J. J. Blanco-Pillado, K. D. Olum and B. Shlaer, Phys. Rev. D 89, no. 2, 023512 (2014) doi:10.1103/PhysRevD.89.023512 [arXiv:1309.6637 [astro-ph.CO]].
  * [135] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 54, 2535-2556 (1996) doi:10.1103/PhysRevD.54.2535 [arXiv:hep-ph/9602271 [hep-ph]].
  * [136] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 65, 043514 (2002) doi:10.1103/PhysRevD.65.043514 [arXiv:hep-ph/0003298 [hep-ph]].
  * [137] P. Amaro-Seoane et al. [LISA Collaboration], arXiv:1702.00786 [astro-ph.IM].
  * [138] W. H. Ruan, Z. K. Guo, R. G. Cai and Y. Z. Zhang, Int. J. Mod. Phys. A 35, no. 17, 2050075 (2020) doi:10.1142/S0217751X2050075X [arXiv:1807.09495 [gr-qc]].
  * [139] J. Luo et al. [TianQin Collaboration], Class. Quant. Grav. 33, no. 3, 035010 (2016) doi:10.1088/0264-9381/33/3/035010 [arXiv:1512.02076 [astro-ph.IM]].
  * [140] V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435 (2006) doi:10.1088/0264-9381/23/7/014 [gr-qc/0512039].
  * [141] N. Seto, S. Kawamura and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001) doi:10.1103/PhysRevLett.87.221103 [astro-ph/0108011].
  * [142] B. Sathyaprakash et al., Class. Quant. Grav. 29, 124013 (2012) Erratum: [Class. Quant. Grav. 30, 079501 (2013)] doi:10.1088/0264-9381/29/12/124013, 10.1088/0264-9381/30/7/079501 [arXiv:1206.0331 [gr-qc]].
  * [143] B. P. Abbott et al. [LIGO Scientific Collaboration], Class. Quant. Grav. 34, no. 4, 044001 (2017) doi:10.1088/1361-6382/aa51f4 [arXiv:1607.08697 [astro-ph.IM]].
  * [144] P. W. Graham, J. M. Hogan, M. A. Kasevich, S. Ra- jendran, and R. W. Romani (MAGIS), (2017), arXiv:1711.02225 [astro-ph.IM].
  * [145] Y. A. El-Neaj et al. (AEDGE), EPJ Quant. Technol. 7, 6 (2020), arXiv:1908.00802 [gr-qc].
  * [146] Y. Gouttenoire, G. Servant and P. Simakachorn, JCAP 07, 032 (2020) doi:10.1088/1475-7516/2020/07/032 [arXiv:1912.02569 [hep-ph]].
  * [147] Y. Cui, M. Lewicki and D. E. Morrissey, Phys. Rev. Lett. 125, no.21, 211302 (2020) doi:10.1103/PhysRevLett.125.211302 [arXiv:1912.08832 [hep-ph]].
  * [148] R. Anantua, R. Easther and J. T. Giblin, Phys. Rev. Lett. 103, 111303 (2009) doi:10.1103/PhysRevLett.103.111303 [arXiv:0812.0825 [astro-ph]].
  * [149] A. D. Dolgov and D. Ejlli, Phys. Rev. D 84, 024028 (2011) doi:10.1103/PhysRevD.84.024028 [arXiv:1105.2303 [astro-ph.CO]].
  * [150] R. Dong, W. H. Kinney and D. Stojkovic, JCAP 1610, 034 (2016) doi:10.1088/1475-7516/2016/10/034 [arXiv:1511.05642 [astro-ph.CO]].
  * [151] F. Li, N. Yang, Z. Fang, R. M. L. Baker, Jr., G. V. Stephenson and H. Wen, Phys. Rev. D 80, 064013 (2009) doi:10.1103/PhysRevD.80.064013 [arXiv:0909.4118 [gr-qc]].
  * [152] A. Arvanitaki and A. A. Geraci, Phys. Rev. Lett. 110, no. 7, 071105 (2013) doi:10.1103/PhysRevLett.110.071105 [arXiv:1207.5320 [gr-qc]].
  * [153] A. Ito, T. Ikeda, K. Miuchi and J. Soda, Eur. Phys. J. C 80, no. 3, 179 (2020) doi:10.1140/epjc/s10052-020-7735-y [arXiv:1903.04843 [gr-qc]].
  * [154] R. w. Hellings and G. s. Downs, Astrophys. J. Lett. 265, L39 (1983). doi:10.1086/183954
  * [155] M. T. Lam et al., Astrophys. J. 834, no. 1, 35 (2017) doi:10.3847/1538-4357/834/1/35 [arXiv:1610.01731 [astro-ph.HE]].
  * [156] G. Hobbs, R. Edwards and R. Manchester, Mon. Not. Roy. Astron. Soc. 369, 655 (2006) doi:10.1111/j.1365-2966.2006.10302.x [astro-ph/0603381].
  * [157] L. Lentati et al., Mon. Not. Roy. Astron. Soc. 453, no. 3, 2576 (2015) doi:10.1093/mnras/stv1538 [arXiv:1504.03692 [astro-ph.CO]].
  * [158] R. M. Shannon et al., Science 349, no. 6255, 1522 (2015) doi:10.1126/science.aab1910 [arXiv:1509.07320 [astro-ph.CO]].
  * [159] Z. Arzoumanian et al. [NANOGRAV Collaboration], Astrophys. J. 859, no. 1, 47 (2018) doi:10.3847/1538-4357/aabd3b [arXiv:1801.02617 [astro-ph.HE]].
  * [160] L. H. Ford, Phys. Rev. D 35, 2955 (1987) doi:10.1103/PhysRevD.35.2955
  * [161] P. J. E. Peebles and A. Vilenkin, Phys. Rev. D 59, 063505 (1999) doi:10.1103/PhysRevD.59.063505 [arXiv:astro-ph/9810509 [astro-ph]].