# Gravitational Wave Sourced by Decay of Massive Particle from Primordial
Black Hole evaporation

Ki-Young Choi  Erdenebulgan Lkhagvadorj  and Satyabrata Mahapatra

###### Abstract

In this article, we investigate the stochastic gravitational waves (GWs)
spectrum, resulting from the emission of gravitons through bremsstrahlung, in
the decay of particles produced by Hawking radiation. Although particle decays
inevitably entail the emission of graviton due to bremsstrahlung, the
associated decay width is notably suppressed due to the Planck scale
suppression in the coupling of matter fields to gravitons. Consequently, the
relic abundance of such GWs constituted of these gravitons undergoes a
corresponding reduction. However, we demonstrate that super-heavy particles,
reaching masses as high as Planck scale, can emerge naturally in the Hawking
radiation of evaporating primordial black holes (PBHs) and can compensate for
this suppression. In addition, we also discuss the stochastic gravitational
waves constituted out of the gravitons directly radiated from such evaporating
PBHs. When the super-heavy particle decays promptly after its production, then
the corresponding GW spectrum remains subdominant to the one arising from
direct PBH evaporation. However, if this particle is long-lived and decays
after PBH evaporation, then the resulting GWs produced in these two processes
have two distinct spectra with their peaks at extremely high frequencies,
providing avenues for proposed ultra-high frequency gravitational wave
detectors. We also show that such gravitational waves contribute significantly
to substantial dark radiation, well within the anticipated sensitivity
thresholds of future experiments like CMB-S4 and EUCLID.

##  1 Introduction

The discovery of Gravitational Waves (GWs) by the LIGO and VIRGO collaboration
[1, 2] from different astrophysical sources as well as the detection of
stochastic gravitational wave background (SGWB), which might be potentially of
cosmological origin, by different Pulsar Timing Array (PTA) collaborations [3,
4, 5, 6], have ignited a great intrigue among the particle physics and
cosmology community. The discovery of GWs hold profound significance as GWs
can traverse the Universe practically unimpeded unlike the electromagnetic
radiation. Thus, GWs can carry undistorted information about their source and
hence they are the purest probes of physical phenomena even at exceptionally
high energy scales. Out of many possible phenomena that can source such GWs,
Hawking evaporation of primordial black holes (PBHs) [7, 8] that may have
formed in the early Universe is of great importance, as the ongoing
elusiveness of any clear, non-gravitational signals from any dark matter(DM)
candidates, has reignited interest in exploring PBHs, which can be potential
DM candidates themselves or play a crucial role in DM production (see e.g. [9,
10] and refrences therein). PBHs have also been explored to address
Baryogenesis [11, 12, 13, 14, 15] as well as in Cogenesis scenarios trying to
address the Baryon-DM coincidence problem (see e.g. [16, 17, 18]).

Black holes are theorized to emit quasi-thermal radiation, known as Hawking
radiation, which leads to their gradual evaporation over time [7, 8]. This
process allows black holes to dissipate energy at predictable rates into all
physical degrees of freedom that have a mass at or below the black hole’s
temperature [19]. The implications of this evaporation extend far beyond the
fate of the black holes themselves. As they radiate away their mass, PBHs
could potentially unleash a cascade of particles, of both visible sector and
the dark sector. Moreover, the gravitons emitted in the Hawking evaporation
can not equilibrate and thus survive until the present day, generating a
stochastic background of gravitational waves of high frequency distinct from
the PTA observations [3, 4, 5, 6]. This offers a unique signature that, if
detected, could provide unprecedented insights into the early Universe. The
prospect of detecting such gravitational waves is particularly fascinating as
it can potentially be sourced from light PBHs with masses below $\sim 5\times
10^{8}$ g. which would have evaporated before the era of big bang
nucleosynthesis (BBN), leaving no other observable traces [20, 21].

In addition, there is a possibility that extremely heavy particles with sub-
Planckian mass could be produced during Hawking evaporation, particularly
towards the end of PBH evaporation when the PBH temperature is extremely high.
These super-heavy particles could have significant implications and unique
detection prospects, which are the central focus of our discussion in this
article. Motivations for such super-heavy particles in various particle
physics models include their potential role as dark matter candidates,
involvement in baryogenesis processes, association with inflationary models,
addressing the hierarchy problem etc. These particles, if produced in Hawking
evaporation, can decay into SM particles or particles of the dark sector,
potentially generating gravitons through bremsstrahlung.

Gravitational waves originating from such quantum gravitational processes have
been studied in some recent works involving scattering of particles [22, 23],
particle annihilation to gravitons [24, 25] and bremsstrahlung in particle
decays [26, 27, 28, 29, 30, 31]. Given that particle decays are invariably
accompanied by graviton emission, gravitational waves stemming from graviton
bremsstrahlung hold particular significance among these quantum gravitational
phenomena. However production of gravitons through such decays is typically
suppressed due to the coupling of matter fields to gravitons by the Planck
mass (the effective interaction can be written as $\mathcal{L}_{\rm
eff}\supset\lambda h_{\mu\nu}T^{\mu\nu}$ with $\lambda\sim M^{-1}_{p}$ where
$h_{\mu\nu}$ is the graviton field and $T^{\mu\nu}$ is the energy momentum
tensor of other matter fields). But in the case of extremely heavy sub-
Planckian mass mother particles, this suppression can be compensated. This
compensation could lead to a sufficient energy transfer to gravitons,
resulting in observable ultra high frequency gravitational wave signatures.

Here, we study the GW stemming from the graviton production in Hawking
evaporation as well as GW production through bremsstrahlung process in the
decay of a long-lived super-heavy scalar produced in PBH evaporation. We
illustrate that GW production through bremsstrahlung can have a remarkable
enhancement over a wide frequency range, without any suppression, being
abundantly sourced from the super-heavy particle. We find that resulting GWs
produced in these two processes have two distinct spectra, thus offering a
potential avenue for probing Planck-scale physics, an endeavor that is
otherwise exceedingly challenging. These spectra of GWs from both sources peak
at extremely high frequencies. We demonstrate that the peak frequency for the
spectra of the GW from PBH evaporation is decided by the PBH initial mass,
whereas that for the GWs from graviton bremsstrahlung in particle decay is
decided by the mass of the super-heavy particle and the corresponding coupling
involved. We show that the resonant cavity detectors [32] have the potential
to detect these GW signals. Also such GWs can contribute to testable dark
radiation that can be constrained by the CMB observations [33] and can be
detected in future CMB experiments with enhanced sensitivities [34, 35, 36].

The manuscript is built up as follows: we begin by discussing the formation
and evaporation characteristics of PBHs in Section 2, followed by the
discussion of gravitational wave generation resulting from Hawking radiation
of gravitons in Section 3. Subsequently, in Section 4, we calculate the
abundance of long-lived super-heavy particles generated during PBH evaporation
and explore the production of ultra-high-frequency gravitational waves from
the graviton bremsstrahlung in their decay and finally conclude in Section 5.

##  2 Primordial Black Hole: Generation and Evaporation

Primordial black holes might have formed in the early universe through various
mechanisms, distinct from those that create stellar black holes. These
mechanisms often involve the gravitational collapse of matter overdensities,
which could have been seeded by inflation or resulted from topological
defects. The concept of PBH formation was first proposed by Zeldovich and
Novikov [37] in the late 1960s and was further developed by Hawking and Carr
in the early 1970s [38]. PBHs are hypothesized to have formed during the
radiation-dominated epoch immediately following the inflationary era [39, 40,
41]. In the radiation-dominated Universe, the initial mass of black hole
denoted as $M_{\rm in}$, at the moment of its formation can be linked to the
background energy density $\rho$ and the Hubble parameter $H$ at the
corresponding temperature $T_{\rm in}$ as [42, 43, 44]

$M_{\rm in}\equiv M_{\rm
BH}(T_{\text{in}})=\gamma\frac{4\pi}{3}\frac{\rho(T_{\rm in})}{H^{3}(T_{\rm
in})},$ (2.1)

where $\gamma\simeq 0.2$ is a numerical correction factor, $\rho(T_{\rm
in})=3M_{p}^{2}H^{2}(T_{\rm in})$ is the total energy density. Here,
$M_{p}=1/\sqrt{8\pi G}\simeq 2.4\times 10^{18}\ \rm{GeV}$ is the reduced
Planck mass. As a result, the temperature of the plasma, denoted as $T_{\rm
in}$ can be expressed in relation to $M_{\rm in}$ as

$T_{\rm in}=\left(\frac{1440\ \gamma^{2}}{g_{*}(T_{\rm
in})}\right)^{1/4}M_{p}\sqrt{\frac{M_{p}}{M_{\rm in}}}\simeq 4.36\times
10^{15}\ \text{GeV}\left(\frac{1\ \rm{g}}{M_{\rm in}}\right)^{1/2},$ (2.2)

where $g_{*}(T_{\rm in})$ denotes the effective degrees of freedom of the
energy density at the formation temperature and ${1\ {\rm g}\simeq 2.34\times
10^{5}M_{p}}$. Once formed, PBH can evaporate by emitting Hawking radiation,
transitioning into particles, if their masses are smaller than the
instantaneous Hawking temperature of PBH, $T_{\text{BH}}$. For Schwarzschild
black hole, $T_{\text{BH}}$ is given by [45, 46]

$T_{\text{BH}}=\frac{M_{p}^{2}}{M_{\text{BH}}}\simeq 10^{13}\ \rm
GeV\left(\frac{1\ \rm g}{M_{\rm BH}}\right).$ (2.3)

In case of Schwarzschild black holes, the energy spectrum of the i-th species
evaporated from the PBH can be written as [47, 48, 49]:

$\frac{d^{2}u_{i}(E,t)}{dEdt}=\frac{g_{i}}{2\pi^{2}}\frac{\sigma_{s_{i}}(M_{\text{BH}},\mu_{i},E_{i})}{e^{E_{i}/T_{\text{BH}}}-(-1)^{2s_{i}}}E_{i}^{3},$
(2.4)

where $g_{i}$ denotes the number of degrees of freedom of $i$-th species,
$E_{i}=\sqrt{\mu_{i}^{2}+p^{2}}$ represents the total energy of particles with
mass $\mu_{i}$ and spin $s_{i}$, and $\sigma_{s_{i}}$ signifies the absorption
cross section of a state with energy $E$ by a single PBH. Thus, the rate of
mass loss for the PBH can be obtained by summing the energy spectrum equation
(Eq. (2.4)), over all possible particle species, and integrating over the
phase space, which can be written as follows [50]:

$\displaystyle\frac{dM_{\text{BH}}}{dt}$
$\displaystyle=-\sum_{i}\int_{0}^{\infty}\frac{d^{2}u_{i}(E,t)}{dEdt}dE=-\varepsilon(M_{\text{BH}})\frac{M^{4}_{p}}{M_{\text{BH}}^{2}},$
(2.5)

where $\varepsilon(M_{\text{BH}})\equiv\sum_{i}g_{i}\varepsilon_{i}(z_{i})$ is
defined as an evaporation function with the function $\varepsilon_{i}(z_{i})$
given by [50]

$\varepsilon_{i}(z_{i})=\frac{27}{128\pi^{3}}\int_{z_{i}}^{\infty}\frac{\psi_{s_{i}}(x)(x^{2}-z_{i}^{2})}{e^{x}-(-1)^{2s_{i}}}xdx.$
(2.6)

Here $x=E_{i}/T_{\text{BH}}$ and $z_{i}=\mu_{i}/T_{\text{BH}}$ with
$T_{\text{BH}}$ given in Eq. (2.3) and we define
$\psi_{s_{i}}(E)=\sigma_{s_{i}}(E)/\left(\frac{27}{64\pi}\frac{M^{2}_{\rm
BH}}{M^{4}_{p}}\right)$.

In the geometric optics limit, where $\psi_{s_{i}}(E)=1$, and neglecting the
mass of the particles, the mass loss rate [51, 50] can be approximated as

$\frac{dM_{\text{BH}}}{dt}\simeq-\frac{27\pi}{4}\frac{\
g_{*}(T_{\text{BH}})}{480}\frac{M_{p}^{4}}{M_{\text{BH}}^{2}},$ (2.7)

where $g_{*}(T_{\rm BH})$ is the total number of relativistic degrees of
freedom emitted by PBH. By solving Eq. (2.7), the time evolution of PBH mass
is given by

$M_{\text{BH}}(t)=M_{\text{in}}\left(1-\frac{t-t_{i}}{\tau}\right)^{1/3},$
(2.8)

with the lifetime $\tau$ as

$\displaystyle\tau=\frac{4}{27}\frac{160\ M_{\text{in}}^{3}}{\pi\
g_{*}(T_{\text{BH}})M_{p}^{4}}\simeq 2.66\times 10^{-28}\
\rm{s}\left(\frac{100}{g_{*}(T_{\text{BH}})}\right)\left(\frac{M_{\rm in}}{1\
\rm{g}}\right)^{3}.$ (2.9)

The majority of the energy loss from PBH occurs during the final stages of the
evaporation process, when its mass has diminished significantly. The
evaporation temperature, denoted as $T_{\text{ev}}$, of the background plasma
immediately following the evaporation of the PBH, can be characterized at the
cosmic time $\tau$ by setting
$\rho_{\text{rad}}(T_{\text{ev}})=3M_{p}^{2}H^{2}(\tau)$. When the PBH
dominates the energy density of the Universe prior to its evaporation, the
Universe enters a matter-dominated phase. In such a scenario,
$H^{2}=\frac{4}{9\tau^{2}}$, resulting in the evaporation temperature given by
[43]:

$\begin{split}T_{\text{ev}}|_{\text{MD}}&\equiv\left(\frac{40M_{p}^{2}}{\pi^{2}g_{*}(T_{\text{ev}})\tau^{2}}\right)^{1/4}\simeq{0.5}\left(\frac{g_{*}^{2}(T_{\rm
BH})}{g_{*}(T_{\rm ev})}\right)^{1/4}\left(\frac{M_{p}^{5}}{M_{\rm
in}^{3}}\right)^{1/2}\simeq{3.55\times 10^{10}}\ \text{GeV}\left(\frac{1\
\text{g}}{M_{\text{in}}}\right)^{3/2},\end{split}$ (2.10)

where $g_{*}(T_{\text{BH}})\simeq g_{*}(T_{\text{ev}})\simeq 100$ is used. In
scenarios where PBH are subdominant before evaporation,
$H^{2}(T_{\text{ev}})=\frac{1}{4\tau^{2}}$ can be utilized leading to an
evaporation temperature which is slightly smaller than that in the matter-
dominated phase [43]

$T_{\text{ev}}|_{\text{RD}}\simeq\frac{\sqrt{3}}{2}T_{\text{ev}}|_{\text{MD}}.$
(2.11)

A lower limit on the mass of PBH is imposed by the maximum Hubble size of the
Universe following inflation, as the tensor-to-scalar ratio is constrained by
observations of cosmic microwave background anisotropy. This constraint is
expressed as $H<2.5\times 10^{-5}M_{p}$ [52]. Moreover, when PBH dominates the
Universe before its evaporation, it is crucial that the evaporation
temperature is sufficiently high to avoid disrupting the predictions of big
bang nucleosynthesis(BBN), necessitating $T_{\text{ev}}>T_{\text{BBN}}\simeq
4\ \text{MeV}$. This requirement sets an upper bound on the PBH mass. By
considering both these constraints, we determine the available range of PBH
mass as:

$0.4\ \text{g}\lesssim M_{\text{in}}\lesssim 9.7\times 10^{8}\ \text{g}.$
(2.12)

##  3 Gravitational Wave production from PBH evaporation

As explained in the previous section, PBH can produce a wide array of
particles with masses $m\lesssim T_{\rm BH}$, including gravitons, through
Hawking evaporation. For the sake of simplicity, considering a Schwarzschild
black hole, the greybody factor, which is a measure of the efficiency of
particle emission, can be approximated as
$\sigma_{s_{i}}=\left(\frac{27}{64\pi}\frac{M^{2}_{\rm BH}}{M^{4}_{p}}\right)$
in the geometric optics limit. Consequently, taking into account the two
graviton polarizations $g_{i}=2$ and the spin of graviton $s=2$, Eq.(2.4)
simplifies to the following expression for graviton production from the direct
evaporation of PBH;

$\frac{d^{2}u_{\rm GW}}{dtdE}\simeq\frac{27}{64\pi^{3}}\frac{M_{\rm
BH}^{2}}{M_{p}^{4}}\frac{E^{3}}{e^{E/T_{\rm BH}}-1}.$ (3.1)

Therefore, calculating the rate of graviton emission for a collective
population of evaporating PBHs, we can write

$\frac{d\rho_{\rm GW}}{dtdE}\simeq n_{\rm BH}(t)\frac{d^{2}u_{\rm GW}}{dtdE}.$
(3.2)

where $n_{\rm BH}(t)$ is the PBH number density. Here, it should be worth
keeping in mind that, the presence of PBH alongside the standard radiation
bath contributes to the overall density and alters the Hubble expansion rate
$H$ of the early universe expressed as:

$H^{2}=\frac{1}{3M_{p}^{2}}(\rho_{r}+\rho_{\text{BH}}).$ (3.3)

Using the comoving energy density of radiation and PBH defined as
$\tilde{\rho}_{r}\equiv a^{4}\times\rho_{r}$ and $\tilde{\rho}_{\rm BH}\equiv
a^{3}\times\rho_{\rm BH}$, where $a$ represents the scale factor, the
Boltzmann equations governing the evolution of the Universe are given by [43,
53, 10, 54, 55]:

$\displaystyle\frac{dM_{\text{BH}}}{d\ln(a)}=-\frac{\varepsilon(M_{\text{BH}})}{H}\frac{M_{p}^{4}}{M_{\text{BH}}^{2}},$
(3.4)
$\displaystyle\frac{d\tilde{\rho}_{\text{BH}}}{d\ln(a)}=\frac{\tilde{\rho}_{\text{BH}}}{M_{\text{BH}}}\frac{dM_{\text{BH}}}{d\ln(a)},$
$\displaystyle\frac{d\tilde{\rho}_{r}}{d\ln(a)}=-\frac{\varepsilon_{\text{SM}}(M_{\text{BH}})}{\varepsilon(M_{\text{BH}})}\frac{a\
\tilde{\rho}_{\text{BH}}}{M_{\text{BH}}}\frac{dM_{\text{BH}}}{d\ln(a)}\,.$

Here $\varepsilon_{\text{SM}}\equiv g_{\rm SM}\sum_{i}\varepsilon_{i}(z_{i})$
represents the contribution solely from the SM to the evaporation excluding
the contribution from gravitons. In this paper, we assume a monochromatic mass
spectrum for the PBH and thus the comoving number density of PBH remain
conserved, i.e. $\tilde{n}_{\rm BH}$=constant. Hence, the PBH number density
at any time can be related to the initial number density $n_{{\rm BH},i}$ as:

$n_{\rm BH}(t)=n_{\text{BH},i}\left(\frac{a_{i}}{a}\right)^{3}.$ (3.5)

Therefore, the energy density of the PBHs can be expressed as

$\rho_{\rm BH}(t)=n_{\rm BH}(t)M_{\rm BH}(t)=\rho_{\rm BH}(t_{\rm
in})\left(\frac{a_{i}}{a}\right)^{3}\left(1-\frac{t-t_{i}}{\tau}\right)^{1/3}.$
(3.6)

Using sudden decay approximation with $\rho_{\rm BH}(t_{\rm
ev})=\rho_{r}(t_{\rm ev})$, we find the relation,

$\begin{split}\rho_{\rm BH}(t_{\rm ev})=\rho_{\rm
BH}(t_{i})\left(\frac{a_{i}}{a_{\rm ev}}\right)^{3}=\rho_{r}(t_{\rm
ev})=\frac{\pi^{2}}{30}g_{*}(T_{\rm ev})T_{\rm ev}^{4}\end{split}$ (3.7)

Therefore,

$\begin{split}n_{\text{BH},i}a_{i}^{3}=\left(\frac{\rho_{r}(t_{\rm ev})a_{\rm
ev}^{3}}{M_{\rm in}}\right)=\left(\frac{\frac{\pi^{2}}{30}g_{*}(T_{\rm
ev})T_{\rm ev}^{4}a_{\rm ev}^{3}}{M_{\rm in}}\right).\end{split}$ (3.8)

Following their production, the energy density and frequency of GW undergo
redshifting, and their subsequent values can be related to those at the time
of complete PBH evaporation by the equations:

$\rho_{\rm GW}=\rho_{\rm GW,ev}\left(\frac{a_{\rm
ev}}{a}\right)^{4},\quad\omega=\omega_{\rm ev}\left(\frac{a_{\rm
ev}}{a}\right).$ (3.9)

Utilizing Eqs. (3.5) and (3.9) for the scaling as well as Eq.(2.8) for the PBH
mass evolution, the energy density for gravitational waves at the time of PBH
evaporation can be estimated as:

$\frac{d\rho_{\rm GW,ev}}{d\ln\omega_{\rm
ev}}=\frac{27}{64\pi^{3}}\frac{M_{\rm in}^{2}}{M_{p}^{4}}n_{\rm
BH}(t_{i})\omega_{\rm ev}^{4}\int_{t_{i}}^{t_{\rm
ev}=t_{i}+\tau}dt\left(1-\frac{t-t_{i}}{\tau}\right)^{2/3}\frac{(a_{i}/a)^{3}}{e^{\omega_{\rm
ev}a_{\rm ev}/aT_{\rm BH}}-1},$ (3.10)

The relationship between the energy density of gravitational waves today and
that at the time of evaporation is expressed as:

$\frac{d\rho_{\rm GW,0}}{d\ln\omega_{0}}=\frac{d\rho_{\rm
GW,ev}}{d\ln\omega_{\rm ev}}\left(\frac{a_{\rm ev}}{a_{0}}\right)^{4},$ (3.11)

with

$\frac{a_{\rm ev}}{a_{0}}\simeq\frac{T_{0}}{T_{\rm
ev}}\left(\frac{g_{*,s}(T_{0})}{g_{*,s}(T_{\rm
ev})}\right)^{1/3}\simeq{6.2\times
10^{-32}}\left(\frac{g_{*,s}(T_{0})}{g_{*,s}(T_{\rm
ev})}\right)^{1/3}\left(\frac{M_{\rm in}}{M_{p}}\right)^{3/2}\simeq 2.3\times
10^{-24}\left(\frac{M_{\rm in}}{1\ \rm{g}}\right)^{3/2},$ (3.12)

where we take $T_{0}\simeq 2.3\times 10^{-13}\ \rm{GeV}$,
$g_{*,s}(T_{0})=3.91$ and $g_{*}(T_{\text{ev}})\simeq
g_{*,s}(T_{\text{ev}})\simeq 100$. Denoting the redshifted frequency values
today as $\omega_{0}=\omega_{\rm ev}a_{\rm ev}$, with $a_{0}=1$, the energy
density of gravitational waves at present can be expressed as

$\displaystyle\frac{d\rho_{\rm GW,0}}{d\ln\omega_{0}}$ $\displaystyle=$
$\displaystyle\frac{27}{64\pi^{3}}\frac{M_{\rm in}^{2}}{M_{p}^{4}}n_{\rm
BH}(t_{i})\omega_{0}^{4}\int_{t_{i}}^{t_{\rm
ev}}dt\left(1-\frac{t-t_{i}}{\tau}\right)^{2/3}\frac{(a_{i}/a)^{3}}{e^{\omega_{0}/aT_{\rm
BH}}-1}$ (3.13) $\displaystyle\simeq$ $\displaystyle{2.8}\times
10^{-4}g_{*}^{2}(T_{\rm BH})\frac{M_{p}^{6}}{M_{\rm
in}^{5}}\omega_{0}^{4}\int_{t_{i}}^{t_{\rm
ev}}dt\left(1-\frac{t-t_{i}}{\tau}\right)^{2/3}\frac{(a_{\rm
ev}/a)^{3}}{e^{\omega_{0}/aT_{\rm BH}}-1},$

where in the last equality, we utilize the fact that the number density of PBH
remains constant over time intervals evident from Eq. (3.8). Using Eq.(3.12)
in Eq.(3.13), we obtain:

$\displaystyle\frac{d\rho_{\rm GW,0}}{d\ln\omega_{0}}$ $\displaystyle\simeq$
$\displaystyle{6.2\times 10^{-77}}\ {\rm GeV}\left(\frac{M_{p}}{M_{\rm
in}}\right)^{1/2}\omega_{0}^{4}\ I(\omega_{0}),$ (3.14)

where the frequency dependent integral is defined as

$I(\omega_{0})=\int_{t_{i}}^{t_{\rm
ev}}dt\left(1-\frac{t-t_{i}}{\tau}\right)^{2/3}\frac{a^{-3}}{e^{\omega_{0}/aT_{\rm
BH}}-1}.$ (3.15)

Therefore, the final relic abundance for gravitational wave can be estimated
as

$\displaystyle h^{2}\Omega_{\rm GW}$ $\displaystyle=$
$\displaystyle\frac{1}{\rho_{\rm cr,0}h^{-2}}\frac{d\rho_{\rm
GW,0}}{d\ln\omega_{0}}\simeq 1.5\times 10^{-33}\ \text{GeV}^{-4}\left(\frac{1\
\rm g}{M_{\rm in}}\right)^{1/2}\omega_{0}^{4}\ I(\omega_{0}),$ (3.16)

where $\rho_{\rm cr,0}\simeq 8h^{2}\times 10^{-47}\ \rm{GeV}^{4}$ with
$h\simeq 0.7$.

We can analytically find the peak frequency at present by extremizing
$h^{2}\Omega_{\rm GW}$ with respect to $\omega_{0}$. Using the blackbody
approximation for the Schwarzschild blackholes, one obtains $\omega_{\rm
peak}\simeq 2.8a_{\rm ev}T_{\rm BH}$ or more explicitly

$f_{\rm peak}\simeq 1.6\times 10^{13}\ {\rm Hz}\left(\frac{M_{\rm in}}{1\ \rm
g}\right)^{1/2}.$ (3.17)

Figure 1: Gravitational wave spectrum from direct evaporation of PBH for
different initial mass of PBH, $M_{\rm in}=1\ \rm{g}$, $10^{4}\ \rm{g}$, and
$10^{8}\ \rm{g}$. The solid lines denote the full numerical integration, while
the dashed lines represent the analytical approximation given in Eq. (3.16).

To evaluate the integration $I(\omega_{0})$, we turn to the matter-dominated
or PBH-dominated epoch, as gravitons are generated during PBH evaporation and
after PBH evaporation, no further gravitons are produced. During the matter-
dominated epoch, the scale factor is related to time as $a=At^{2/3}$, where
$A$ is the coefficient determined by $a=a_{\rm
ev}\left(\frac{t}{\tau}\right)^{2/3}=At^{2/3}$. Consequently using Eqs. (2.9)
and (3.12), $A$ is given by:

$\displaystyle A=a_{\rm ev}\left(\frac{1}{\tau}\right)^{2/3}\simeq{1.17\times
10^{-31}}\left(\frac{M_{p}^{7/6}}{M_{\rm in}^{1/2}}\right),$ (3.18)

Thus the integral Eq. (3.15) can be written as

$\displaystyle I(\omega_{0})$ $\displaystyle=$ $\displaystyle
A^{-3}\int_{t_{i}}^{t_{\rm
ev}}dt\left(1-\frac{t-t_{i}}{\tau}\right)^{2/3}\frac{t^{-2}}{\exp\left[\alpha
t^{-2/3}\left(1-\frac{t-t_{i}}{\tau}\right)^{1/3}\right]-1},$ (3.19)

where $\alpha\equiv\frac{\omega_{0}M_{\rm in}}{AM_{p}^{2}}$. Typically, the
frequency-dependent integral necessitates numerical evaluation for each
$\omega_{0}$, which can then be converted to linear frequency as
$f=\omega/2\pi$. For this, we evaluate the formation time $t_{i}$ by using
Eq.(2.1), and the total energy density at the formation time of the PBH;
$\rho(T_{\rm in})=3M_{p}^{2}H^{2}(T_{\rm in})$, as:

$t_{i}=\frac{M_{\rm in}}{8\pi\gamma M^{2}_{p}},$ (3.20)

and the evaporation time is approximated as $t_{\rm ev}=t_{i}+\tau\simeq\tau$.

The gravitational wave spectra from direct evaporation of PBH is shown in Fig.
1 as a function of frequency $f$ for the different initial mass of PBH,
$M_{\rm in}=1,10^{4},10^{8}\ \rm g$. The solid lines represent numerical
result taking into account for the full greybody factors, while the dashed
lines depict the analytic results for which we assume the geometric optics
limit for the evaporation of PBH. Here it is worth noticing that as the PBH
mass increases, the peak frequency of the GW spectrum increases, consequently
shifting the spectrum towards higher frequencies. This is due to the fact that
larger mass of PBH corresponds to formation of PBH at a later time,
consequently leading to a relatively shorter period of cosmological redshift
thereby shifting the spectrum to higher frequencies, however, the peak
amplitude remains unchanged.

##  4 Gravitational Wave production from massive scalar particle

Massive scalar particles of the sub-Planckian mass can be produced from PBHs
at the late stage of its lifetime as a consequence of the increase in Hawking
temperature as the PBH mass decreases. Once this massive scalar particle is
produced from PBH, it can decay into the SM fermions or dark sector particles,
with its lifetime being decided by its mass and couplings. The Boltzmann
equation to track the evolution of such a scalar $\chi$ can be written as

$\displaystyle\frac{d\tilde{n}_{\chi}}{d\ln(a)}=\frac{\tilde{\rho}_{\text{BH}}}{M_{\text{BH}}}\frac{\Gamma_{\text{BH}\rightarrow\chi}}{H}-\frac{\Gamma_{\chi\rightarrow
f\bar{f}}}{H}\tilde{n}_{\chi},$ (4.1)

which has to be solved along with Eq.(3.4) simultaneously. Here
$\tilde{n}_{\chi}\equiv n_{\chi}a^{3}$ is the comoving number density of
$\chi$ and $\Gamma_{\chi}$ is the decay rate of scalar particle $\chi$
decaying into a pair of fermions. This decay width is given by

$\Gamma_{\chi\rightarrow f\bar{f}}=\frac{y_{f}^{2}}{8\pi}m_{\chi},$ (4.2)

with Yukawa coupling $y_{f}$ and mass of scalar $\chi$ denoted as $m_{\chi}$.
Thus the lifetime of $\chi$ can be estimated to be

$\tau_{\chi}\simeq 7\times 10^{-26}\ {\rm
s}\left(\frac{10^{-7}}{y_{f}}\right)^{2}\left(\frac{10^{-2}M_{p}}{m_{\chi}}\right).$
(4.3)

Clearly, we can see that, by an appropriate choice of the Yukawa coupling
$y_{f}$, the scalar $\chi$ can be long lived enough to decay after PBH
evaporation if $\tau_{\chi}>\tau_{\rm BH}$, while still being safe from
cosmological and phenomenological constraints. This requirement imposes an
upper limit on the value of $y_{f}$ for our analysis.

Furthermore, $\Gamma_{\text{BH}\rightarrow\chi}$ is the momentum integrated
emission rate of scalar particle from PBH evaporation given by

$\displaystyle\Gamma_{\text{BH}\rightarrow\chi}(t)=\int
dp\frac{d^{2}N_{\chi}}{dp\
dt}=\frac{27g_{\chi}}{128\pi^{3}}\frac{M_{p}^{2}}{M_{\text{BH}}(t)}\int_{z}^{\infty}\frac{\psi_{\chi}(x)(x^{2}-z^{2})}{e^{x}-1}dx,$
(4.4)

where $x=E/T_{\text{BH}}$, $z=m_{\chi}/T_{\rm BH}$. Similar to the previous
section, again resorting to the geometric optics limit for the greybody
factors, we put $\psi_{\chi}(x)=1$, which thus allows the integration to be
analytically solved. Utilizing Eq. (2.8) for PBH mass evolution, the emission
rate is given by [50, 49]:

$\Gamma_{\text{BH}\rightarrow\chi}(t)=\frac{27g_{\chi}}{64\pi^{3}}\frac{M_{p}^{2}}{M_{\text{in}}}\left(1-\frac{t-t_{i}}{\tau}\right)^{-1/3}\mathcal{G}(z).$
(4.5)

Here $\mathcal{G}(z)=[z{\rm Li}_{2}(e^{-z})+{\rm Li}_{3}(e^{-z})]$ and
$\rm{Li}_{n}$ is the polylog function of order $n$. It is worth noting here
that for temperatures significantly higher than $m_{\chi}$,(i.e. $z\ll 1$),
$\mathcal{G}(z)$ approaches the value $\zeta(3)=1.20206$ where $\zeta(n)$ is
the Riemann zeta function. However, for our numerical calculation, we utilise
the "ULYSSES" package which incorporates the modified graybody factors [56].

In Fig. 2, the left panel shows the numerical results for the evolution of the
comoving energy density of radiation (cyan), PBH (black) and the massive
scalar particle (red) as a function of the scale factor $a$. One can see that
PBHs dominate the energy density of the universe before they evaporate
completely. The right panel shows the number density of massive scalar
particle in terms of the scale factor $a$. An enhancement in the number
density for $\chi$ can be noticed at the scale factor of approximately $a\sim
10^{-23}$, which corresponds to the fact that most of the massive scalar
particle can be produced towards the end of PBH lifetime when $T_{\rm BH}$
increases beyond $m_{\chi}$ due to the PBH mass loss.

Figure 2: Evolution of the comoving energy densities for radiation, PBH and
scalar particle (left) and the number density of scalar particle (right) as a
function of scale factor. Here, we take the initial mass of PBH as $M_{\rm
in}=1\ \rm{g}$.

If this produced scalar particle $\chi$ has a short lifetime such that it
decays promptly into SM fermions or dark sector particles, after its
production even before the end of complete evaporation of PBH, then its number
density diminishes extremely quickly such that the production of the scalar
$\chi$ loses its phenomenological significance as there will not be any
interesting consequences. Conversely, we focus on the scenario where this
scalar has a long lifespan and decays well after the PBH evaporation. In such
a case, we can disregard the second term on the right-hand side of Eq. (4.1)
for the production of scalars from PBH. Therefore, the number density of
massive scalar particle, $n_{\chi}$, right after PBH evaporation can be
written as

$(a^{3}n_{\chi})|_{\rm ev}=n_{\rm BH}(t_{i})a_{i}^{3}\int_{t_{i}}^{t_{\rm
ev}}\Gamma_{BH\rightarrow\chi}(t)\,dt,$ (4.6)

where we used the constancy of $n_{\rm BH}(t_{i})a_{i}^{3}$. With Eq.(4.5), it
becomes

$\displaystyle n_{\chi}(a_{\rm ev})$ $\displaystyle=$
$\displaystyle\frac{27g_{\chi}}{64\pi^{3}}\frac{M_{p}^{2}}{M_{\rm{in}}}n_{\rm
BH}(t_{i})\left(\frac{a_{i}}{a_{\rm ev}}\right)^{3}\int_{t_{i}}^{t_{\rm
ev}}dt\left(1-\frac{t-t_{i}}{\tau}\right)^{-1/3}~{}\mathcal{G}(z),$ (4.7)
$\displaystyle\simeq$
$\displaystyle\frac{15~{}g_{\chi}~{}\zeta(3)}{\pi^{4}g_{*}(T_{\rm
BH})}\frac{M^{2}_{p}}{m^{2}_{\chi}}n_{\rm BH}(t_{i})\left(\frac{a_{i}}{a_{\rm
ev}}\right)^{3}\simeq 1.54\times 10^{21}\left(\frac{1g}{M_{\rm
in}}\right)^{7}\left({\frac{10^{-2}M_{p}}{m_{\chi}}}\right)^{2},$

where we have used Eq.(2.8) for PBH mass evolution and an approximation
$\mathcal{G}(z)\simeq\zeta(3)$ for $t>t_{1}$ in the second line. Here $t_{1}$
is characterized as the point when the production of $\chi$ through Hawking
evaporation becomes viable i.e. when $T_{\rm BH}=m_{\chi}$ or $M_{\rm
BH}(t_{1})={M_{p}^{2}}/{m_{\chi}}$ and consequently we have:

$\begin{split}1-\frac{t_{1}}{\tau}=\left(\frac{M_{p}^{2}}{M_{\rm
in}m_{\chi}}\right)^{3}.\end{split}$ (4.8)

Once these super heavy scalars are produced from the PBH evaporation, their
decay into some fermions, also produces a graviton through bremsstrahlung
process as depicted in Fig. 3. Here it is worth mentioning that, if the $\chi$
is short lived and decays immediately after its production during PBH
evaporation process, then the associated energy density of the gravitons
produced in the bremsstrahlung will always be subdominant to the energy
density of gravitons produced directly in PBH evaporation. Thus in such a
scenario, it is not possible to have a distinguishable GW signature for the
scalar $\chi$. Hence, we focus on the case where $\chi$ decays much later
after PBH evaporation.

Figure 3: Feynman diagrams representing the decay of a scalar particle
($\chi$) into fermion ($f$) along with a graviton ($h_{\mu\nu}$)
bremsstrahlung.

The evolution of the energy density of gravitons $\rho_{\rm GW}$ produced in
the bremsstrahlung of $\chi$, can be tracked by the Boltzmann equation

$\frac{d}{dt}\left(\frac{d\rho_{\rm GW}}{dE_{\rm
GW}}\right)+4H\frac{d\rho_{\rm GW}}{dE_{\rm GW}}=n_{\chi}(a_{\rm
ev})\left(\frac{a_{\rm ev}}{a}\right)^{3}\frac{d\Gamma_{\chi\rightarrow\rm
GW}}{dE_{\rm GW}}E_{\rm GW},$ (4.9)

or in terms of the scale factor $a$ as

$\displaystyle\frac{d}{da}\left(a^{4}\frac{d\rho_{\rm GW}}{d\ln E_{\rm
GW}}\right)=\frac{n_{\chi}(a_{\rm ev})a_{\rm
ev}^{3}}{H}\frac{d\Gamma_{\chi\rightarrow\rm GW}}{dE_{\rm GW}}E_{\rm GW}^{2},$
(4.10)

where $d\Gamma_{\chi\rightarrow\rm GW}/dE_{\rm GW}$ is the differential decay
rate of scalar particle. When the decay happens after $\chi$ becomes non-
relativistic, it is given by [30]

$\frac{d\Gamma_{\chi\rightarrow\rm GW}}{dE_{\rm
GW}}=\frac{y_{f}^{2}m_{\chi}^{3}}{64\pi^{3}M_{p}^{2}E_{\rm GW}}F(E_{\rm
GW}/m_{\chi}),$ (4.11)

with form factor: $F(x)=(1-2x)(1-2x+2x^{2})$. We consider that the generation
of gravitational waves can occur from the period spanning the complete
evaporation of PBH ($a_{\rm ev}$) to the decay of the scalar particle
($a_{\chi}$).

$\displaystyle a_{\chi}^{4}\frac{d\rho_{\rm GW}(a_{\chi})}{d\ln E_{\rm GW}}$
$\displaystyle=$
$\displaystyle\frac{y_{f}^{2}m_{\chi}^{3}}{64\pi^{3}M_{p}^{2}}E_{\rm
GW}F(E_{\rm GW}/m_{\chi})n_{\chi}(a_{\rm ev})a_{\rm ev}^{3}\int_{a_{\rm
ev}}^{a_{\chi}}da\frac{1}{H}$ (4.12) $\displaystyle=$
$\displaystyle\frac{y_{f}^{2}m_{\chi}^{3}}{64\pi^{3}M_{p}^{2}}E_{\rm
GW}F(E_{\rm GW}/m_{\chi})n_{\chi}(a_{\rm ev})a_{\rm
ev}^{3}~{}\frac{2}{3A^{2}}(a_{\chi}^{3}-a_{\rm ev}^{3}),$

where we use $H=\dot{a}/a$ and the relation $a=At^{1/2}$ with $A=a_{\rm
ev}/\tau^{1/2}$ as the Universe is radiation dominated when the particle
$\chi$ decays. Therefore, $A$ can be estimated to be

$A=\frac{a_{\rm ev}}{\tau^{1/2}}\simeq 7.4\times 10^{-32}M_{p}^{1/2}.$ (4.13)

Now using the result of Eq.(4.7), one can find the energy density of
gravitational wave at $a_{\chi}$ as:

$\displaystyle\frac{d\rho_{\rm GW}(a_{\chi})}{d\ln E_{\rm GW}}$
$\displaystyle\simeq$ $\displaystyle 2.3\times
10^{58}~{}\frac{y_{f}^{2}m_{\chi}M_{p}^{9}}{M_{\rm in}^{7}}E_{\rm GW}F(E_{\rm
GW}/m_{\chi})\frac{a_{\rm ev}^{3}}{a_{\chi}}\left[1-\left(\frac{a_{\rm
ev}}{a_{\chi}}\right)^{3}\right]$ (4.14)

Thus the relic abundance of gravitational wave can be estimated as

$h^{2}\Omega_{\rm GW}=\frac{1}{\rho_{\rm cr,0}h^{-2}}\frac{d\rho_{\rm
GW}(a_{\chi})}{d\ln E_{\rm GW}}\left(\frac{a_{\chi}}{a_{0}}\right)^{4},$
(4.15)

where $E_{\rm GW}=2\pi f(a_{0}/a_{\chi})$. Therefore, $h^{2}\Omega_{\rm GW}$
can be written by

$\displaystyle h^{2}\Omega_{\rm GW}$ $\displaystyle\simeq$
$\displaystyle\frac{2.3\times 10^{58}}{\rho_{\rm
cr,0}h^{-2}}\frac{y_{f}^{2}m_{\chi}M_{p}^{9}}{M_{\rm in}^{7}}E_{\rm
GW}F(E_{\rm GW}/m_{\chi})\frac{a_{\rm
ev}^{3}}{a_{\chi}}\left[1-\left(\frac{a_{\rm
ev}}{a_{\chi}}\right)^{3}\right]\left(\frac{a_{\chi}}{a_{0}}\right)^{4}$
(4.16) $\displaystyle\simeq$ $\displaystyle 1.8\times 10^{105}\
{\rm{GeV}^{-4}}\frac{y_{f}^{2}m_{\chi}M_{p}^{9}}{M_{\rm in}^{7}}F(E_{\rm
GW}/m_{\chi})f~{}a_{\rm ev}^{3}\left(\frac{a_{\chi}}{a_{0}}\right)^{2}.$

The temperature at which the massive scalar particle decays into fermions and
thus produces the GW via bremsstrahlung, denoted as $T_{\chi}$, can be
determined from the relation $H=1/2\tau_{\chi}$, during the radiation
dominated epoch as:

$T_{\chi}={\left(\frac{3\sqrt{10}y_{f}^{2}M_{p}m_{\chi}}{16\pi^{2}g_{*}^{1/2}(T_{\chi})}\right)^{1/2}\simeq\frac{3y_{f}}{4\pi}\left(\frac{M_{p}m_{\chi}}{g_{*}^{1/2}(T_{\chi})}\right)^{1/2}}.$
(4.17)

From the entropy conservation, $a_{\chi}/a_{0}$ is given by

$\displaystyle\frac{a_{\chi}}{a_{0}}$ $\displaystyle\simeq$
$\displaystyle\frac{T_{0}}{T_{\chi}}\left(\frac{g_{*,s}(T_{0})}{g_{*,s}(T_{\chi})}\right)^{1/3}\simeq\frac{{10^{-12}}}{y_{f}}\left(\frac{\text{GeV}^{2}}{M_{p}m_{\chi}}\right)^{1/2},$
(4.18)

where we use $g_{*,s}(T_{\chi})\simeq g_{*}(T_{\chi})\simeq 100$ in the last
equality. Finally, using Eqs. (3.12) and (4.18), we can find the relic
abundance for GW as

$h^{2}\Omega_{\rm GW}\simeq 9\times 10^{-17}\left(\frac{f}{10^{17}\ \rm
Hz}\right)\left(\frac{1\ \rm{g}}{M_{\rm in}}\right)^{5/2}F(E_{\rm
GW}/m_{\chi}).$ (4.19)

Figure 4: Gravitational wave spectrums for graviton bremsstrahlung from decay
of massive scalar particle for two different masses $m_{\chi}=10^{-2}M_{p}$
(left) and $m_{\chi}=0.1M_{p}$ (right) are shown for different values of the
coupling $y_{f}$. The numerical result (solid lines) and the analytical
estimation in Eq.(4.19) (dashed lines) are shown together.

Since we consider the GWs generated by the scalar particle decay, the GWs at
production could carry at most energy $E_{\rm GW}\leq m_{\chi}/2$.
Cosequently, the peak frequency is calculated as

$\displaystyle f_{\rm
peak}=\frac{m_{\chi}}{4\pi}\left(\frac{a_{\chi}}{a_{0}}\right)$
$\displaystyle\simeq$ $\displaystyle 1.2\times 10^{17}\
\text{Hz}\left(\frac{10^{-7}}{y_{f}}\right)\left(\frac{m_{\chi}}{10^{-2}M_{p}}\right)^{1/2},$
(4.20)

where we use the result of Eq.(4.18).

In Fig. 4, we illustrate the amplitude of the gravitational wave (GW)
spectrum, $\Omega_{\rm GW}h^{2}$, plotted against the frequency $f$ of the GWs
for various $y_{f}$ values, as indicated in the plot legend. It is evident
that the GW spectrum remains unaffected by different Yukawa couplings $y_{f}$
and scalar particle mass $m_{\chi}$. However, the peak position can be shift
towards the higher frequency region by tuning the coupling $y_{f}$ to smaller
values which in turn makes the scalar $\chi$ to decay at relatively late
times. This is due to the fact that if $\chi$ decays late then it leads to a
relatively shorter period of cosmological redshift. Clearly, we can see that
our numerical result and the analytical estimation are in good agreement with
each other.

The dimensionless characteristic strain, $h_{c}$, can be defined as

$h_{c}=f^{-1}\sqrt{\frac{3H_{0}^{2}}{4\pi^{2}}\Omega_{\rm GW}}\simeq
8.93\times 10^{-19}\sqrt{\Omega_{\rm GW}h^{2}}\left(\frac{\rm Hz}{f}\right),$
(4.21)

where the Hubble rate at present $H_{0}=100h\ \text{km}\ \text{s}^{-1}\
\text{Mpc}^{-1}\simeq 3.24\times 10^{-18}h\ \text{s}^{-1}$ . Our result
Eq.(4.16) is written in terms of the dimensionless characteristic strain as

$h_{c}\simeq 8.5\times 10^{-44}\left(\frac{10^{17}\ \rm
Hz}{f}\right)^{1/2}\left(\frac{1\ \rm{g}}{M_{\rm in}}\right)^{5/4}.$ (4.22)

Figure 5: GW characteristic strain $h_{c}$ for the initial mass of PBH $M_{\rm
in}=1\ \rm{g}$ and several different values of $y_{f}=10^{-7},10^{-9},$ and
$10^{-11}$ with the scalar particle of mass $m_{\chi}=10^{-2}M_{p}$ (left) and
$m_{\chi}=M_{p}$ (right). Sensitivities of some current and future proposed GW
detectors [57, 58, 59, 60, 61] are reflected by shaded areas. The brown line
illustrates the parameter space which can be probed by resonant cavities [32].
The dark yellow lines are CMB bound on $\Delta N_{\rm eff}$ from Planck [33]
with solid line, and future sensitivity from CMB-S4 [34] with dashed line,
Euclid [35] with dot-dashed line, and proposed CMB-CVL experiment [36] with
dotted line, respectively.

In Fig. 5, we present the dimensionless characteristic strain $h_{c}$ plotted
as a function of the GW frequency $f$. Here, we compare the contributions from
both the direct evaporation of PBH with a red solid line and the graviton
bremsstrahlung from the massive scalar particle with different values of the
coupling $y_{f}$. Shaded regions represent the limits from various proposed GW
detectors, including LISA [57], the Big Bang Observer (BBO) [57], CE [58],
DECIGO [59], and the advanced LIGO [60]. Moreover, we also showcase the
existing and aniticipated promising experimental techniques that has the
potential to detect high-frequency gravitational waves such as optically
levitated sensors, enhanced magnetic conversions (EMC), inverse Gertsenshtein
effect etc., including JURA, ALPS II, OSQAR, IAXO and CAST, which we adopt
from [61]. The brown line illustrates the parameter space that can be probed
by resonant cavities [32]. It is evident that the resonant cavities have the
potential to test both the GW spectra when the scale of the scalar particle
mass $m_{\chi}$ is around $1M_{p}$. Thus, it would imply a direct probe of the
Planck-scale physics, which is otherwise challenging to explore.

Energy density of GW background decreases with expansion of the universe
following the behavious of relativistic degrees of freedom $\rho_{\rm
GW}\propto a^{-4}$. Thus it acts as an additional radiation field in the
Universe, contributing to the effective number of neutrino species, $N_{\rm
eff}$, which is defined as

$\rho_{\rm
rad}=\rho_{\gamma}\left(1+\frac{7}{8}\left(\frac{4}{11}\right)^{4/3}N_{\rm
eff}\right),$ (4.23)

where $N_{\rm eff}=N_{\rm eff}^{\rm SM}+\Delta N_{\rm eff}$. Here $N_{\rm
eff}^{\rm SM}=3.046$ [62] represents the contribution from SM neutrinos and
$\Delta N_{\rm eff}$ parameterizes the additional neutrino-like relativistic
component. Therefore, taking $\rho_{\gamma}=\frac{\pi^{2}}{15}T_{0}^{4}$ and
$\rho_{\rm GW}\simeq\Omega_{\rm GW}^{\rm max}\rho_{\rm cr,0}$, GW contribution
to $\Delta N_{\rm eff}$ becomes

$\Delta N_{\rm eff}=\frac{8}{7}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{\rm
GW}}{\rho_{\gamma}}=\frac{120}{7\pi^{2}}\left(\frac{11}{4}\right)^{4/3}\frac{\rho_{\rm
cr,0}}{T_{0}^{4}}\Omega_{\rm GW}^{\rm max}.$ (4.24)

Finally, the GW relic abundance can be written in terms of $\Delta N_{\rm
eff}$ as:

$h^{2}\Omega_{\rm GW}^{\rm max}\lesssim 5.6269\times 10^{-6}\Delta N_{\rm
eff}.$ (4.25)

The Planck collaboration has already put a constraint $\Delta
N_{\text{eff}}<0.30$ at 95$\%$ confidence level [33]. Future experiments such
as CMB-S4 [34], Euclid [35], and potentially cosmic-variance-limited (CVL) CMB
polarization experiments are anticipated to improve the constraint on $\Delta
N_{\text{eff}}$ to $\lesssim 0.06$, $\lesssim 0.013$, and $\simeq 3.1\times
10^{-6}$, respectively [36]. These constraints are depicted with dark yellow
lines in Fig. 5. The GWs originating from the direct evaporation of PBH and
those stemming from the massive scalar particle bremsstrahlung decay align
with the current constraints on $\Delta N_{\rm eff}$. However, they can be
probed in the future CMB observations.

##  5 Conclusion

In this study, we have explored the production of stochastic gravitational
waves not only through the direct evaporation of PBHs but also via the
bremsstrahlung process during the decay of a massive scalar particle $\chi$.
The GWs spectrum resulting from PBH evaporation reaches it peak frequency
approximately at the Hawking temperature which tends to be considerably high
for PBH with masses smaller than $\sim 10^{8}$g as required for evaporation
before the onset of BBN. Additionally, we have examined the scenario of a
long-lived, super-heavy scalar particle emerging from PBH evaporation. This
particle ultimately decays into fermions, simultaneously emitting gravitons
through the bremsstrahlung process. The investigation into the bremsstrahlung
process as a novel source of GW production is especially intriguing for its
implications in high-frequency GW scenarios. While these high-frequency GWs
may be within the detection capabilities of proposed high-frequency detectors,
given their sensitivity projections, they remain way below the sensitivity of
the present detectors. Nonetheless, the distinct spectral signatures from both
PBH direct evaporation and the bremsstrahlung process associated with the
decay of the massive scalar particle offer promising avenues of detection by
the resonant cavity detectors. Moreover, GWs emanating from the decay of
massive scalar particles, characterized by sub-Planckian mass scales and small
Yukawa coupling, could significantly contribute to dark radiation. These
contributions are expected to fall within the detection thresholds of future
experiments, such as CMB-S4 and EUCLID, presenting an exciting frontier for
observational cosmology.

## Acknowledgments

The authors acknowledge the financial support from National Research
Foundation(NRF) grant funded by the Korea government (MEST)
NRF-2022R1A2C1005050.

## References

  * [1] LIGO Scientific, Virgo collaboration, _Observation of Gravitational Waves from a Binary Black Hole Merger_ , _Phys. Rev. Lett._ 116 (2016) 061102 [1602.03837].
  * [2] LIGO Scientific, Virgo collaboration, _GW170817: Observation of Gravitational Waves from a Binary Neutron Star Inspiral_ , _Phys. Rev. Lett._ 119 (2017) 161101 [1710.05832].
  * [3] NANOGrav collaboration, _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8 [2306.16213].
  * [4] EPTA, InPTA: collaboration, _The second data release from the European Pulsar Timing Array - III. Search for gravitational wave signals_ , _Astron. Astrophys._ 678 (2023) A50 [2306.16214].
  * [5] D.J. Reardon et al., _Search for an Isotropic Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 951 (2023) L6 [2306.16215].
  * [6] H. Xu et al., _Searching for the Nano-Hertz Stochastic Gravitational Wave Background with the Chinese Pulsar Timing Array Data Release I_ , _Res. Astron. Astrophys._ 23 (2023) 075024 [2306.16216].
  * [7] S.W. Hawking, _Black hole explosions_ , _Nature_ 248 (1974) 30.
  * [8] S.W. Hawking, _Particle Creation by Black Holes_ , _Commun. Math. Phys._ 43 (1975) 199.
  * [9] P. Gondolo, P. Sandick and B. Shams Es Haghi, _Effects of primordial black holes on dark matter models_ , _Phys. Rev. D_ 102 (2020) 095018 [2009.02424].
  * [10] N. Bernal and O. Zapata, _Dark Matter in the Time of Primordial Black Holes_ , _JCAP_ 03 (2021) 015 [2011.12306].
  * [11] B.J. Carr, _Some cosmological consequences of primordial black-hole evaporations_ , _Astrophys. J._ 206 (1976) 8.
  * [12] D. Baumann, P.J. Steinhardt and N. Turok, _Primordial Black Hole Baryogenesis_ , hep-th/0703250.
  * [13] A. Hook, _Baryogenesis from Hawking Radiation_ , _Phys. Rev. D_ 90 (2014) 083535 [1404.0113].
  * [14] Y. Hamada and S. Iso, _Baryon asymmetry from primordial black holes_ , _PTEP_ 2017 (2017) 033B02 [1610.02586].
  * [15] K.-Y. Choi, J. Kim and E. Lkhagvadorj, _Non-thermal WIMPy baryogenesis with primordial black hole_ , _JCAP_ 02 (2024) 020 [2309.16122].
  * [16] L. Morrison, S. Profumo and Y. Yu, _Melanopogenesis: Dark Matter of (almost) any Mass and Baryonic Matter from the Evaporation of Primordial Black Holes weighing a Ton (or less)_ , _JCAP_ 05 (2019) 005 [1812.10606].
  * [17] N. Smyth, L. Santos-Olmsted and S. Profumo, _Gravitational baryogenesis and dark matter from light black holes_ , _JCAP_ 03 (2022) 013 [2110.14660].
  * [18] B. Barman, D. Borah, S. Jyoti Das and R. Roshan, _Gravitational wave signatures of a PBH-generated baryon-dark matter coincidence_ , _Phys. Rev. D_ 107 (2023) 095002 [2212.00052].
  * [19] D.N. Page, _Particle Emission Rates from a Black Hole. 2. Massless Particles from a Rotating Hole_ , _Phys. Rev. D_ 14 (1976) 3260.
  * [20] B. Carr and F. Kuhnel, _Primordial black holes as dark matter candidates_ , _SciPost Phys. Lect. Notes_ 48 (2022) 1 [2110.02821].
  * [21] A. Escrivà, F. Kuhnel and Y. Tada, _Primordial Black Holes_ , 2211.05767.
  * [22] J. Ghiglieri, G. Jackson, M. Laine and Y. Zhu, _Gravitational wave background from Standard Model physics: Complete leading order_ , _JHEP_ 07 (2020) 092 [2004.11392].
  * [23] A. Ringwald, J. Schütte-Engel and C. Tamarit, _Gravitational Waves as a Big Bang Thermometer_ , _JCAP_ 03 (2021) 054 [2011.04731].
  * [24] J. Ghiglieri, J. Schütte-Engel and E. Speranza, _Freezing-in gravitational waves_ , _Phys. Rev. D_ 109 (2024) 023538 [2211.16513].
  * [25] G. Choi, W. Ke and K.A. Olive, _Minimal Production of Prompt Gravitational Waves during Reheating_ , 2402.04310.
  * [26] K. Nakayama and Y. Tang, _Stochastic Gravitational Waves from Particle Origin_ , _Phys. Lett. B_ 788 (2019) 341 [1810.04975].
  * [27] D. Huang and L. Yin, _Stochastic Gravitational Waves from Inflaton Decays_ , _Phys. Rev. D_ 100 (2019) 043538 [1905.08510].
  * [28] B. Barman, N. Bernal, Y. Xu and O. Zapata, _Gravitational wave from graviton Bremsstrahlung during reheating_ , _JCAP_ 05 (2023) 019 [2301.11345].
  * [29] N. Bernal, S. Cléry, Y. Mambrini and Y. Xu, _Probing reheating with graviton bremsstrahlung_ , _JCAP_ 01 (2024) 065 [2311.12694].
  * [30] S. Kanemura and K. Kaneta, _Gravitational Waves from Particle Decays during Reheating_ , 2310.12023.
  * [31] W. Hu, K. Nakayama, V. Takhistov and Y. Tang, _Gravitational Wave Probe of Planck-scale Physics After Inflation_ , 2403.13882.
  * [32] N. Herman, L. Lehoucq and A. Fúzfa, _Electromagnetic antennas for the resonant detection of the stochastic gravitational wave background_ , _Phys. Rev. D_ 108 (2023) 124009 [2203.15668].
  * [33] Planck collaboration, _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6 [1807.06209].
  * [34] CMB-S4 collaboration, _CMB-S4 Science Book, First Edition_ , 1610.02743.
  * [35] R. Laureijs, J. Amiaux, S. Arduini, J.L. Auguères, J. Brinchmann, R. Cole et al., _Euclid definition study report_ , 2011.
  * [36] I. Ben-Dayan, B. Keating, D. Leon and I. Wolfson, _Constraints on scalar and tensor spectra from $N_{eff}$_, _JCAP_ 06 (2019) 007 [1903.11843].
  * [37] Y.B. Zel’dovich and I.D. Novikov, _The Hypothesis of Cores Retarded during Expansion and the Hot Cosmological Model_ , _Sov. Astron._ 10 (1967) 602.
  * [38] B.J. Carr and S.W. Hawking, _Black holes in the early Universe_ , _Mon. Not. Roy. Astron. Soc._ 168 (1974) 399.
  * [39] B.J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _New cosmological constraints on primordial black holes_ , _Phys. Rev. D_ 81 (2010) 104019.
  * [40] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on primordial black holes_ , _Reports on Progress in Physics_ 84 (2021) 116902.
  * [41] B. Carr and F. Kuhnel, _Primordial Black Holes as Dark Matter: Recent Developments_ , _Ann. Rev. Nucl. Part. Sci._ 70 (2020) 355 [2006.02838].
  * [42] T. Fujita, K. Harigaya, M. Kawasaki and R. Matsuda, _Baryon asymmetry, dark matter, and density perturbation from primordial black holes_ , _Phys. Rev. D_ 89 (2014) 103501.
  * [43] I. Masina, _Dark matter and dark radiation from evaporating primordial black holes_ , _Eur. Phys. J. Plus_ 135 (2020) 552 [2004.04740].
  * [44] P. Gondolo, P. Sandick and B. Shams Es Haghi, _Effects of primordial black holes on dark matter models_ , _Phys. Rev. D_ 102 (2020) 095018.
  * [45] S.W. Hawking, _Black hole explosions?_ , _Nature_ 248 (1974) 30.
  * [46] S.W. Hawking, _Particle creation by black holes_ , _Communications in Mathematical Physics_ 43 (1975) 199.
  * [47] T. Ukwatta, D. Stump, J. Linnemann, J. MacGibbon, S. Marinelli, T. Yapici et al., _Primordial black holes: Observational characteristics of the final evaporation_ , _Astroparticle Physics_ 80 (2016) 90.
  * [48] C. Lunardini and Y.F. Perez-Gonzalez, _Dirac and Majorana neutrino signatures of primordial black holes_ , _JCAP_ 08 (2020) 014 [1910.07864].
  * [49] Y.F. Perez-Gonzalez and J. Turner, _Assessing the tension between a black hole dominated early universe and leptogenesis_ , _Phys. Rev. D_ 104 (2021) 103021 [2010.03565].
  * [50] A. Cheek, L. Heurtier, Y.F. Perez-Gonzalez and J. Turner, _Primordial black hole evaporation and dark matter production. i. solely hawking radiation_ , _Phys. Rev. D_ 105 (2022) 015022.
  * [51] I. Baldes, Q. Decant, D.C. Hooper and L. Lopez-Honorez, _Non-cold dark matter from primordial black hole evaporation_ , _Journal of Cosmology and Astroparticle Physics_ 2020 (2020) 045.
  * [52] Planck collaboration, _Planck 2018 results. X. Constraints on inflation_ , _Astron. Astrophys._ 641 (2020) A10 [1807.06211].
  * [53] G.F. Giudice, E.W. Kolb and A. Riotto, _Largest temperature of the radiation era and its cosmological implications_ , _Phys. Rev. D_ 64 (2001) 023508 [hep-ph/0005123].
  * [54] S. Jyoti Das, D. Mahanta and D. Borah, _Low scale leptogenesis and dark matter in the presence of primordial black holes_ , _JCAP_ 11 (2021) 019 [2104.14496].
  * [55] B. Barman, D. Borah, S.J. Das and R. Roshan, _Non-thermal origin of asymmetric dark matter from inflaton and primordial black holes_ , _JCAP_ 03 (2022) 031 [2111.08034].
  * [56] A. Granelli, K. Moffat, Y.F. Perez-Gonzalez, H. Schulz and J. Turner, _ULYSSES: Universal LeptogeneSiS Equation Solver_ , _Comput. Phys. Commun._ 262 (2021) 107813 [2007.09150].
  * [57] E. Thrane and J.D. Romano, _Sensitivity curves for searches for gravitational-wave backgrounds_ , _Phys. Rev. D_ 88 (2013) 124032.
  * [58] D. Reitze et al., _Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO_ , _Bull. Am. Astron. Soc._ 51 (2019) 035 [1907.04833].
  * [59] N. Seto, S. Kawamura and T. Nakamura, _Possibility of direct measurement of the acceleration of the universe using 0.1 hz band laser interferometer gravitational wave antenna in space_ , _Phys. Rev. Lett._ 87 (2001) 221103.
  * [60] KAGRA, LIGO Scientific, Virgo, VIRGO collaboration, _Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA_ , _Living Rev. Rel._ 21 (2018) 3 [1304.0670].
  * [61] N. Aggarwal et al., _Challenges and opportunities of gravitational-wave searches at MHz to GHz frequencies_ , _Living Rev. Rel._ 24 (2021) 4 [2011.12414].
  * [62] G. Mangano, G. Miele, S. Pastor, T. Pinto, O. Pisanti and P.D. Serpico, _Relic neutrino decoupling including flavor oscillations_ , _Nucl. Phys. B_ 729 (2005) 221 [hep-ph/0506164].