# Explaining PTA Data with Inflationary GWs in a PBH-Dominated Universe

Satyabrata Datta<EMAIL_ADDRESS>Saha Institute of Nuclear
Physics, 1/AF, Bidhannagar, Kolkata 700064, India

###### Abstract

We show that an ultralight primordial black hole (PBH) dominated phase makes
blue-tilted inflationary gravitational waves (BGW) compatible with the recent
detection of an nHz stochastic GW background by pulsar-timing arrays (PTAs),
for high reheating temperatures. This PBH-dominated phase suppresses the BGW
spectrum via entropy dilution and generates a new GW spectrum from PBH density
fluctuations. This combined spectrum is detectable at ongoing and planned
near-future GW detectors and exhibits a unique shape with a low-frequency peak
explaining PTA data, a mid-range dip, and a sharp peak followed by a third
peak at high-frequency. This distinctive shape sets it apart from spectra
generated by other matter dominations or exotic physics. Therefore, while
important for studying GWs in the nHz range, the recent PTA result also sets
the stage for testing and constraining various well-studied mechanisms
following a PBH domination, using low-frequency measurements and correlated
observations of unique high-frequency GW spectral features.

##  I  Introduction

Recently, collaborations of pulsar-timing arrays (PTAs) such as NANOGrav,
EPTA, and PPTA, along with InPTA and CPTA, have published their latest data.
This data presents substantial evidence for a stochastic gravitational wave
background (SGWB) at nHz frequencies [1, 2, 3, 4]. A similar discovery,
although with less statistical significance, has been present for the past two
years, generating considerable interest within the scientific community [5, 6,
7]. Intriguingly, this time, the signal displays the distinct angular
correlations of pulsars, referred to as the quadrupolar Hellings-Downs curve
[8], a feature that is specific to an SGWB. While the origins of such GWs are
still uncertain, the favored power-law $\Omega_{\rm GW}\propto f^{1.8\pm
0.6}$, for instance, in the new NANOGrav data, does not rule out the
straightforward GW-driven models of supermassive black hole binaries (SMBHBs)
at $3\sigma$. Yet another intriguing prospect is to explore the GWs that
originate from cosmological sources. In a companion theory paper [9], the
NANOGrav collaboration (we shall specifically focus on NANOGrav 15 yrs data
[1]; the results of other PTAs align well) has compiled a comprehensive list
discussing numerous cosmological sources that are consistent with the data
111Contrary to the 12.5 years of NANOGrav data, which was well-matched by
stable cosmic strings [10, 11, 12, 13, 14, 15], the latest data seems to
challenge this fit, suggesting that stable cosmic strings may not be the best
explanation [1].. Following this, various articles have explored these sources
either within the framework of different cosmological models or by re-
analyzing the fit to the new data, incorporating the results from other PTAs
[16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 22, 23, 30, 31, 32,
33, 34, 35, 36, 37, 38, 39, 40, 41, 17, 42, 43, 44, 21, 45, 46, 47, 48, 49,
50, 51, 52, 53, 48, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63]. Inflationary GWs
with a significant tensor blue-tilt, also known as Blue Tilted Gravitational
Waves (BGWs), are a notable aspect that aligns well with both past and recent
data [2, 1, 64, 65, 66, 67, 68]; although, It’s important to mention that such
BGWs can be generated in models that may not necessarily align with the
conventional slow-roll inflation paradigm, see, e.g., [69, 70, 71, 72, 73, 74,
75, 76]. However, the range of parameters that allow for such a fit is, in
fact, limited. This is primarily because GWs with large blue-tilt saturate big
bang nucleosynthesis (BBN) bound on the effective number of neutrino species,
disfavoring any post-inflationary cosmology founded on high reheating
temperature ($T_{\rm RH}\gtrsim 10$ GeV) after inflation [9, 68]. Nonetheless,
if a non-standard matter epoch leads to entropy production between the
reheating after inflation and the most recent radiation domination before the
BBN [77], BGWs get suppressed and provide a good fit to PTA data for high
reheating temperatures. Now, contrary to the standard case, such a scenario
allows the overall GW spectrum to span decades of frequencies with
characteristic spectral features testable by high-frequency detectors, e.g.,
by the future LIGO run [78, 79, 80].

Building on the methodology outlined in the ref.[81, 82], we carried out a
tomographic analysis of early matter domination (EMD) resulting from PBHs and
their imprint on BGWs. While previous research has explored the impact of EMD
induced by PBHs on the GW spectrum generated by cosmic strings [13, 14, 15,
83, 84], no studies have specifically focused on BGWs and GWs originating from
PBHs. Such ultralight PBH domination has been extensively studied within the
context of Beyond Standard Model (BSM) physics scenarios, such as superheavy
dark matter (DM) [14] and leptogenesis [13, 15].

Motivated by these, we have undertaken a systematic investigation into the
impact of a PBH-dominated epoch on BGWs, with a particular focus on explaining
the results of the NANOGrav at low frequencies. We have demonstrated that the
characteristic features of PBHs, specifically their mass and initial energy
fraction, leave a mark on BGWs. This necessitates a larger initial fraction,
which in turn leads to a higher entropy injection. This higher entropy
injection suppresses the BGWs at high frequencies to ensure compliance with
the constraints imposed by the LIGO and BBN. Among the various mechanisms
through which PBHs can generate GWs, we focus on the one induced by PBH
density fluctuations [85, 86, 87, 88, 89, 90]. While a requirement for a
higher initial energy fraction of PBHs appears to conflict directly with GWs
from density fluctuations saturating the constraints imposed by the LIGO and
BBN, we have successfully identified a complementary region. This region not
only addresses PTA data at low frequency and a peak at high frequency, but it
also presents a unique, sharply peaked GW spectrum at mid-frequencies due to
PBH density fluctuations. Thus, the overall spectrum introduces a rich
phenomenology with characteristic features that can be explored in a
complementary manner by future gravitational wave detectors.

The structure of the paper is as follows. In Section II, we provide a brief
discussion on GWs from inflation with tensor blue tilt. Section III is
dedicated to the dynamics of PBHs and the necessary components for studying
their effect on BGWs. In Section IV, we conduct a numerical analysis,
prioritizing the satisfaction of NANOGrav data and examining PBH archaeology
with BGWs. Section V discusses a unique possibility of detecting a
complementary signal from BGWs and GWs resulting from PBH density
fluctuations. Finally, in Section VI, we conclude our findings.

##  II  Blue-tilted GWs from inflation

One of the most plausible explanations for the origin of primordial GWs is
cosmic inflation [91, 92]. In this section, we will briefly discuss how GWs
are produced during inflation and how they travel through different cosmic
eras until they reach the present day. GWs are described as a perturbation in
the FLRW line element:

$\displaystyle
ds^{2}=a(\tau)\left[-d\tau^{2}+(\delta_{ij}+h_{ij})dx^{i}dx^{j})\right],$
(II.1)

with $\tau$ being the conformal time, and $a(\tau)$ the scale factor. The GWs
are interpreted by the transverse and traceless ($\partial_{i}h^{ij}=0$,
$\delta^{ij}h_{ij}=0$) part of the of the $3\times 3$ symmetric matrix
$h_{ij}$. Since the GWs are too feeble, $|h_{ij}|\ll 1$, the linearized
evolution equation

$\displaystyle\partial_{\mu}(\sqrt{-g}\partial^{\mu}h_{ij})=16\pi
a^{2}(\tau)\mathcal{\pi}_{ij}$ (II.2)

is sufficient for studying their propagation. The tensor part of the
anisotropy stress, $\pi_{ij}$, coupled with $h_{ij}$ acts as an external
source. It is useful to express $h_{ij}$ in the Fourier space:

$\displaystyle
h_{ij}(\tau,\vec{x})=\sum_{\lambda}\int\frac{d^{3}\vec{k}}{(2\pi)^{3/2}}e^{i\vec{k}.\vec{x}}\epsilon_{ij}^{\lambda}(\vec{k})h_{\vec{k}}^{\lambda}(\tau),$
(II.3)

where the index $\lambda=``+/-"$ denotes that the GWs have two polarisation
states. The polarization tensors, in addition to being transverse and
traceless, also fulfill the conditions:

$\displaystyle\text{(i)}\>\epsilon^{(\lambda)ij}(\vec{k})\epsilon_{ij}^{(\lambda^{\prime})}(\vec{k})=2\delta_{\lambda\lambda^{\prime}}$
(II.4)
$\displaystyle\text{(ii)}\>\epsilon^{(\lambda)}_{ij}(-\vec{k})=\epsilon^{(\lambda)}_{ij}(\vec{k}).$

Assuming that each polarization state evolves identically and isotropically,
we can rename the notation by letting $h_{\vec{k}}^{\lambda}(\tau)$ as
$h_{k}(\tau)$, where $k=|\vec{k}|=2\pi f$ with $f$ being the frequency of the
GWs today at $a_{0}=1$. Taking into account the sub-dominant contribution from
$\pi_{ij}$, the equation for the propagation of GWs in Fourier space can be
written as

$\displaystyle\ddot{h}_{k}+2\frac{\dot{a}}{a}\dot{h}_{k}+k^{2}h_{k}=0,$ (II.5)

where the dot represents a derivative with respect to the conformal time. By
utilizing Eq.(II.3) and Eq.(II.5), we can compute the energy density of the
GWs as [93]

$\displaystyle\rho_{\rm GW}=\frac{1}{32\pi
G}\int\frac{dk}{k}\left(\frac{k}{a}\right)^{2}T_{T}^{2}(\tau,k)P_{T}(k),$
(II.6)

where $T_{T}^{2}(\tau,k)=|h_{k}(\tau)|^{2}/|h_{k}(\tau_{i})|^{2}$ is a
transfer function that is derived from Eq.(II.5), with $\tau_{i}$ being an
initial conformal time. The primordial power spectrum,
$P_{T}(k)=\frac{k^{3}}{\pi^{2}}|h_{k}(\tau_{i})|^{2}$ is linked to the
inflation models with specific forms and is parametrised as a power-law given
by

$\displaystyle P_{T}(k)=rA_{s}(k_{*})\left(\frac{k}{k_{*}}\right)^{n_{T}},$
(II.7)

where $r\lesssim 0.06$ [94] is the tensor-to-scalar-ratio, $A_{s}\simeq
2\times 10^{-9}$ represents the scalar perturbation amplitude at the pivot
scale $k_{*}=0.01\rm Mpc^{-1}$ and the tensor spectral index is denoted by
$n_{T}$. Interestingly, the standard slow-roll inflation models satisfy a
consistency relation $n_{T}=-r/8$ [95], which results in slightly red-tilted
GWs with $n_{T}\lesssim 0$. However, in this work, we have considered GWs with
a significant blue tilt, where $n_{T}>0$, and we have assumed it to be
constant throughout. The GW energy density, which is crucial for detection
purposes, can be expressed as

$\displaystyle\Omega_{\rm GW}(k)=\frac{k}{\rho_{c}}\frac{d\rho_{\rm GW}}{dk},$
(II.8)

where the quantity $\rho_{c}=3H_{0}^{2}/8\pi G$ with $H_{0}\simeq 2.2\times
10^{-4}~{}\rm Mpc^{-1}$ being the present-day Hubble constant. From Eq.(II.6),
the quantity $\Omega_{\rm GW}(k)$ is computed as

$\displaystyle\Omega_{\rm
GW}(k)=\frac{1}{12H_{0}^{2}}\left(\frac{k}{a_{0}}\right)^{2}T_{T}^{2}(\tau_{0},k)P_{T}(k),\text{
with}~{}~{}\tau_{0}=1.4\times 10^{4}{\rm~{}Mpc}.$ (II.9)

There have been various attempts to compute the transfer function analytically
[96, 97, 98, 99], and one of the commonly used one for standard reheating is
given by[100, 101]

$\displaystyle T_{T}^{2}(\tau_{0},k)=F(k)T_{1}^{2}(\zeta_{\rm
eq})T_{2}^{2}(\zeta_{R}),$ (II.10)

where $F(k)$ reads

$\displaystyle F(k)=\Omega_{m}^{2}\left(\frac{g_{*}(T_{k,\rm
in})}{g_{*0}}\right)\left(\frac{g_{*s0}}{g_{*s}(T_{k,\rm
in})}\right)^{4/3}\left(\frac{3j_{1}(k\tau_{0})}{k\tau_{0}}\right)^{2}.$
(II.11)

Here $j_{1}(k\tau_{0})$ is the spherical Bessel function, $\Omega_{m}=0.31$,
$g_{*0}=3.36$, and $g_{*0s}=3.91$. The scale-dependent $g_{*0(s)}(T_{k,\rm
in})$, used in Eq.(II.11) can be approximated analytically as [101, 102, 103]

$\displaystyle g_{*0(s)}(T_{k,\rm in})=g_{*0}\left(\frac{A+{\rm
tanh~{}k_{1}}}{A+1}\right)\left(\frac{B+{\rm tanh~{}k_{2}}}{B+1}\right),$
(II.12)

where

$\displaystyle
A=\frac{-1-10.75/g_{*0(s)}}{-1+10.75/g_{*0(s)}},~{}~{}B=\frac{-1-g_{max}/10.75}{-1+g_{max}/10.75},$
(II.13)

and

$\displaystyle k_{1}=-2.5~{}{\rm log}_{10}\left(\frac{k/2\pi}{2.5\times
10^{-12}{\rm Hz}}\right),$ (II.14) $\displaystyle k_{2}=-2.0~{}{\rm
log}_{10}\left(\frac{k/2\pi}{6.0\times 10^{-9}{\rm Hz}}\right).$ (II.15)

The transfer functions are given by

$\displaystyle T_{1}^{2}(\zeta)=1+1.57\zeta+3.42\zeta^{2},$ (II.16)
$\displaystyle
T_{2}^{2}(\zeta)=\left(1-0.22\zeta^{1.5}+0.65\zeta^{2}\right)^{-1},$ (II.17)

where $\zeta_{i}\equiv k/k_{i}$, with $k_{i}$s being the wave number of the
modes entering the horizon at different epochs and are derived as

$\displaystyle k_{\rm eq}$ $\displaystyle=$ $\displaystyle 7.1\times
10^{-2}\Omega_{m}h^{2}{\rm Mpc^{-1}}.$ (II.18)

In this paper, we shall use $k_{*}=0.01$ Mpc-1 and $h=0.7$. It’s essential to
highlight that a significant limitation to consider regarding the potential
for BGWs is the $\Delta N_{\rm eff}$ bound from BBN, as well as the absence of
any SGWB detection by LIGO [79, 80]. In the following section, we will
demonstrate that if any late-time entropy production occurs through the
evaporation of ultralight PBHs after significant PBH domination before
reheating, it can considerably alter the transfer function during the PBH-
dominated era and suppress the GW spectrum for modes that entered the horizon
during the PBH-dominated phase. This could potentially ameliorate the
constraints from BBN and LIGO.

##  III  Diluting BGWs through Entropy injection in PBH Domination

The dynamical evolution of energy densities of the black holes ($\rho_{\rm
BH}$), and radiation ($\rho_{\rm R}$) is governed by the following Friedmann
equations [13, 104]:

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

where $H$ is the Hubble parameter and $z=T_{\rm Bf}/T$. The scale factor $a$
and the quantity $\tilde{H}$ evolve as

$\displaystyle\tilde{H}=\left(H+\mathcal{K}\right),~{}\frac{da}{dz}=\left(1-\frac{\mathcal{K}}{\tilde{H}}\right)\frac{a}{z},$
(III.3)

where $\mathcal{K}=\frac{\dot{M}_{\rm BH}}{M_{\rm BH}}\frac{\rho_{\rm
BH}}{4\rho_{\rm R}}$. In the derivation of Eq.(III.1)-Eq.(III.3), we have made
the assumption that entropy ($g_{*s}$) and the energy ($g_{*\rho}$) degrees of
freedom are equal and constant. If the initial energy fraction of PBHs,
denoted as $\beta\equiv\frac{\rho_{\rm BH}(T_{\rm Bf})}{\rho_{\rm R}(T_{\rm
Bf})}$, surpasses a critical value of
$\beta_{c}\equiv\gamma^{-1/2}\sqrt{\left({\mathcal{G}g_{*B}(T_{\rm
BH})}/{10240\pi}\right)}\frac{M_{\rm Pl}}{M_{\rm BH}}$, it can result in early
matter domination. For a given $\beta$, the above equations can be solved
numerically to determine the temperatures of PBH domination and evaporation,
as well as the entropy production $\Delta_{\rm
PBH}=\tilde{S}_{2}/\tilde{S}_{1}$, where $\tilde{S}_{1,\>(2)}\propto
a_{1,\>(2)}^{3}/z_{1,\>(2)}^{3}$ is the total entropy before (after) the PBH
evaporation. We will ultimately utilize analytical expressions that are in
close agreement with numerical results and can be approximated accordingly
[15]

$\displaystyle\Delta_{\rm PBH}\simeq 233\beta\left(\frac{M_{\rm BH}}{M_{\rm
Pl}}\right)\left(\frac{\gamma}{g_{*B}(T_{\rm BH})\mathcal{G}}\right)^{1/2},$
(III.4)

where $\gamma\simeq 0.2$ is the formation efficiency of PBHs, $g_{*B}\simeq
100$ is the no of relativistic d.o.f. below the Hawking temperature $T_{\rm
BH}$, and $\mathcal{G}\simeq 3.8$ is the greybody factor.

In the present scenario where late-time entropy production takes place after
reheating through PBH evaporation, the background evolution follows the
sequence of MD (inflation-dominated) $\rightarrow$ RD $\rightarrow$ MD (PBH-
dominated) $\rightarrow$ RD. In such a case, the transfer function is defined
by

$\displaystyle T_{T}^{2}(\tau_{0},k)=F(k)T_{1}^{2}(\zeta_{\rm
eq})T_{2}^{2}(\zeta_{\text{ev}})T_{3}^{2}(\zeta_{\text{dom}})T_{2}^{2}(\zeta_{R}),$
(III.5)

where $T_{3}^{2}(\zeta)$ is the transfer function corresponding to the
transition from the first RD phase to the PBH-dominated phase and is given by

$\displaystyle T_{3}^{2}(\zeta)=1+0.59\zeta+0.65\zeta^{2},$ (III.6)

and the modes $k_{i}$’s given by

$\displaystyle k_{\text{ev}}$ $\displaystyle=$ $\displaystyle 1.7\times
10^{14}\left(\frac{g_{*s}(T_{\text{ev}})}{106.75}\right)^{1/6}\left(\frac{T_{\text{ev}}}{10^{7}\rm
GeV}\right){\rm Mpc^{-1}},$ (III.7) $\displaystyle k_{\text{dom}}$
$\displaystyle=$ $\displaystyle 1.7\times
10^{14}\left(\frac{g_{*s}(T_{\text{dom}})}{106.75}\right)^{1/6}\left(\frac{T_{\text{dom}}}{10^{7}\rm
GeV}\right){\rm Mpc^{-1}},$ (III.8)

and

$\displaystyle k_{\rm R}=1.7\times 10^{14}\Delta_{\rm
PBH}^{-1/3}\left(\frac{g_{*s}(T_{\rm
R})}{106.75}\right)^{1/6}\left(\frac{T_{\rm R}}{10^{7}\rm GeV}\right){\rm
Mpc^{-1}}$ (III.9)

which reenters the horizon at PBH evaporation temperature ($T_{\rm ev}$), PBH
domination temperature ($T_{\rm dom}$), and inflationary reheating temperature
($T_{\rm R}$), respectively where [13]

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

and

$\displaystyle T_{\rm dom}=\beta T_{\rm Bf},$ (III.11)

where $\tau$, and $T_{\rm Bf}$ correspond to the BH lifetime and BH formation
temperature respectively, and are given by

$\displaystyle\tau=\frac{10240\pi M_{\rm BH}^{3}}{\mathcal{G}g_{*B}(T_{\rm
BH})M_{\rm Pl}^{4}},$ (III.12)

and

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

Note also that, by construction, the formation of BHs occurs after the initial
inflationary reheating, i.e. $T_{\rm Bf}\lesssim T_{\rm R}$. For simplicity,
we shall use $T_{\rm R}=T_{\rm Bf}$ throughout the rest of our discussion.

##  IV A fit to the NANOGrav data in a PBH-dominated scenario and PBH
archaeology with BGWs

Taking into account the aforementioned equations and utilizing $\Delta_{\rm
PBH}$ from Eq.(III.4), we evaluate Eq.(II.9) to get the GW spectrum and to fit
the NANOGrav data, while considering two additional constraints. I) The LIGO
O3 constraint on SGWB, which can be approximately expressed as $\Omega_{\rm
GW}(25\text{\>Hz})h^{2}\leq 2.2\times 10^{-9}$ [78], and II) BBN bound on the
effective number of neutrino species: $\int_{f_{\rm low}}^{f_{\rm
high}}f^{-1}df\Omega_{\rm GW}(f)h^{2}\lesssim 5.6\times 10^{-6}\Delta N_{\rm
eff}$, where $\Delta N_{\rm eff}\lesssim 0.2$ [105]. The lower limit of the
integration is associated with the frequency that entered the horizon during
the BBN epoch. Conversely, the Hubble rate at the end of inflation sets the
upper limit: $f_{\rm high}=a_{\rm end}H_{\rm end}/2\pi$. For numerical
calculations, $f_{\rm high}\simeq 10^{5}$ Hz is sufficient as the spectrum
falls and the integration saturates.

Figure 1: Top: The red $\bullet$ represents the amplitude $A$ and the spectral
index $\gamma$ corresponding to a benchmark ($n_{T}$, $r$)$\equiv$ ($1.4$,
$3\times 10^{-6}$). NANOGrav $1,2,3\sigma$ regions are shown by the cyan
shades.

We follow the NANOGrav parametrization for the GW energy density to conduct a
power-law fit to the new data within the frequency range
$f\in\left[10^{-9}~{}{\rm Hz},f_{\rm yr}\right]$. The parametrization reads

$\Omega_{\rm GW}(f)=\Omega_{\rm yr}\left(\frac{f}{f_{\rm
yr}}\right)^{(5-\gamma)}$ (IV.1)

with $\Omega_{\rm yr}=\frac{2\pi^{2}}{3H_{0}^{2}}A^{2}f_{\rm yr}^{2}$ and
$f_{\rm yr}=1{\rm yr}^{-1}\simeq 32$ nHz. To fit the data, one must compare
Eq.(II.9) and Eq.(IV.1), then extract the values of the amplitude $A$ and the
spectral index $\gamma$ that fall within the $1,\>2,\>3\sigma$ contours as
reported by the NANOGrav [1]. In Fig.1, we present a particular benchmark
where $n_{T}=1.4$ and $r=3\times 10^{-6}$ for which the BGW spectrum
intersects with the 1$\sigma$ range of NANOGrav. For this specific benchmark
point, the parameters that are vital to PBH dynamics, namely $\beta$ and
$M_{\rm BH}$, have been adjusted freely. Due to this flexibility, for a wide
range of $n_{T}$ and $r$, the BGW spectrum, which was previously ruled out by
stringent constraints from LIGO and BBN at higher frequencies, can now be
permitted by choosing appropriate values for $\beta$ and $M_{\rm BH}$. To
provide a quantitative understanding of the impact of PBH domination, we
consider the same benchmark point $n_{T}=1.4$ and $r=3\times 10^{-6}$ as a
reference.

In the top/bottom left plot of Fig.2, we varied $\beta$ for three distinct
values: $10^{-6}$ (depicted in red), $10^{-5}$ (depicted in brown), and
$10^{-4}$ (depicted in purple), while keeping $M_{\rm BH}=10^{7}$ grams fixed.
Conversely, in the top/bottom right plot of Fig.2, we varied $M_{\rm BH}$ for
three different values: $10^{5.5}$ grams (shown in brown), $10^{6}$ grams
(shown in red), and $10^{7}$ grams (shown in purple) with $\beta=3\times
10^{-5}$ held constant.

Interestingly, it is possible to find the analytical expressions of three
turning point frequencies as a function of $\beta$ and $M_{\rm BH}$ (c.f.
Eq.(III.7), (III.8), (III.9)) and these are shown as follows:

$\displaystyle f_{\rm ev}=4.88\times 10^{-10}M_{\rm
Pl}\left(\frac{g_{\rm*s}(T_{\rm
ev})}{106.75}\right)^{-1/12}\left(\frac{g_{\rm*B}}{100}\right)^{1/2}\left(\frac{M_{\rm
Pl}}{M_{\rm BH}}\right)^{3/2},$ (IV.2) $\displaystyle f_{\rm dom}=2.04\times
10^{-17}M_{\rm
Pl}\left(\frac{\beta}{10^{-8}}\right)\left(\frac{g_{\rm*B}}{100}\right)^{-1/4}\left(\frac{g_{\rm*s}(T_{\rm
dom})}{106.75}\right)^{1/6}\left(\frac{M_{\rm Pl}}{M_{\rm BH}}\right)^{1/2},$
(IV.3) $\displaystyle f_{\rm R}=5.42\times 10^{-7}M_{\rm
Pl}\left(\frac{\beta}{10^{-8}}\right)^{-1/3}\left(\frac{g_{\rm*B}}{100}\right)^{-1/12}\left(\frac{g_{\rm*s}(T_{\rm
R})}{106.75}\right)^{1/6}\left(\frac{M_{\rm Pl}}{M_{\rm BH}}\right)^{5/6}.$
(IV.4)

Figure 2: Top-left: Evolution of the radiation energy density (solid), PBH
energy density (dashed), and the total entropy (dot-dashed) with the inverse
temperature. The benchmark points are represented in different colors: red for
$\beta=10^{-6}$, brown for $\beta=10^{-5}$, and purple for $\beta=10^{-4}$. In
all these cases, the PBH mass, $M_{\rm BH}$ is fixed at $10^{7}$ grams. Top-
right: Evolution of the same quantities, but for different PBH masses and
fixed initial energy fraction. Bottom-left/right: The BGW spectra comply with
the quantities that determine the nature of the corresponding plots on the
top.

From the inspection of the subplots in Fig.2 and Eqs.(II.9), (III.5), (IV.2),
(IV.3), (IV.4), we can extract the following pieces of information, I) The
variation of $\beta$ has a significant impact on the second peak and the mid-
frequency dip, i.e., Eq.(IV.4) and (IV.3), respectively. For larger values of
$\beta$, the high-frequency peak shifts towards lower frequency according to
Eq.(IV.4). Additionally, the dip shifts towards higher frequency as per
Eq.(IV.3). Furthermore, with increasing $\beta$, the BGW amplitude at high
frequency gets suppressed while the amplitude at low frequency remains
unaffected. The rationale for this can be readily comprehended by examining
the top-left plot of Fig.2. As $\beta$ increases, the strength of matter
domination will be much stronger, leading to a larger entropy injection that
dilutes BGWs at high frequency. The evaporation temperature (see Eq.(III.10)),
which is not influenced by $\beta$, indicates that the BGW amplitude at low-
frequency peak (c.f. Eq.(III.5)), which depends on the evaporation
temperature, remains unchanged. II) Altering $M_{\rm BH}$ impacts the entire
BGW spectrum and all the turning point frequencies. As can be observed from
the top-right plot of Fig.2, an increase in $M_{\rm BH}$ extends the lifetime
of PBHs (See Eq.(III.12)). Consequently, this leads to a larger entropy
injection, resulting in a more diluted BGW spectrum.

We end this section by highlighting a potential challenge: in the future, if
GW detectors identify such inflationary double-peak GWs, we may not be able to
distinguish between ultralight PBH domination and other intermediate matter
domination scenarios. However, there is an intriguing solution to this issue.
For certain PBH model parameters, as we will demonstrate next, they can
generate their own GWs from PBH density fluctuations.

##  V Detection Prospects for the Unique GW Spectrum(BGW+PBH density
fluctuations) and Realistic BSM Physics Scenarios

It’s noteworthy to state that ultralight PBHs play multiple roles in the
generation of GWs. For instance, the initial curvature perturbations that lead
to the formation of PBHs also give rise to GWs (see, e.g., [106, 107, 108]),
PBHs emit gravitons, which constitute high-frequency GWs [109], PBHs also
merge, releasing GWs [110], Lastly, the inhomogeneous distribution of PBHs,
which results in density fluctuations, triggers the production of GWs [87, 88,
89]. In this section, we will concentrate on the last point mentioned above.
As recently highlighted in ref.[87] and further elaborated in refs. [88, 89],
it’s observed that PBHs, immediately after their formation, are distributed
randomly in space following Poisson distribution. Hence, even though the PBH
gas, on average, behaves like pressure-less dust, the uneven spatial
distribution results in density fluctuations, which are inherently
isocurvature. When PBHs start to dominate the energy density of the Universe,
the isocurvature component transitions into curvature perturbations, which
subsequently give rise to secondary GWs. Due to the significant density
fluctuations at smaller scales (equivalent to the average separation of PBHs
at $T_{\rm Bf}$), substantial GWs are generated.

These GWs are further amplified due to the nearly instantaneous evaporation of
PBHs. The amplitude of such induced GWs in the present day is given by222Note
that the amplitude of the induced GWs is inherently independent of the history
of PBH formation and, by design, does not rely on any non-standard inflation
dynamics [88, 89].

$\displaystyle\Omega_{\rm GW}^{\rm PBH}(t_{0},f)\simeq\Omega^{p}_{\rm
GW}\left(\frac{f}{f_{p}}\right)^{11/3}\Theta(f_{p}-f),$ (V.1)

where the peak amplitude is given by

$\displaystyle\Omega^{p}_{\rm GW}\simeq 2\times
10^{-6}\left(\frac{\beta}{10^{-8}}\right)^{16/3}\left(\frac{M_{\rm
BH}}{10^{7}g}\right)^{34/9},$ (V.2)

and the peak frequency is given by

$\displaystyle f_{p}\simeq 1.7\times 10^{3}\text{\>Hz}\left(\frac{M_{\rm
BH}}{10^{4}g}\right)^{-5/6}.$ (V.3)

The ultraviolet cutoff represented by the $\Theta$ function, which is on par
with the frequency that corresponds to the comoving scale, signifies the
average separation of PBHs at their time of formation.

Figure 3: Left: Allowed region of combined GW spectra from BGW and PBH density
fluctuations for $n_{T}=0.9$ and $r=0.06$, Right: Combined GW spectrum
subjected to the benchmark point ($M_{\rm BH}=7.08\times 10^{6}$gram,
$\beta=10^{-8}$)

As previously demonstrated, PBHs with $\beta>\beta_{c}$ have the potential to
cause spectral distortion in the BGWs from inflation due to the domination and
subsequent evaporation of PBHs. Concurrently, this also gives rise to a new
source of GWs originating from density fluctuations leading to a complex GW
spectral shape characterized by three distinct peak features. The left plot in
Fig.3 shows a red-shaded region that is excluded due to the condition
$\Omega_{\rm BGW}(f^{p})>\Omega_{\rm GW}^{\rm PBH}(f^{p})$. This implies that
the characteristic peak from PBH density fluctuations will always be
suppressed below the BGW spectrum, rendering the PBH scenario
indistinguishable from any other intermediate matter domination. We present an
allowed parameter space on the $\beta$-$M_{\rm BH}$ plane for a benchmark
point $n_{T}=0.9$ and $r=0.06$. Although this does not provide the best fit to
the NANOGrav data, it fits the data at $2\sigma$. The green region is allowed,
satisfying both the required frequency and amplitude by NANOGrav and the bound
on SGWB from LIGO and BBN. Interestingly, throughout the entire parameter
space, the condition $\Omega_{\rm BGW}(f^{p})<\Omega_{\rm GW}^{\rm
PBH}(f^{p})$ holds true. This allows for the possibility of a characteristic
GW signal from PBHs that exhibits three peaks. For this specific benchmark
point, we have also ensured that the inequality $f_{\rm dom}<f^{p}<f_{\rm R}$
holds true. An important point to note is that while the parameters $n_{T}$
and $r$ can be varied to identify new allowed regions that fit the NANOGrav
data more exhaustively, there is no qualitative difference. The allowed values
of $n_{T}$ and $r$ are close to the presented benchmark. This is because, as
seen in the bottom left plot of Fig.2, higher $\beta$ values are needed to
suppress the BGW spectrum at higher frequencies to satisfy the BBN and LIGO
constraints. However, from Eq.(V.2), it can be seen that the GW peak amplitude
from PBH density fluctuations increases with $\beta$. These two features, in
conjunction with the need for consistency with the NANOGrav data, strongly
constrain the parameter space.

Considering the broad overview provided earlier, we can highlight how such a
spectrum could serve as an indicator for many BSM models in particle physics.
Firstly, while such light PBHs may ultimately evaporate, they could generate
relics that are either stable or unstable. A stable relic could potentially
represent dark matter, while an unstable particle, such as right-handed
neutrinos in the seesaw mechanism, that emerges from PBH could potentially
trigger baryogenesis via leptogenesis (see, e.g., refs.[111, 112, 113, 114]).
Moreover, based on the mass and initial energy fraction, a universe dominated
by PBHs can modify the standard parameter space of numerous particle physics
models, such as dark matter models [115]. Therefore, PBHs (their mass and
initial energy fraction) serve as a gateway linking GWs to the parameters of
high-energy particle physics models. This further demands a unified
investigation of these models using both GWs and conventional particle physics
experiments which is beyond the scope of this work. To conclude, we would like
to emphasize that one of the main motivations of this paper is to show that
the effects of an EMD caused by the standard long-lived fields or some
additional new physics [116, 81, 117] and PBHs could be very different.
Specifically, distorting BGWs with EMD leads to a low-frequency peak followed
by a dip and a high-frequency peak in both scenarios. However, in the case of
EMD from PBH, depending on the initial energy fraction of PBHs, a third
distinct sharp peak might appear at mid-frequency. This sets the PBH scenarios
apart. Furthermore, the impacts of ultra-light PBHs, which are highly
intriguing to investigate, for instance, in the context of dark matter
production, baryogenesis, and axions, could make these spectral features in
the GWs a unique indicator of these models, regardless of the strength of the
particle physics coupling.

##  VI Conclusion

We have discussed a unique framework that uses NANOGrav PTA data, interpreted
as a stochastic gravitational wave background (SGWB) from inflation, to
investigate an epoch dominated by primordial black holes (PBH). With an
appropriate choice of PBH parameters, a PBH-dominated epoch can accommodate
almost any value of the inflationary parameters, fitting the GW spectrum at
the NANOGrav frequency while also satisfying the stringent high-frequency
constraints from LIGO and BBN. The presence of any significant intermediate
matter domination (here PBH domination) after inflation results in a double-
peaked GW spectrum with a dip in between, while the introduction of ultralight
PBHs contributes an additional GW spectrum from density fluctuations. The
combined GW spectrum exhibits a unique shape, characterized by a low-frequency
peak that explains the PTA data, a dip in the middle, and a sharp tilted peak,
followed by a third high-frequency peak. Such a combined spectrum, along with
the PTA fit and constraints from LIGO and BBN, predicts a highly constrained
PBH parameter space that remains distinguishable from any other intermediate
matter domination scenario. In addition to the distinct features that can be
verified in GW detectors, the framework discussed in our work can also lead to
a plethora of beyond the Standard model phenomenological implications, such as
the production of dark matter from PBH evaporation and high-scale
leptogenesis, among others.

## acknowledgement

The author would like to thank Rome Samanta for useful discussions and careful
reading of the manuscript.

## References

  * [1] NANOGrav collaboration, G. Agazie et al., _The NANOGrav 15 yr Data Set: Evidence for a Gravitational-wave Background_ , _Astrophys. J. Lett._ 951 (2023) L8, [2306.16213].
  * [2] EPTA collaboration, J. Antoniadis et al., _The second data release from the European Pulsar Timing Array III. Search for gravitational wave signals_ , 2306.16214.
  * [3] 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].
  * [4] 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].
  * [5] EPTA collaboration, S. Chen et al., _Common-red-signal analysis with 24-yr high-precision timing of the European Pulsar Timing Array: inferences in the stochastic gravitational-wave background search_ , _Mon. Not. Roy. Astron. Soc._ 508 (2021) 4970–4993, [2110.13184].
  * [6] NANOGrav collaboration, Z. Arzoumanian et al., _The NANOGrav 12.5 yr Data Set: Search for an Isotropic Stochastic Gravitational-wave Background_ , _Astrophys. J. Lett._ 905 (2020) L34, [2009.04496].
  * [7] B. Goncharov et al., _On the Evidence for a Common-spectrum Process in the Search for the Nanohertz Gravitational-wave Background with the Parkes Pulsar Timing Array_ , _Astrophys. J. Lett._ 917 (2021) L19, [2107.12112].
  * [8] R. w. Hellings and G. s. Downs, _UPPER LIMITS ON THE ISOTROPIC GRAVITATIONAL RADIATION BACKGROUND FROM PULSAR TIMING ANALYSIS_ , _Astrophys. J. Lett._ 265 (1983) L39–L42.
  * [9] NANOGrav collaboration, A. Afzal et al., _The NANOGrav 15 yr Data Set: Search for Signals from New Physics_ , _Astrophys. J. Lett._ 951 (2023) L11, [2306.16219].
  * [10] S. Blasi, V. Brdar and K. Schmitz, _Has NANOGrav found first evidence for cosmic strings?_ , _Phys. Rev. Lett._ 126 (2021) 041305, [2009.06607].
  * [11] J. Ellis and M. Lewicki, _Cosmic String Interpretation of NANOGrav Pulsar Timing Data_ , _Phys. Rev. Lett._ 126 (2021) 041304, [2009.06555].
  * [12] R. Samanta and S. Datta, _Gravitational wave complementarity and impact of NANOGrav data on gravitational leptogenesis_ , _JHEP_ 05 (2021) 211, [2009.13452].
  * [13] S. Datta, A. Ghosal and R. Samanta, _Baryogenesis from ultralight primordial black holes and strong gravitational waves from cosmic strings_ , _JCAP_ 08 (2021) 021, [2012.14981].
  * [14] R. Samanta and F. R. Urban, _Testing super heavy dark matter from primordial black holes with gravitational waves_ , _JCAP_ 06 (2022) 017, [2112.04836].
  * [15] D. Borah, S. Jyoti Das, R. Samanta and F. R. Urban, _PBH-infused seesaw origin of matter and unique gravitational waves_ , _JHEP_ 03 (2023) 127, [2211.15726].
  * [16] J. Ellis, M. Lewicki, C. Lin and V. Vaskonen, _Cosmic Superstrings Revisited in Light of NANOGrav 15-Year Data_ , 2306.17147.
  * [17] Z. Wang, L. Lei, H. Jiao, L. Feng and Y.-Z. Fan, _The nanohertz stochastic gravitational-wave background from cosmic string Loops and the abundant high redshift massive galaxies_ , 2306.17150.
  * [18] N. Kitajima, J. Lee, K. Murai, F. Takahashi and W. Yin, _Nanohertz Gravitational Waves from Axion Domain Walls Coupled to QCD_ , 2306.17146.
  * [19] G. Franciolini, A. Iovino, Junior., V. Vaskonen and H. Veermae, _The recent gravitational wave observation by pulsar timing arrays and primordial black holes: the importance of non-gaussianities_ , 2306.17149.
  * [20] E. Megias, G. Nardini and M. Quiros, _Pulsar Timing Array Stochastic Background from light Kaluza-Klein resonances_ , 2306.17071.
  * [21] K. Fujikura, S. Girmohanta, Y. Nakai and M. Suzuki, _NANOGrav Signal from a Dark Conformal Phase Transition_ , 2306.17086.
  * [22] C. Han, K.-P. Xie, J. M. Yang and M. Zhang, _Self-interacting dark matter implied by nano-Hertz gravitational waves_ , 2306.16966.
  * [23] L. Zu, C. Zhang, Y.-Y. Li, Y.-C. Gu, Y.-L. S. Tsai and Y.-Z. Fan, _Mirror QCD phase transition as the origin of the nanohertz Stochastic Gravitational-Wave Background_ , 2306.16769.
  * [24] J. Yang, N. Xie and F. P. Huang, _Implication of nano-Hertz stochastic gravitational wave background on ultralight axion particles_ , 2306.17113.
  * [25] S.-Y. Guo, M. Khlopov, X. Liu, L. Wu, Y. Wu and B. Zhu, _Footprints of Axion-Like Particle in Pulsar Timing Array Data and JWST Observations_ , 2306.17022.
  * [26] Z.-Q. Shen, G.-W. Yuan, Y.-Y. Wang and Y.-Z. Wang, _Dark Matter Spike surrounding Supermassive Black Holes Binary and the nanohertz Stochastic Gravitational Wave Background_ , 2306.17143.
  * [27] G. Franciolini, D. Racco and F. Rompineve, _Footprints of the QCD Crossover on Cosmological Gravitational Waves at Pulsar Timing Arrays_ , 2306.17136.
  * [28] G. Lambiase, L. Mastrototaro and L. Visinelli, _Astrophysical neutrino oscillations after pulsar timing array analyses_ , 2306.16977.
  * [29] Y.-Y. Li, C. Zhang, Z. Wang, M.-Y. Cui, Y.-L. S. Tsai, Q. Yuan et al., _Primordial magnetic field as a common solution of nanohertz gravitational waves and Hubble tension_ , 2306.17124.
  * [30] P. Athron, A. Fowlie, C.-T. Lu, L. Morris, L. Wu, Y. Wu et al., _Can supercooled phase transitions explain the gravitational wave background observed by pulsar timing arrays?_ , 2306.17239.
  * [31] E. Babichev, D. Gorbunov, S. Ramazanov, R. Samanta and A. Vikman, _NANOGrav spectral index $\gamma=3$ from melting domain walls_, 2307.04582.
  * [32] G. Lazarides, R. Maji and Q. Shafi, _Superheavy quasi-stable strings and walls bounded by strings in the light of NANOGrav 15 year data_ , 2306.17788.
  * [33] R. Maji and W.-I. Park, _Supersymmetric $U(1)_{B-L}$ flat direction and NANOGrav 15 year data_, 2308.11439.
  * [34] H.-L. Huang, Y. Cai, J.-Q. Jiang, J. Zhang and Y.-S. Piao, _Supermassive primordial black holes in multiverse: for nano-Hertz gravitational wave and high-redshift JWST galaxies_ , 2306.17577.
  * [35] J.-Q. Jiang, Y. Cai, G. Ye and Y.-S. Piao, _Broken blue-tilted inflationary gravitational waves: a joint analysis of NANOGrav 15-year and BICEP/Keck 2018 data_ , 2307.15547.
  * [36] M. R. Gangopadhyay, V. V. Godithi, K. Ichiki, R. Inui, T. Kajino, A. Manusankar et al., _Is NanoGRAV signals pointing towards resonant particle creation during inflation?_ , 2309.03101.
  * [37] A. Nikandish, S. R. Zadeh, R. Naderi, F. Elahi and H. Mehrabpour, _Signature of Non-Minimal Scalar-Gravity Coupling with an Early Matter Domination on the Power Spectrum of Gravitational Waves_ , 2309.10554.
  * [38] S. Kawai and J. Kim, _Probing inflationary moduli space with gravitational waves_ , 2308.13272.
  * [39] L. Frosina and A. Urbano, _On the inflationary interpretation of the nHz gravitational-wave background_ , 2308.06915.
  * [40] Z. Yi, Z.-Q. You and Y. Wu, _Model-independent reconstruction of the primordial curvature power spectrum from PTA data_ , 2308.05632.
  * [41] N. Bhaumik, R. K. Jain and M. Lewicki, _Ultra-low mass PBHs in the early universe can explain the PTA signal_ , 2308.07912.
  * [42] S. A. Hosseini Mansoori, F. Felegray, A. Talebian and M. Sami, _PBHs and GWs from $\mathbb{T}$2-inflation and NANOGrav 15-year data_, _JCAP_ 08 (2023) 067, [2307.06757].
  * [43] M. Yamada and K. Yonekura, _Dark baryon from pure Yang-Mills theory and its GW signature from cosmic strings_ , 2307.06586.
  * [44] G. Ye and A. Silvestri, _Can the gravitational wave background feel wiggles in spacetime?_ , 2307.05455.
  * [45] E. Madge, E. Morgante, C. Puchades-Ibáñez, N. Ramberg, W. Ratzinger, S. Schenk et al., _Primordial gravitational waves in the nano-Hertz regime and PTA data – towards solving the GW inverse problem_ , 2306.14856.
  * [46] T. Bringmann, P. F. Depta, T. Konstandin, K. Schmidt-Hoberg and C. Tasillo, _Does NANOGrav observe a dark sector phase transition?_ , 2306.09411.
  * [47] P. F. Depta, K. Schmidt-Hoberg, P. Schwaller and C. Tasillo, _Do pulsar timing arrays observe merging primordial black holes?_ , 2306.17836.
  * [48] T. Broadhurst, C. Chen, T. Liu and K.-F. Zheng, _Binary Supermassive Black Holes Orbiting Dark Matter Solitons: From the Dual AGN in UGC4211 to NanoHertz Gravitational Waves_ , 2306.17821.
  * [49] S. Ge, _Stochastic gravitational wave background: birth from axionic string-wall death_ , 2307.08185.
  * [50] V. K. Oikonomou, _Flat energy spectrum of primordial gravitational waves versus peaks and the NANOGrav 2023 observation_ , _Phys. Rev. D_ 108 (2023) 043516, [2306.17351].
  * [51] V. K. Oikonomou, _A Stiff Pre-CMB Era with a Mildly Blue-tilted Tensor Inflationary Era can Explain the 2023 NANOGrav Signal_ , 2309.04850.
  * [52] N. Kitajima and K. Nakayama, _Nanohertz gravitational waves from cosmic strings and dark photon dark matter_ , 2306.17390.
  * [53] A. Eichhorn, R. R. Lino dos Santos and J. a. L. Miqueleto, _From quantum gravity to gravitational waves through cosmic strings_ , 2306.17718.
  * [54] C. Unal, A. Papageorgiou and I. Obata, _Axion-Gauge Dynamics During Inflation as the Origin of Pulsar Timing Array Signals and Primordial Black Holes_ , 2307.02322.
  * [55] P. Di Bari and M. H. Rahat, _The split majoron model confronts the NANOGrav signal_ , 2307.03184.
  * [56] X. K. Du, M. X. Huang, F. Wang and Y. K. Zhang, _Did the nHZ Gravitational Waves Signatures Observed By NANOGrav Indicate Multiple Sector SUSY Breaking?_ , 2307.02938.
  * [57] S. Antusch, K. Hinze, S. Saad and J. Steiner, _Singling out SO(10) GUT models using recent PTA results_ , 2307.04595.
  * [58] Z. Zhang, C. Cai, Y.-H. Su, S. Wang, Z.-H. Yu and H.-H. Zhang, _Nano-Hertz gravitational waves from collapsing domain walls associated with freeze-in dark matter in light of pulsar timing array observations_ , 2307.11495.
  * [59] H. An, B. Su, H. Tai, L.-T. Wang and C. Yang, _Phase transition during inflation and the gravitational wave signal at pulsar timing arrays_ , 2308.00070.
  * [60] S. He, L. Li, S. Wang and S.-J. Wang, _Constraints on holographic QCD phase transitions from PTA observations_ , 2308.07257.
  * [61] J. Ellis, M. Fairbairn, G. Franciolini, G. Hütsi, A. Iovino, M. Lewicki et al., _What is the source of the PTA GW signal?_ , 2308.08546.
  * [62] M. Kawasaki and K. Murai, _Enhancement of gravitational waves at Q-ball decay including non-linear density perturbations_ , 2308.13134.
  * [63] M. X. Huang, F. Wang and Y. K. Zhang, _The Interplay Between the Muon $g-2$ Anomaly and the PTA nHZ Gravitational Waves from Domain Walls in NMSSM_, 2309.06378.
  * [64] S. Vagnozzi, _Implications of the NANOGrav results for inflation_ , _Mon. Not. Roy. Astron. Soc._ 502 (2021) L11–L15, [2009.13432].
  * [65] S. Bhattacharya, S. Mohanty and P. Parashari, _Implications of the NANOGrav result on primordial gravitational waves in nonstandard cosmologies_ , _Phys. Rev. D_ 103 (2021) 063532, [2010.05071].
  * [66] S. Kuroyanagi, T. Takahashi and S. Yokoyama, _Blue-tilted inflationary tensor spectrum and reheating in the light of NANOGrav results_ , _JCAP_ 01 (2021) 071, [2011.03323].
  * [67] M. Benetti, L. L. Graef and S. Vagnozzi, _Primordial gravitational waves from NANOGrav: A broken power-law approach_ , _Phys. Rev. D_ 105 (2022) 043520, [2111.04758].
  * [68] S. Vagnozzi, _Inflationary interpretation of the stochastic gravitational wave background signal detected by pulsar timing array experiments_ , _JHEAp_ 39 (2023) 81–98, [2306.16912].
  * [69] A. Gruzinov, _Elastic inflation_ , _Phys. Rev. D_ 70 (2004) 063518, [astro-ph/0404548].
  * [70] T. Kobayashi, M. Yamaguchi and J. Yokoyama, _G-inflation: Inflation driven by the Galileon field_ , _Phys. Rev. Lett._ 105 (2010) 231302, [1008.0603].
  * [71] S. Endlich, A. Nicolis and J. Wang, _Solid Inflation_ , _JCAP_ 10 (2013) 011, [1210.0569].
  * [72] D. Cannone, G. Tasinato and D. Wands, _Generalised tensor fluctuations and inflation_ , _JCAP_ 01 (2015) 029, [1409.6568].
  * [73] A. Ricciardone and G. Tasinato, _Primordial gravitational waves in supersolid inflation_ , _Phys. Rev. D_ 96 (2017) 023508, [1611.04516].
  * [74] Y.-F. Cai, J.-O. Gong, S. Pi, E. N. Saridakis and S.-Y. Wu, _On the possibility of blue tensor spectrum within single field inflation_ , _Nucl. Phys. B_ 900 (2015) 517–532, [1412.7241].
  * [75] T. Fujita, S. Kuroyanagi, S. Mizuno and S. Mukohyama, _Blue-tilted Primordial Gravitational Waves from Massive Gravity_ , _Phys. Lett. B_ 789 (2019) 215–219, [1808.02381].
  * [76] Y. Mishima and T. Kobayashi, _Revisiting slow-roll dynamics and the tensor tilt in general single-field inflation_ , _Phys. Rev. D_ 101 (2020) 043536, [1911.02143].
  * [77] R. H. Cyburt, B. D. Fields, K. A. Olive and T.-H. Yeh, _Big Bang Nucleosynthesis: 2015_ , _Rev. Mod. Phys._ 88 (2016) 015004, [1505.01076].
  * [78] KAGRA, Virgo, LIGO Scientific collaboration, R. Abbott et al., _Upper limits on the isotropic gravitational-wave background from Advanced LIGO and Advanced Virgo’s third observing run_ , _Phys. Rev. D_ 104 (2021) 022004, [2101.12130].
  * [79] A. Peimbert, M. Peimbert and V. Luridiana, _The primordial helium abundance and the number of neutrino families_ , _Rev. Mex. Astron. Astrofis._ 52 (2016) 419–424, [1608.02062].
  * [80] LIGO Scientific, Virgo collaboration, B. P. Abbott et al., _Upper Limits on the Stochastic Gravitational-Wave Background from Advanced LIGO’s First Observing Run_ , _Phys. Rev. Lett._ 118 (2017) 121101, [1612.02029].
  * [81] S. Datta and R. Samanta, _Gravitational waves-tomography of Low-Scale-Leptogenesis_ , _JHEP_ 11 (2022) 159, [2208.09949].
  * [82] S. Datta, _Inflationary gravitational waves, pulsar timing data and low-scale-leptogenesis_ , 2307.00646.
  * [83] D. Borah, S. Jyoti Das and R. Roshan, _Probing high scale seesaw and PBH generated dark matter via gravitational waves with multiple tilts_ , 2208.04965.
  * [84] D. Borah, S. Jyoti Das, R. Roshan and R. Samanta, _Imprint of PBH domination on gravitational waves generated by cosmic strings_ , _Phys. Rev. D_ 108 (2023) 023531, [2304.11844].
  * [85] B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, _Constraints on primordial black holes_ , _Rept. Prog. Phys._ 84 (2021) 116902, [2002.12778].
  * [86] K. Inomata, K. Kohri, T. Nakama and T. Terada, _Enhancement of Gravitational Waves Induced by Scalar Perturbations due to a Sudden Transition from an Early Matter Era to the Radiation Era_ , _Phys. Rev. D_ 100 (2019) 043532, [1904.12879].
  * [87] T. Papanikolaou, V. Vennin and D. Langlois, _Gravitational waves from a universe filled with primordial black holes_ , _JCAP_ 03 (2021) 053, [2010.11573].
  * [88] G. Domènech, C. Lin and M. Sasaki, _Gravitational wave constraints on the primordial black hole dominated early universe_ , _JCAP_ 04 (2021) 062, [2012.08151].
  * [89] G. Domènech, V. Takhistov and M. Sasaki, _Exploring evaporating primordial black holes with gravitational waves_ , _Phys. Lett. B_ 823 (2021) 136722, [2105.06816].
  * [90] G. Domènech, _Scalar Induced Gravitational Waves Review_ , _Universe_ 7 (2021) 398, [2109.01398].
  * [91] A. H. Guth, _The Inflationary Universe: A Possible Solution to the Horizon and Flatness Problems_ , _Phys. Rev. D_ 23 (1981) 347–356.
  * [92] A. D. Linde, _A New Inflationary Universe Scenario: A Possible Solution of the Horizon, Flatness, Homogeneity, Isotropy and Primordial Monopole Problems_ , _Phys. Lett. B_ 108 (1982) 389–393.
  * [93] WMAP collaboration, L. Page et al., _Three year Wilkinson Microwave Anisotropy Probe (WMAP) observations: polarization analysis_ , _Astrophys. J. Suppl._ 170 (2007) 335, [astro-ph/0603450].
  * [94] BICEP2, Keck Array collaboration, P. A. R. Ade et al., _BICEP2 / Keck Array x: Constraints on Primordial Gravitational Waves using Planck, WMAP, and New BICEP2/Keck Observations through the 2015 Season_ , _Phys. Rev. Lett._ 121 (2018) 221301, [1810.05216].
  * [95] A. R. Liddle and D. H. Lyth, _The Cold dark matter density perturbation_ , _Phys. Rept._ 231 (1993) 1–105, [astro-ph/9303019].
  * [96] N. Seto and J. Yokoyama, _Probing the equation of state of the early universe with a space laser interferometer_ , _J. Phys. Soc. Jap._ 72 (2003) 3082–3086, [gr-qc/0305096].
  * [97] L. A. Boyle and P. J. Steinhardt, _Probing the early universe with inflationary gravitational waves_ , _Phys. Rev. D_ 77 (2008) 063504, [astro-ph/0512014].
  * [98] K. Nakayama, S. Saito, Y. Suwa and J. Yokoyama, _Probing reheating temperature of the universe with gravitational wave background_ , _JCAP_ 06 (2008) 020, [0804.1827].
  * [99] S. Kuroyanagi, T. Chiba and N. Sugiyama, _Precision calculations of the gravitational wave background spectrum from inflation_ , _Phys. Rev. D_ 79 (2009) 103501, [0804.3249].
  * [100] K. Nakayama and J. Yokoyama, _Gravitational Wave Background and Non-Gaussianity as a Probe of the Curvaton Scenario_ , _JCAP_ 01 (2010) 010, [0910.0715].
  * [101] S. Kuroyanagi, T. Takahashi and S. Yokoyama, _Blue-tilted Tensor Spectrum and Thermal History of the Universe_ , _JCAP_ 02 (2015) 003, [1407.4785].
  * [102] Y. Watanabe and E. Komatsu, _Improved Calculation of the Primordial Gravitational Wave Spectrum in the Standard Model_ , _Phys. Rev. D_ 73 (2006) 123515, [astro-ph/0604176].
  * [103] K. Saikawa and S. Shirai, _Primordial gravitational waves, precisely: The role of thermodynamics in the Standard Model_ , _JCAP_ 05 (2018) 035, [1803.01038].
  * [104] I. Masina, _Dark matter and dark radiation from evaporating primordial black holes_ , _Eur. Phys. J. Plus_ 135 (2020) 552, [2004.04740].
  * [105] Planck collaboration, N. Aghanim et al., _Planck 2018 results. VI. Cosmological parameters_ , _Astron. Astrophys._ 641 (2020) A6, [1807.06209].
  * [106] R. Saito and J. Yokoyama, _Gravitational wave background as a probe of the primordial black hole abundance_ , _Phys. Rev. Lett._ 102 (2009) 161101, [0812.4339].
  * [107] R.-g. Cai, S. Pi and M. Sasaki, _Gravitational Waves Induced by non-Gaussian Scalar Perturbations_ , _Phys. Rev. Lett._ 122 (2019) 201101, [1810.11000].
  * [108] K. Inomata, M. Kawasaki, K. Mukaida, T. Terada and T. T. Yanagida, _Gravitational Wave Production right after a Primordial Black Hole Evaporation_ , _Phys. Rev. D_ 101 (2020) 123533, [2003.10455].
  * [109] R. Anantua, R. Easther and J. T. Giblin, _GUT-Scale Primordial Black Holes: Consequences and Constraints_ , _Phys. Rev. Lett._ 103 (2009) 111303, [0812.0825].
  * [110] D. Hooper, G. Krnjaic, J. March-Russell, S. D. McDermott and R. Petrossian-Byrne, _Hot Gravitons and Gravitational Waves From Kerr Black Holes in the Early Universe_ , 2004.00618.
  * [111] T. Fujita, M. Kawasaki, K. Harigaya and R. Matsuda, _Baryon asymmetry, dark matter, and density perturbation from primordial black holes_ , _Phys. Rev. D_ 89 (2014) 103501, [1401.1909].
  * [112] O. Lennon, J. March-Russell, R. Petrossian-Byrne and H. Tillim, _Black Hole Genesis of Dark Matter_ , _JCAP_ 04 (2018) 009, [1712.07664].
  * [113] G. Domènech and M. Sasaki, _Gravitational wave hints black hole remnants as dark matter_ , _Class. Quant. Grav._ 40 (2023) 177001, [2303.07661].
  * [114] G. Franciolini and P. Pani, _Stochastic gravitational-wave background at 3G detectors as a smoking gun for microscopic dark matter relics_ , 2304.13576.
  * [115] 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].
  * [116] T. Asaka and H. Okui, _Neutrino masses and gravitational wave background_ , _Phys. Lett. B_ 814 (2021) 136074, [2012.13527].
  * [117] D. Borah, S. Jyoti Das and R. Samanta, _Inflationary origin of gravitational waves with \textitMiracle-less WIMP dark matter in the light of recent PTA results_ , 2307.00537.