#  Exploring the Spectral Shape of Gravitational Waves Induced by Primordial
Scalar Perturbations and Connection with the Primordial Black Hole Scenarios

Ioannis Dalianis<EMAIL_ADDRESS>Physics Division, National Technical
University of Athens, 15780 Zografou campus, Athens, Greece  Konstantinos
Kritos<EMAIL_ADDRESS>Physics Division, National Technical
University of Athens, 15780 Zografou campus, Athens, Greece

###### Abstract

There is a growing expectation that the gravitational wave detectors will
start probing the stochastic gravitational wave backgrounds in the following
years. We explore the spectral shapes of gravitational waves induced to second
order by scalar perturbations and presumably have been produced in the early
universe. We calculate the gravitational wave spectra generated during
radiation and kination eras together with the associated primordial black hole
counterpart. We employ power spectra for the primordial curvature perturbation
generated by $\alpha$-attractors and nonminimal derivative coupling inflation
models as well as Gaussian and delta-type shapes. We demonstrate the ability
of the tensor modes to constrain the spectrum of the primordial curvature
perturbations and discriminate among inflationary models. Gravitational wave
production during kination and radiation eras can also be distinguished by
their spectral shapes and amplitudes.

##  I Introduction

A network of operating and designed gravitational wave detectors raise
expectations that, sooner or later, the relic gravitational radiation
background will be detected. Gravitational waves (GWs) are thought to be
generated by several processes in the early universe Kuroyanagi  _et al._
(2018); Christensen (2019); Caprini and Figueroa (2018). Among them there is a
model independent source, the primordial density scalar perturbations, that
unambiguously source tensor modes at second order perturbation theory
Matarrese  _et al._ (1993, 1998); Noh and Hwang (2004); Carbone and Matarrese
(2005); Nakamura (2007), the so-called induced gravitational waves (IGW).

Tensor modes are induced to second order by scalar perturbations despite their
decoupling at first order. The cosmological implications of the scalar induced
GWs to second order were first discussed in the pioneering works Mollerach
_et al._ (2004), where the effects of IGWs on cosmic microwave background
(CMB) polarization was computed, and Ananda  _et al._ (2007), where the
spectrum of the IGWs on small scales was studied. Further notable results
followed that extended our understanding and formulation of IGW physics
Assadullahi and Wands (2009); Baumann  _et al._ (2007); Saito and Yokoyama
(2009) and gave us insights of how to probe the primordial power spectrum,
${\cal P_{R}}(k)$, at small scales.

Primordial scalar density perturbations are directly measured only at the very
large scales of the CMB with wave number $k_{\text{cmb}}\sim 10^{-2}$ Mpc-1
and with amplitude $\delta\rho/\rho\sim 10^{-5}$. At smaller scales $k\gg
k_{\text{cmb}}$ there are weak constraints coming mainly from bounds on the
abundance of primordial black holes (PBHs) Dalianis (2019); Kalaja  _et al._
(2019); Carr  _et al._ (2020). Additional constraints on the scalar density
amplitude can come from gravitational wave experiments. Already, at the
frequency range $f\sim 10^{-9}$ Hz and $f\sim 10^{2}$ Hz the PTA and
LIGO/VIRGO experiments, respectively, have provided upper bounds.

The weak bounds on the primordial ${\cal P_{R}}(k)$ has prompted the building
of several scenarios that exhibit an enhanced power at small scales. A
subclass of these scenarios predicts a significant PBH abundance, see Ref.
Sasaki  _et al._ (2018) for a review. The new observational window to the very
early universe through the IGWs can test these scenarios Bugaev and Klimai
(2010a); Saito and Yokoyama (2010); Bugaev and Klimai (2010b, 2011); Nakama
_et al._ (2017); Orlofsky  _et al._ (2017); Garcia-Bellido  _et al._ (2017);
Cai  _et al._ (2019a); Inomata and Nakama (2019); Espinosa  _et al._ (2018);
Bartolo  _et al._ (2019); Byrnes  _et al._ (2019); Bartolo  _et al._ (2019);
Cai  _et al._ (2019b); Wang  _et al._ (2019); Gong (2019); Mahbub (2020);
Özsoy and Tasinato (2020); Bartolo  _et al._ (2020); Braglia  _et al._ (2020);
Ballesteros  _et al._ (2020). The interpretation of a GW signal is
nevertheless a nontrivial task and one needs to take into account several
parameters, such as the reheating temperature and the equation of state of the
early universe, as well as whether the primordial fluctuations are Gaussian or
a non-negligible non-Gaussianity exists Nakama  _et al._ (2017); Cai  _et al._
(2019a); Bartolo  _et al._ (2019); Ota (2020). From a different point of view,
detection of the relic GWs stochastic background is a direct probe of the very
early cosmic history, which is unknown for $t\lesssim 1$s Allahverdi  _et al._
(2020). In this context, the recent works Inomata  _et al._ (2020a, b);
Domènech (2019); Domènech  _et al._ (2020) have studied the impact of
different early cosmological evolution on the IGWs.

In this work we explore the shape of the IGWs produced in the early universe
assuming different scenarios and models. Prior to BBN there are no direct
tests of the thermodynamical state of the universe and we choose to presume
either the standard radiation equation of state or that of the stiff fluid, as
benchmark cases. We make a step forward compared to previous works taking into
account the transition details from kination to radiation. We focus on a
${\cal P_{R}}(k)$ shape with a peak. ${\cal P_{R}}(k)$ with a peak is a rather
interesting scenario because it predicts IGWs and might also generate PBHs
that can be part of the dark matter in the galaxies. We target at cosmological
parameter values that can be tested by the existing and designed GW
experiments, such as LIGO and Laser Interferometer Space Antenna (LISA). We
utilize Dirac $\delta$ and Gaussian-type ${\cal P_{R}}(k)$ shapes and give a
simple analytic expression for the energy density parameter of the IGWs,
$\Omega_{\text{IGW}}(t_{0})$, produced during radiation and kination regime.

Furthermore, we derive the scalar-induced GW spectral shapes for particular
models of inflation: two $\alpha$-attractors inflation models Kallosh  _et
al._ (2013) and a generalized Galileon-type Horndeski (1974); Kobayashi
(2019); Deffayet  _et al._ (2011) model of inflation. We focus on particular
three models Dalianis  _et al._ (2019); Dalianis and Tringas (2019); Dalianis
_et al._ (2020) because each model yields a ${\cal P_{R}}(k)$ with a different
shape form: different width and scaling both in the IR and UV part of the
scalar spectrum. We compare features of the ${\cal P_{R}}(k)$ spectral shape
with the corresponding $\Omega_{\text{IGW}}(f)$, examining the amplitudes and
the scaling. This approach makes possible a demonstrative exploration of the
IGW spectra with different shapes. A similar study in this direction has been
performed in Refs. Cai  _et al._ (2019c); Pi and Sasaki (2020) where,
following analytic steps, general features of the IGW spectrum have been also
described. Our analysis provides observational tests for this sort of viable
inflation models as well.

In addition to the IGW spectra we examine the abundance and distribution of
the PBHs associated with the ${\cal P_{R}}(k)$. We examine scenarios that
predict a sizable amount of both PBHs and IGWs, and scenarios that predict
significant IGWs with the PBH counterpart being either negligibly small or
promptly evaporating. We present analytic and numerical results and we
manifest the connection between the PBH mass and the GW counterpart making use
of figures and tables.

The paper is organized as follows. In Sec. II, we review the method that is
usually followed in the literature to calculate the spectra of scalar-induced
GWs. In Sec. III we focus on particular inflationary models that predict
${\cal P_{R}}(k)$ with a peak as well as general Gaussian ${\cal P_{R}}(k)$
types, and we specialize into radiation and kination domination early universe
scenarios. In Sec. IV we turn into the PBH counterpart and calculate mass
distributions and abundances. Our results are presented in Sec. V, and
finally, in Sec. VI we draw our conclusions.

##  II Methodology for calculating Induced Gravitational Waves

In this section we review the formulas required to calculate the energy
density of the induced gravitational waves (IGWs) by first order scalar
perturbations and adjust the notation in our context. We follow Ref. Kohri and
Terada (2018) and for a further study in a relevant context we refer the
reader to Kohri and Terada (2018); Domènech (2019); Pi and Sasaki (2020).

We consider a flat Friedmann-Robertson-Walker model with first order scalar
perturbations $\Phi$, $\Psi$ and second order (induced) tensor perturbations
$h_{ij}$, and ignore vector perturbations and first order tensor
perturbations. We describe our perturbations in the Newtonian gauge. The line
element of the perturbed metric is written as,

$ds^{2}=a^{2}(\eta)\ \big{[}-(1+2\Phi)\ d\eta^{2}\\\ +\big{(}(1-2\Psi)\
\delta_{ij}+{1\over 2}h_{ij}\big{)}\ dx^{i}\ dx^{j}\big{]},$ (1)

where $a(\eta)$ is the scale factor and $\delta_{ij}$ is the Kronecker tensor.
For future reference we also denote the conformal Hubble function by
$\mathcal{H}(\eta)\equiv{a{{}^{\prime}}(\eta)\over a(\eta)}$. Here, and later,
prime (′) denotes conformal time differentiation. We will also ignore
anisotropies111For the effects of the anisotropic stress due to ${\cal
O}(10)\%$ difference between $\Phi$ and $\Psi$ see Ref. Baumann  _et al._
(2007). and set $\Phi=\Psi$. In the Newtonian gauge the scalar perturbations
play the role of the gravitational potential and the tensor modes describe
gravitational waves.

We define $x\equiv k\eta$ and express the evolution of the scalar
perturbations in terms of the scalar transfer function, $\Phi$, defined by
$\Phi_{\textbf{k}}(\eta)\equiv\Phi(x)\ \phi_{\textbf{k}}$ in momentum space,
where $\phi_{\textbf{k}}$ the Fourier mode of the primordial scalar
perturbations. Considering the field equations of general relativity for a
single perfect cosmological fluid with the standard density-pressure equation
of state (EOS) relation $p=w\rho$, the evolution equations for the
gravitational potential is obtained by

$\displaystyle\Phi{{}^{\prime\prime}}(x)+3(1+w)\
\mathcal{H}(\eta)\Phi{{}^{\prime}}(x)+w\ k^{2}\ \Phi(x)=0\,.$ (2)

We consider adiabatic perturbations (zero entropy perturbations) and initial
conditions $\Phi(0)=1$ and $\Phi{{}^{\prime}}(0)=0$ for the scalar transfer
function. The general solution222 Here, we neglect the contribution from the
Bessel function of the second kind, $Y_{\alpha}(\sqrt{w}x)$, to avoid the
singularity at $\eta=0$. of Eq. (2) is given in terms of the Bessel function
of the first kind, $J$, and the gamma function, $\Gamma$, Baumann  _et al._
(2007),

$\displaystyle\Phi(x)={2^{\alpha}\
\Gamma(\alpha+1)\over\left(\sqrt{w}x\right)^{\alpha}}\,J_{\alpha}(\sqrt{w}x),$
(3)

where $\alpha={5+3w\over 2(1+3w)}$.

The evolution of the induced tensor modes is given by the equation,

$h_{\textbf{k}}{{}^{\prime\prime}}(\eta)+2\mathcal{H}(\eta)h_{\textbf{k}}{{}^{\prime}}(\eta)+k^{2}h_{\textbf{k}}(\eta)=\mathcal{S}_{\textbf{k}}(\eta),$
(4)

where the source function ${\cal S}_{\textbf{k}}$ at the left-hand side plays
a critical role. It is a convolution of scalar perturbations at different wave
numbers given by

$\displaystyle\mathcal{S}_{\mathbf{k}}(\eta)=4\int\ {d^{3}q\over(2\pi)^{3/2}}\
e^{ij}(\textbf{k})\ q_{i}\ q_{j}\ \phi_{\textbf{k}}\
\phi_{\textbf{k}-\textbf{q}}$ $\displaystyle\times\
f(|\textbf{q}|/k,|\textbf{k}-\textbf{q}|/k,\eta,k).$ (5)

Here, $e^{ij}$ is the polarization tensor, and $f$ is an auxiliary function
defined by

$f(u,v,\eta,k)=2\ \Phi(ux)\Phi(vx)+{4\over
3(1+w)}\left[\Phi(ux)+{uk\over\mathcal{H}(\eta)}\
\Phi{{}^{\prime}}(ux)\right]\left[\Phi(vx)+{vk\over\mathcal{H}(\eta)}\
\Phi{{}^{\prime}}(vx)\right].$ (6)

Gaussian fluctuations are best described with the power spectral density. For
scalar and tensor perturbations the power spectral densities, called
$\mathcal{P}_{\Phi}$ and $\mathcal{P}_{h}$, respectively, are defined in terms
of two point correlation functions as

$\langle\phi_{\textbf{q}}\ \phi_{\textbf{k}}\rangle\equiv{2\pi\over k^{3}}\
\mathcal{P}_{\Phi}(k)\ \delta^{(3)}(\textbf{q}+\textbf{k}),$ (7)

and

$\langle h^{\oplus,\otimes}_{\textbf{q}}(\eta)\
h^{\oplus,\otimes}_{\textbf{k}}(\eta)\rangle\equiv{2\pi\over k^{3}}\
\mathcal{P}_{h}(\eta,k)\ \delta^{(3)}(\textbf{q}+\textbf{k}),$ (8)

where $\delta^{(3)}$ is the three-dimensional delta distribution and
$\oplus,\otimes$ denote the two polarization modes. The power spectral density
of the comoving curvature perturbation, $\mathcal{P}_{\mathcal{R}}$, is
related to $\mathcal{P}_{\Phi}$ via

$\mathcal{P}_{\mathcal{R}}(k)=\left({5+3w\over 3+3w}\right)^{2}\
\mathcal{P}_{\Phi}(k).$ (9)

The density parameter of IGWs is defined to be the ratio of the energy density
in IGWs over the total energy density, taking into account both $\oplus$ and
$\otimes$ GW polarizations. A useful definition for the energy density of IGWs
per unit logarithmic frequency interval is given in terms of the tensor power
spectrum as

$\Omega_{\textrm{IGW}}(\eta,k)\equiv{1\over 24}\
\left({k\over\mathcal{H}(\eta)}\right)^{2}\
\overline{\mathcal{P}_{h}(\eta,k)},$ (10)

where the overline denotes the oscillation average. The power spectral density
of the tensor perturbations is expressed as a double integral involving the
power spectrum of the curvature perturbations,

$\overline{\mathcal{P}_{h}(\eta,k)}=\int_{0}^{\infty}dv\int_{|1-v|}^{1+v}du\
\mathcal{T}(u,v,\eta,k)\ \mathcal{P}_{\mathcal{R}}(uk)\
\mathcal{P}_{\mathcal{R}}(vk).$ (11)

The $\mathcal{T}$ is the tensor transfer function given by

$\mathcal{T}(u,v,\eta,k)=4\ \left(4v^{2}-(1+v^{2}-u^{2})^{2}\over
4uv\right)^{2}\\\ \times\left({3+3w\over 5+3w}\right)^{2}\
\overline{I^{2}(u,v,\eta,k)},$ (12)

where we have defined the kernel function by

$I(u,v,\eta,k)=\int_{0}^{x}dy\ {a(y/k)\over a(\eta)}\ kG_{k}(x,y)\
f(u,v,y/k,k),$ (13)

and have used the Green’s function $kG_{k}$ given in terms of Bessel
functions, as $kG_{k}(x,y)={\pi\over 2}\sqrt{x\ y}\cdot\left[Y_{\nu}(x)\
J_{\nu}(y)-Y_{\nu}(y)\ J_{\nu}(x)\right]$, with $\nu\equiv{3(1-w)\over
2(1+3w)}$, Tomikawa and Kobayashi (2020). The oscillation average is given by

$\overline{I^{2}(u,v,\eta,k)}\equiv{1\over 2\pi}\ \int_{x}^{x+2\pi}dy\
I^{2}(u,v,\eta,y/\eta).$ (14)

The IGWs are sourced gravitational waves. At the time $t_{\text{c}}$ the
production of IGWs ceases and afterward the GWs propagate freely. In a
radiation dominated background the energy density parameter of the GWs also
remains constant. The energy density parameter of the IGW today, $t_{0}$, is
therefore given by Eq. (10) at the moment $t_{\text{c}}$ times the current
radiation density parameter, $\Omega_{\gamma,0}h^{2}=4.2\times 10^{-5}$,
modulo changes in the number of the relativistic degrees of freedom $g_{*}$ in
the radiation fluid,

$\Omega_{\textrm{IGW}}(t_{0},f)h^{2}=0.39\times\left({g_{\ast}\over
106.75}\right)^{-{1\over 3}}\ \Omega_{\gamma,0}h^{2}\\\
\times\Omega_{\textrm{IGW}}(t_{\text{c}},f).$ (15)

For the total IGW density parameter, we integrate Eq. (15) over the
logarithmic interval of frequency,

$\Omega_{\textrm{IGW}}(t_{0})=\int_{f_{\text{min}}}^{f_{\text{max}}}d\ln f\
\Omega_{\textrm{IGW}}(t_{0},f).$ (16)

Since our spectra will be peaked at some frequency $f_{\text{IGW}}^{\text{p}}$
we may take $f_{\text{min}}=f_{\text{IGW}}^{\text{p}}/100$ and
$f_{\text{max}}=10f_{\text{IGW}}^{\text{p}}$. Here
$\Omega_{\text{IGW}}(t_{0})$ should not be confused with
$\Omega_{\text{IGW}}(t_{0},f)$ where the later function is frequency dependent
and denotes the density parameter of IGWs per logarithmic frequency interval.

In the following we will examine particular ${\cal P_{R}}(k)$ types and we
will focus on two benchmark early universe cosmological scenarios: the
radiation and kination dominated eras.

##  III Study of the IGWs produced by explicit models and in different
cosmological eras

It is well known that the equation of state of the Universe has not been
directly probed for times prior to BBN $t\sim 1$ s. IGWs have amplitude and
spectral shape that depend on the details of cosmological era at the time of
their production. Indeed, the $\Phi$ that sources the GWs evolves in a
different way for different equation of state $w$ and, additionally, the
energy density of the produced GWs scales differently with respect to the
background.

In order to study the IGW spectrum and demonstrate its properties and
behaviour we assume that the curvature power spectrum ${\cal P_{R}}(k)$
features a peak at the wave number $k_{\text{p}}$. If the amplitude of the
${\cal P_{R}}(k)$ peak is significant it produces a strong IGW signal, that is
possibly detectable, and additionally, enhances exponentially the production
probability for PBHs. Therefore, explicit models can be tested and upper
bounds on the IGWs can be derived.

We will study ${\cal P_{R}}(k)$ peaks generated by complete and explicit
inflationary models, namely the $\alpha$-attractors and Horndeski general
nonminimal derivative coupling (GNMDC) inflation introduced in Refs. Dalianis
_et al._ (2019) and Dalianis  _et al._ (2020) respectively. The qualitative
element of the first model is that features a near-inflection point Garcia-
Bellido and Ruiz Morales (2017); Germani and Prokopec (2017); Motohashi and Hu
(2017); Ballesteros and Taoso (2018) and the second model features a so-called
high friction regime Germani and Kehagias (2010). From the structural side the
two inflationary models Dalianis  _et al._ (2019, 2020) have a noncanonical
kinetic term,

$\displaystyle\frac{{\cal
L}^{\alpha-\text{attractors}}}{\sqrt{-g}}=\frac{R}{2}-\frac{(\partial_{\mu}\tilde{\varphi})^{2}}{2(1-\frac{\tilde{\varphi}^{2}}{6\alpha})^{2}}+V(\tilde{\varphi})$
(17) $\displaystyle\frac{{\cal
L}^{\text{GNMDC}}}{\sqrt{-g}}=\frac{R}{2}-\hat{f}(\varphi)G^{\mu\nu}\partial_{\mu}\varphi\partial_{\nu}\varphi+V(\varphi)$
(18)

We note that we examine the second order induced tensor perturbations, not the
first order tensor modes produced by the inflationary stage itself; for the
latter, see e.g., Kuroyanagi  _et al._ (2015); Bernal and Hajkarim (2019) for
relevant studies. What is of interest for our analysis is that, respectively,
the form of the $\alpha$-attractor potential and the form of the GNMDC
coupling can generate an enhanced ${\cal P_{R}}(k)$ spectrum at $k_{\text{p}}$
and with a variety of shape forms. The $\alpha$-attractors potentials that we
utilize generate ${\cal P_{R}}(k)$ slopes that increase like $k^{3}$ and
$k^{4}$ and decrease like $k^{-1}$ or $k^{-4}$. On the other hand the ${\cal
P_{R}}(k)$ slopes of the Galileon GNMDC model increase like $k^{2}$ and
$k^{4}$ and decrease like $k^{-4}$. For a comprehensive description and
details about the inflationary Lagrangians (17) and (18), as well as the
generated ${\cal P}_{R}(k)$ shapes we refer the reader to the original works
Dalianis  _et al._ (2019, 2020). A brief description can be found in Appendix
B.

Moreover, irrespective of the inflationary theory, we assume general shapes
for the ${\cal P_{R}}(k)$ spectra. We utilize Gaussian and Dirac $\delta$-type
distributions,

$\displaystyle{\cal P_{R}}^{\text{G}}(z)=A_{\cal
R}\,e^{-(z/\epsilon)^{2}}=\frac{A_{0}}{\epsilon\sqrt{\pi}}e^{-(z/\epsilon)^{2}}$
(19) $\displaystyle{\cal P_{R}}^{\text{D}}(z)=A_{0}\,\delta(z)$ (20)

where $z\equiv\ln(k/k_{\text{p}})$ and $A_{\cal R}={\cal P_{R}}(k_{p})$. In
the limit $\epsilon\rightarrow 0$ the Gaussian distribution approaches the
$\delta$-distribution. The advantage of these distributions is that analytic
or semianalytic results can be obtained and, in addition, realistic models can
be parametrized. Gaussian spectra have been also studied recently in Pi and
Sasaki (2020).

We are interested in three characteristics for the power spectra ${\cal
P_{R}}(k)$: i) the amplitude $A_{\cal R}$; ii) the wave number of the peak,
$k_{\text{p}}$, where ${\cal P_{R}}(k_{\text{p}})=A_{\cal R}$; and iii) the
width $\epsilon$. We choose a large amplitude $A_{\cal R}$ so that the IGW
signal is significant. We slightly modulate the $A_{\cal R}$ in order to
maximize or minimize to a negligible amount the PBH abundance. The wave number
$k_{\text{p}}$ is chosen so that either the PBHs constitute a significant
fraction of the dark matter in the galaxies, or the frequency of the IGWs lays
in the sensitivity range of the gravitational detectors.

Let us now examine two benchmark cosmological scenarios separately: the
radiation and the kination.

###  III.1 The radiation domination scenario

Let us assume that the ${\cal P_{R}}(k)$ peak enters the horizon during the
radiation era. This is the standard early universe cosmological scenario. The
IGWs produced during radiation era have been thoroughly studied, thus, we will
not repeat known results and technical details. We remind the reader that the
potential $\Phi(x)$ oscillates with a $x^{-2}$ decaying amplitude in
subhorizon scales. The formalism of the Sec. II is applied for $w=1/3$. We
quote the analytic expression that we use to obtain the IGWs results for the
models (17), (18) and (19). Relying on analytical methods, we obtain the
tensor transfer function in late times as was done in Kohri and Terada (2018);
the resulting expression is

$\displaystyle\lim_{x\to\infty}x\ \mathcal{T}_{\textrm{RD}}(u,v,x)=\
2\left({4v^{2}-(1+v^{2}-u^{2})^{2}\over
4uv}\right)^{2}\left({3(u^{2}+v^{2}-3)\over
4u^{3}v^{3}}\right)^{2}\left[\left(-4uv+(u^{2}+v^{2}-3)\ln\left|{3-(u+v)^{2}\over
3-(u-v)^{2}}\right|\right)^{2}\right.$
$\displaystyle\left.+\pi^{2}(u^{2}+v^{2}-3)^{2}\Theta\left(u+v\sqrt{3}\right)\right]\
\ \ \ \ $ (21)

where $\Theta$ is the unit step function.

We then perform a numerical integration to calculate the ${\cal
P}_{h}(\eta,k)$ as Eq. (11) dictates. In order for this to be done, an input
for the power spectrum shape ${\cal P_{R}}(k)$ is required. We assume two
Gaussian power spectra with medium and narrow widths, $\epsilon=1$ and
$\epsilon=0.1$, respectively, as well as power spectra produced by
$\alpha$-attractor and Horndeski-type inflationary theories; see Fig. 4. The
${\cal P_{R}}(k)$ produced by these inflationary theories is found after
solving the Mukhanov-Sasaki equation numerically.

###  III.2 The kination domination scenario

In this subsection, let us assume the cosmological scenario that the very
early universe has been dominated by a phase whose equation of state is
stiffer than radiation. Kination is a regime where the kinetic energy of a
scalar field is dominant against its potential energy. In the limiting case
$w=1$, where the sound velocity is equal to the speed of light, the energy
density of the stiff fluid scales as $\rho\sim a^{-6}$. Such a scenario is
natural in theories with runaway potentials, e.g. theories with moduli fields.

We examine the kination scenario because it has striking implications both for
the PBH formation and the IGW produced. Our numerical results complement
previous studies Domènech (2019) and we show that an early universe kination
domination (eKD) era shapes in a different way the spectrum of the tensor
perturbations. This fact renders the eKD era testable.

We comment that the attribute "early" to kination might sound redundant since
there is no kination era observed in the universe. Nevertheless we use it to
emphasize that the kination era assumed precedes the standard radiation era.

The redshift of the stiff fluid energy density is the fastest and any ambient
radiation will sooner or latter dominate the early universe. Hence, a kination
era ends naturally if there is an extra fluid with softer equation of state.
Besides this gradual transition, an eKD era can end suddenly. The later is
simpler to examine because the scalar transfer function can be computed
explicitly due to the discrete $w$ values and single fluid analysis.333 A
mechanism that can implement a sudden transition can be easily conceived: an
extra field direction that ends the kination stage via a waterfall transition
to a global minimum where the field decays and reheats the universe.

Let us consider an early kination domination scenario that transits suddenly
into the RD phase. The transition takes place at the reheating temperature
$T_{\text{rh}}$. The transfer function in the kination domination cannot be
written in closed analytical expressions. One can only rely on approximate
methods, since the Bessel functions of order one involved can only be written
as infinite series and not further simplified. Hence, we treat the scenario
eKD to RD numerically. Further, in order to simplify our calculations we
utilize a monochromatic power spectrum of scalar perturbations modeled by a
delta distribution peaked at $k_{\text{p}}$, given by Eq. (20). Such a choice
models a narrow and sharply peaked $\mathcal{P}_{\mathcal{R}}$. In the
conclusions, we will discuss a correspondence between $\delta$-type and
Gaussian distributions, and anticipate features for the IGWs produced by a
wider ${\cal P_{R}}(k)$.

In a sudden eKD$\rightarrow$RD scenario, the scale factor is a piecewise
defined function, that is continuous at the point of transition
$\eta=\eta_{\text{rh}}$,

$a(\eta)=\begin{cases}\displaystyle\sqrt{\eta\over\eta_{\text{rh}}},&\eta<\eta_{\text{rh}}\\\
\displaystyle{\eta+\eta_{\text{rh}}\over
2\eta_{\text{rh}}},&\eta\geq\eta_{\text{rh}}\end{cases}$ (22)

A similar expression can be obtained for the conformal Hubble parameter,
${\cal H}$. The scalar transfer function, Eq. (2), is given by

$\Phi(x)=\begin{cases}\displaystyle{2\over x}\ J_{1}(x),&x<x_{\text{rh}}\\\
\displaystyle{3\over x^{2}}\left[C_{1}\left({\sin(x/\sqrt{3})\over
x/\sqrt{3}}-\cos(x/\sqrt{3})\right)+C_{2}\left({\cos(x/\sqrt{3})\over
x/\sqrt{3}}+\sin(x/\sqrt{3})\right)\right],&x\geq x_{\text{rh}}\end{cases}$
(23)

After the transition, into the RD era, we write the solution to Eq. (2) as a
linear combination of $(x/\sqrt{3})^{3/2}\cdot J_{3/2}(x/\sqrt{3})$ and
$(x/\sqrt{3})^{3/2}\cdot Y_{3/2}(x/\sqrt{3})$ with two constant coefficients
$C_{1}$ and $C_{2}$. We determine the $C_{1}$ and $C_{2}$ by the continuity of
the potential and its derivative at the point of the transition, after
expressing $J_{3/2}$ and $Y_{3/2}$ in terms of spherical Bessel functions. The
plot of the absolute value of Eq. (23) is depicted in Fig. 1. The $\Phi(x)$
remains constant as long as $x\ll 1$. At the horizon crossing,
$\eta=\eta_{\text{entry}}$, the gravitational potential starts decaying as Eq.
(23) dictates: roughly as $x^{-3/2}$ during eKD and after the transition to
the RD phase, as $x^{-2}$.

Figure 1: The transfer function in the case of the eKD$\rightarrow$RD sudden
transition in the Newtonian gauge. We plot the absolute value of the
gravitational potential, Eq. (23), for $x_{\text{rh}}=10,\ 50,\ 100$. The
horizontal axis is normalized so that the transition occurs at 1.

The reheating conformal time $\eta_{\text{rh}}$ is inversely proportional to
the reheating temperature, $\eta_{\text{rh}}\propto T_{\text{rh}}^{-1}$. This
means that higher reheating temperature scenarios predict a faster reheating
transition into RD. During eKD the growth of IGWs goes proportional to
conformal time for a given wave mode $k$; i.e.,
$\Omega_{\textrm{IGW}}^{\textrm{(KD)}}\sim(k_{\text{p}}\eta)^{2}\overline{I^{2}_{\textrm{KD}}}\propto\eta\propto
a^{2}$. This growth stalls once the transition is completed and radiation
takes over. In the cases that $\eta_{\text{entry}}\ll\eta_{\text{rh}}$, we
neglect the effect from RD era, because by that time the gravitational
potential sourcing the tensor modes is negligibly small. Numerically, in those
cases we calculate the IGWs at the time of transition444After the transition,
a sharp peak appears at $k={2\over\sqrt{3}}k_{\text{p}}$ superimposed on the
spectral shape obtained at $\eta=\eta_{\text{rh}}$. This corresponds to a
logarithmic divergence manifest in the radiation domination as can be seen in
Eq. (III.1). We neglect such a resonance peak as it has no physical meaning,
but it is rather a mathematical singularity associated to the delta power
spectrum we used, Eq. (20). Therefore, sudden eKD to RD scenarios may show up
in the experiments with an enhanced IGW spectrum.

When there is a fast reheating, i.e., right after the horizon entry of the
perturbation the universe transits into the radiation regime, the growth,
$\Omega_{\text{IGW}}\propto a^{2}$, of IGWs ceases early. However, the
gravitational potential has not been suppressed significantly, and the
contribution from RD era is important. In such a case of fast reheating we
calculate the IGWs at a conformal time $\eta_{\text{c}}=50\
\eta_{\text{rh}}\gg\eta_{\text{rh}}$, thus taking into account the non-
negligible contribution from the RD era. In Fig. 2, we present the evolution
of the time-dependent piece of the IGW spectral shape with respect to the
conformal time for the transition eKD$\rightarrow$RD.

Figure 2: The behavior of the time-dependent piece of the IGW spectrum under a
sudden eKD$\rightarrow$RD transition. The transition occurs when
$\eta/\eta_{\text{rh}}=1$. Here we plot the combination
$\left(\eta/\eta_{\text{rh}}\right)|I|$ for $k_{\text{p}}/k=1$ and for four
choices of the ratio $k_{\text{rh}}/k_{\text{p}}$ corresponding to different
colors in this figure. The growth is proportional to $\sqrt{\eta}$ during eKD.
After the transition, this function always approaches an oscillatory state
with a constant amplitude in the late times.

The final formula for the density parameter of IGWs in the eKD$\rightarrow$RD
transition is given by

$\displaystyle\Omega_{\textrm{IGW}}^{\textrm{eKD}\rightarrow\textrm{RD}}(\eta_{\text{c}},k)={1\over
6}A_{\mathcal{R}}^{2}\
\left({k_{\text{p}}\over\mathcal{H}(\eta_{\text{c}})}\right)^{2}\
\left[1-\left(k\over 2k_{\text{p}}\right)^{2}\right]^{2}$
$\displaystyle\times\
\overline{I_{\textrm{eKD}\rightarrow\textrm{RD}}^{2}\left({k_{\text{p}}\over
k},{k_{\text{p}}\over k},\eta_{\text{c}},k\right)}\ \Theta\left(1-{k\over
2k_{\text{p}}}\right).\ \ \ \ \ \ $ (24)

Here the unit step function $\Theta$ is included by conservation of momentum,
so that tensor modes with $k>2k_{\text{p}}$ are cut off. We remind the reader
that we have assumed a delta-type power spectrum for the curvature
perturbation. Later, in our results in Sec. V and conclusions in Sec. VI, we
estimate the amplitude of the $\Omega_{\text{IGW}}$ for a broad curvature
power spectrum, of Gaussian type, at the peak frequency.

Finally, let us mention that the existence of an early kination phase implies
that the energy density of the produced IGW gets enhanced and might backreact
on the geometry. This backreaction must be limited in order that the energy
density of the IGWs during the BBN does not increase the expansion rate to a
disturbing level for the observed abundances of the relic nuclei. In order to
comply with this constraint we find out a new bound on the reheating
temperature that must be satisfied,

$\displaystyle{T_{\text{rh}}}\gtrsim 10^{7}\,{\text{GeV}}\,A_{\cal
R}^{3/2}\left(\frac{M_{\text{}}/\gamma}{10^{20}\text{g}}\right)^{-1/2}.$ (25)

This is a rough, conservative bound that connects horizon mass $M/\gamma$ and
the $A_{\cal R}$ with the $T_{\text{rh}}$, and guarantees the success of the
BBN predictions. Further discussion and derivation details can be found in
Appendix A.

##  IV Primordial black holes and the associated Gravitational Wave signal

PBHs form from the collapse of large-amplitude primordial inhomogeneities Carr
and Hawking (1974); Meszaros (1974); Carr (1975); Lindley (1980). An
inhomogeneity decouples from the background expansion if the power spectrum
${\cal P_{R}}(k)$ is enhanced at a scale $k^{-1}$, characteristic of the PBH
mass. During RD, typical values for the ${\cal P_{R}}(k)\sim{\cal O}(10^{-2})$
are required. If PBHs have mass $M<10^{15}$ g evaporates at time scales less
than the age of the universe, whereas PBHs with $M>10^{15}$ g survive till
today being dynamically cold component of the dark matter in galactic
structures.

The PBH dark matter scenario is constrained in a wide range for the mass
parameter $M$ by several observational experiments, labeled with extra
galacric gamma-ray background (EGB), V, GC, HSC, O, EM, LIGO/VIRGO, GW2,
supernova (SN), wide binaries (WB), E, x-ray binaries (XB), Planck (PA) in
Fig. 6. In the present universe, the abundance of light PBH is constrained
from the extra galactic gamma-ray background Page and Hawking (1976);
MacGibbon and Carr (1991); Carr and MacGibbon (1998); Barrau  _et al._ (2003);
Carr  _et al._ (2016a), as well as from positron constraints (V) Dasgupta  _et
al._ (2019). Black holes of mass above $10^{17}$g are subject to gravitational
lensing constraints Barnacka  _et al._ (2012); Tisserand  _et al._ (2007);
Niikura  _et al._ (2019), given by Subaru (HSC), Ogle (O), EROS (E), and MACHO
(M), microlensing of SN and others. The CMB anisotropies measured by PA
constrains the PBH with mass above $10^{33}$g Ricotti  _et al._ (2008); Carr
_et al._ (2016b); Clesse and García-Bellido (2017); Bird  _et al._ (2016);
Poulin  _et al._ (2017); Serpico  _et al._ (2020). At the large mass region
there are also constraints from accretion limits in X-ray and radio
observations Gaggero  _et al._ (2017) and XB Carr  _et al._ (2020). There are
also dynamical limits from disruption of WB, survival of star clusters in
Eridanus II. Second order tensor perturbations (GW2) generated by scalar
perturbations already constrain the abundance of PBH masses approximately in
the range $10^{30}-10^{33}$ g. Also, the expected GWs generated by PBH
binaries with mass ${\cal O}(1-10)M_{\odot}$ coalescing at the present epoch
is constrained by LIGO/Virgo Sasaki  _et al._ (2016); Eroshenko (2018); Raidal
_et al._ (2017) as well as compact binary systems with component masses in the
range $0.2-1M_{\odot}$ Abbott  _et al._ (2018a) shown in Fig. 6. For a recent
update on the PBH constraints see Ref. Carr  _et al._ (2020).

###  IV.1 PBH abundance for general EOS

The present relic energy density parameter of primordial black holes with mass
$M$ produced at the cosmic time $t$ is

$\Omega_{\text{PBH}}(M)={\Omega_{\text{m}}}\gamma\beta(M)\left(\frac{M_{\text{rh}}}{M/\gamma}\right)^{\frac{2w}{1+w}}\left(\frac{M_{\text{eq}}}{M_{\text{rh}}}\right)^{1/2}\tilde{g}(g_{*})$
(26)

where
$\tilde{g}(g_{*})=2^{1/4}(g_{*}(t)/g_{*}(t_{\text{rh}}))^{-s/4}(g_{*}(t_{\text{rh}})/g_{*}(t_{\text{eq}}))^{-1/4}$
and $g_{*}$ the thermalized degrees of freedom. The parameter $s$ is equal to
1 for $t>t_{\text{rh}}$ and 0 for $t<t_{\text{rh}}$. $\Omega_{\text{m}}$ is
the total matter density parameter today, $M_{\text{eq}}$ is the horizon mass
at the moment of matter radiation equality, and $M_{\text{rh}}$ the horizon
mass at the moment of reheating.

The $\beta(M)$ is the mass fraction of the universe with horizon mass
$M/\gamma$ that collapsed and formed PBHs; it can be interpreted as the black
hole formation probability. Assuming Gaussian statistics, for a spherically
symmetric region it is

$\beta_{\text{}}(M)=\int_{\delta_{c}}d\delta\frac{1}{\sqrt{2\pi\sigma^{2}(M)}}e^{-\frac{\delta^{2}}{2\sigma^{2}(M)}}\,,$
(27)

where $\sigma(M)$ is the variance of the density perturbations and $\delta$
the density contrast. The PBH abundance has an exponential sensitivity to the
variance of the perturbations $\sigma(M)$ and the threshold value
$\delta_{c}$, which is $w$ dependent. In this work, we assign values to
$\delta_{c}$ following the findings of Ref. Harada  _et al._ (2013).

The horizon mass increases with a different rate for different expansion
rates. Using the relation $f=k/(2\pi)$ we find the frequency-PBH mass
correspondence for general equation of state $w$ and reheating temperature
$T_{\text{rh}}$,

$\displaystyle f_{\text{hor}}(M,T_{\text{rh}},$ $\displaystyle w)\simeq
2.7\times
10^{2}\text{Hz}\left(\frac{T_{\text{rh}}}{10^{10}\text{GeV}}\right)^{\frac{1-3w}{3(1+w)}}$
$\displaystyle\left(\frac{M/\gamma}{10^{12}\text{g}}\right)^{-\frac{3w+1}{3(1+w)}}\left(\frac{g_{*}}{106.75}\right)^{-\frac{2w}{6(1+w)}}\,.$
(28)

In the above relation, 555For the case of kination domination the minor
correction of replacing $T_{\text{rh}}$ by $2^{1/4}T_{\text{rh}}$ should be
considered if there is equipartition between radiation and scalar field energy
density at the time of reheating. we have assumed a one-to-one correspondence
between $k$ (or$f$) and the PBH mass $M$. This is true for the approximation
of a monochromatic PBH mass spectrum, which is practically the case in many
models. We note that the PBH mass distribution peaks at a value $M(k)$ that is
in a slight offset with the value $M(k_{\text{p}})$, as our following
numerical analysis shows. In the $k$-space the peak of the $\sigma(k)$ is
about at $0.7k_{\text{p}}$ Wang  _et al._ (2019). The $f(M,T_{\text{rh}},w)$
relation is depicted in Fig. 3. The $f_{\text{hor}}(M)$ frequency and
$f_{\text{p,IGW}}$ frequency, where the IGW spectrum maximizes, have similar
size but do not coincide.

Therefore, searching for GWs with frequency $f_{\text{IGW}}$, one can probe
PBH scenarios and the reheating temperature of the universe.

###  IV.2 GW detectors and PBH scenarios

Currently, there is a network of operating ground-based GW detectors that
focus on the high and low frequency regime, roughly at $10^{2}$ and $10^{-9}$
Hz, of the GW spectrum. Additionally, there are several designed and proposed
experiments sensitive enough to detect or constrain the stochastic GW
background at a large range of frequencies. In this subsection, we outline the
GW experiments labeled in Fig. 5 that can test the predictions of the models
described in this work.

The operating ground-based GW detector system is the advanced Laser
Interferometer Gravitational wave Observatory (aLIGO). In Fig. 5 we show for
reference O1 and O5 which correspond to the first observing run and the design
sensitivity respectively Abbott  _et al._ (2018b). The future ground-based
laser third generation interferometers is the Einstein Telescope (ET)
Sathyaprakash  _et al._ (2012), that probes the range of frequencies near 100
Hz. The space-based proposed experiments are DECIGO Seto  _et al._ (2001);
Sato  _et al._ (2017) and BBO Crowder and Cornish (2005), that will probe in
the decihertz frequency window. In the millihertz band of the spectrum there
is the scheduled LISA, Amaro-Seoane  _et al._ (2017), and proposed space-based
experiments such as the TianQin, Luo  _et al._ (2016). Pulsar Timing Arrays,
with the international PTA (IPTA) Hobbs  _et al._ (2010), and the planned
Square Kilometer Array (SKA) Janssen  _et al._ (2015) are projects that probe
and already pose constrains on the nanohertz regime.

Fitting functions regarding the sensitivity curves for LISA, aLIGO and ET are
provided by Sathyaprakash and Schutz (2009); for DECIGO and BBO, we consulted
Yagi and Seto (2011) and adjusted the sensitivities of these experiments to
current values and expectations; for TianQin, we used the fitting function
provided in Ref. Huang  _et al._ (2020).

IGWs can constrain a large window of the PBHs masses. Currently, there are
indirect constraints from the pulsar timing array experiments on IGWs
associated with the formation of relatively massive PBHs at the epoch of
horizon entry. Notably, a very severe constraint exists, $\beta(M)\lesssim
10^{-52}$, in the solar mass range, from pulsar timing data Saito and Yokoyama
(2009); Chen  _et al._ (2020). We note that Ref. Cai  _et al._ (2019d) pointed
out that the tension with PTA constraints can be relieved if the perturbations
are locally non-Gaussian.

Figure 3: The figure depicts the relation between the PBH mass
$M_{\text{PBH}}$ and the corresponding frequency of horizon with size
$k^{-1}=(2\pi f_{\text{hor}})^{-1}$ and mass $M_{\text{PBH}}/\gamma$ for two
equation of states $w=1/3$ (orange line) and $w=1$ (green lines); see Eq.
(IV.1). The reheating temperature is $T_{\text{rh}}=10^{4}$ GeV ($5\times
10^{2}$ GeV) for the dashed (dot-dashed) lines. The gray strip highlights the
frequency band of LISA optimal sensitivity.

##  V Results

In this section we outline our findings. We assume ${\cal P_{R}}(k)$
amplitudes and wave numbers having as a ruler the PBH abundance, and then we
calculate IGWs aiming at constraining early universe cosmological scenarios
and inflationary models.

We consider ${\cal P_{R}}(k)$ peaks generated by the $\alpha$-attractors
Dalianis  _et al._ (2019) and nonminimal derivative coupling Horndeski
inflation Dalianis  _et al._ (2020), as well as, irrespectively of the
inflationary theory, Gaussian and Dirac $\delta$-type distributions. The peak
choice of the $\mathcal{P}_{\mathcal{R}}$ is spanning a wide range of
frequencies in the GW spectrum; see Fig. 4. We compare the resulting IGW shape
with the sensitivity of current, planned and proposed GW experiments. We
choose amplitudes and shapes for the ${\cal P_{R}}(k)$ that induce GWs and PBH
abundances of cosmologically significant amount. We always make sure that our
input parameters respect the current constraints for the PBH abundance and the
GWs. We also take care the position and the width of the ${\cal P_{R}}(k)$
peak to respect constraints coming from the Hawking radiation Dalianis (2019).

We assemble and present our results considering first ${\cal P_{R}}(k)$ that
generate a sizable amount of both PBHs and IGWs and secondly ${\cal P_{R}}(k)$
that generate only IGWs with the PBH counterpart being either negligibly
small, $\Omega_{\text{PBH}}\lesssim 10^{-10}$, or promptly evaporating.

###  V.1 IGWs in scenarios with dominant and subdominant PBH dark matter
component

####  V.1.1 Significant PBH abundance

In the mass window $M_{\text{}}=10^{17}-5\times 10^{22}$ g, the PBH abundance
can reach its maximum value,666In that range there are feeble constraints
coming from white dwarfs and neutron stars Capela  _et al._ (2013a, b); Graham
_et al._ (2015); Brandt (2016) but these constraints, together with others in
the same mass range such as femtolensing and picolensing, are seen as insecure
Katz  _et al._ (2018); Montero-Camacho  _et al._ (2019); Carr  _et al._
(2020). $\Omega_{\text{PBH}}/\Omega_{\text{DM}}=1$; see Ref. Carr  _et al._
(2020). If PBHs constitute a significant fraction of the total dark matter
then an IGW counterpart must exist that can be tested by LISA and by proposed
experiments such as BBO and DECIGO.

Another motivated PBH mass range is $M={\cal O}(1-100)M_{\odot}$ where LIGO
detected several coalescence events the last years with the most recently
published the intriguing event Abbott  _et al._ (2020). For that PBH mass
range the IGWs are expected to produce GWs with nanohertz frequency detectable
by PTA experiments Chen  _et al._ (2020).

The exact PBH mass parameters that we choose can be found in Tables 2 and 3
and in Fig. 6. In the tables, the order of magnitude value for PBH abundance
is quoted due to the exponential sensitivity to the model-dependent threshold
value $\delta_{c}$.

Figure 4: The power spectra of curvature perturbations sourcing IGWs during
RD. We consider a medium and a narrow Gaussian shapes for ${\cal P_{R}}(k)$ as
well as ${\cal P_{R}}(k)$ produced by $\alpha$-attractors (red curves) and
Horndeski nonminimal derivative coupling inflation models (green curve). The
scaling of each ${\cal P_{R}}(k)$ slope is given in Table 1. Oscillations in
the red curves correspond to numerical effects. Left panel: ${\cal P_{R}}(k)$
that generates significant PBH abundances, $\Omega_{\text{PBH}}\sim{\cal
O}(0.01-1)$. Right panel: ${\cal P_{R}}(k)$ that generates negligible PBH
abundances, $\Omega_{\text{PBH}}\lesssim 10^{-10}$.

Figure 5: The IGW spectral shapes for two early universe scenarios, the RD
(top panels) and the eKD$\rightarrow$RD (bottom panels) sourced by ${\cal
P_{R}}(k)$ depicted in Fig. 4. In the left panels we consider the scenario of
abundant PBH production and in the right panels negligible PBH abundances. In
the background of the IGW spectral curves the sensitivity curves of
current/planned/proposed GW detectors are shown as described in the subsection
IV.2. The parameters for the ${\cal P_{R}}(k)$ and the derived values for the
IGWs are listed in the Tables 2-5.

Figure 6: The fractional abundance of the PBHs produced during RD (left panel)
and eKD (right panel) for our ${\cal P_{R}}(k)$ scenarios. The colored upper
bounds are the present experimental constraints as described in Sec. IV. In
the left panel, the red curve corresponds to the $\alpha$-attractor
inflationary model, the green to the Horndeski nonminimal derivative coupling
model, the orange to a Gaussian of width $\epsilon=1$, and the dashed to a
Gaussian of width $\epsilon=0.1$. The ${\cal P_{R}}(k)$ for each curve is
depicted in Fig.4. In the right panel, we take a $\delta$-distribution for the
$\cal P_{R}$; also, for the orange curves the reheating temperature is higher
than for the dashed curves. Detailed values can be found in Tables 2 and 3,
respectively; the associated IGWs are depicted in Fig. 5. $M_{\odot}$ stands
for the sun mass, roughly $2\times 10^{33}$ g.

####  V.1.2 Negligible PBH abundance

Contrary to the exponential sensitivity of the $\Omega_{\text{PBH}}$ on the
$A_{\cal R}$, the sensitivity of the IGWs amplitude on the $A_{\cal R}$ is of
a power-law type. Therefore, a minute decrease in the ${\cal P_{R}}(k)$
amplitude can nearly disappear the PBHs abundance while decrease the IGW
amplitude only by a minor amount. This is rather fascinating because
primordial density perturbations washed out completely by the RD era have left
a relic GW behind that reveals their presence.

In addition, there are PBH masses that evaporate promptly in the early
universe. These PBHs do not survive until today and, unless a stable remnant
is left behind, the corresponding $\Omega_{\text{PBH}}$ is zero and cannot
account for the dark matter in the galaxies. Nevertheless, the associated IGWs
might be strong enough to be detectable. Such a possibility is rather
interesting because the PBH remnants can constitute a significant fraction of
the total $\Omega_{\text{DM}}$ and, moreover, the inflationary scenarios that
predict mini-PBHs can have a spectral index $n_{s}$ in full accordance with
the Planck-2018 data Dalianis and Tringas (2019), see Eq. (B.1) for the
$\alpha$-attractors potential that we use.

Our results for zero or negligible PBH abundance are listed in Tables 4 and 5.
We note that, instead of $\Omega_{\text{PBH}}$, we write the $\beta(M)$
values, see Eq. (27), since PBH with mass $M\lesssim 10^{15}$ g has already
evaporated.

###  V.2 Cosmological eras

####  V.2.1 Radiation

The spectral shapes of the IGWs produced during the radiation era are depicted
in the upper panels of Fig. 5. For narrow Gaussian ${\cal P_{R}}(k)$ peaks,
$\epsilon=0.1$, the spectral shape obtained has a double peak structure with a
sharp peak at $f\simeq 2f_{\text{p}}/\sqrt{3}$, resembling the IGW spectrum
generated by a $\delta$-distribution for ${\cal P_{R}}(k)$. For a broader
${\cal P_{R}}(k)$ peak a single wide peak forms in the IGWs at approximately
$2f_{\textrm{p}}/\sqrt{3}$, see top panels of Fig. 5.

The $\alpha$-attractor inflationary model, described by the potential (38),
generates an enhanced ${\cal P_{R}}(k)$ as depicted in Fig. 4. The predicted
IGW spectrum resembles more that of the broad Gaussian, $\epsilon=1$, with a
peak shifted towards $k_{\text{p}}$, i.e. $f_{\text{p,IGW}}\simeq
f_{\text{p}}$. The general nonminimal derivative coupling model, see Eqs.
(18), (40) and (41), features an IGW spectrum that resembles the broad
Gaussian; hence it also resembles the $\alpha$-attractor model around the
peak. However, the IGW spectra of the two inflationary models deviate at large
frequencies following the different scaling at the high-$k$ tale of the ${\cal
P_{R}}(k)$, as the inspection of Figs. 4 and 5 shows.

For the case the $A_{\cal R}$ amplitude is decreased to a level that the PBH
abundance becomes negligible, $\Omega_{\text{PBH}}\lesssim 10^{-10}$, the
amplitude of the IGWs decreases a little and the spectral shape
characteristics remain unchanged. For the case of evaporating mini-PBHs, we
examine a different $\alpha$-attractors model, described by the potential
(B.1). This model generates a ${\cal P_{R}}(k)$ peak that is relatively
narrow. The spectrum of the IGWs is found to feature a double peak structure,
manifest in the Figs. 4 and 5.

Regarding the energy density, we find that the growth of the IGWs amplitude
occurs rapidly, for a period $\eta_{\text{c}}\sim O(10)k_{\text{p}}^{-1}$, and
then the growth ceases. The growth is realized from the time of entry until
the potential decreases significantly, since $\Phi(x)$ oscillates with a
$x^{-2}$ decaying amplitude in subhorizon scales. Afterward the source term
becomes negligible and the GW propagates freely. At that time $\eta_{c}$ the
$\Omega_{\text{IGW}}$ produced by a Gaussian ${\cal P_{R}}(k)$ peak with width
$\epsilon$ and at the frequency $f_{\text{p}}$ is

$\Omega^{(\text{RD})}_{\text{IGW}}(\eta_{\text{c}},f_{\text{p}})\sim\left(A_{\cal
R}\,\epsilon\sqrt{\pi}\,\right)^{2}.$ (29)

A similar relation has been found following analytic steps in Ref. Pi and
Sasaki (2020). We further comment on the scaling of the IGW spectrum, at
scales beyond the peak, in the next subsection. Note that in Tables 2-5 the
present total energy density value, $\Omega_{\textrm{IGW}}(t_{0})=\int d\ln f\
\Omega_{\textrm{IGW}}(t_{0},f)$, is listed.

####  V.2.2 Early kination

The IGWs produced during the early kination era and sourced by
$\delta$-distributions for ${\cal P_{R}}(k)$ are depicted in the lower panels
of Fig. 5. The increase of the IGW amplitude is depicted in Fig. 2. What is of
interest is that the shape of the IGW spectrum depends on the reheating
temperature, always assuming a sudden change from eKD to RD at the temperature
$T_{\text{rh}}$.

If the $\delta$ peak enters well before the reheating moment,
$\eta_{\text{entry}}\ll\eta_{\text{rh}}$, the induced tensor power spectrum
has a distinctive shape, characteristic of the kination stage. It is
distinctive because the scalar perturbations decay as $x^{3/2}$ experiencing a
maximal pressure. No sharp peak in the IGWs appears even for the monochromatic
${\cal P_{R}}(k)$. This feature has been also noted in Domènech (2019).
Additionally, the energy density of the IGWs gets enhanced due to the redshift
of the stiffer background. We call this transition "slow eKD$\rightarrow$ RD."

If, on the other hand, the reheating happens fast after the horizon entry of
the $\delta$-peak, the scalar perturbations experience the maximal pressure of
the stiff fluid only for a while. It is actually the radiation phase that
mostly forms the IGW spectrum. We call this transition "fast eKD$\rightarrow$
RD". This is readily seen at the IGW spectrum in the PTA frequency range, as
depicted in Fig. 5, where the BBN bound for the $T_{\text{rh}}$ limits the
duration of the kination era. The parameters chosen for this case are
$k_{\text{p}}\eta_{\text{rh}}\approx 1.5$ which means that transition
essentially occurs at the moment of entry of the curvature perturbation.
Apparently, the distinction between the pure RD and the very fast
eKD$\rightarrow$RD transition is hard to be made. By decreasing the
$\eta_{\text{entry}}$, the differences in the amplitude and the $f$-scaling
become manifest.

We note that, in our examples, the amplitude $A_{0}$ for the
$\delta$-distribution together with the $T_{\text{rh}}$ value has been chosen
so that either the $\Omega_{\text{PBH}}$ or the $\Omega_{\text{IGW}}$ is
maximized, given the observational constraints. Thus, these eKD scenarios can
be probed in the near future by gravitational wave observatories. Also, a
change in the $T_{\text{rh}}$, apart from uplifting/downlifting the
$\Omega_{\text{IGW}}$ shifts the position of the spectrum changing the
frequency $f_{\text{p,IGW}}$.

Regarding the energy density, the IGWs at the time $\eta_{\text{c}}$ in the
radiation era have an energy density parameter

$\Omega^{(\text{KD})}_{\text{IGW}}(\eta_{\text{c}})\approx\frac{\eta_{\text{rh}}}{\eta_{\text{entry}}}A_{0}^{2}\,.$
(30)

This result is found for a monochromatic scalar power spectrum. Utilizing the
correspondence between the $\delta$ and the Gaussian distribution with width
$\epsilon$,

$A_{0}\longleftrightarrow A_{\cal R}\equiv\frac{A_{0}}{\epsilon\sqrt{\pi}},$
(31)

one can find the energy density parameter of the IGWs produced during the
kination era for wide ${\cal P_{R}}(k)$ distributions as well.

###  V.3 Fitting broken power laws for IGW spectra

In this subsection we make a few notes on the power laws that describe the IGW
spectral curves in Fig. 5, and facilitate the detectability of our results.
The power laws for each case are listed in Table 1 together with the power-law
scaling of the "mother" curvature power spectrum. We separately describe the
scaling in the IR part of the spectrum for $f<f_{\text{p, IGW}}$ and in the UV
part for $f>f_{\text{p, IGW}}$. For a recent independent investigation
regarding the infrared scaling of the IGWs spectra see also Cai  _et al._
(2019c), and Kuroyanagi  _et al._ (2018) for a broken power law analysis for
general stochastic GW backgrounds.

  * •

The Gaussian models produce IGWs that scale as a power law $f^{3}$ for
$f<f_{\text{p,IGW}}$, despite that the Gaussian ${\cal P_{R}}(k)$ demonstrates
very steep slopes such as $k^{8}$ ($\epsilon=1$) or $k^{40}$ ($\epsilon=0.1$),
estimated at about the half of the curve’s peak height. In the UV frequency
band, $f>f_{\text{p,IGW}}$, the $\epsilon=0.1$ Gaussian is almost cut off
whereas the $\epsilon=1$ falls off roughly as $f^{-10}$.

  * •

The first $\alpha$-attractors inflationary model, given by Eq. (38), produces
a ${\cal P_{R}}(k)$ that increases like $k^{4}$ in the beginning, like $k^{3}$
before the peak and then falls like $k^{-1.2}$. The spectrum of the IGWs
scales mainly $f^{3}$ in the IR frequencies, similarly to the power-law
scaling of the GNMDC model. In the UV it scales like a $f^{-2.7}$ . The IGW
spectrum is found to be broader.

  * •

The second $\alpha$-attractors inflationary model, given by Eq. (B.1),
triggers evaporating mini-PBH formations and induced GWs near the LIGO
frequency band. The produced ${\cal P_{R}}(k)$ increases like $k^{4}$ and
falls like $k^{-4}$. The associated IGWs scale like $f^{3}$ in the IR band and
like $f^{-9}$ in the UV band.

  * •

The Horndeski general nonminimal derivative coupling (GNMC) inflation model
produces a ${\cal P_{R}}(k)$ that scales like $k^{2/5}$ at small wave numbers,
after like $k^{2}$ and then falls like $k^{-4}$. The IGW spectrum has mainly a
$f^{3}$ scaling in the IR band and $f^{-7}$ power law scaling in the UV.

Apparently, the distinction between $\alpha$-attractors and the Horndeski
GNMDC inflationary models can be made in the large-$f$ region.

  * •

For the eKD case, the scaling of IGWs follows a power-law behavior,
proportional to $f$ in the IR band. The UV part drops abruptly, consistent
with our choice of delta distribution for $\mathcal{P}_{\mathcal{R}}$. In the
case of a fast transition into RD, the scaling, as expected, goes like $f^{2}$
in the IR.

We summarize the scaling for IGWs in Table 1.

era | $\bm{{\cal P_{R}}(k)}$ | $\bm{k_{\text{IR}}}$ | $\bm{k_{\text{UV}}}$ | $\bm{f_{\text{IR}}}$ | $\bm{f_{\text{UV}}}$
---|---|---|---|---|---
|  Gauss. $\\#1$ |  |  | $f^{3}$ | cut-off
|  Gauss. $\\#2$ |  |  | $f^{3}$ | $f^{-10}$
RD |  $\alpha$-attr. $\\#1$ | $k^{4},k^{3}$ | $k^{-1.2}$ | $f^{3}$ | $f^{-2.7}$
|  $\alpha$-attr. $\\#2$ | $k^{4}$ | $k^{-4}$ | $f^{3}$ | $f^{-9}$
| GNMDC |  $k^{2/5}$, $k^{2}$ | $k^{-4}$ | $f^{3}$ | $f^{-7}$
| $\delta$ | - | - | $f^{2}$ | cut-off
fast eKD | $\delta$ | - | - | $f^{2}$ | cut-off
slow eKD | $\delta$ | - | - | $f$ | cut-off
Table 1: We list the approximate power-law scaling that describes the IGW
spectral shapes for the ${\cal P_{R}}(k)$ models we studied, from IR to UV
bands. $k_{\text{IR}}$ $(k_{\text{UV}})$ represents the wave numbers
$k<k_{\text{p}}$ ($k>k_{\text{p}}$), and $f_{\text{IR}}$ $(f_{\text{UV}})$ the
frequencies $f<f_{\text{p,IGW}}$ ($f>f_{\text{p,IGW}}$). The "cutoff" means
that the IGW spectrum falls off very abruptly in the UV. The Gauss. $\\#$1 and
Gauss. $\\#$2 models stand for the narrow and medium width Gaussians
$\epsilon=0.1$, $\epsilon=1$, respectively. The $\alpha$-attr. $\\#1$ and
$\alpha$-attr. $\\#2$ stand for the $\alpha$-attractors inflation models with
potentials (38) and (B.1), respectively.

In the IR band, as a general rule, GWs induced during the radiation era have a
spectrum that follows a $f^{3}$ power law under for a broad $\cal P_{R}$ and a
$f^{2}$ power-law for a monochromatic $\cal P_{R}$.777We comment that our IGW
spectra are found to be broad enough and are not expected to experience a
deformation, mentioned recently in Domcke  _et al._ (2020), via the Sachs-
Wolfe and integrated Sachs-Wolfe effect as the GWs travel toward the
detectors. Our result agrees with the conclusions of Cai  _et al._ (2019c) as
well as with the GW spectrum described in Pi and Sasaki (2020).888Note that
our $\epsilon$ is $\sqrt{2}\Delta$ in Pi and Sasaki (2020). In the slow
eKD$\rightarrow$RD scenario, where most of the GWs are induced during the
kination era, the spectrum increases as $f$ in IR frequency band.

In the UV frequency band, for almost all cases, the spectrum falls off very
abruptly and it is cutoff in the monochromatic cases. The exception is the
first $\alpha-$attractors inflation model (38) that features a broader ${\cal
P_{R}}(k)$ that falls off slowly, ${\cal P_{R}}(k)\sim k^{-1.2}$. Hence, broad
scalar power spectra can be distinguished from the narrow ones via the
spectrum of the induced GWs.

##  VI Conclusions

In this work we have explored the features of the scalar-induced GW spectrum
produced by different types of ${\cal P_{R}}(k)$ peaks and for two different
early universe cosmological scenarios: radiation and kination domination. In
addition, we examined two explicit inflationary models that generate PBHs,
$\alpha$-attractors, and Horndeski general nonminimal derivative coupling
models, and tested the predicted IGW signals against the observational
constraints.

Assuming a Gaussian ${\cal P_{R}}(k)$ with amplitude $A_{\cal R}$ and width
$\epsilon$ we find that that the IGWs produced during radiation era have a
spectral energy density parameter today, at the frequency $f_{\text{p}}$ where
the ${\cal P_{R}}$ maximizes, given by

$\displaystyle\Omega^{(\text{RD})}_{\text{IGW}}(t_{0},f_{\text{p}})\sim\,5.2\times
10^{-9}\,\epsilon^{2}\,\left(\frac{g_{*}}{106.75}\right)^{-1/3}\left(\frac{A_{\cal
R}}{10^{-2}}\right)^{2}.$ (32)

For IGWs produced during kination domination, with $T_{\text{rh}}$ the
reheating temperature, the corresponding expression for the energy density
today is

$\displaystyle\Omega^{(\text{KD})}_{\text{IGW}}(t_{0},f_{\text{p}})\sim\,\pi\,\Omega_{\text{IGW}}^{\text{(RD)}}(t_{0},f_{\text{p}})\left(\frac{10^{7}\,\text{GeV}}{T_{\text{rh}}}\right)\left(\frac{f_{\text{p}}}{\text{Hz}}\right)$
(33)

where $\Omega_{\text{IGW}}^{\text{(RD)}}(t_{0},f_{\text{p}})$ is given by Eq.
(32).

The RD and an eKD era may be distinguished by their power-law scaling in the
small-$f$ band, $f<f_{\text{p,IGW}}$. An eKD scenario with a slow reheating
and sudden transition to RD predicts spectral shapes with large amplitudes for
IGWs. Over the next years, aLIGO, reaching its design sensitivity, will put
constraints on the eKD scenario in the high frequency band of the spectrum.
Hints for the reheating temperature can be found if the IGW spectrum has been
modified due to the transition into the radiation era.

A calculation involving a $\delta$-distribution is simpler to implement,
compared to the Gaussian or any other realistic distribution such as those of
$\alpha$-attractors, since the integrals (11) can be computed analytically. We
used $\delta$, i.e, monochromatic, distributions only for the eKD case,
because a broad distribution is computationally more costly in that case.
Utilizing the correspondence, Eq. (31), between $\delta$ and Gaussian
distributions, $A_{0}\longleftrightarrow A_{\cal R}$, the maximum
$\Omega_{\text{IGW}}$ can be found following analytic steps and in a good
approximation, either for the RD or the eKD case.

The spectral shape for the $\Omega_{\text{IGW}}(t_{0},f)$ depends on the
features of the source, the scalar spectrum ${\cal P_{R}}(k)$. It maximizes at
a frequency $f_{\text{p,IGW}}$, in a little offset from $f_{\text{p}}$,
depending on the width of the ${\cal P_{R}}(k)$, as the results listed in the
Tables 4 and 5 demonstrate. Although the power-law scaling of the IGWs
increases roughly universally as $f^{3}$ for RD, it falls off with different
scaling. It is interesting to mention that, for the radiation domination case
at least, the ${\cal P_{R}}(k)$ shape is projected in a much more informative
manner onto the $\Omega_{\text{IGW}}(f)$ spectrum than on the PBH mass
distribution, which is predominantly monochromatic. By observing the IGW
spectral shape and the power law scaling in the UV frequency band one can
infer the width and the amplitude of the scalar spectrum, the generator of the
IGWs.

Consequently, we can say that the detection of the $\Omega_{\text{IGW}}(f)$
spectrum is a portal to the primordial power spectrum of curvature
perturbations, ${\cal P_{R}}(k)$. It can be used to discriminate inflationary
models, and our analysis aimed at contributing to this direction.

## Acknowledgments

We would like to thank A. Kehagias, C. Kouvaris and A. Riotto for discussions.
The work of I.D. is supported by IKY Scholarship, cofinanced by Greece and the
European Union (European Social Fund), through the Operational Program "Human
Resources Development, Education and Lifelong Learning" in the context of the
project “Reinforcement of Postdoctoral Researchers - 2nd Cycle” (MIS-5033021),
implemented by the State Scholarships Foundation.

Ag Al

PBH Mass $\bm{{\cal P_{R}}(k)}$ type $\bm{A_{\cal R}}$
$\bm{{\Omega_{\text{PBH}}}/{\Omega_{\text{DM}}}}$
$\bm{h^{2}\,\Omega_{\text{IGW}}}$ $\bm{f_{\text{p,IGW}}}$ Experi-
ment  lrt Gaussian 0.1 $1.3\times 10^{-1}$ lrt $9.3\times 10^{-10}$ $1.4\times
10^{-1}$ Hz lrt lr Gaussian 1 $4.1\times 10^{-2}$ lr $3.3\times 10^{-9}$
$1.4\times 10^{-1}$ Hz lr lr $\alpha$-attractors $\\#1$ $2.1\times 10^{-2}$ lr
$1.1\times 10^{-9}$ $9.1\times 10^{-2}$ Hz lr lrb$\bm{10^{18}}$ g Galileon
$2.5\times 10^{-2}$ lrb${\cal O}(1)$ $1.1\times 10^{-9}$ $1.4\times 10^{-1}$
Hz lrbDECIGO +
WD/Lensing lrt Gaussian 0.1 $1.5\times 10^{-1}$ lrt $1.4\times 10^{-9}$
$1.2\times 10^{-3}$ Hz lrt lrb$\bm{10^{22}}$ g Gaussian 1 $4.8\times 10^{-2}$
lrb${\cal O}(1)$ $4.4\times 10^{-9}$ $1.5\times 10^{-3}$ Hz lrbLISA +
Lensing lrt Gaussian 0.1 $2.2\times 10^{-1}$ lrt $8.9\times 10^{-9}$
$4.0\times 10^{-10}$ Hz lrt lrb$\bm{10^{35}}$ g Gaussian 1 $6.2\times 10^{-2}$
lrb${\cal O}(10^{-5})$ $2.6\times 10^{-8}$ $4.1\times 10^{-10}$ Hz lrbPTA +
Lensing, X-rays

Table 2:  The values for the PBHs and the associated IGWs produced during
radiation domination are listed. The types of ${\cal P_{R}}(k)$ are listed
together with the $A_{\cal R}$, the amplitude of the ${\cal P_{R}}(k)$ peak,
and an estimation of the produced PBH abundance.
$h^{2}\Omega_{\textrm{IGW}}=h^{2}\int d\ln f\ \Omega_{\textrm{IGW}}(t_{0},f)$
is the present total energy density value. The $f_{\text{p,IGW}}$ is the
frequency that the IGW energy density peaks. In the last column we quote the
experiments that probe each scenario.

Ag Al

PBH Mass $\bm{{\cal P_{R}}(k)}$ type $\bm{A_{0}}$
$\bm{{\Omega_{\text{PBH}}}/{\Omega_{\text{DM}}}}$ $\bm{T_{\text{rh}}}$
$\bm{h^{2}\,\Omega_{\text{IGW}}}$ $\bm{f_{\text{p,IGW}}}$ Experi-
ment  lrt lrt $1.6\times 10^{-2}$ ${\cal O}(1)$ $10^{6}$ GeV $6.2\times
10^{-8}$ $1.1\times 10^{-1}$ Hz lrt lrb$\bm{10^{18}}$ g lrt $1.2\times
10^{-2}$ ${\cal O}(10^{-5})$ $5\times 10^{4}$ GeV $1.3\times 10^{-6}$
$5.4\times 10^{-1}$ Hz lrb DECIGO
\+
WD/Lensing lrt lrb$\delta$-distribution $1.7\times 10^{-2}$ lrt $10^{4}$ GeV
$1.1\times 10^{-6}$ $1.0\times 10^{-3}$ Hz lrt lrb$\bm{10^{22}}$ g lr
$1.5\times 10^{-2}$ lrb${\cal O}(10^{-1})$ $5\times 10^{2}$ GeV $2.2\times
10^{-6}$ $4.2\times 10^{-3}$ Hz lrb LISA
\+
Lensing $\bm{10^{35}}$ g $2.3\times 10^{-2}$ ${\cal O}(10^{-5})$ $0.01$ GeV
$3.1\times 10^{-8}$ $5.5\times 10^{-10}$ Hz PTA +
Lensing/X-rays

Table 3:  As in Table 2, the values for the PBHs and the associated IGWs
produced during kination domination are listed. The $A_{0}$ is the coefficient
of the $\delta$ distribution and $T_{\text{rh}}$ the reheating temperature for
each scenario.

Ag Al

$\bm{f_{\text{hor}}}$ $\bm{{\cal P_{R}}(k)}$ type $\bm{A_{\cal R}}$
$\bm{M_{\text{PBH}}}$ $\bm{\beta_{\text{PBH}}}$
$\bm{h^{2}\,\Omega_{\text{IGW}}}$ $\bm{f_{\text{p,IGW}}}$ Experi-
ment  lrt Gaussian 0.1 $7.8\times 10^{-2}$ lrt lrt $3.4\times 10^{-10}$
$1.1\times 10^{2}$ Hz lrt lr Gaussian 1 $2.5\times 10^{-2}$ lr lr $1.2\times
10^{-9}$ $1.0\times 10^{2}$ Hz lr lrb$\bm{10^{2}}$ Hz $\alpha$-attractors
$\\#2$ $2.0\times 10^{-2}$ lrb$7\times 10^{12}$ g lrb${\cal O}(10^{-24})$
$1.9\times 10^{-10}$ $3.7\times 10^{1}$ Hz lrbLIGO lrt Gaussian 0.1 $9.0\times
10^{-2}$ lrt lrt $4.5\times 10^{-10}$ $1.1\times 10^{-1}$ Hz lrt lr Gaussian 1
$2.7\times 10^{-2}$ lr lr $1.4\times 10^{-9}$ $1.0\times 10^{-1}$ Hz lr
lrb$\bm{10^{-1}}$ Hz $\alpha$-attractors $\\#2$ $1.4\times 10^{-2}$
lrb$7\times 10^{18}$ g lrb${\cal O}(10^{-22})$ $5.2\times 10^{-10}$ $3.4\times
10^{-2}$ Hz lrbDECIGO lrt Gaussian 0.1 $8.7\times 10^{-2}$ lrt lrt $4.2\times
10^{-10}$ $1.1\times 10^{-3}$ Hz lrt lrb$\bm{10^{-3}}$ Hz Gaussian 1
$3.0\times 10^{-2}$ lrb$7\times 10^{22}$ g lrb${\cal O}(10^{-20})$ $1.7\times
10^{-9}$ $1.0\times 10^{-3}$ Hz lrbLISA lrt Gaussian 0.1 $1.2\times 10^{-1}$
lrt lrt $1.7\times 10^{-9}$ $1.1\times 10^{-9}$ Hz lrt lrb$\bm{10^{-9}}$ Hz
Gaussian 1 $3.8\times 10^{-2}$ lrb$7\times 10^{34}$ g lrb${\cal O}(10^{-16})$
$5.9\times 10^{-9}$ $1.0\times 10^{-9}$ Hz lrbPTA

Table 4:  The values for IGWs produced during radiation domination with zero
or negligible PBH abundance, $\Omega_{\text{PBH}}/\Omega_{\text{DM}}\sim
10^{-10}$. The parameters are as in Table 2. The $\beta(M)$ values are listed,
Eq. (27), since PBHs with mass $M\lesssim 10^{15}$ g have already evaporated.

Ag Al

$\bm{f_{\text{hor}}}$ $\bm{{\cal P_{R}}(k)}$ type $\bm{A_{0}}$
$\bm{M_{\text{PBH}}}$ $\bm{\beta_{\text{PBH}}}$ $\bm{T_{\text{rh}}}$
$\bm{h^{2}\,\Omega_{\text{IGW}}}$ $\bm{f_{\text{p,IGW}}}$ Experi-
ment  lrt lrt $9\times 10^{-3}$ $10^{13}$ g ${\cal O}(10^{-25})$ $10^{9}$ GeV
$5.3\times 10^{-9}$ $7.1\times 10^{1}$ Hz lrt lrb$\bm{10^{2}}$ Hz lrt $4\times
10^{-4}$ $5\times 10^{14}$ g ${\cal O}(10^{-550})$ $10^{6}$ GeV $1.7\times
10^{-8}$ $6.2\times 10^{1}$ Hz lrbLIGO $\bm{10^{-1}}$ Hz
lrb$\delta$-distribution $8.8\times 10^{-3}$ $5\times 10^{20}$ g ${\cal
O}(10^{-26})$ $10^{3}$ GeV $2.3\times 10^{-6}$ $3.1\times 10^{-2}$ Hz DECIGO
$\bm{10^{-3}}$ Hz lr $10^{-2}$ $6\times 10^{24}$ g ${\cal O}(10^{-24})$ 10 GeV
$2.2\times 10^{-6}$ $3.1\times 10^{-4}$ Hz LISA $\bm{10^{-9}}$ Hz $1.5\times
10^{-2}$ $10^{35}$ g ${\cal O}(10^{-16})$ 0.02 GeV $7.7\times 10^{-9}$
$1.6\times 10^{-9}$ Hz PTA

Table 5:  As in Table 4, the values for IGWs produced during kination
domination for different reheating temperatures and with zero or negligible
PBH abundance, $\Omega_{\text{PBH}}/\Omega_{\text{DM}}\sim 10^{-10}$, are
listed.

##  Appendix A GWs during kination era and constraints on the ${\cal
P_{R}}(k)$ and the reheating temperature.

The gravitational wave energy density gets enhanced during the kination regime
Giovannini (1999); Riazuelo and Uzan (2000); Yahiro  _et al._ (2002); Boyle
and Buonanno (2008). The energy density of the GWs does not alter BBN
predictions if

$\displaystyle I\equiv
h^{2}\int^{k_{\text{max}}}_{k_{\text{BBN}}}\Omega_{\text{GW}}(k,\eta_{0})d\ln
k\leq 2\times 10^{-6}\,.$ (34)

Equivalently, the above constraint can be written in terms of the
$\Omega_{\text{GW}}(k,\eta_{\text{c}})$. The duration of the kination era is
constrained due to energy density of the runaway field $\varphi$ and the GWs,
see Dimopoulos and Owen (2017); Akrami  _et al._ (2018); Artymowski  _et al._
(2018) for a recent discussion in this context. Here, we derive a bound on the
duration of the kination era coming from the amplitude of the IGWs, thus the
${\cal P_{R}}(k)$ at small scales, see also Kohri and Terada (2018); Domènech
(2019). Let us assume that the energy of GWs is stored mainly in a narrow wave
band $(k_{1},k_{2})$ with central wave number $k_{\text{p}}$ that enters the
horizon at $\eta_{\text{entry}}$. During kination regime the GW energy density
parameter scales as $\Omega_{\text{GW}}\propto a^{2}$. At the time of
reheating it is

$\displaystyle\Omega_{\text{GW}}(\eta_{\text{rh}})\simeq\frac{1}{2}\frac{\rho_{\text{GW}}(\eta_{\text{entry}})}{\rho_{\text{tot}}(\eta_{\text{entry}})}\left(\frac{a(\eta_{\text{rh}})}{a(\eta_{\text{entry}})}\right)^{2}$

where we assumed equipartition of energy densities at $T_{\text{rh}}$; in the
case of sudden transition, the 1/2 factor should be dropped. Therefore,
$\Omega_{\text{GW}}(\eta_{\text{c}})\simeq\Omega_{\text{GW}}(\eta_{\text{rh}})\simeq\frac{1}{2}\Omega_{\text{GW}}(\eta_{\text{entry}})\,({k_{\text{p}}}/{k_{\text{rh}}})$.
Let us make the top-hat approximation for the GW spectrum,
$\Omega_{\text{GW}}=A_{\text{GW}}$ in the interval $(k_{1},k_{2})$. Then, the
integral (34) can be estimated. Plugging in numbers we find approximately the
constraint

$\displaystyle\frac{k_{\text{p}}}{k_{\text{rh}}}\,A_{\text{GW}}\,\ln\left(\frac{k_{2}}{k_{1}}\right)\lesssim
0.4\,.$ (35)

For a ballpark analytic estimation, let us assume the GW spectrum width
$k_{2}/k_{1}=10$ and consider the scaling
$\Omega_{\text{GW}}(\eta_{\text{c}})\sim A_{\cal R}^{2}$ during radiation
domination era.

From the relation $k_{\text{rh}}=2\times
10^{7}(T_{\text{rh}}/\text{GeV})\,\text{Mpc}^{-1}$ and the kination era
relation $k_{\text{p}}=5.4\times
10^{28}g_{*}^{-1/6}(T_{\text{rh}}/\text{GeV})^{-1/3}(M/\text{g})^{-2/3}\gamma^{2/3}\text{Mpc}^{-1}$,
see Eq. (IV.1), a rough conservative lower bound on the reheating temperature
for the kination regime is found,

$\displaystyle{T_{\text{rh}}}\gtrsim 10^{7}\,{\text{GeV}}\,A_{\cal
R}^{3/2}\left(\frac{M_{\text{hor}}}{10^{20}\text{g}}\right)^{-1/2},$ (36)

where $A_{\cal R}\equiv{\cal P_{R}}(k_{\text{p}})$. If this bound is violated,
the BBN predictions are endangered.

##  Appendix B Explicit PBH generating inflation models

In inflationary cosmology the primordial perturbations are produced from
quantum fluctuations with a vast range of wave lenghts. Here, we briefly
present two sort of inflationary models that generate an enhanced amplitude
for the ${\cal P_{R}}(k)$ at small scales: the $\alpha$-attractors Dalianis
_et al._ (2019); Dalianis and Tringas (2019) and the GNMDC inflation models
Dalianis  _et al._ (2020); Fu  _et al._ (2019).

###  B.1 $\alpha$-attractors inflation

In the $\alpha$-attractors inflation scenario Kallosh  _et al._ (2013), the
effective Lagrangian for the inflaton field $\varphi$ turns out to be

$\frac{\cal
L}{\sqrt{-g}}=\frac{1}{2}R-\frac{1}{2}\Big{(}\partial_{\mu}\varphi\Big{)}^{2}-f^{2}\Big{(}\tanh\frac{\varphi}{\sqrt{6\alpha}}\Big{)}\,,$
(37)

where Re$\Phi=\sqrt{3}\tanh({\varphi}/{\sqrt{6\alpha}})$ is a chiral
superfield, and $R$ is the Ricci scalar. We also took $M_{\text{Pl}}=1$. The
canonically normalized inflaton potential that determines the inflationary
trajectory is $V(\varphi)=f^{2}(\tanh\varphi/\sqrt{6\alpha})$. Polynomial,
trigonometric and exponential forms for the function $f(\text{Re}\Phi)$ can
feature an inflection point plateau sufficient to generate a significant dark
matter abundance in accordance with the observational constraints Dalianis
_et al._ (2019). In this work we presented tensor power spectra induced by
scalar power spectra ${\cal P_{R}}(k)$ predicted by the inflationary
potentials

$V(\varphi)=f^{2}_{0}\left[\sum_{n=0}^{3}c_{n}\left(\tanh({\varphi}/{\sqrt{6\alpha}}\right)^{n}\right]^{2}$
(38)

and

$\displaystyle V(\varphi)=f^{2}_{0}\,\big{[}c_{0}+$ $\displaystyle
c_{1}e^{\lambda_{1}\tanh\varphi/\sqrt{6}}\,+$ $\displaystyle
c_{2}e^{\lambda_{2}\left(\tanh(\varphi/\sqrt{6})-\tanh(\varphi_{\text{P}}/\sqrt{6})\right)}\big{]}^{2}\,.$
(39)

The form of the potentials (38) and (B.1) is depicted in Fig. 7. The
inflection point plateau is the field region where the curvature perturbations
get enhanced. The potential (38) implements an RD cosmological
postinflationary phase with high reheating temperature together with large
amplitudes for ${\cal P_{R}}(k)$ that induce GWs. The potential (B.1)
implements an inflationary phase followed by a kination phase. The kination
phase can end gradually, as in the original proposal Dalianis and Tringas
(2019), or via a sudden transition due to e.g. an extra field direction toward
a global minimum. The power spectra generated are depicted in Fig. 4. Details
about the range of the parameter values for the inflationary models (38) and
(B.1) as well details about the reheating temperature and the ${\cal
P_{R}}(k)$ amplitude and the PBH abundances can be found in the works Dalianis
_et al._ (2019) and Dalianis and Tringas (2019).

Figure 7: The figure depicts the $a$-attractors inflationary potentials that
generate CMB anisotropies, enhanced small-scale perturbations that trigger PBH
formation, and induced GWs. A reheating phase follows inflation for the
potential in red and a kination phase for the potential in green. The
amplitude of the potential in green has been adjusted to fit in the same plot.
Figure 8: The main plot depicts the evolution of the inflaton field with
GNMDC. The inner plots depict a Higgs-like potential that together with the
field-dependent GNMDC, $A(\varphi)=1+3H^{2}f(\varphi)$ generates a peak in the
${\cal P_{R}}(k)$, a PBH abundance, and induced GWs.

###  B.2 Horndeski GNMDC inflation

Horndeski theories are scalar-tensor theories of gravity Horndeski (1974);
Kobayashi (2019). The Horndeski inflation model that we examine is the one
with the nonminimal derivative coupling of the scalar field to the Einstein
tensor, with a general $\varphi$-dependent form for the function
$G_{5}(\varphi,X)=\hat{f}(\varphi)$ that is part of the the generalized
Galileon action Deffayet  _et al._ (2011). This term was named general non-
minimal derivative coupling (GNMDC) in Dalianis  _et al._ (2020). The inflaton
action is given by Eq. (18) and the potential that we use is

$\displaystyle V(\varphi)=\frac{\lambda}{4}\varphi^{4}$ (40)

that can be identified with the Higgs Dalianis  _et al._ (2020). The form of
the GNMDC is

$\displaystyle\hat{f}(\varphi)=\frac{\alpha\varphi^{\alpha-1}}{M^{\alpha+1}}\left(1+f_{II}(\varphi)\right)~{},$
(41)

with $f_{II}=d(\left((\varphi-\varphi_{0})/s\right)^{2}+1)^{-1/2}$ that
features a sharp peak at $\varphi_{0}$ Fu  _et al._ (2019). The parameters $s$
and $d$ determine, respectively, the width and the amplitude of the $f_{II}$
term. The GNMDC inflation features a noncanonical kinetic term that acts like
friction, decelerating the inflaton field about the $\varphi_{0}$ value, while
the field rolls down the $\lambda\varphi^{4}$ potential. The form of the GNMDC
enhances significantly the ${\cal P_{R}}(k)$ and triggers PBH formation
Dalianis  _et al._ (2020). The evolution of the inflaton scalar field with a
Higgs-like potential is depicted in Fig. 8 when a field-dependent GNMDC,
called ${\cal A}(\varphi)$, acts.

## References

  * Kuroyanagi  _et al._ (2018) Sachiko Kuroyanagi, Takeshi Chiba, and Tomo Takahashi, “Probing the Universe through the Stochastic Gravitational Wave Background,” JCAP 11, 038 (2018), arXiv:1807.00786 [astro-ph.CO] .
  * Christensen (2019) Nelson Christensen, “Stochastic Gravitational Wave Backgrounds,” Rept. Prog. Phys. 82, 016903 (2019), arXiv:1811.08797 [gr-qc] .
  * Caprini and Figueroa (2018) Chiara Caprini and Daniel G. Figueroa, “Cosmological Backgrounds of Gravitational Waves,” Class. Quant. Grav. 35, 163001 (2018), arXiv:1801.04268 [astro-ph.CO] .
  * Matarrese  _et al._ (1993) Sabino Matarrese, Ornella Pantano, and Diego Saez, “A General relativistic approach to the nonlinear evolution of collisionless matter,” Phys. Rev. D 47, 1311–1323 (1993).
  * Matarrese  _et al._ (1998) Sabino Matarrese, Silvia Mollerach, and Marco Bruni, “Second order perturbations of the Einstein-de Sitter universe,” Phys. Rev. D 58, 043504 (1998), arXiv:astro-ph/9707278 .
  * Noh and Hwang (2004) Hyerim Noh and Jai-chan Hwang, “Second-order perturbations of the Friedmann world model,” Phys. Rev. D 69, 104011 (2004).
  * Carbone and Matarrese (2005) Carmelita Carbone and Sabino Matarrese, “A Unified treatment of cosmological perturbations from super-horizon to small scales,” Phys. Rev. D 71, 043508 (2005), arXiv:astro-ph/0407611 .
  * Nakamura (2007) Kouji Nakamura, “Second-order gauge invariant cosmological perturbation theory: Einstein equations in terms of gauge invariant variables,” Prog. Theor. Phys. 117, 17–74 (2007), arXiv:gr-qc/0605108 .
  * Mollerach  _et al._ (2004) Silvia Mollerach, Diego Harari, and Sabino Matarrese, “CMB polarization from secondary vector and tensor modes,” Phys. Rev. D 69, 063002 (2004), arXiv:astro-ph/0310711 .
  * Ananda  _et al._ (2007) Kishore N. Ananda, Chris Clarkson, and David Wands, “The Cosmological gravitational wave background from primordial density perturbations,” Phys. Rev. D 75, 123518 (2007), arXiv:gr-qc/0612013 .
  * Assadullahi and Wands (2009) Hooshyar Assadullahi and David Wands, “Gravitational waves from an early matter era,” Phys. Rev. D 79, 083511 (2009), arXiv:0901.0989 [astro-ph.CO] .
  * Baumann  _et al._ (2007) Daniel Baumann, Paul J. Steinhardt, Keitaro Takahashi, and Kiyotomo Ichiki, “Gravitational Wave Spectrum Induced by Primordial Scalar Perturbations,” Phys. Rev. D 76, 084019 (2007), arXiv:hep-th/0703290 .
  * Saito and Yokoyama (2009) Ryo Saito and Jun’ichi Yokoyama, “Gravitational wave background as a probe of the primordial black hole abundance,” Phys. Rev. Lett. 102, 161101 (2009), [Erratum: Phys.Rev.Lett. 107, 069901 (2011)], arXiv:0812.4339 [astro-ph] .
  * Dalianis (2019) Ioannis Dalianis, “Constraints on the curvature power spectrum from primordial black hole evaporation,” JCAP 08, 032 (2019), arXiv:1812.09807 [astro-ph.CO] .
  * Kalaja  _et al._ (2019) Alba Kalaja, Nicola Bellomo, Nicola Bartolo, Daniele Bertacca, Sabino Matarrese, Ilia Musco, Alvise Raccanelli, and Licia Verde, “From Primordial Black Holes Abundance to Primordial Curvature Power Spectrum (and back),” JCAP 10, 031 (2019), arXiv:1908.03596 [astro-ph.CO] .
  * Carr  _et al._ (2020) Bernard Carr, Kazunori Kohri, Yuuiti Sendouda, and Jun’ichi Yokoyama, “Constraints on Primordial Black Holes,” (2020), arXiv:2002.12778 [astro-ph.CO] .
  * Sasaki  _et al._ (2018) Misao Sasaki, Teruaki Suyama, Takahiro Tanaka, and Shuichiro Yokoyama, “Primordial black holes—perspectives in gravitational wave astronomy,” Class. Quant. Grav. 35, 063001 (2018), arXiv:1801.05235 [astro-ph.CO] .
  * Bugaev and Klimai (2010a) E.V. Bugaev and P.A. Klimai, “Bound on induced gravitational wave background from primordial black holes,” JETP Lett. 91, 1–5 (2010a), arXiv:0911.0611 [astro-ph.CO] .
  * Saito and Yokoyama (2010) R. Saito and J. Yokoyama, “Gravitational-wave constraints on the abundance of primordial black holes,” Progress of Theoretical Physics 123, 867–886 (2010).
  * Bugaev and Klimai (2010b) Edgar Bugaev and Peter Klimai, “Induced gravitational wave background and primordial black holes,” Phys. Rev. D 81, 023517 (2010b), arXiv:0908.0664 [astro-ph.CO] .
  * Bugaev and Klimai (2011) Edgar Bugaev and Peter Klimai, “Constraints on the induced gravitational wave background from primordial black holes,” Phys. Rev. D 83, 083521 (2011), arXiv:1012.4697 [astro-ph.CO] .
  * Nakama  _et al._ (2017) Tomohiro Nakama, Joseph Silk, and Marc Kamionkowski, “Stochastic gravitational waves associated with the formation of primordial black holes,” Phys. Rev. D 95, 043511 (2017), arXiv:1612.06264 [astro-ph.CO] .
  * Orlofsky  _et al._ (2017) Nicholas Orlofsky, Aaron Pierce, and James D. Wells, “Inflationary theory and pulsar timing investigations of primordial black holes and gravitational waves,” Phys. Rev. D 95, 063518 (2017), arXiv:1612.05279 [astro-ph.CO] .
  * Garcia-Bellido  _et al._ (2017) Juan Garcia-Bellido, Marco Peloso, and Caner Unal, “Gravitational Wave signatures of inflationary models from Primordial Black Hole Dark Matter,” JCAP 09, 013 (2017), arXiv:1707.02441 [astro-ph.CO] .
  * Cai  _et al._ (2019a) Rong-gen Cai, Shi Pi, and Misao Sasaki, “Gravitational Waves Induced by non-Gaussian Scalar Perturbations,” Phys. Rev. Lett. 122, 201101 (2019a), arXiv:1810.11000 [astro-ph.CO] .
  * Inomata and Nakama (2019) Keisuke Inomata and Tomohiro Nakama, “Gravitational waves induced by scalar perturbations as probes of the small-scale primordial spectrum,” Phys. Rev. D 99, 043511 (2019), arXiv:1812.00674 [astro-ph.CO] .
  * Espinosa  _et al._ (2018) José Ramón Espinosa, Davide Racco, and Antonio Riotto, “A Cosmological Signature of the SM Higgs Instability: Gravitational Waves,” JCAP 09, 012 (2018), arXiv:1804.07732 [hep-ph] .
  * Bartolo  _et al._ (2019) N. Bartolo, V. De Luca, G. Franciolini, M. Peloso, D. Racco, and A. Riotto, “Testing primordial black holes as dark matter with LISA,” Phys. Rev. D 99, 103521 (2019), arXiv:1810.12224 [astro-ph.CO] .
  * Byrnes  _et al._ (2019) Christian T. Byrnes, Philippa S. Cole, and Subodh P. Patil, “Steepest growth of the power spectrum and primordial black holes,” JCAP 06, 028 (2019), arXiv:1811.11158 [astro-ph.CO] .
  * Cai  _et al._ (2019b) Yi-Fu Cai, Chao Chen, Xi Tong, Dong-Gang Wang, and Sheng-Feng Yan, “When Primordial Black Holes from Sound Speed Resonance Meet a Stochastic Background of Gravitational Waves,” Phys. Rev. D 100, 043518 (2019b), arXiv:1902.08187 [astro-ph.CO] .
  * Wang  _et al._ (2019) Sai Wang, Takahiro Terada, and Kazunori Kohri, “Prospective constraints on the primordial black hole abundance from the stochastic gravitational-wave backgrounds produced by coalescing events and curvature perturbations,” Phys. Rev. D 99, 103531 (2019), [Erratum: Phys.Rev.D 101, 069901 (2020)], arXiv:1903.05924 [astro-ph.CO] .
  * Gong (2019) Jinn-Ouk Gong, “Analytic integral solutions for induced gravitational waves,” (2019), arXiv:1909.12708 [gr-qc] .
  * Mahbub (2020) Rafid Mahbub, “Primordial black hole formation in inflationary $\alpha$-attractor models,” Phys. Rev. D 101, 023533 (2020), arXiv:1910.10602 [astro-ph.CO] .
  * Özsoy and Tasinato (2020) Ogan Özsoy and Gianmassimo Tasinato, “On the slope of curvature power spectrum in non-attractor inflation,” JCAP 04, 048 (2020), arXiv:1912.01061 [astro-ph.CO] .
  * Bartolo  _et al._ (2020) N. Bartolo, D. Bertacca, V. De Luca, G. Franciolini, S. Matarrese, M. Peloso, A. Ricciardone, A. Riotto, and G. Tasinato, “Gravitational wave anisotropies from primordial black holes,” JCAP 02, 028 (2020), arXiv:1909.12619 [astro-ph.CO] .
  * Braglia  _et al._ (2020) Matteo Braglia, Dhiraj Kumar Hazra, Fabio Finelli, George F. Smoot, L. Sriramkumar, and Alexei A. Starobinsky, “Generating PBHs and small-scale GWs in two-field models of inflation,” (2020), arXiv:2005.02895 [astro-ph.CO] .
  * Ballesteros  _et al._ (2020) Guillermo Ballesteros, Julián Rey, Marco Taoso, and Alfredo Urbano, “Primordial black holes as dark matter and gravitational waves from single-field polynomial inflation,” (2020), arXiv:2001.08220 [astro-ph.CO] .
  * Ota (2020) Atsuhisa Ota, “Induced superhorizon tensor perturbations from anisotropic non-Gaussianity,” Phys. Rev. D 101, 103511 (2020), arXiv:2001.00409 [astro-ph.CO] .
  * Allahverdi  _et al._ (2020) Rouzbeh Allahverdi _et al._ , “The First Three Seconds: a Review of Possible Expansion Histories of the Early Universe,” (2020), arXiv:2006.16182 [astro-ph.CO] .
  * Inomata  _et al._ (2020a) Keisuke Inomata, Kazunori Kohri, Tomohiro Nakama, and Takahiro Terada, “Enhancement of gravitational waves induced by scalar perturbations due to a sudden transition from an early matter era to the radiation era,”  _Proceedings, 16th International Conference on Topics in Astroparticle and Underground Physics (TAUP 2019): Toyama, Japan_ , J. Phys. Conf. Ser. 1468, 012002 (2020a).
  * Inomata  _et al._ (2020b) Keisuke Inomata, Kazunori Kohri, Tomohiro Nakama, and Takahiro Terada, “Gravitational waves induced by scalar perturbations during a gradual transition from an early matter era to the radiation era,”  _Proceedings, 16th International Conference on Topics in Astroparticle and Underground Physics (TAUP 2019): Toyama, Japan_ , J. Phys. Conf. Ser. 1468, 012001 (2020b).
  * Domènech (2019) Guillem Domènech, “Induced gravitational waves in a general cosmological background,”  (2019), 10.1142/S0218271820500285, arXiv:1912.05583 [gr-qc] .
  * Domènech  _et al._ (2020) Guillem Domènech, Shi Pi, and Misao Sasaki, “Induced gravitational waves as a probe of thermal history of the universe,” (2020), arXiv:2005.12314 [gr-qc] .
  * Kallosh  _et al._ (2013) Renata Kallosh, Andrei Linde, and Diederik Roest, “Superconformal Inflationary $\alpha$-Attractors,” JHEP 11, 198 (2013), arXiv:1311.0472 [hep-th] .
  * Horndeski (1974) Gregory Walter Horndeski, “Second-order scalar-tensor field equations in a four-dimensional space,” Int. J. Theor. Phys. 10, 363–384 (1974).
  * Kobayashi (2019) Tsutomu Kobayashi, “Horndeski theory and beyond: a review,” Rept. Prog. Phys. 82, 086901 (2019), arXiv:1901.07183 [gr-qc] .
  * Deffayet  _et al._ (2011) C. Deffayet, Xian Gao, D.A. Steer, and G. Zahariade, “From k-essence to generalised Galileons,” Phys. Rev. D 84, 064039 (2011), arXiv:1103.3260 [hep-th] .
  * Dalianis  _et al._ (2019) Ioannis Dalianis, Alex Kehagias, and George Tringas, “Primordial black holes from $\alpha$-attractors,” JCAP 01, 037 (2019), arXiv:1805.09483 [astro-ph.CO] .
  * Dalianis and Tringas (2019) Ioannis Dalianis and George Tringas, “Primordial black hole remnants as dark matter produced in thermal, matter, and runaway-quintessence postinflationary scenarios,” Phys. Rev. D 100, 083512 (2019), arXiv:1905.01741 [astro-ph.CO] .
  * Dalianis  _et al._ (2020) Ioannis Dalianis, Stelios Karydas, and Eleftherios Papantonopoulos, “Generalized Non-Minimal Derivative Coupling: Application to Inflation and Primordial Black Hole Production,” JCAP 06, 040 (2020), arXiv:1910.00622 [astro-ph.CO] .
  * Cai  _et al._ (2019c) Rong-Gen Cai, Shi Pi, and Misao Sasaki, “Universal infrared scaling of gravitational wave background spectra,” (2019c), arXiv:1909.13728 [astro-ph.CO] .
  * Pi and Sasaki (2020) Shi Pi and Misao Sasaki, “Gravitational Waves Induced by Scalar Perturbations with a Lognormal Peak,” (2020), arXiv:2005.12306 [gr-qc] .
  * Kohri and Terada (2018) Kazunori Kohri and Takahiro Terada, “Semianalytic calculation of gravitational wave spectrum nonlinearly induced from primordial curvature perturbations,” Phys. Rev. D97, 123532 (2018), arXiv:1804.08577 [gr-qc] .
  * Tomikawa and Kobayashi (2020) Keitaro Tomikawa and Tsutomu Kobayashi, “Gauge dependence of gravitational waves generated at second order from scalar perturbations,” Phys. Rev. D 101, 083529 (2020), arXiv:1910.01880 [gr-qc] .
  * Garcia-Bellido and Ruiz Morales (2017) Juan Garcia-Bellido and Ester Ruiz Morales, “Primordial black holes from single field models of inflation,” Phys. Dark Univ. 18, 47–54 (2017), arXiv:1702.03901 [astro-ph.CO] .
  * Germani and Prokopec (2017) Cristiano Germani and Tomislav Prokopec, “On primordial black holes from an inflection point,” Phys. Dark Univ. 18, 6–10 (2017), arXiv:1706.04226 [astro-ph.CO] .
  * Motohashi and Hu (2017) Hayato Motohashi and Wayne Hu, “Primordial Black Holes and Slow-Roll Violation,” Phys. Rev. D 96, 063503 (2017), arXiv:1706.06784 [astro-ph.CO] .
  * Ballesteros and Taoso (2018) Guillermo Ballesteros and Marco Taoso, “Primordial black hole dark matter from single field inflation,” Phys. Rev. D 97, 023501 (2018), arXiv:1709.05565 [hep-ph] .
  * Germani and Kehagias (2010) Cristiano Germani and Alex Kehagias, “New Model of Inflation with Non-minimal Derivative Coupling of Standard Model Higgs Boson to Gravity,” Phys. Rev. Lett. 105, 011302 (2010), arXiv:1003.2635 [hep-ph] .
  * Kuroyanagi  _et al._ (2015) Sachiko Kuroyanagi, Kazunori Nakayama, and Jun’ichi Yokoyama, “Prospects of determination of reheating temperature after inflation by DECIGO,” PTEP 2015, 013E02 (2015), arXiv:1410.6618 [astro-ph.CO] .
  * Bernal and Hajkarim (2019) Nicolás Bernal and Fazlollah Hajkarim, “Primordial Gravitational Waves in Nonstandard Cosmologies,” Phys. Rev. D 100, 063502 (2019), arXiv:1905.10410 [astro-ph.CO] .
  * Carr and Hawking (1974) Bernard J. Carr and S.W. Hawking, “Black holes in the early Universe,” Mon. Not. Roy. Astron. Soc. 168, 399–415 (1974).
  * Meszaros (1974) P. Meszaros, “The behaviour of point masses in an expanding cosmological substratum,” Astron. Astrophys. 37, 225–228 (1974).
  * Carr (1975) Bernard J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201, 1–19 (1975).
  * Lindley (1980) D. Lindley, “A COLD BIG BANG WITH SMALL BLACK HOLES.” in  _International School of Nuclear Physics: Nuclear Astrophysics_ (1980) pp. 279–283.
  * Page and Hawking (1976) Don N. Page and S.W. Hawking, “Gamma rays from primordial black holes,” Astrophys. J. 206, 1–7 (1976).
  * MacGibbon and Carr (1991) Jane H. MacGibbon and Bernard J. Carr, “Cosmic rays from primordial black holes,” Astrophys. J. 371, 447–469 (1991).
  * Carr and MacGibbon (1998) Bernard J. Carr and J.H. MacGibbon, “Cosmic rays from primordial black holes and constraints on the early universe,” Phys. Rept. 307, 141–154 (1998).
  * Barrau  _et al._ (2003) A. Barrau, G. Boudoul, and L. Derome, “An improved gamma-ray limit on the density of pbhs,” in  _28th International Cosmic Ray Conference_ (2003) pp. 1697–1699, arXiv:astro-ph/0304528 .
  * Carr  _et al._ (2016a) B.J. Carr, Kazunori Kohri, Yuuiti Sendouda, and Jun’ichi Yokoyama, “Constraints on primordial black holes from the Galactic gamma-ray background,” Phys. Rev. D 94, 044029 (2016a), arXiv:1604.05349 [astro-ph.CO] .
  * Dasgupta  _et al._ (2019) Basudeb Dasgupta, Ranjan Laha, and Anupam Ray, “Neutrino and positron constraints on spinning primordial black hole dark matter,” (2019), arXiv:1912.01014 [hep-ph] .
  * Barnacka  _et al._ (2012) A. Barnacka, J.F. Glicenstein, and R. Moderski, “New constraints on primordial black holes abundance from femtolensing of gamma-ray bursts,” Phys. Rev. D 86, 043001 (2012), arXiv:1204.2056 [astro-ph.CO] .
  * Tisserand  _et al._ (2007) P. Tisserand _et al._ (EROS-2), “Limits on the Macho Content of the Galactic Halo from the EROS-2 Survey of the Magellanic Clouds,” Astron. Astrophys. 469, 387–404 (2007), arXiv:astro-ph/0607207 .
  * Niikura  _et al._ (2019) Hiroko Niikura _et al._ , “Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations,” Nature Astron. 3, 524–534 (2019), arXiv:1701.02151 [astro-ph.CO] .
  * Ricotti  _et al._ (2008) Massimo Ricotti, Jeremiah P. Ostriker, and Katherine J. Mack, “Effect of Primordial Black Holes on the Cosmic Microwave Background and Cosmological Parameter Estimates,” Astrophys. J. 680, 829 (2008), arXiv:0709.0524 [astro-ph] .
  * Carr  _et al._ (2016b) Bernard Carr, Florian Kuhnel, and Marit Sandstad, “Primordial Black Holes as Dark Matter,” Phys. Rev. D 94, 083504 (2016b), arXiv:1607.06077 [astro-ph.CO] .
  * Clesse and García-Bellido (2017) Sebastien Clesse and Juan García-Bellido, “The clustering of massive Primordial Black Holes as Dark Matter: measuring their mass distribution with Advanced LIGO,” Phys. Dark Univ. 15, 142–147 (2017), arXiv:1603.05234 [astro-ph.CO] .
  * Bird  _et al._ (2016) Simeon Bird, Ilias Cholis, Julian B. Muñoz, Yacine Ali-Haïmoud, Marc Kamionkowski, Ely D. Kovetz, Alvise Raccanelli, and Adam G. Riess, “Did LIGO detect dark matter?” Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Poulin  _et al._ (2017) Vivian Poulin, Pasquale D. Serpico, Francesca Calore, Sebastien Clesse, and Kazunori Kohri, “CMB bounds on disk-accreting massive primordial black holes,” Phys. Rev. D 96, 083524 (2017), arXiv:1707.04206 [astro-ph.CO] .
  * Serpico  _et al._ (2020) Pasquale D. Serpico, Vivian Poulin, Derek Inman, and Kazunori Kohri, “Cosmic microwave background bounds on primordial black holes including dark matter halo accretion,” Phys. Rev. Res. 2, 023204 (2020), arXiv:2002.10771 [astro-ph.CO] .
  * Gaggero  _et al._ (2017) Daniele Gaggero, Gianfranco Bertone, Francesca Calore, Riley M. T. Connors, Mark Lovell, Sera Markoff, and Emma Storm, “Searching for Primordial Black Holes in the radio and X-ray sky,” Phys. Rev. Lett. 118, 241101 (2017), arXiv:1612.00457 [astro-ph.HE] .
  * Sasaki  _et al._ (2016) Misao Sasaki, Teruaki Suyama, Takahiro Tanaka, and Shuichiro Yokoyama, “Primordial Black Hole Scenario for the Gravitational-Wave Event GW150914,” Phys. Rev. Lett. 117, 061101 (2016), [Erratum: Phys.Rev.Lett. 121, 059901 (2018)], arXiv:1603.08338 [astro-ph.CO] .
  * Eroshenko (2018) Yu N. Eroshenko, “Gravitational waves from primordial black holes collisions in binary systems,” J. Phys. Conf. Ser. 1051, 012010 (2018), arXiv:1604.04932 [astro-ph.CO] .
  * Raidal  _et al._ (2017) Martti Raidal, Ville Vaskonen, and Hardi Veermäe, “Gravitational Waves from Primordial Black Hole Mergers,” JCAP 09, 037 (2017), arXiv:1707.01480 [astro-ph.CO] .
  * Abbott  _et al._ (2018a) B.P. Abbott _et al._ (LIGO Scientific, Virgo), “Search for Subsolar-Mass Ultracompact Binaries in Advanced LIGO’s First Observing Run,” Phys. Rev. Lett. 121, 231103 (2018a), arXiv:1808.04771 [astro-ph.CO] .
  * Harada  _et al._ (2013) Tomohiro Harada, Chul-Moon Yoo, and Kazunori Kohri, “Threshold of primordial black hole formation,” Phys. Rev. D 88, 084051 (2013), [Erratum: Phys.Rev.D 89, 029903 (2014)], arXiv:1309.4201 [astro-ph.CO] .
  * Abbott  _et al._ (2018b) B.P. Abbott _et al._ (KAGRA, LIGO Scientific, VIRGO), “Prospects for Observing and Localizing Gravitational-Wave Transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel. 21, 3 (2018b), arXiv:1304.0670 [gr-qc] .
  * Sathyaprakash  _et al._ (2012) B. Sathyaprakash _et al._ , “Scientific Objectives of Einstein Telescope,” Class. Quant. Grav. 29, 124013 (2012), [Erratum: Class.Quant.Grav. 30, 079501 (2013)], arXiv:1206.0331 [gr-qc] .
  * Seto  _et al._ (2001) Naoki Seto, Seiji Kawamura, and Takashi 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, 221103 (2001), arXiv:astro-ph/0108011 .
  * Sato  _et al._ (2017) Shuichi Sato _et al._ , “The status of DECIGO,” J. Phys. Conf. Ser. 840, 012010 (2017).
  * Crowder and Cornish (2005) Jeff Crowder and Neil J. Cornish, “Beyond LISA: Exploring future gravitational wave missions,” Phys. Rev. D 72, 083005 (2005), arXiv:gr-qc/0506015 .
  * Amaro-Seoane  _et al._ (2017) Pau Amaro-Seoane _et al._ (LISA), “Laser Interferometer Space Antenna,” (2017), arXiv:1702.00786 [astro-ph.IM] .
  * Luo  _et al._ (2016) Jun Luo _et al._ (TianQin), “TianQin: a space-borne gravitational wave detector,” Class. Quant. Grav. 33, 035010 (2016), arXiv:1512.02076 [astro-ph.IM] .
  * Hobbs  _et al._ (2010) G. Hobbs _et al._ , “The international pulsar timing array project: using pulsars as a gravitational wave detector,” Class. Quant. Grav. 27, 084013 (2010), arXiv:0911.5206 [astro-ph.SR] .
  * Janssen  _et al._ (2015) Gemma Janssen _et al._ , “Gravitational wave astronomy with the SKA,” PoS AASKA14, 037 (2015), arXiv:1501.00127 [astro-ph.IM] .
  * Sathyaprakash and Schutz (2009) B.S. Sathyaprakash and B.F. Schutz, “Physics, Astrophysics and Cosmology with Gravitational Waves,” Living Rev. Rel. 12, 2 (2009), arXiv:0903.0338 [gr-qc] .
  * Yagi and Seto (2011) Kent Yagi and Naoki Seto, “Detector configuration of DECIGO/BBO and identification of cosmological neutron-star binaries,” Phys. Rev. D 83, 044011 (2011), [Erratum: Phys.Rev.D 95, 109901 (2017)], arXiv:1101.3940 [astro-ph.CO] .
  * Huang  _et al._ (2020) Shun-Jia Huang, Yi-Ming Hu, Valeriya Korol, Peng-Cheng Li, Zheng-Cheng Liang, Yang Lu, Hai-Tian Wang, Shenghua Yu, and Jianwei Mei, “Science with the TianQin Observatory: Preliminary Results on Galactic Double White Dwarf Binaries,” (2020), arXiv:2005.07889 [astro-ph.HE] .
  * Chen  _et al._ (2020) Zu-Cheng Chen, Chen Yuan, and Qing-Guo Huang, “Pulsar Timing Array Constraints on Primordial Black Holes with NANOGrav 11-Year Data Set,” Phys. Rev. Lett. 124, 251101 (2020), arXiv:1910.12239 [astro-ph.CO] .
  * Cai  _et al._ (2019d) Rong-Gen Cai, Shi Pi, Shao-Jiang Wang, and Xing-Yu Yang, “Pulsar Timing Array Constraints on the Induced Gravitational Waves,” JCAP 10, 059 (2019d), arXiv:1907.06372 [astro-ph.CO] .
  * Capela  _et al._ (2013a) Fabio Capela, Maxim Pshirkov, and Peter Tinyakov, “Constraints on Primordial Black Holes as Dark Matter Candidates from Star Formation,” Phys. Rev. D 87, 023507 (2013a), arXiv:1209.6021 [astro-ph.CO] .
  * Capela  _et al._ (2013b) Fabio Capela, Maxim Pshirkov, and Peter Tinyakov, “Constraints on primordial black holes as dark matter candidates from capture by neutron stars,” Phys. Rev. D 87, 123524 (2013b), arXiv:1301.4984 [astro-ph.CO] .
  * Graham  _et al._ (2015) Peter W. Graham, Surjeet Rajendran, and Jaime Varela, “Dark Matter Triggers of Supernovae,” Phys. Rev. D 92, 063007 (2015), arXiv:1505.04444 [hep-ph] .
  * Brandt (2016) Timothy D. Brandt, “Constraints on MACHO Dark Matter from Compact Stellar Systems in Ultra-Faint Dwarf Galaxies,” Astrophys. J. Lett. 824, L31 (2016), arXiv:1605.03665 [astro-ph.GA] .
  * Katz  _et al._ (2018) Andrey Katz, Joachim Kopp, Sergey Sibiryakov, and Wei Xue, “Femtolensing by Dark Matter Revisited,” JCAP 12, 005 (2018), arXiv:1807.11495 [astro-ph.CO] .
  * Montero-Camacho  _et al._ (2019) Paulo Montero-Camacho, Xiao Fang, Gabriel Vasquez, Makana Silva, and Christopher M. Hirata, “Revisiting constraints on asteroid-mass primordial black holes as dark matter candidates,” JCAP 08, 031 (2019), arXiv:1906.05950 [astro-ph.CO] .
  * Abbott  _et al._ (2020) R. Abbott _et al._ (LIGO Scientific, Virgo), “GW190814: Gravitational Waves from the Coalescence of a 23 Solar Mass Black Hole with a 2.6 Solar Mass Compact Object,” Astrophys. J. 896, L44 (2020), arXiv:2006.12611 [astro-ph.HE] .
  * Domcke  _et al._ (2020) Valerie Domcke, Ryusuke Jinno, and Henrique Rubira, “Deformation of the gravitational wave spectrum by density perturbations,” (2020), arXiv:2002.11083 [astro-ph.CO] .
  * Giovannini (1999) Massimo Giovannini, “Production and detection of relic gravitons in quintessential inflationary models,” Phys. Rev. D 60, 123511 (1999), arXiv:astro-ph/9903004 .
  * Riazuelo and Uzan (2000) Alain Riazuelo and Jean-Philippe Uzan, “Quintessence and gravitational waves,” Phys. Rev. D 62, 083506 (2000), arXiv:astro-ph/0004156 .
  * Yahiro  _et al._ (2002) M. Yahiro, G.J. Mathews, K. Ichiki, T. Kajino, and M. Orito, “Constraints on cosmic quintessence and quintessential inflation,” Phys. Rev. D 65, 063502 (2002), arXiv:astro-ph/0106349 .
  * Boyle and Buonanno (2008) Latham A. Boyle and Alessandra Buonanno, “Relating gravitational wave constraints from primordial nucleosynthesis, pulsar timing, laser interferometers, and the CMB: Implications for the early Universe,” Phys. Rev. D 78, 043531 (2008), arXiv:0708.2279 [astro-ph] .
  * Dimopoulos and Owen (2017) Konstantinos Dimopoulos and Charlotte Owen, “Quintessential Inflation with $\alpha$-attractors,” JCAP 06, 027 (2017), arXiv:1703.00305 [gr-qc] .
  * Akrami  _et al._ (2018) Yashar Akrami, Renata Kallosh, Andrei Linde, and Valeri Vardanyan, “Dark energy, $\alpha$-attractors, and large-scale structure surveys,” JCAP 06, 041 (2018), arXiv:1712.09693 [hep-th] .
  * Artymowski  _et al._ (2018) Michal Artymowski, Olga Czerwinska, Zygmunt Lalak, and Marek Lewicki, “Gravitational wave signals and cosmological consequences of gravitational reheating,” JCAP 04, 046 (2018), arXiv:1711.08473 [astro-ph.CO] .
  * Fu  _et al._ (2019) Chengjie Fu, Puxun Wu, and Hongwei Yu, “Primordial Black Holes from Inflation with Nonminimal Derivative Coupling,” Phys. Rev. D 100, 063532 (2019), arXiv:1907.05042 [astro-ph.CO] .