# Testing Super-Heavy Dark Matter from Primordial Black Holes with
Gravitational Waves

Rome Samanta<EMAIL_ADDRESS>Federico R. Urban<EMAIL_ADDRESS>CEICO, Institute of Physics of the Czech Academy of Sciences, Na Slovance
1999/2, 182 21 Prague 8, Czech Republic

###### Abstract

Ultra-light primordial black holes with masses $M_{BH}<10^{9}$ g evaporate
before big-bang nucleosynthesis producing all matter fields, including dark
matter, in particular super-heavy dark matter: $M_{DM}\gtrsim 10^{10}$ GeV. If
the dark matter gets its mass via $U(1)$ symmetry-breaking, the phase
transition that gives a mass to the dark matter also produces cosmic strings
which radiate gravitational waves. Because the symmetry-breaking scale
$\Lambda_{CS}$ is of the same order as $M_{DM}$, the gravitational waves
radiated by the cosmic strings have a large enough amplitude to be detectable
across all frequencies accessible with current and planned experimental
facilities. Moreover, an epoch of early primordial black hole domination
introduces a unique spectral break in the gravitational wave spectrum whose
frequency is related to the super-heavy dark matter mass. Hence, the features
of a stochastic background of primordial gravitational waves could indicate
that super-heavy dark matter originated from primordial black holes. In this
perspective, the recent finding of a stochastic common-spectrum process across
many pulsars by two nano-frequency pulsar timing arrays would fix the dark
matter mass to be $3\times 10^{13}~{}\text{GeV}\lesssim M_{DM}\lesssim
10^{14}~{}\text{GeV}$. The (non-)detection of a spectral break at
$0.2~{}\text{Hz}\lesssim f_{*}\lesssim 0.4~{}\text{Hz}$ would (exclude)
substantiate this interpretation of the signal.

##  I introduction

Ultra-light primordial black holes (PBHs), that is, black holes with mass
$M_{BH}\lesssim 10^{9}$ g) [1, 2, 3] that exist and evaporate prior to Big
Bang Nucleosynthesis (BBN) via Hawking radiation [4], can leave observable
imprints by producing Gravitational Waves (GWs) [5, 6, 7, 8, 9, 10, 11, 13],
generating the Baryon asymmetry of the Universe (BAU) [15, 16, 14, 17, 18, 19,
20, 21] and cosmologically stable relics [22, 23]. Most importantly, owing to
the fact that PBHs must be agnostic about Standard Model (SM) quantum numbers,
as they evaporate they must also produce Dark Matter (DM) [14, 24, 25, 26, 27,
28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41].

The existence of DM in the Universe is strongly corroborated by experimental
data [42, 43]. However, thus far only the gravitational effects of DM have
been observed, and, despite extensive dedicated experimental programmes, any
attempt to detect possible DM interactions with the SM have yielded null
results111In fact, this is a good reason not to discount the possibility that
what looks like a new particle is instead a deviation from General Relativity
[44, 45].. Therefore, it might well be that DM does not interact with the SM
_at all_. In this case, we need to rethink about (a) how DM could be generated
in the early Universe, and (b) how to further constrain and hopefully detect
and test DM via its gravitational effects—possibly the only way to do so. In
this article we show how, if ultra-light PBHs existed in the early Universe,
(a) DM, in fact super-heavy dark matter (SHDM) with masses above
$M_{DM}\gtrsim 10^{10}$ GeV, can be easily generated, and (b) the PBH origin
of SHDM is testable through its peculiar features in the spectrum of
primordial Gravitational Waves that this mechanism would produce.

Owing primarily to its very large energy scale, SHDM is very difficult to
make. Existing mechanisms are, e.g., gravitational production at the end of
inflation [46, 47, 48, 49], supersymmetry breaking [50], freeze-in [51, 52],
thermal freeze-out [53, 54, 55, 56] and phase transitions [57, 58]. The
mechanism of production of SHDM via PBHs, which we adopt in this work, has two
distinct advantages. First, with PBHs one can access a large range of SHDM
masses, potentially up to the Planck scale, without having to specify the
details of the couplings of the dark sector to the SM or gravity [59, 60].
Second, for a large enough initial PBH energy density, PBHs come to dominate
the energy density of the Universe before evaporation; when they finally
evaporate they inject entropy in the system and significantly dilute any
preexisting model-dependent SHDM abundance, thereby making the ${\rm
PBH\rightarrow SHDM}$ channel unambiguously dominant.

SHDM is also difficult to test because it is too heavy to be produced directly
in colliders, even assuming that SHDM interacts with the SM in the first
place. If SHDM does interact, it is expected to decay into ultra-high-energy
cosmic rays (UHECRs), either charged nuclei or neutral photons or neutrinos
[61]. The flux of UHECRs observed by experiments is very low, and there are
large uncertainties about the energy spectrum and chemical composition.
Nonetheless, from UHECR data we know that SHDM can only make up for a small
fraction of the UHECR flux below $E\approx 10^{20}$ eV, but could come to
dominate at higher energies (see, e.g., [62]). UHECR data constrain the
lifetime of SHDM to be larger than about $10^{22}$ yr to $10^{23}$ yr [63]. In
this work we aim to show how ${\rm PBH\rightarrow SHDM}$ can in fact be tested
with GWs in a distinctive manner.

Our idea is based on the following arguments: if, like all known elementary
particles, the SHDM acquires a mass via a phase transition (PT) mechanism, the
scale of the transition (typically defined at a temperature $T\sim v$, where
$v$ is the vacuum expectation value of the scalar field responsible for the
PT) should be at least of the order of the SHDM mass. The same PT, if it is
strongly first-order [64] or if it produces topological defects [65], can also
be a source of GWs. It is well known that cosmic strings [65, 66, 67] which
originate as topological defects due to the breaking of gauge symmetries such
as $U(1)$, produce detectable GWs if the symmetry breaking scale $T\sim
v\sim\Lambda_{CS}\gtrsim 10^{10}$ GeV [68, 69, 70, 71, 72, 73, 74, 75]. Thus,
cosmic strings provide an outstanding opportunity to probe the physics
involving super-heavy particles [76, 77, 78, 79] that acquire a mass via PTs,
because detectable GWs from cosmic strings can be naturally associated with
the SHDM mass generation via the $U(1)$ breaking at a scale
$\Lambda_{CS}\gtrsim M_{DM}$222From a more fundamental theory perspective, it
is quite easy to naturally obtain a $U(1)$-symmetric SHDM mass term at the
intermediate stage of a Grand Unified Theory (GUT) to SM breaking, wherein the
$U(1)$ can be identified as $U(1)_{B-L}$ [76, 85, 86]..

The most distinct feature of GWs from cosmic strings is a scale invariant
spectrum spanning a wide range of frequencies, which makes them an ideal
candidate for multi-frequency GW studies [73, 74, 75]. If ultralight PBHs with
large enough initial energy density existed in the early universe, apart from
efficiently producing SHDM and diluting any preexisting SHDM density via
entropy injection, they would introduce one further signature in the GW
spectrum, namely a spectral break at a frequency $f_{*}$, beyond which the
spectrum turns red. The frequency $f_{*}$ for a given SHDM mass is unique
provided that the SHDM generated by the evaporation of PBHs makes up all of
the observed DM. We find that, for the allowed range of SHDM masses, $f_{*}$
lies in a very constrained and testable range. This peculiar feature in the
GWs spectrum makes the ${\rm PBH\rightarrow SHDM}$ scenario quite distinct
from standard scenarios of GWs from cosmic strings333GWs can also arise
directly from either PBH production mechanisms or their evolution [9, 10, 11,
12]. Specifically, in Sec.V, we discuss the induced GWs [11, 12] from the
density fluctuations of PBHs which potentially makes our scenario unique in
terms of its GW signatures..

Note also that in this set-up heavier SHDM requires larger symmetry breaking
scales $\Lambda_{CS}$. Since the amplitude of GWs from cosmic strings
increases with $\Lambda_{CS}$, the detectability improves with increasing SHDM
mass. Adopting this point of view, it is tempting to interpret the recent
finding of a high-amplitude stochastic common-spectrum process across 45
pulsars by the NANOGrav pulsar-timing array [80] (confirmed also by PPTA in an
updated analysis [81]) in terms of a background of GWs radiated by CSs [82,
83, 84]. In this case, we show that the data can be fitted with a SHDM mass
$3\times 10^{13}~{}\text{GeV}\lesssim M_{DM}\lesssim 10^{14}~{}\text{GeV}$,
and that our model can be distinguished from other cosmic strings fits owing
to its unique prediction of a spectral break at a turning-point frequency
$0.2~{}\text{Hz}\lesssim f_{*}\lesssim 0.4~{}\text{Hz}$.

Let us briefly summarise the main ideas here. We assume that ultra-light PBHs
existed sometime in the very early Universe, before BBN. These PBHs have a
large initial energy density and eventually come to dominate the energy budget
of the Universe before evaporating. As they evaporate into all fields, they
must also produce DM, which, as it turns out, has to be super-heavy. If, and
this is our second assumption, SHDM gets its mass via $U(1)$ symmetry-
breaking, the phase transition that gives a mass to the DM also produces
cosmics strings which radiate GWs. The symmetry-breaking scale $\Lambda_{CS}$
has to be very large because SHDM is super-heavy. Therefore, the GWs radiated
by cosmic strings have a large enough amplitude to be detectable. Lastly, the
PBH domination epoch introduces a unique spectral break in the GW spectrum
whose frequency is related to the SHDM mass. Therefore, these marked features
of a primordial GWs signal, potentially detectable across all frequencies
accessible with current and planned experimental facilities, would indicate
that SHDM originated from PBHs. A possible timeline for the proposed scenario
is shown in Fig.1.

The rest of the paper is organised as follows: in Sec.II we obtain an
expression for the GW background we expect from the cosmic string network. In
Sec.III we explain how PBHs produce SHDM. In Sec.IV we present our numerical
results including the fit to NANOGrav data. In Sec. V we discuss how GWs from
PBH density perturbation could be a complementary signature of our model. We
also discuss a possible extension of our scenario to address baryogenesis.
Finally, in Sec.VI, we conclude and give an outlook for further developments
of this mechanism.

Figure 1: A possible timeline for the proposed scenario.

##  II Gravitational waves from cosmic strings and their spectral features

We consider a $U(1)$-charged complex scalar field $\Phi$ with the tree level
potential
$V(\Phi)=-\frac{m_{\Phi}^{2}}{2}|\Phi|^{2}+\frac{\lambda}{4}|\Phi|^{4}$, where
$m_{\Phi}$ and $\lambda$ are the field’s mass and self-interacting coupling
constant, respectively. After the spontaneous breaking of the $U(1)$ symmetry,
the field $\Phi$ gets a vacuum expectation value ($\Lambda_{CS}\equiv
v_{\Phi}=m_{\Phi}/\sqrt{\lambda}$) and generates the mass of the DM as well as
gives rise to cosmic strings. As a prototype, let the DM be a Majorana
fermion444This, although not necessary for our purposes, is nonetheless useful
in order to embed the scenario in a more realistic model which also includes
baryogenesis, as we discuss in Sec.V. $\chi$ having $U(1)$ charge
$q_{\chi}=-1$; therefore, to make the mass term $\mathcal{L}_{DM}\sim
y_{\chi}\chi\chi\Phi$ symmetric under $U(1)$, we set $q_{\Phi}=2$. Once the
symmetry is broken the DM becomes massive ($M_{DM}=y_{\chi}v_{\Phi}$) and
cosmic strings are formed. After their formation, the strings, which are
randomly distributed in space, form close loops and a network of horizon-size
long strings [87, 88]. When two segments of long strings cross each other they
inter-commute and form loops. Long strings are characterised by a correlation
length $L=\sqrt{\mu/\rho_{\infty}}$, where $\rho_{\infty}$ is the long string
energy density and $\mu$ is the string tension defined as $\mu=\pi
v_{\Phi}^{2}h\left(\lambda/2g^{2}\right)$. The quantity $h$ is a slowly
varying function of its arguments ($\lambda$ and the gauge coupling $g$) with
$h(1)\simeq 1$. In numerical simulations of cosmic strings, usually both the
couplings are taken to be $\mathcal{O}(1)$ and therefore, the string tension
$\mu$ is defined as $\mu=\pi v_{\Phi}^{2}$. We shall work with
$y_{\chi}=\lambda=1$ and $g=1/\sqrt{2}$ so that the string tension becomes
$\mu=\pi M_{DM}^{2}$. The choice $M_{DM}=v_{\Phi}$ is well motivated to make
the model free from DM relics that may originate due to the pair production of
SHDM from the scalar ($M_{\Phi}=\sqrt{2}v_{\Phi}$) and the heavy gauge boson
($M_{Z}=v_{\Phi}$). In any case, at the end we will see that any pre-existing
SHDM relic before PBH evaporation will be diluted significantly.

A string network interacts strongly with the thermal plasma, thereby its
motion gets damped [89]. Once the damping phase is over, the strings oscillate
and enter a phase of scaling evolution in which two competing dynamics
coexist, namely, the stretching of the correlation length owing to the cosmic
expansion and the fragmentation of the long strings into close loops which
oscillate independently and produce GWs or particle radiation [68, 69, 70].
Between these two competing dynamics, there exists an attractor solution
called the scaling regime [90, 91, 92] in which the characteristic length $L$
scales as cosmic time $t$. In this regime, for a constant string tension, we
have $\rho_{\infty}\propto t^{-2}$. Therefore, the network tracks any
cosmological background energy density $\rho_{bg}\propto t^{-2}$ with a small
constant number proportional to $G\mu$, where $G$ is the Newton constant. This
scaling behaviour prevents the cosmic string network from dominating the
energy density of the Universe unlike all other cosmological defects.

Motivated by the results of cosmic string network simulations (see Refs. [71,
72, 90, 91]) we always assume that the network is in the scaling regime while
contributing to the GWs. The loops radiate GWs at a constant rate $\Gamma$
which determines the time evolution of a loop of initial (i.e., at creation
time $t_{i}$) size $l_{i}=\alpha t_{i}$ as $l(t)=\alpha t_{i}-\Gamma
G\mu(t-t_{i})$, where $\Gamma\simeq 50$ [68, 70] and $\alpha\simeq 0.1$ are
determined by numerical simulations [73, 74]. The total energy loss from a
loop can be decomposed into a set of normal-mode oscillations with frequencies
$f_{k}=2k/l=a(t_{0})/a(t)f$, where $k=1,2,3...k_{max}$
($k_{max}\rightarrow\infty$), $f$ is the frequency observed at present time
$t_{0}$ and $a$ is the scale factor of the Universe. The GW density parameter
is given by $\Omega_{GW}(t_{0},f)\equiv
f\rho_{c}^{-1}d\rho_{GW}/df=\sum_{k}\Omega_{GW}^{(k)}(t_{0},f)$, with the
$k$-th mode amplitude $\Omega_{GW}^{(k)}(t_{0},f)$ as [73]

$\displaystyle\Omega_{GW}^{(k)}(f)=\frac{2kG\mu^{2}\Gamma_{k}}{f\rho_{c}}\int_{t_{osc}}^{t_{0}}\left[\frac{a(t)}{a(t_{0})}\right]^{5}n\left(t,l_{k}\right)dt,$
(II.1)

where $n\left(t,l_{k}\right)$ is the loop number density, $\rho_{c}$ is the
critical energy density of the Universe, $\rho_{GW}$ is the GW energy density
and $l_{k}$ is loop length corresponding to the frequency $f_{k}$. The small
scale structures of the loop sets the quantity $\Gamma_{k}$ which is given by
$\Gamma_{k}=\frac{\Gamma k^{-\delta}}{\zeta(\delta)}$, with $\delta=4/3$ and
$5/3$ for loops containing cusps and kinks, respectively [93]. Eq.II.1 is
valid only for $t_{i}>t_{osc}={\rm Max}\,\left[{\rm
network~{}formation~{}time}~{}(t_{F}),{\rm
end~{}of~{}damping}~{}(t_{\text{fric}})\right]$ and $t_{i}>l_{crit}/\alpha$,
with $l_{crit}$ the critical length above which GWs dominate over the massive
particle radiation [94, 95]. As discussed in Sec.V, both these bounds set a
high-frequency cut-off in the spectrum (a comprehensive analysis can be found
in Ref. [96]).

With the standard cosmological evolution, Eq.II.1 leads to the following
spectral features: (I) a peak at low frequency due to the GWs from the loops
that are produced in the radiation era and decay in the following matter era
and (II) a plateau at high frequency [73, 74, 99] given by (see Appendix A)

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

where $\epsilon_{r}=\alpha/\Gamma G\mu\gg 1$, $A_{r}\simeq 5.4$ being the loop
production efficiency in radiation domination and $\Omega_{r}\sim 9\times
10^{-5}$ is the present radiation energy density fraction—this behaviour
arises from the loop dynamics in the radiation era. Since
$\mu\sim\Lambda_{CS}^{2}$, the plateau amplitude is
$\Omega_{GW}^{plt}\sim\Lambda_{CS}$: this is a property that makes cosmic
strings an outstanding probe of super-high-scale physics—in our context this
scale is related to the SHDM mass scale.

For a very long period of radiation domination, the plateau remains intact
across all frequencies accessible by current and planned GW detectors such as
LISA [100] and LIGO [101, 102]. However, if a non-standard cosmological
evolution is introduced before the onset (at $T_{*}$) of the radiation
domination era that leads to BBN ($T_{\text{BBN}}\sim 5$ MeV) [103, 104, 105],
the spectrum deviates from the constant plateau at high frequencies: this
happens at a turning-point frequency $f_{*}$ which depends on the end of the
non-standard evolution epoch. In case of an early matter domination, the
spectral slope changes from $\Omega_{GW}\sim f^{0}$ to $\Omega_{GW}\sim
f^{-1}$ beyond $f_{*}$. The sensitivity reach of GW detectors to probe $f_{*}$
and $T_{*}$ can be found, e.g., in Fig. 4 and 5 of Ref. [75]. This scenario is
relevant in our case, since the PBHs introduce an early matter domination era
before BBN. For the $k=1$ mode, the $f_{*}$ can be calculated in a simple
analytical way that gives [78, 75]

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

where $z_{\rm eq}\simeq 3387$ is the red-shift at the usual matter-radiation
equality, which takes place at time $t_{\rm eq}$.

Note that even though the GW spectrum beyond the turning-point frequency
behaves as $\Omega(f>f_{*})\sim f^{-1}$ for $k=1$, when the contributions from
all the modes are taken into account one finds $\Omega(f>f_{*})\sim f^{-1/3}$
(see Appendix B). Moreover, even without a non-standard cosmological epoch,
GWs from cosmic strings may exhibit a non-flat spectrum in the case of global
strings [106, 107], inflation-diluted strings [108, 109] and melting strings
[110].

##  III  generalities of PBH dynamics and Dark matter production

A useful parameter to track the evolution of the energy density of the black
holes is

$\displaystyle\beta\equiv\frac{\rho_{BH}(t_{Bf})}{\rho_{\rm R}(t_{Bf})}\,,$
(III.1)

where $\rho_{\rm R}(t_{Bf})$ and $\rho_{BH}(t_{Bf})$ are the energy densities
of PBHs and radiation at the PBH formation time $t_{Bf}$, respectively. If the
PBHs do not dominate the energy density of the Universe at their evaporation,
the parameter $\beta$ is bounded from above as [9] (Appendix C)

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

where $\gamma\simeq 0.2$ is the black hole formation efficiency,
$T_{BH}=M_{Pl}^{2}/8\pi M_{BH}$ [4] is the PBH temperature, $\mathcal{G}\simeq
3.8$ being the graybody factor, $g_{*B}\simeq 100$ is the number of
relativistic particle species below $T_{BH}$ in the SM with three left-handed
light neutrinos [113] and we approximate non-rotating PBHs with a
monochromatic mass spectrum. In Fig. 2 (left panel) we show the $\beta$
parameter as a function of the PBH mass. The region above (below) the dashed
black line corresponds to PBH (radiation) domination. The region between the
dashed black and solid blue lines is allowed for the induced GWs by PBH
density fluctuation not to saturate the BBN bound[11] (see discussion in
Sec.V). For future reference, we also indicate three benchmark values
BP$(1,2,3)=10^{-6,-9,-12}$ for $M_{BH}=10^{5}$ g. Note that BP3 is in the
radiation domination region, i.e., for $\beta=10^{-12}$, PBHs of mass $10^{5}$
g can never dominate the energy density of the Universe. On the contrary, PBHs
come to dominate for some time for the other two BPs.

Figure 2: Left panel: the parameter $\beta$ of Eq.III.1 as a function of
$M_{BH}$. The region above (below) the dashed black line corresponds to PBH
(radiation) domination. The region between the dashed black and solid blue
lines represents the allowed ranges of $\beta$ (cf. Ref.[11]). The horizontal
red region is excluded by the CMB constraint on Hubble scale of inflation
$H_{\rm inf}\lesssim 3\times 10^{14}$ GeV at 95$\%$ CL [112]. Right panel:
evolution of the energy density fractions of the PBHs (dashed lines) and
radiation (solid lines) for $M_{BH}=10^{5}$ g. In both plots the benchmark
points correspond to $\beta=10^{-6,-9,-12}$. The dashed vertical line
corresponds to the PBH evaporation temperature calculated analytically in
Eq.C.11 for radiation domination. Eq.D.6 would correspond to a very similar
line.

Larger values of $\beta$ correspond to longer period of PBH domination: As we
discuss later in this section, a large value of $\beta$ is preferred in our
model to make ${\rm PBH\rightarrow SHDM}$ the dominant channel. The time
evolution of the energy densities can be understood by solving the following
Friedmann equations555These two equations are derived from the standard
Friedmann equations in presence of PBHs and the entropy non-conservation
equation due the evaporation of PBHs, see Appendix C.  [30, 114]

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

where $z=M_{0}/T$ with $M_{0}$ being an arbitrary mass scale which we fix as
$M_{0}=T_{Bf}$ (PBH formation temperature) and the parameter $\mathbb{K}$
defined 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)\equiv-\mathbb{K}$
(III.5)

with $g_{*s}(T)$ being the number of entropy degrees of freedom. The solutions
of Eq.III.3 and Eq.III.4 are shown in the right panel of Fig. 2 for the same
BPs as in the left panel. It can be seen that, the larger $\beta$, the longer
the period of PBH domination. Even though the non-thermal DM production from
PBHs is independent of $\beta$ provided PBHs dominate the energy density of
the Universe, large values of $\beta$ make this scenario more predictive as
any preexisting relic DM density will be mostly washed out, making the ${\rm
PBHs\rightarrow DM}$ channel the dominant one.

Let us now turn to the production of DM from PBHs. A general expression for
the DM relic energy density fraction is given by

$\displaystyle\Omega_{\rm DM}h^{2}=\frac{M_{\rm
DM}n_{\gamma}^{0}}{10.54f({T_{\rm ev}},T_{0}){\rm
GeVm^{-3}}}\left(\frac{N_{{\rm DM}}}{N_{\gamma}}\right)_{T_{\rm ev}}\simeq
1.45\times 10^{6}\left(\frac{N_{{\rm DM}}}{N_{\gamma}}\right)_{T_{\rm
ev}}\left(\frac{M_{\rm DM}}{\rm GeV}\right),$ (III.6)

where $N_{\rm DM}$ and $N_{\gamma}$ are the DM and photon number densities,
respectively, normalised to the ultra-relativistic equilibrium number density
$n_{f,\rm eq}^{\rm ur}=g_{f}T^{3}/\pi^{2}$ of a spin-$1/2$ fermion666This is a
standard ‘$T^{3}$’ normalisation. One can also normalise it with respect to
the entropy density $s\sim T^{3}$ [30]., $n_{\gamma}^{0}\simeq 410.7\times
10^{6}\rm m^{-3}$ and $f(T_{\rm ev},T_{0})\simeq 27.3$ are the relic photon
number density at the present time and photon dilution factor respectively.
Given that $\Omega_{\rm DM}h^{2}\simeq 0.12$ [42], from Eq.III.6 we find

$\displaystyle N_{\rm DM}^{\rm ev}=N_{\rm DM}^{\rm Obs}\simeq 1.1\times
10^{-7}\left(\frac{\rm GeV}{M_{\rm DM}}\right).$ (III.7)

One has to compare Eq.III.7 with the number density of DM produced by PBHs,
i.e.,

$\displaystyle N_{\rm DM}^{\rm ev}=N_{\rm BH}^{\rm ev}\bar{n}_{\rm DM},$
(III.8)

where $\bar{n}_{\rm DM}$ is the number of DM particles produced by a black
hole and $N_{\rm BH}^{\rm ev}$ is the normalised number density of PBHs at
evaporation.

Computing $N_{\rm BH}^{\rm ev}\bar{n}_{\rm DM}$ in a PBH dominated Universe
(see Appendix D) and noting that $N_{DM}^{\rm ev}\lesssim N_{DM}^{Obs}$ we
obtain

$\displaystyle M_{DM}\lesssim 3\times
10^{-7}\left(\frac{M_{BH}}{M_{Pl}}\right)^{1/2}{\rm GeV}~{}~{}{\rm
for}~{}~{}T_{BH}>M_{DM},$ (III.9)

and

$\displaystyle M_{DM}\gtrsim 4.5\times
10^{3}\left(\frac{M_{BH}}{M_{Pl}}\right)^{-5/2}M_{Pl}^{2}~{}~{}{\rm
GeV^{-1}}~{}~{}{\rm for}~{}~{}T_{BH}<M_{DM}$ (III.10)

where the sign $\simeq$ corresponds to $\Omega_{\rm DM}h^{2}\simeq 0.12$ and
in Eq.III.10, dark matter is produced when the PBHs becomes light enough so
that the corresponding $T_{BH}$ becomes larger than $M_{DM}$ (see Appendix D).
Eq.III.9 and Eq.III.10 represent to the light and the heavy DM scenarios,
respectively. However, there is another constraint that needs to be taken into
account in the case of non-thermal light DM. The constraint comes from the
fact that the free-streaming length of light DM should not be too large at
$z_{\rm eq}$ (Ly$\alpha$ cloud constraint) or it will interfere with structure
formation [115, 116, 117, 118]. This has been studied extensively in [25, 30,
31], and the result is that the DM mass is constrained from below as

$\displaystyle M_{DM}>4.4\times
10^{-6}\left(\frac{M_{BH}}{M_{pl}}\right)^{1/2}.$ (III.11)

Figure 3: Parameter space in the $M_{DM}-M_{BH}$ plane: the light green region
in the top-right corner is the allowed region ($\Omega_{DM}h^{2}<0.12$) and
the thick solid line marks the values for which we obtain the observed dark
matter abundance $\Omega_{DM}h^{2}\simeq 0.12$. We exclude $M_{DM}>10^{15}$
GeV assuming that the symmetry breaking that generates the dark matter mass
takes place after inflation, which also implies that the cosmic string loop
number density will be undiluted. The gray region is excluded due to the too
large free-streaming length of the DM. Figure 4: The evolution of the PBH
energy density fraction (dashed lines), a preexisting dark matter abundance
(solid lines) and the total entropy density (dot-dashed lines) for a long PBH-
domination (red) or short PBH-domination (blue) period, or in radiation
domination throughout (green).

The constraint in Eq.III.11 is incompatible with that in Eq.III.9 and rules
out the ${\rm PBHs\rightarrow DM}$ scenario for light DM. Therefore the only
viable option is Eq.III.10, i.e., the heavy (in fact, super-heavy) DM
scenario. All the constraints on the $M_{DM}-M_{BH}$ parameter space are shown
in Fig. 3. In general, the light green region in the top-right corner is
allowed, i.e., for this region one has $\Omega_{DM}h^{2}<0.12$. Note that, the
produced DM is super-heavy and in principle the mass $M_{DM}$ spans the range
$10^{10}{\rm GeV}\lesssim M_{DM}\lesssim M_{Pl}$. This range is much wider
than what can be achieved in scenarios of gravitational production [46, 47,
48]. However, conservatively, we will limit ourselves to a maximum mass of
$M_{DM}^{max}\simeq 10^{15}\sim T_{RH}^{max}$ GeV, i.e., a post-inflationary
symmetry breaking that generates the DM mass as well as an undiluted cosmic
strings loop number density [108, 109]. The upper limit on $M_{DM}$
corresponds to a lower bound on the PBH mass: $M_{BH}\gtrsim 10^{6}$ g.

Let us conclude this section by analysing the fate of a preexisting DM relic
($N_{DM}^{i}$) for the three benchmark values of $\beta$ and for
$M_{BH}=10^{7}$ g, which corresponds to $M_{DM}\simeq 5\times 10^{12}$ GeV
(cf. Fig. 3). First of all, in presence of PBHs, the Boltzmann equation for a
preexisting DM relic can be derived as

$\displaystyle\frac{dN_{DM}}{dz}-\frac{3}{\mathbb{K}z}\left(\mathbb{K}-H\right)N_{DM}=0$
(III.12)

which in the case of $\mathbb{K}\rightarrow H$, i.e., when the PBHs are absent
or do not dominate the energy density of the universe, leads to
$N_{DM}(z\rightarrow\infty)=N_{DM}^{i}$, as expected. However, when PBHs
dominate for some period the situation is very different: due the appearance
of the second term in the LHS of Eq.III.12, any preexisting relic will be
diluted by an amount that depends on $\beta$ and hence on the duration of the
PBH domination epoch. The final value of $N_{DM}^{i}$ is then given by

$\displaystyle N_{DM}^{f}=N_{DM}^{i}{\rm
exp}\left[\int_{1}^{\infty}\frac{3}{\mathbb{K}z}\frac{\dot{M}_{BH}}{M_{BH}}\frac{\rho_{BH}}{4\rho_{R}}\right].$
(III.13)

This dilution follows from the entropy non-conservation equation (Eq.C.19)
that leads to $\mathbb{K}\neq H$ in Eq.III.5. Neglecting the time variation of
$g_{*s}$, the amount of total entropy $\tilde{S}\sim a^{3}/z^{3}$ production
can be tracked using

$\displaystyle\frac{da}{dz}=\left(1-\frac{\dot{M}_{BH}}{M_{BH}}\frac{\rho_{BH}}{4\rho_{R}\mathbb{K}}\right)\frac{a}{z}$
(III.14)

that can be easily derived from Eq.C.19.

In Fig. 4, we show the evolution of $N_{DM}^{i}(z)$, $\Omega_{BH}(z)$ and
$\tilde{S}(z)$ in solid, dashed and dot-dashed lines respectively; we show the
cases $\beta=10^{-4}$ (red, strong PBH domination), $\beta=10^{-9}$ (blue,
mild PBH domination), $\beta=10^{-14}$ (green, radiation domination). Note
that the longer the PBH domination, the more the produced entropy, and the
stronger is the dilution of $N_{DM}^{i}$. For example, a PBH of mass $10^{7}$
g with $\beta=10^{-4}$ can dilute any preexisting DM density by a factor of
$10^{8}$. However, the amount of a preexisting relic is strongly model
dependent, e.g., in our case, if kinematically allowed, the $U(1)$ gauge boson
or scalar can produce the initial $N_{DM}^{i}$. Nonetheless, as shown already,
in case of a strong PBH domination such a relic density is strongly suppressed
at the time of PBH evaporation. In scenarios like those of Refs. [10, 11],
where the GWs that come from the density perturbations of PBHs can potentially
give an upper bound on $\beta$, a large value of $N_{DM}^{i}$ may be a cause
of concern.

##  IV Imprints of SHDM on the GWs from DM mass origin/cosmic strings

We start from the scenario in which the SHDM becomes massive before the PBHs
are formed (Fig.1). The DM mass scale corresponding to correct DM density is
larger than the PBH formation temperature for $M_{DM}\gtrsim M_{DM}^{c}\simeq
10^{12}$ GeV (Fig. 3, thick black line). This leads to a symmetry breaking
scale $v_{\Phi}\equiv\Lambda_{CS}\simeq{\rm Max}\left[T_{Bf},M_{DM}\right]$.
To evaluate $G\mu$, one can therefore set $\mu\sim M_{DM}^{2}$ for
$M_{DM}>M_{DM}^{c}$ and $\mu\sim T_{Bf}^{2}$ for $M_{DM}<M_{DM}^{c}$777Note
that $\mu\sim T_{Bf}^{2}$ is a possibility (Fig.1), in which case, the
condition $M_{DM}\lesssim T_{Bf}$ requires $y_{\chi}\lesssim 1$. Because we
discuss only the $\mu\sim M_{DM}^{2}$ scenario, we shall not explore this
possibility further.. As mentioned earlier, in our set-up, the GW amplitude as
well as the spectral features are sensitive to the details of the ${\rm
PBHs\rightarrow SHDM}$ mechanism. Let us look at the amplitude first. Consider
the case $v_{\Phi}=M_{DM}>10^{12}$ GeV. The spectra of GWs for
$M_{DM}=10^{15,13,12}$ GeV are shown in Fig. 5 (left panel) in solid, dashed,
and dotted blue lines, respectively. The spectrum is in accordance with the
standard expectation that with large $v_{\Phi}$ the amplitude of GWs
increases, i.e., in our case, the prospects for detection improve as the SHDM
becomes more massive. The most interesting fact is that even if the SHDM mass
is as ‘light’ as $10^{11}$ GeV, mid-band detectors such as LISA [100], DECIGO
[119] and BBO [120] are able to test the parameter space fully.

In principle we can relax the assumption of $y_{\chi}=1$, e.g., a SHDM of mass
$M_{DM}\sim 10^{14}$ GeV can be obtained with a symmetry breaking scale
$v_{\Phi}\sim 10^{15}$ GeV when $y_{\chi}=0.1$. In this case, the string
tension $G\mu$ ranges from $\sim Gv_{\Phi}^{2}$ to $\sim GM_{DM}^{2}$.
Therefore, as we go to lower DM masses (compared to $v_{\Phi}$), the allowed
range of $G\mu$ increases. We indicate this in Fig. 5 (left panel). For
$M_{DM}=10^{13}(10^{12})$ GeV the GW spectrum can be anywhere within the
region between the solid and the dashed (dotted) lines. In the case of
$y_{\chi}\ll 1$, the SHDM could be produced by the scalar and the gauge boson
having abundant initial number densities. In that case, one has to consider
the period of PBH domination accordingly to dilute the generated $N_{DM}^{i}$.
In this article, we do not present any explicit computation in this direction.

We turn now to the spectral shape by computing the turning-point frequency
$f_{*}$. First of all, the PBH evaporation temperature in the case where they
dominate the energy density, is derived as (see Appendix D):

$\displaystyle
T_{ev}=\left(\frac{5M_{Pl}^{2}}{\pi^{3}g_{*}(T_{ev})\tau^{2}}\right)^{1/4},$
(IV.1)

with $\tau$ being the life-time of a PBH. Therefore, from Eq.IV.1 and
Eq.III.10 we obtain the evaporation temperature in terms of the SHDM mass as

$\displaystyle T_{ev}=2.1\times 10^{-8}\left(\frac{M_{DM}}{{\rm
GeV}}\right)^{3/5}.$ (IV.2)

From Eq.II.3, the turning point frequency $f_{*}$ can then be expressed in
terms of $T_{ev}$ and hence the SHDM mass as

$\displaystyle f_{*}\simeq 2.1\times 10^{-8}\sqrt{\frac{50}{z_{\rm
eq}\alpha\Gamma
G\mu}}\left(\frac{M_{DM}}{T_{0}}\right)^{3/5}T_{0}^{-2/5}t_{0}^{-1},$ (IV.3)

where $T_{0}=2.7$ K is the photon temperature today. The expression IV.3 is of
immense interest and it is one of the most important results of this work. For
a given value of the DM mass and $G\mu$ (for $y_{\chi}=1$, $G\mu\sim
GM_{DM}^{2}$), one can determine approximately at which frequency the spectrum
changes slope from a plateau described by $f^{0}$ to $f^{-1/3}$ (see the
discussion in Sec.II). For example, a SHDM of mass $10^{15}$ GeV would
correspond to a break at $f_{*}\sim 0.1$ Hz, as shown in Fig. 5 (left panel).
A robust prediction of our ${\rm PBHs\rightarrow SHDM}$ set-up is that, given
the allowed range of $M_{DM}$ (Fig. 3) , the turning-point frequencies always
lie within the sensitivity range of mid-band detectors:
$0.1~{}\text{Hz}\lesssim f_{*}\lesssim 7~{}\text{Hz}$. This window will
enlarge for $y_{\chi}\ll 1$. Reconstruction techniques for such spectral
shapes can be found in [121].

Figure 5: Left panel: the gravitational wave spectrum for our
$\rm{PBHs\rightarrow SHDM}$ scenario. The solid blue line represents the
spectrum for $M_{DM}\simeq 10^{15}$ GeV. As the SHDM mass decreases the
parameter space on the $\Omega_{GW}h^{2}-f$ plane enlarges, the allowed region
for $M_{DM}\simeq 10^{15}{\rm GeV}-10^{13}$ GeV lies between the solid and the
dashed lines, whereas the region between the solid and the dotted lines is
allowed for $M_{DM}\simeq 10^{15}{\rm GeV}-10^{12}$ GeV. Right panel: for a
given value of $G\mu$ we show the potential of various detectors to probe the
PBH evaporation temperature and hence $M_{DM}$. The blue regions are the
allowed parameter space in our ${\rm PBHs\rightarrow SHDM}$ scenario. The red
shades are excluded due particle production from the loops and thermal
friction that freezes the network, see main text.

In Fig. 5 (right panel), we show the allowed range in the $G\mu-M_{DM}$
parameter space that corresponds to such spectral breaks. The dark-blue region
corresponds to the case when $v_{\Phi}>T_{Bf}$ (see Fig.1), whereas the
parameter space enlarges to include the light-blue region when the condition
$v_{\Phi}>T_{Bf}$ is relaxed. As mentioned earlier in Sec.II, the loops
effectively contribute to the GWs for $t>t_{i}$ which is bounded from below in
presence of thermal friction and particle production from the loops. A lower
bound on $t_{i}$ (upper bound on temperature) also sets a cut-off on the GW
spectrum at a high frequency beyond which the spectrum falls. The effects of
thermal friction is sub-dominant when $H\sim T^{2}/M_{Pl}>\Gamma_{fric}\sim
T^{3}/\mu$ [89]. This implies that the motion of the network gets damped until
$T_{fric}\sim\mu/M_{Pl}$. For the GW spectrum to fall due to PBH domination,
which is our requirement, one should have $T_{fric}>T_{ev}$. This corresponds
to the constraint

$\displaystyle G\mu>2.1\times 10^{-8}\left(\frac{M_{DM}}{{\rm
GeV}}\right)^{3/5}\left(\frac{M_{Pl}}{\rm GeV}\right)^{-1}.$ (IV.4)

A very recent numerical simulation reported that there is a critical size
$l<l_{crit}$ for which particle production is the dominant radiation mode and
the loop number density does not scale [94, 95]. This corresponds to a lower
bound $t_{i}>l_{crit}/\alpha$ and therefore a cut-off on the GW spectrum at a
high frequency. Now proceeding in the same way as in the friction case, a
bound on $G\mu$ is obtained as [78]

$\displaystyle G\mu>1.7\times 10^{-22}\left(\frac{M_{DM}}{{\rm
GeV}}\right)^{12/25},$ (IV.5)

where we have considered particle production from the cusps with
$l_{crit}=\frac{\mu^{-1/2}}{(\Gamma G\mu)^{2}}$ [94, 95] and $\alpha=0.1$. The
constraints in Eq.IV.4 and Eq.IV.5 are shown in Fig.5 (right) by dark-red and
light-red shades respectively. It is evident that the particle production cut-
off is much stronger and rules out a portion of the allowed parameter space.

We now discuss how, in this set-up, the recent finding of a stochastic common
spectrum process by the NANOGrav pulsar timing array would constrain the range
of SHDM mass and how the turning point frequencies potentially can be a
complementary probe of our ${\rm PBHs\rightarrow SHDM}$ scenario. With their
recently released 12.5 yrs data set, the NANOGrav collaboration has reported
strong evidence for a stochastic common-spectrum process across 45 pulsars
[80]. Nonetheless, since the time residuals do not show the characteristic
angular correlation described by the Hellings–Downs curve [122], the detection
has not been claimed as GWs. Moreover, systematics such as solar system
effects [123] and pulsar spin noise [124] may affect the signal.

Figure 6: A fit to the NANOGrav data for $2\times 10^{13}~{}{\rm GeV}\lesssim
M_{DM}\lesssim 2\times 10^{14}~{}{\rm GeV}$ (red line). Points marked with
different shapes correspond to the different values of the SHDM mass reported
in the legenda to the right of the plot. The NANOGrav$@2\sigma$ range is well
fitted with a SHDM mass within $3\times 10^{13}~{}{\rm GeV}\lesssim
M_{DM}\lesssim 1\times 10^{14}~{}{\rm GeV}$.

Propitiously, if the signal is interpreted as GWs, cosmic strings provide an
excellent explanation of such a finding [82, 83, 84]. Let us mention that,
while the NANOGrav 12.5 yr data [80] is nearly consistent with previous EPTA
data [125] and with a very recent analysis of the PPTA data set [81], the
older NANOGrav 11 yrs data set [126] is in tension with the new one. Ref. [80]
mentions that, with the older data, the tension would be reduced once an
improved prior for the intrinsic pulsar red noise will be used. The NANOGrav
common-spectrum process is expressed in terms of a power-law signal with
characteristic strain given by

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

with $f_{yr}=1~{}{\rm yr}^{-1}$, $A$ and $\gamma_{N}$ being the characteristic
strain amplitude and the timing-residual cross-power spectral index
($\gamma_{N}=13/3$ for super-massive black hole mergers), respectively. The
abundance of GWs can be recast in the standard form as:

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

The power law approximation was fitted to 5 bins covering approximately the
frequency range $f\in\left[2.5\times 10^{-9},1.2\times 10^{-8}\right]$ Hz,
while higher-frequency bins in this data set are dominated by noise. The
NANOGrav $1\sigma$ and $2\sigma$ contours are given as the light and dark
green regions in Fig. 6. In the relevant frequency range, approximating the
GWs from cosmic strings as a power law, we perform a logarithmic fit
(represented by the red line in Fig. 6) similar to that of Refs. [82, 83, 84]
for $G\mu\sim Gv^{2}_{\Phi}\sim GM_{DM}^{2}$ within the DM mass range
$M_{DM}\in\left[2\times 10^{14}~{}\text{GeV},2\times
10^{13}~{}\text{GeV}\right]$. Each point, marked a different shape,
corresponds to a different benchmark value of the DM mass, as shown in the
legenda. The interesting fact is that at $2\sigma$, the data constrains the
SHDM mass to be within the range $3\times 10^{13}~{}{\rm GeV}\lesssim
M_{DM}\lesssim 10^{14}~{}{\rm GeV}$; this implies that the turning-point
frequency is bound to be within $\sim 0.2-0.4$ Hz (from Eq.IV.2 and Eq.IV.3).
Thus, whether the GWs possibly detected by NANOGrav come from the phase
transition that gave a mass to the SHDM will be confirmed (or excluded) by the
upcoming detectors LISA and DECIGO. Note also that such a spectral behaviour
also distinguishes this scenario from those of Refs. [82, 83, 84] where at
high-frequency the spectrum is flat ($\Omega_{GW}\sim f^{0}$).

##  V  Discussion and extensions

The predicted range of $f_{*}$ in this ${\rm PBH\rightarrow SHDM}$ mechanism
could be obtained with a high-scale $U(1)$ breaking and a well motivated early
matter dominated phase. For instance, in Refs. [77, 127], with a suitable
choice of $G\mu$ and the decay width ($\Gamma_{D}$) of the long-lived field,
the signal shown in Fig.5 (left) can be reproduced. Therefore to more robustly
differentiate our scenario from models such as Refs. [77, 127], we should look
for additional signatures. In the simplest scenario, i.e., without going into
the detail of a larger symmetry group such as GUT, one such signature could
arise both from the DM and the PBH sectors of our model. Firstly, for certain
interaction properties of DM, ultra-high-energy cosmic rays can play a crucial
role to test or constrain the model (see appendix E for a constraint on the
SHDM mass and life-time). On the other hand, there could be further GWs from
the PBH sector which could make the signatures of the model unique. The most
relevant example is GWs from the PBH density fluctuations that are very strong
in amplitude for large $\beta$, which is a value that is preferred in our set-
up as discussed in Sec.III.

Gravitational waves from PBH density fluctuations In a PBH dominated universe,
strong and sharply peaked gravitational waves arise due to the inhomogeneous
distribution of PBHs, as discussed recently in Ref. [10] and further developed
in Refs. [11, 12]. The present-day energy density of such GWs is given by [11,
12]

$\displaystyle\Omega_{GW}(t_{0},f)^{I}\simeq\Omega_{GW}^{\rm
peak}\left(\frac{f}{f_{\rm peak}}\right)^{5}\Theta\left(f_{\rm
peak}-f\right),$ (V.1)

where

$\displaystyle\Omega_{GW}^{\rm peak}\simeq 2\times
10^{-6}\left(\frac{\beta}{10^{-8}}\right)^{16/3}\left(\frac{M_{BH}}{10^{7}\rm
g}\right)^{34/9},~{}f_{\rm peak}\simeq 1.7\times 10^{3}{\rm
Hz}\left(\frac{M_{BH}}{10^{4}\rm g}\right)^{-5/6}.$ (V.2)

Two important points should be noted about Eq.V.2. First, not to saturate the
BBN bound on the effective number of neutrino species [103], the initial PBH
energy density fraction parameter $\beta$ can not be arbitrarily large as
shown in Fig.2, where the solid blue line corresponds to $\beta_{\rm
max}(M_{BH})$ [11]. Second, for the PBH mass range
$M_{BH}~{}\epsilon~{}\left[10^{6}{\rm g}-10^{8}{\rm g}\right]$ — an allowed
window corresponding to SHDM production (see Fig.3), the peak frequency
$f_{\rm peak}$ lies within the testable range $f~{}\epsilon~{}\left[36{\rm
Hz}\sim 0.7{\rm Hz}\right]$. Since, once the PBHs dominate, the number density
of SHDM is independent of $\beta$, the prediction of the exact peak amplitude
of $\Omega_{GW}(t_{0},f)^{I}$ requires additional input. Nonetheless, as the
location of the turning point frequency $f_{*}$ is derived considering the
fact that the entire DM relic density is in the form of SHDMs from PBHs, any
other SHDM relics that may originate e.g., from gravitational production
mechanisms [46, 47] should be diluted by considering large $\beta<\beta_{\rm
max}$. Therefore, in this set-up, particularly for light SHDMs (e.g.,
$M_{DM}<10^{12}$ GeV), a spectral bump on the flat-plateau is strongly
expected. Conversely, the non-observation of such a peak would put further
upper bounds on $N_{DM}^{i}$ in the $\rm PBH\rightarrow SHDM$ scenario. As a
matter of fact, as discussed briefly in the next section, such peaked
$\beta$-dependent GWs are not only relevant to constrain $N_{DM}^{i}$, but
also useful to test or constrain leptogenesis mechanisms along with SHDM in a
more well-motivated model.

Lastly, let us point out a remarkable coincidence that in this scenario there
is a significant overlap of the predicted ranges of $f_{*}$ and $f_{\rm
peak}$. In general $f_{*}$ and $f_{\rm peak}$ do not coincide and vary in
opposite directions as $M_{DM}$ increases. This is shown in Fig.7 with the
blue curves. Interestingly, at $M_{DM}\simeq 2.2\times 10^{11}$ GeV and
$f\simeq 2$Hz, $f_{*}$ and $f_{\rm peak}$ coincide and the overall
($\Omega_{GW}^{CS}[f]+\Omega_{GW}^{I}[f]$) GW spectral shape becomes unique.

Figure 7: $M_{DM}$ vs. $f_{*,~{}{\rm peak}}$. The blue solid line represents
the turning point frequencies of the GWs from cosmic strings whereas the blue
dashed line represents the peak frequency of the GWs induced by PBH density
fluctuation. The black dashed lines correspond to the cases where the scale of
symmetry breaking and the SHDM mass are different. The red circle marks the
frequency $f=f_{*}=f_{\rm peak}$.

Leptogenesis Our scenario can be embedded in the well-motivated anomaly-free
$U(1)_{B-L}$ extension of the SM (${\rm SM_{U(1)_{B-L}}}$) [128, 129, 130]
with three right handed (RH) neutrinos. A comprehensive summary of the
appearance of $U(1)_{B-L}$ in GUT models can be found in [76]. The inclusion
of RH neutrinos in the SM solves the problem of non-vanishing neutrino mass
via the seesaw mechanism [131, 132]; this model can also explain the baryon
asymmetry of the Universe via leptogenesis [133].

The recent articles [84, 76, 77, 78] consider the ${\rm SM_{U(1)_{B-L}}}$
model where all the three RH neutrinos decay quickly and show how the GWs from
cosmic strings from the $U(1)$ symmetry breaking can be a probe of
leptogenesis. However, this RH neutrino extension can be made more interesting
by taking one of them to be rigorously stable due to an extra parity symmetry
or stable on cosmological scale because of its very feeble interactions with
the SM: such RH neutrino then qualifies as a DM candidate. The model can still
explain the neutrino oscillation data as well as baryogenesis via leptogenesis
with the other two RH neutrinos [134, 135]. A number of works (e.g., see [136]
and references therein) are dedicated to this possibility because a three-
fermion SM extension solves three issues, namely, the neutrino masses, the
baryon asymmetry of the Universe and the DM.

Whereas the majority of these works consider freeze-out or freeze-in
production mechanisms for light DM particles, in our set-up the RH neutrino DM
has to be super-heavy. As shown in Ref. [79], if the RH neutrinos that produce
the lepton/baryon asymmetry are of PBH origin, then for a successful
leptogenesis one finds an upper bound on the PBH mass that reads
$M_{BH}\lesssim 13$ g; this is obviously not possible in our ${\rm
PBHs\rightarrow SHDM}$ scenario (cf. Fig. 3). Therefore, we need to consider
thermal leptogeneis. In the simplest scenario in which
$\Lambda_{CS}=M_{DM}=M_{R}^{2}\simeq M_{R}^{1}$, where $M_{R}^{1,2}$ are the
RH neutrino masses, leptogensis will be realised in the resonant regime [137].
Defining the next-to-the-lightest RH neutrino mass
$M_{R}^{2}=M_{R}^{1}\left(1+\delta\right)$ with $\delta\ll 1$, the baryon
asymmetry can be obtained as $N_{B-L}\simeq 6\times
10^{-17}\kappa\frac{\Lambda_{CS}}{\delta}$, where $\kappa\sim 10^{-2}$ is the
efficiency of lepton asymmetry production [137]. Let us consider a benchmark
value $\Lambda_{CS}=2.2\times 10^{11}$ as indicated with the red circle in
Fig.7. With this choice, $N_{B-L}$ becomes $N_{B-L}\simeq 1.4\times
10^{-7}/\delta$. The final baryon asymmetry is then $N_{B-L}^{f}\simeq
1.4\times 10^{-7}\Delta\left(\beta\right)/\delta$, where
$\Delta\left(\beta\right)$ is an entropy dilution factor that originates due
to PBH domination. One has to compare $N_{B-L}^{f}$ with the observed value
$N_{B-L}^{Obs}\simeq 10^{-8}$ [42]. Since $\delta\ll 1$, the overproduced
$N_{B-L}^{f}$ should be brought down to the observed one considering large
$\beta$ and therefore small $\Delta\left(\beta\right)$ (see Fig.4). Recall
that, as $\beta$ increases the peak amplitude of $\Omega_{GW}(t_{0},f)^{I}$
also increases. Therefore, in this model there is an explicit relation between
baryon asymmetry and $\Omega_{GW}(t_{0},f)^{I}$. From the observed baryon
asymmetry, as the quasi-degeneracy between the RH neutrinos becomes strong
($\delta$ becomes smaller), the amplitude of $\Omega_{GW}(t_{0},f)^{I}$ grows.
A detailed discussion of leptogenesis and SHDM with $U(1)_{B-L}$ in a PBH
dominated universe will be provided in our future publication.

##  VI Summary and outlook

In this paper we consider a novel way to test the possibility that dark
matter, specifically super-heavy dark matter with masses $M_{DM}\gtrsim
10^{10}$ GeV, is generated in the decay of primordial black holes
$M_{BH}\lesssim 10^{9}$ g that evaporate before Big Bang Nucleosynthesis. We
show how this possibility can be tested with gravitational waves, which have a
strong amplitude because the dark matter is super-heavy. We explicitly build
an example case in which the dark matter mass is generated by the breaking of
a gauged $U(1)$ symmetry, which brings about a network of cosmic strings; it
is these cosmic strings that are the sources of the strong gravitational waves
that are within the sensitivity reach of current and planned gravitational
waves detectors at nearly all frequencies.

This model possesses a further, unique feature that allows it to be
distinguished from other scenarios. The primordial black holes, in addition to
acting as a very natural source of super-heavy dark matter, cause a break in
the gravitational wave spectrum at higher frequencies wherein the spectral
slope changes from $\Omega_{GW}\sim f^{0}$ to $\Omega_{GW}\sim f^{-1/3}$
beyond a turning-point frequency $f_{*}$. Given a measured amplitude
($\Omega_{GW}\sim M_{DM}$) at low frequencies, the model forecasts a spectral
break $f_{*}\sim M_{DM}^{-2/5}$ at high frequency. Interestingly, if the dark
matter mass is in the range $10^{11}$ GeV $\lesssim M_{DM}\lesssim 10^{15}$
GeV, the spectral break $f_{*}$ would be precisely within the sensitivity
ranges of mid-band detectors such as BBO and DECIGO (and possibly LISA). Thus
the model exhibits a unique way to detect SHDM indirectly via the low
frequency-high frequency complementary measurements of GW spectral properties.

We speculate that the recent finding of a stochastic common-spectrum process
by NANOGrav and PPTA could be explained by the gravitational wave signal
expected by a super-heavy dark matter with mass of about $3\times
10^{13}~{}\text{GeV}\lesssim M_{DM}\lesssim 10^{14}~{}\text{GeV}$. Espousing
this interpretation, we predict that the gravitational wave spectrum exhibits
a spectral break $f_{*}$ in the $\mathcal{O}(0.1)$ Hz range, and that this
would be detectable by BBO and DECIGO. Therefore, whether dark matter is
super-heavy and produced by the evaporation of primordial black holes would be
tested by gravitational waves across nano-Hertz to Hertz frequencies, by
probing both the low-frequency extended plateau and cross-checking this with
the mid-frequency spectral break. We discuss how further GWs originating from
the density fluctuations of PBHs makes our scenario even more unique.

As possible extensions of this work, we briefly discussed a simple embedding
into a $U(1)_{B-L}$ model with three right-handed neutrinos, detail its
prospects for also including a realistic explanation for the neutrino masses
and baryogenesis. Furthermore, in this work we consider the effect of the PBHs
on the GW spectrum from cosmic strings at the level of the background
expansion of the Universe. However, the cosmic strings could potentially
interact with the PBHs and form a black hole-string network that may lead to
further spectral distortions in the GW spectrum [138]. A precise evolution of
such a network would require numerical simulations, see [139].

## acknowledgement

RS is supported by the MSCA-IF IV FZU - CZ.02.2.69/0.0/0.0/20 079/0017754
project and acknowledges European Structural and Investment Fund and the Czech
Ministry of Education, Youth and Sports. FU is supported by the European
Regional Development Fund (ESIF/ERDF) and the Czech Ministry of Education,
Youth and Sports (MEYS) through Project CoGraDS -
CZ.02.1.01/0.0/0.0/$15\\_003$/0000437.

##  Appendix A Scaling loop number density and the plateau integral

The energy density parameter of gravitational waves at $t_{0}$ 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)

where [73]

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

The quantity $F=f\left[\frac{a(t_{0})}{a(t)}\right]$ accounts for the red-
shifting of the frequency and $P_{GW}(t,f_{k})$ corresponds to the power
emitted by the loops, defined as [73]

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

Integrating Eq.A.3 over the loop lengths gives

$\displaystyle
P_{GW}(t,f_{k})=\frac{2kG\mu^{2}\Gamma_{k}}{f_{k}^{2}}n(t,l_{k})=\frac{2kG\mu^{2}\Gamma_{k}}{f^{2}\left[\frac{a(t_{0})}{a(t)}\right]^{2}}n\left(t,\frac{2k}{f}\left[\frac{a(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_{osc}}^{t_{0}}\left[\frac{a(t)}{a(t_{0})}\right]^{5}n\left(t,\frac{2k}{f}\left[\frac{a(t)}{a(t_{0})}\right]\right)dt.$
(A.5)

Therefore the energy density for the $k$-th mode is given by

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

The loop number density
$n\left(t,l_{k}=\frac{2k}{f}\left[\frac{a(t)}{a(t_{0})}\right]\right)$ can be
calculated with the analytical Velocity-dependent-One-Scale (VOS) model [97,
98, 99] that assumes the loop production function to be a delta function,
i.e., all the loops are created with the same size. For a general equation of
state parameter $\omega$, the scaling number density
$n_{\omega}\left(t,l_{k}\right)$ can be computed as [99, 78]

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

where $\beta=2/3(1+\omega)$. The parameter
$A_{\beta}=5.4~{}(\text{for}~{}w=1/3)$ or
$A_{\beta}=0.39~{}(\text{for}~{}\omega=0)$ is constant and represents the loop
production efficiency in different cosmological epochs [99]. The constant
values of $A_{\beta}$ are subjected to the assumption that the scaling regime
is reached instantaneously when the network goes from one cosmological epoch
to another [96].

The analytic expression for the $\Omega_{GW}^{plt}$ can be derived by
analysing the fundamental mode ($k=1$) and considering loop production as well
as decay during radiation domination. A convenient way to obtain an expression
for $\Omega_{GW}^{plt}$ is to convert the $t$ integral in Eq.II.1 into a scale
factor integral as

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

where

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

and the loop number density
$n\left(t,l_{1}(t)\equiv\frac{2}{f}\left[\frac{a(t)}{a(t_{0})}\right]\right)$
in Eq.A.7 (in radiation domination) is expressed as

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

Let us mention also that the VOS model overestimates the number density of the
loops approximately by an order of magnitude compared to the numerical
simulations [73]. This is due to the fact that VOS model assumes that all the
loops are of same size (a fraction of the horizon) at their production. In
reality, there could be a distribution of $\alpha$. Numerical simulations find
that only 10$\%$ of the energy from the long-string network goes to large
loops ($\alpha\simeq 0.1$) while the remaining $90\%$ goes to highly boosted
smaller loops which do not radiate GWs. Therefore, to be consistent with
numerical simulations, we add the so-called calibration or normalisation
factor $\mathcal{F}_{\alpha}\sim 0.1$ in Eq.A.7 [99].

##  Appendix B Summing over all the modes

The contributions from all the Fourier $k$ modes can easily be accounted for
owing to the fact that

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

Specifically, the spectral shape beyond the turning-point frequency $f_{*}$
changes in a significant way. For example, expanding the RHS of Eq.B.1 for the
first few modes, i.e.,

$\displaystyle\Omega_{GW}(f)$ $\displaystyle=$
$\displaystyle\sum_{k}k^{-\delta}\Omega^{(1)}(f/k)$ (B.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 $m$, $n$ and $r$ obey $1<m<n<r$, we see that, should we
keep on increasing the number of modes, there will be a critical point
$k\equiv k_{\Delta}$ for which the amplitude
$\Omega_{GW}^{(1)}(f_{\Delta}=f/k_{\Delta})$ would contribute to the
$\Omega_{GW}(f)$. Therefore, the sum can be performed in two parts:

$\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})$
(B.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}$ (B.4)

That gives

$\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}$ (B.5)

upon using the asymptotic expansion of the Euler-Maclaurin series for the
first term and the expansion of the Hurwitz zeta function for the second term.

##  Appendix C PBH factsheet

Life-time, evaporation and condition for PBH domination: For the radiation
energy density to dominate at the time of PBH evaporation $t_{ev}$, one has

$\displaystyle r(t_{ev})\equiv\frac{\rho_{BH}(t_{ev})}{\rho_{\rm
R}(t_{ev})}<1.$ (C.1)

The ratio of the $r$ parameters at the time of black hole formation and
evaporation is therefore

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

Using the first 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}$ (C.3)

along with the expression for the Hubble parameter $H(t)$ in radiation
domination and the radiation energy density

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

where $g_{*}(T)$ ($\simeq$ 106.75 in SM) is the effective degrees of freedom
that contribute to the radiation. Eq.C.2 can be rewritten 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}}.$
(C.5)

Eq.C.1 now 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}}.$
(C.6)

The above condition on $\beta$ can be expressed fully in terms a single free
parameter $M_{BH}$—we work in the approximation that the PBH mass spectrum is
monochromatic. Assuming radiation domination at black hole formation, the mass
of a PBH originating from gravitational collapse is approximately given by the
energy density enclosed in a post-inflationary particle horizon [111], i.e.,

$\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}}.$
(C.7)

The quantity $\gamma\simeq 0.2$ depends on the details of the gravitational
collapse mechanism. From Eq.C.7, the PBH formation time $t_{Bf}$ is calculated
as

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

Now using Eq.C.3 and Eq.C.4 one obtains the PBH 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}.$
(C.9)

Let us, for now, assume that radiation dominates throughout the evolution of
the PBH—we will discuss the case in which PBHs come to dominate for a period
of time in Appendix D. The Hubble parameter is given by

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

where the $\tau$ is the lifetime of the PBH. Therefore, again, using Eq.C.3
and Eq.C.4, $T_{ev}$ can be obtained as

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

The lifetime of the PBH is obtained from the dynamics of the mass loss of a
PBH via Hawking evaporation [4]. The rate at which PBH loses mass is given by

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

where the $\frac{dE}{dt}$ is obtained 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}.$ (C.13)

The quantity $f_{ev}$ is the efficiency of PBH evaporation and $r_{BH}$ is the
Schwarzschild radius. The quantity $g_{BH}(T_{BH})$ counts the bosonic and
fermionic degrees of freedom for $T<T_{BH}$. Recalling the PBH temperature
($T_{BH}=M_{Pl}^{2}/8\pi M_{BH}$ [4]) and using $r_{BH}=2GM_{BH}$, Eq.C.12 is
recast as

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

where one uses $f_{ev}g_{BH}(T_{BH})=\mathcal{G}g_{*B}(T_{BH})$. The lifetime
$\tau$ is then obtained 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}}.$ (C.15)

Now combining Eq.C.15, Eq.C.11 and Eq.C.9, we can recast Eq.C.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}}.$
(C.16)

The Friedmann equations:

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

where the Eq.C.19 represents the non-conservation of entropy $\tilde{S}\sim
sa^{3}$ due to PBH evaporation.

##  Appendix D DM and PBH number densities

In a simplified approach (neglecting the momentum distribution of the emitted
particles), the differential number of a particle species $`X$’ emitted by a
black hole can be computed as [16]

$\displaystyle
d\bar{n}=dE/3T_{BH}=\frac{M_{Pl}^{2}}{24\pi}\frac{1}{T_{BH}^{3}}dT_{BH},$
(D.1)

where we have used

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

and the mean energy of the radiated particles is $\bar{E}=3T$. The total
number of $X$ particles emitted in the case of a complete evaporation is given
by

$\displaystyle\bar{n}_{X}=\frac{g_{X}}{g_{*B}}\int_{T_{BH}}^{\infty}d\bar{n}=\frac{4\pi}{3}\frac{g_{X}}{g_{*B}}\left(\frac{M_{BH}}{M_{Pl}}\right)^{2}~{}~{}{\rm
for}~{}~{}T_{BH}>M_{X},$ (D.3)
$\displaystyle\bar{n}_{X}=\frac{g_{X}}{g_{*B}}\int_{M_{X}}^{\infty}d\bar{n}=\frac{1}{48\pi}\frac{g_{X}}{g_{*B}}\left(\frac{M_{Pl}}{M_{X}}\right)^{2}~{}~{}{\rm
for}~{}~{}T_{BH}<M_{X},$ (D.4)

where $M_{X}$ and $g_{X}$ are the mass and internal degrees of freedom of the
radiated particles. The quantity $N_{\rm BH}^{\rm ev}$ is obtained as

$\displaystyle N_{\rm BH}^{\rm ev}=\left(\frac{n_{\rm BH}}{n_{f,\rm eq}^{\rm
ur}}\right)_{T_{ev}}=\frac{M_{pl}^{2}}{6\pi\tau^{2}M_{BH}}\left(\frac{g_{f}T_{ev}^{3}}{\pi^{2}}\right)^{-1},$
(D.5)

where we have used $n_{\rm BH}=\rho_{\rm BH}/M_{BH}$ with $\rho_{\rm BH}$
calculated from the Friedmann equation when PBHs dominate: $H\simeq 2/3\tau$.
Proceeding in the same way as in Eq.C.11, the evaporation temperature in the
case of PBH domination can be calculated as

$\displaystyle
T_{ev}=\left(\frac{5M_{Pl}^{2}}{\pi^{3}g_{*}(T_{ev})\tau^{2}}\right)^{1/4}.$
(D.6)

Combining Eq.D.5 and Eq.D.6 we find

$\displaystyle N_{\rm BH}^{\rm
ev}=\left(\frac{\pi}{6g_{f}}\right)\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}.$
(D.7)

##  Appendix E UHECR constraints

If the DM is indeed super-heavy and it is not exactly stable, it is possible
to test and constrain its properties using UHECRs data, see, e.g., [62, 61]
and references therein. Indeed, when SHDM particles decay, they initiate a
cascade which produces a slew of very high energy photons, neutrinos, and
protons, which can be detected in UHECR experimental facilities. The most
stringent constraints come from the non-detection of ultra-high-energy
photons, as they are more numerous than protons and much easier to detect than
neutrinos. For the range of masses we are interested in, the most significant
result comes from the non-observation of any UHECR, and therefore also
photons, at primary energies above $E=10^{11.3}~{}\text{GeV}$ at the Pierre
Auger Observatory [63], which translates in an upper limit on the integral
flux above this energy of

$\displaystyle\Phi_{\text{exp}}(E>10^{11.3}~{}\text{GeV})\lesssim\frac{3.6\cdot
10^{-5}}{\text{km}^{2}\,\text{sr}\,\text{yr}}\,.$ (E.1)

At lower energies the limits from the Pierre Auger Observatory and the
Telescope Array collaboration are similar and approximately two orders of
magnitude lower [140, 141].

The predicted flux of photons from SHDM decay instead can be written as [62,
61, 142]

$\displaystyle J(E)=\frac{1}{4\pi
M_{DM}\tau_{DM}}\frac{\mathrm{d}N_{\gamma}}{\mathrm{d}E}\int\mathrm{d}s\,\rho_{DM}(s)\,,$
(E.2)

where the integral of the SHDM density $\rho_{DM}(s)$ over the line of sight
$s$ can be written as an integral over Galactic radius $r$ as

$\displaystyle\int\mathrm{d}s\,\rho_{DM}(s)=\left[2\int_{r_{\odot}\sin\theta}^{r_{\odot}}+\int_{r_{\odot}}^{r_{H}}\right]\mathrm{d}r\,r\,\frac{\rho_{DM}(r)}{\sqrt{r^{2}-r_{\odot}^{2}\sin^{2}\theta}}\,,$
(E.3)

with the SHDM density itself—we adopt the Einasto profile here—being given by

$\displaystyle\rho_{DM}(r)\equiv\rho_{s}\exp\left\\{-\frac{2}{\alpha}\left[\left(\frac{r}{r_{s}}\right)^{\alpha}-1\right]\right\\}\,.$
(E.4)

In these expressions $\rho_{s}=0.033~{}\text{GeV}/\text{cm}^{3}$,
$r_{s}=28.44~{}\text{kpc}$, $\alpha=0.11$, $r_{H}=260~{}\text{kpc}$ and
$r_{\odot}=8.33~{}\text{kpc}$. The SHDM density integral Eq.E.3 depends on the
angle $\theta$ between the line of sight and the Sun-Galactic Centre axis; in
our numerical estimates we will average the photon flux over the whole sky for
simplicity.

The flux Eq.E.2 depends on the SHDM lifetime $\tau_{DM}$ and mass $M_{DM}$, as
well as on the differential number of photons produced by each SHDM decay

$\displaystyle\frac{\mathrm{d}N_{\gamma}}{\mathrm{d}E}\equiv\frac{2}{M_{DM}}\frac{\mathrm{d}N_{\gamma}}{\mathrm{d}x}\equiv\frac{2\hat{N}_{\gamma}}{M_{DM}}\left(\frac{x}{\hat{x}}\right)^{-\gamma}\,,$
(E.5)

with $\gamma=1.9$, $\hat{N}_{\gamma}\approx 10^{8}$ and $\hat{x}=10^{-5}$.
This simplification is within 10% (or better) of the true numerical result
provided the lower limit of the integral flux satisfies $x\equiv
2E/M_{DM}\lesssim 0.1$—this is certainly true for the UHECR energies and SHDM
masses we are considering in this work.

Integrating the flux Eq.E.2 above the energy cut
$E_{\text{cut}}=10^{11.3}~{}\text{GeV}$ gives us the predicted integral flux
$\Phi_{\gamma}$ which needs to satisfy the constraint Eq.E.1. For the SHDM
that are of most interest in this work we therefore obtain

$\displaystyle\tau_{DM}\gtrsim
10^{24}\left(\frac{M_{DM}}{\text{GeV}}\right)^{-0.1}~{}\text{yr}\,.$ (E.6)

For lower values of the mass $M_{DM}\in[10^{10},10^{12}]~{}\text{GeV}$ the
limit reads $\tau_{DM}\gtrsim 10^{21}~{}\text{yr}$.

## 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] S. W. Hawking, Commun. Math. Phys. 43, 199 (1975) Erratum: [Commun. Math. Phys. 46, 206 (1976)]. doi:10.1007/BF02345020, 10.1007/BF01608497
  * [5] 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]].
  * [6] A. D. Dolgov and D. Ejlli, Phys. Rev. D 84, 024028 (2011) doi:10.1103/PhysRevD.84.024028 [arXiv:1105.2303 [astro-ph.CO]].
  * [7] 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]].
  * [8] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott and R. Petrossian-Byrne, [arXiv:2004.00618 [astro-ph.CO]].
  * [9] 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]].
  * [10] T. Papanikolaou, V. Vennin and D. Langlois, arXiv:2010.11573 [astro-ph.CO].
  * [11] G. Domènech, C. Lin and M. Sasaki, arXiv:2012.08151 [gr-qc].
  * [12] G. Domènech, V. Takhistov and M. Sasaki, Phys. Lett. B 823, 136722 (2021) doi:10.1016/j.physletb.2021.136722 [arXiv:2105.06816 [astro-ph.CO]].
  * [13] G. Domènech, [arXiv:2109.01398 [gr-qc]].
  * [14] 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]].
  * [15] 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
  * [16] D. Baumann, P. J. Steinhardt and N. Turok, hep-th/0703250 [HEP-TH].
  * [17] A. Hook, Phys. Rev. D 90, no. 8, 083535 (2014) doi:10.1103/PhysRevD.90.083535 [arXiv:1404.0113 [hep-ph]].
  * [18] Y. Hamada and S. Iso, PTEP 2017, no. 3, 033B02 (2017) doi:10.1093/ptep/ptx011 [arXiv:1610.02586 [hep-ph]].
  * [19] D. Hooper and G. Krnjaic, arXiv:2010.01134 [hep-ph].
  * [20] Y. F. Perez-Gonzalez and J. Turner, arXiv:2010.03565 [hep-ph].
  * [21] G. Aliferis and V. Zarikas, Phys. Rev. D 103, no.2, 023509 (2021) doi:10.1103/PhysRevD.103.023509 [arXiv:2006.13621 [gr-qc]].
  * [22] A. G. Polnarev and M. Y. Khlopov, Sov. Phys. Usp. 28, 213 (1985) [Usp. Fiz. Nauk 145, 369 (1985)]. doi:10.1070/PU1985v028n03ABEH003858
  * [23] 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].
  * [24] 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]].
  * [25] 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]].
  * [26] 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]].
  * [27] D. Hooper, G. Krnjaic and S. D. McDermott, JHEP 1908, 001 (2019) doi:10.1007/JHEP08(2019)001 [arXiv:1905.01301 [hep-ph]].
  * [28] C. Lunardini and Y. F. Perez-Gonzalez, JCAP 2008, 014 (2020) doi:10.1088/1475-7516/2020/08/014 [arXiv:1910.07864 [hep-ph]].
  * [29] A. Chaudhuri and A. Dolgov, arXiv:2001.11219 [astro-ph.CO].
  * [30] I. Masina, Eur. Phys. J. Plus 135, no. 7, 552 (2020) doi:10.1140/epjp/s13360-020-00564-9 [arXiv:2004.04740 [hep-ph]].
  * [31] 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]].
  * [32] 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]].
  * [33] N. Bernal and Ó. Zapata, arXiv:2011.12306 [astro-ph.CO].
  * [34] J. Auffinger, I. Masina and G. Orlando, arXiv:2012.09867 [hep-ph].
  * [35] T. Kitabayashi, [arXiv:2101.01921 [hep-ph]].
  * [36] S. Jyoti Das, D. Mahanta and D. Borah, [arXiv:2104.14496 [hep-ph]].
  * [37] A. Cheek, L. Heurtier, Y. F. Perez-Gonzalez and J. Turner, [arXiv:2107.00013 [hep-ph]].
  * [38] T. Kitabayashi, [arXiv:2107.11692 [hep-ph]].
  * [39] N. Bernal, F. Hajkarim and Y. Xu, [arXiv:2107.13575 [hep-ph]].
  * [40] P. Sandick, B. S. Es Haghi and K. Sinha, [arXiv:2108.08329 [astro-ph.CO]].
  * [41] B. Barman, D. Borah, S. J. Das and R. Roshan, [arXiv:2111.08034 [hep-ph]].
  * [42] P. A. R. Ade et al. [Planck], Astron. Astrophys. 594, A13 (2016) doi:10.1051/0004-6361/201525830 [arXiv:1502.01589 [astro-ph.CO]].
  * [43] P. Bull, Y. Akrami, J. Adamek, T. Baker, E. Bellini, J. Beltran Jimenez, E. Bentivegna, S. Camera, S. Clesse and J. H. Davis, et al. Phys. Dark Univ. 12, 56-99 (2016) doi:10.1016/j.dark.2016.02.001 [arXiv:1512.05356 [astro-ph.CO]].
  * [44] C. Skordis and T. Zlosnik, Phys. Rev. Lett. 127, no.16, 161302 (2021) doi:10.1103/PhysRevLett.127.161302 [arXiv:2007.00082 [astro-ph.CO]].
  * [45] C. Skordis and T. Zlosnik, [arXiv:2109.13287 [gr-qc]].
  * [46] D. J. H. Chung, E. W. Kolb and A. Riotto, Phys. Rev. D 59, 023501 (1998) doi:10.1103/PhysRevD.59.023501 [arXiv:hep-ph/9802238 [hep-ph]].
  * [47] V. Kuzmin and I. Tkachev, JETP Lett. 68, 271-275 (1998) doi:10.1134/1.567858 [arXiv:hep-ph/9802304 [hep-ph]].
  * [48] E. W. Kolb, D. J. H. Chung and A. Riotto, AIP Conf. Proc. 484, no.1, 91-105 (1999) doi:10.1063/1.59655 [arXiv:hep-ph/9810361 [hep-ph]].
  * [49] D. J. H. Chung, P. Crotty, E. W. Kolb and A. Riotto, Phys. Rev. D 64, 043503 (2001) doi:10.1103/PhysRevD.64.043503 [arXiv:hep-ph/0104100 [hep-ph]].
  * [50] S. Chang, C. Coriano and A. E. Faraggi, Nucl. Phys. B 477, 65-104 (1996) doi:10.1016/0550-3213(96)00371-9 [arXiv:hep-ph/9605325 [hep-ph]].
  * [51] M. Chianese and S. F. King, JCAP 09, 027 (2018) doi:10.1088/1475-7516/2018/09/027 [arXiv:1806.10606 [hep-ph]].
  * [52] M. Chianese, B. Fu and S. F. King, JCAP 03, 030 (2020) doi:10.1088/1475-7516/2020/03/030 [arXiv:1910.12916 [hep-ph]].
  * [53] A. Berlin, Phys. Rev. Lett. 119, 121801 (2017) doi:10.1103/PhysRevLett.119.121801 [arXiv:1704.08256 [hep-ph]].
  * [54] H. Kim and E. Kuflik, Phys. Rev. Lett. 123, no.19, 191801 (2019) doi:10.1103/PhysRevLett.123.191801 [arXiv:1906.00981 [hep-ph]].
  * [55] E. D. Kramer, E. Kuflik, N. Levi, N. J. Outmezguine and J. T. Ruderman, Phys. Rev. Lett. 126, no.8, 081802 (2021) doi:10.1103/PhysRevLett.126.081802 [arXiv:2003.04900 [hep-ph]].
  * [56] L. Bian, X. Liu and K. P. Xie, JHEP 11, 175 (2021) doi:10.1007/JHEP11(2021)175 [arXiv:2107.13112 [hep-ph]].
  * [57] A. Azatov, M. Vanvlasselaer and W. Yin, JHEP 03, 288 (2021) doi:10.1007/JHEP03(2021)288 [arXiv:2101.05721 [hep-ph]].
  * [58] M. Ahmadvand, JHEP 10, 109 (2021) doi:10.1007/JHEP10(2021)109 [arXiv:2108.00958 [hep-ph]].
  * [59] E. Babichev, D. Gorbunov and S. Ramazanov, Phys. Lett. B 794, 69-76 (2019) doi:10.1016/j.physletb.2019.05.030 [arXiv:1812.03516 [hep-ph]].
  * [60] E. Babichev, D. Gorbunov, S. Ramazanov and L. Reverberi, JCAP 09, 059 (2020) doi:10.1088/1475-7516/2020/09/059 [arXiv:2006.02225 [hep-ph]].
  * [61] O. K. Kalashev and M. Y. Kuznetsov, Phys. Rev. D 94 (2016) no.6, 063535 doi:10.1103/PhysRevD.94.063535 [arXiv:1606.07354 [astro-ph.HE]].
  * [62] L. Marzola and F. R. Urban, Astropart. Phys. 93 (2017), 56-69 doi:10.1016/j.astropartphys.2017.04.005 [arXiv:1611.07180 [astro-ph.HE]].
  * [63] E. Alcantara, L. A. Anchordoqui and J. F. Soriano, Phys. Rev. D 99 (2019) no.10, 103016 doi:10.1103/PhysRevD.99.103016 [arXiv:1903.05429 [hep-ph]].
  * [64] A. Mazumdar and G. White, Rept. Prog. Phys. 82, no.7, 076901 (2019) doi:10.1088/1361-6633/ab1f55 [arXiv:1811.01948 [hep-ph]].
  * [65] T. W. B. Kibble, J. Phys. A 9, 1387-1398 (1976) doi:10.1088/0305-4470/9/8/029
  * [66] H. B. Nielsen and P. Olesen, Nucl. Phys. B 61, 45-61 (1973) doi:10.1016/0550-3213(73)90350-7
  * [67] R. Jeannerot, J. Rocher and M. Sakellariadou, Phys. Rev. D 68, 103514 (2003) doi:10.1103/PhysRevD.68.103514 [arXiv:hep-ph/0308134 [hep-ph]].
  * [68] A. Vilenkin, Phys. Lett. B 107, 47-50 (1981) doi:10.1016/0370-2693(81)91144-8
  * [69] N. Turok, Nucl. Phys. B 242, 520-541 (1984) doi:10.1016/0550-3213(84)90407-3
  * [70] T. Vachaspati and A. Vilenkin, Phys. Rev. D 31, 3052 (1985) doi:10.1103/PhysRevD.31.3052
  * [71] C. Ringeval, M. Sakellariadou and F. Bouchet, JCAP 02, 023 (2007) doi:10.1088/1475-7516/2007/02/023 [arXiv:astro-ph/0511646 [astro-ph]].
  * [72] 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]].
  * [73] 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]].
  * [74] 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]].
  * [75] Y. Cui, M. Lewicki, D. E. Morrissey and J. D. Wells, JHEP 01, 081 (2019) doi:10.1007/JHEP01(2019)081 [arXiv:1808.08968 [hep-ph]].
  * [76] 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]].
  * [77] S. Blasi, V. Brdar and K. Schmitz, Phys. Rev. Res. 2, no.4, 043321 (2020) doi:10.1103/PhysRevResearch.2.043321 [arXiv:2004.02889 [hep-ph]].
  * [78] R. Samanta and S. Datta, JHEP 11, 017 (2021) doi:10.1007/JHEP11(2021)017 [arXiv:2108.08359 [hep-ph]].
  * [79] S. Datta, A. Ghosal and R. Samanta, JCAP 08, 021 (2021) doi:10.1088/1475-7516/2021/08/021 [arXiv:2012.14981 [hep-ph]].
  * [80] Z. Arzoumanian et al. [NANOGrav], Astrophys. J. Lett. 905, no.2, L34 (2020) doi:10.3847/2041-8213/abd401 [arXiv:2009.04496 [astro-ph.HE]].
  * [81] B. Goncharov, R. M. Shannon, D. J. Reardon, G. Hobbs, A. Zic, M. Bailes, M. Curylo, S. Dai, M. Kerr and M. E. Lower, et al. doi:10.3847/2041-8213/ac17f4 [arXiv:2107.12112 [astro-ph.HE]].
  * [82] J. Ellis and M. Lewicki, Phys. Rev. Lett. 126, no.4, 041304 (2021) doi:10.1103/PhysRevLett.126.041304 [arXiv:2009.06555 [astro-ph.CO]].
  * [83] S. Blasi, V. Brdar and K. Schmitz, Phys. Rev. Lett. 126, no.4, 041305 (2021) doi:10.1103/PhysRevLett.126.041305 [arXiv:2009.06607 [astro-ph.CO]].
  * [84] R. Samanta and S. Datta, JHEP 05, 211 (2021) doi:10.1007/JHEP05(2021)211 [arXiv:2009.13452 [hep-ph]]. l
  * [85] W. Buchmüller, V. Domcke, K. Kamada and K. Schmitz, JCAP 10, 003 (2013) doi:10.1088/1475-7516/2013/10/003 [arXiv:1305.3392 [hep-ph]].
  * [86] 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]].
  * [87] M. B. Hindmarsh and T. W. B. Kibble, Rept. Prog. Phys. 58, 477-562 (1995) doi:10.1088/0034-4885/58/5/001 [arXiv:hep-ph/9411342 [hep-ph]].
  * [88] A.VilenkinandE.P.S.Shellard,CosmicStringsandOtherTopologicalDefects(CambridgeUniversity Press,2000).
  * [89] A. Vilenkin, Phys. Rev. D 43, 1060-1062 (1991) doi:10.1103/PhysRevD.43.1060
  * [90] D. P. Bennett and F. R. Bouchet, Phys. Rev. Lett. 60, 257 (1988).
  * [91] D. P. Bennett and F. R. Bouchet, Phys. Rev. Lett. 63, 2776 (1989).
  * [92] A. Albrecht and N. Turok, Phys. Rev. D 40, 973-1001 (1989).
  * [93] T. Damour and A. Vilenkin, Phys. Rev. D 64, 064008 (2001).
  * [94] D. Matsunami, L. Pogosian, A. Saurabh and T. Vachaspati, Phys. Rev. Lett. 122, no.20, 201301 (2019).
  * [95] P. Auclair, D. A. Steer and T. Vachaspati, Phys. Rev. D 101, no.8, 083511 (2020).
  * [96] Y. Gouttenoire, G. Servant and P. Simakachorn, JCAP 07, 032 (2020).
  * [97] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 54, 2535-2556 (1996).
  * [98] C. J. A. P. Martins and E. P. S. Shellard, Phys. Rev. D 65, 043514 (2002).
  * [99] P. Auclair, J. J. Blanco-Pillado, D. G. Figueroa, A. C. Jenkins, M. Lewicki, M. Sakellariadou, S. Sanidas, L. Sousa, D. A. Steer and J. M. Wachter, et al. JCAP 04, 034 (2020).
  * [100] P. Amaro-Seoane et al. [LISA], [arXiv:1702.00786 [astro-ph.IM]].
  * [101] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 6, 061102 (2016).
  * [102] B. P. Abbott et al. [LIGO Scientific and Virgo Collaborations], Phys. Rev. Lett. 116, no. 24, 241103 (2016).
  * [103] R. H. Cyburt, B. D. Fields, K. A. Olive and T. H. Yeh, Rev. Mod. Phys. 88, 015004 (2016)
  * [104] M. Kawasaki, K. Kohri and N. Sugiyama, Phys. Rev. Lett. 82, 4168 (1999).
  * [105] T. Hasegawa, N. Hiroshima, K. Kohri, R. S. L. Hansen, T. Tram and S. Hannestad, JCAP 12, 012 (2019).
  * [106] M. Gorghetto, E. Hardy and H. Nicolaescu, JCAP 06, 034 (2021) doi:10.1088/1475-7516/2021/06/034 [arXiv:2101.11007 [hep-ph]].
  * [107] C. F. Chang and Y. Cui, [arXiv:2106.09746 [hep-ph]].
  * [108] G. S. F. Guedes, P. P. Avelino and L. Sousa, Phys. Rev. D 98, no.12, 123505 (2018) doi:10.1103/PhysRevD.98.123505 [arXiv:1809.10802 [astro-ph.CO]].
  * [109] 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]].
  * [110] W. T. Emond, S. Ramazanov and R. Samanta, [arXiv:2108.05377 [hep-ph]].
  * [111] 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]].
  * [112] Y. Akrami et al. [Planck Collaboration], Astron. Astrophys. 641, A10 (2020) doi:10.1051/0004-6361/201833887 [arXiv:1807.06211 [astro-ph.CO]].
  * [113] J. H. MacGibbon, Phys. Rev. D 44, 376 (1991). doi:10.1103/PhysRevD.44.376
  * [114] G. F. Giudice, E. W. Kolb and A. Riotto, Phys. Rev. D 64, 023508 (2001) doi:10.1103/PhysRevD.64.023508 [hep-ph/0005123].
  * [115] J. Hisano, K. Kohri and M. M. Nojiri, Phys. Lett. B 505, 169 (2001) doi:10.1016/S0370-2693(01)00395-1 [hep-ph/0011216].
  * [116] F. Takahashi, Phys. Lett. B 660, 100 (2008) doi:10.1016/j.physletb.2007.12.048 [arXiv:0705.0579 [hep-ph]].
  * [117] R. Samanta, A. Biswas and S. Bhattacharya, JCAP 01, 055 (2021) doi:10.1088/1475-7516/2021/01/055 [arXiv:2006.02960 [hep-ph]].
  * [118] M. Viel, G. D. Becker, J. S. Bolton and M. G. Haehnelt, Phys. Rev. D 88, 043502 (2013) doi:10.1103/PhysRevD.88.043502 [arXiv:1306.2314 [astro-ph.CO]].
  * [119] N. Seto, S. Kawamura and T. Nakamura, Phys. Rev. Lett. 87, 221103 (2001) doi:10.1103/PhysRevLett.87.221103 [arXiv:astro-ph/0108011 [astro-ph]].
  * [120] V. Corbin and N. J. Cornish, Class. Quant. Grav. 23, 2435-2446 (2006) doi:10.1088/0264-9381/23/7/014 [arXiv:gr-qc/0512039 [gr-qc]].
  * [121] C. Caprini, D. G. Figueroa, R. Flauger, G. Nardini, M. Peloso, M. Pieroni, A. Ricciardone and G. Tasinato, JCAP 11, 017 (2019) doi:10.1088/1475-7516/2019/11/017 [arXiv:1906.09244 [astro-ph.CO]].
  * [122] R. w. Hellings and G. s. Downs, Astrophys. J. Lett. 265, L39 (1983). doi:10.1086/183954
  * [123] 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].
  * [124] 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]].
  * [125] 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]].
  * [126] 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]].
  * [127] D. Borah, S. J. Das, A. K. Saha and R. Samanta, [arXiv:2202.10474 [hep-ph]].
  * [128] A. Davidson, Phys. Rev. D 20, 776 (1979).
  * [129] R. E. Marshak and R. N. Mohapatra, Phys. Lett. 91B, 222 (1980).
  * [130] R. N. Mohapatra and R. E. Marshak, Phys. Rev. Lett. 44, 1316 (1980) Erratum: [Phys. Rev. Lett. 44, 1643 (1980)].
  * [131] P. Minkowski, Phys. Lett. B 67, 421-428 (1977) doi:10.1016/0370-2693(77)90435-X
  * [132] T. Yanagida, Prog. Theor. Phys. 64, 1103 (1980) doi:10.1143/PTP.64.1103
  * [133] M. Fukugita and T. Yanagida, Phys. Lett. B 174, 45-47 (1986) doi:10.1016/0370-2693(86)91126-3
  * [134] S. F. King, Nucl. Phys. B 576, 85-105 (2000) doi:10.1016/S0550-3213(00)00109-7 [arXiv:hep-ph/9912492 [hep-ph]].
  * [135] R. Samanta and M. Sen, JHEP 01, 193 (2020) doi:10.1007/JHEP01(2020)193 [arXiv:1908.08126 [hep-ph]].
  * [136] N. Okada and O. Seto, Phys. Rev. D 82, 023507 (2010) doi:10.1103/PhysRevD.82.023507 [arXiv:1002.2525 [hep-ph]].
  * [137] A. Pilaftsis and T. E. J. Underwood, Nucl. Phys. B 692, 303-345 (2004) doi:10.1016/j.nuclphysb.2004.05.029 [arXiv:hep-ph/0309342 [hep-ph]].
  * [138] A. Vilenkin, Y. Levin and A. Gruzinov, JCAP 11, 008 (2018) doi:10.1088/1475-7516/2018/11/008 [arXiv:1808.00670 [astro-ph.CO]].
  * [139] `https://www.youtube.com/watch?v=QWI5F8--jLM`
  * [140] P. Abreu et al. [Pierre Auger], PoS ICRC2021 (2021), 373 doi:10.22323/1.395.0373
  * [141] R. Abbasi et al. [Telescope Array], PoS ICRC2021 (2021), 864 doi:10.22323/1.395.0864
  * [142] M. Chianese, D. F. G. Fiorillo, R. Hajjar, G. Miele and N. Saviano, JCAP 11, 035 (2021) doi:10.1088/1475-7516/2021/11/035 [arXiv:2108.01678 [hep-ph]].