# On the Anisotropy of the Stochastic Gravitational Wave Background
from Sub-Horizon-Collapsed Primordial Black Hole Mergers

Stefano Profumo1,2 Fengwei Yang3 1Department of Physics, University of
California, Santa Cruz Santa Cruz, CA 95064, USA  2Santa Cruz Institute for
Particle Physics, Santa Cruz, CA 95064, USA  3Department of Physics and
Astronomy, University of Utah, Salt Lake City, UT 84112, USA

###### Abstract

We study the properties of the stochastic gravitational wave background (SGWB)
resulting from the mergers of primordial black holes (PBH) that formed from
the collapse of sub-horizon regions in the early universe. We adopt a model-
independent approach, where we parameterize the fraction $f_{H}$ of the
wavelength of the perturbation mode in units of the horizon radius when the
patch starts to gravitationally collapse. Assuming a monochromatic spectrum of
isocurvature perturbations and spherically-symmetric density perturbations, we
investigate the isotropic SGWB energy density and angular power spectrum at
various frequencies, PBH masses, and horizon size fractions. The key effect of
sub-horizon formation is a change in the PBH mass function and formation
redshift, which, in turn, affects gravitational wave (GW) observables. We find
that sub-horizon PBH formation in general enhances the isotropic SGWB energy
density and the absolute angular power spectrum. However, the quasi-monotonic
increases in both quantities as $f_{H}$ decreases cease when the chirp mass of
the binary PBHs reaches a mass threshold determined by the frequency of
observation; the isotropic SGWB energy density spectrum significantly drops
above the corresponding cutoff frequency.

##  I Introduction

The intriguing possibility that black holes form from large density
fluctuations in the early universe instead of from the gravitational collapse
of large astrophysical objects in the late universe was proposed long ago
Hawking (1971). Such “primordial black holes” (PBHs) might be related to the
cosmological dark matter Carr and Kuhnel (2020), to the matter-antimatter
asymmetry (see e.g. Smyth  _et al._ (2022) and references therein), or both
Morrison  _et al._ (2019); they could play a role in seeding supermassive
black holes Bean and Magueijo (2002), or even explain part (or perhaps all) of
the black hole-black hole mergers observed by gravitational wave (GW)
interferometers Bird  _et al._ (2016). While the direct detection of even a
single black hole with a mass below the Tolman-Oppenheimer-Volkoff limit
Belczynski  _et al._ (2010) would indicate the existence of PBHs (or of
unexpected new physics Giffin  _et al._ (2022)), other observables may point
towards a non-astrophysical origin for black holes. For instance, the spin
distribution expected from PBHs, whether produced in matter- or radiation-
domination in the early universe, is expected to be markedly, or at least
statistically significantly different than that expected for astrophysical
black holes Fernandez and Profumo (2019).

While no evidence for a stochastic background of GWs has been conclusively
detected yet (see however Arzoumanian  _et al._ (2020)), there are several
reasons to believe that it will soon be. A perfectly plausible source of such
diffuse background are unresolved binary black hole mergers at all redshifts,
including of PBHs Shannon  _et al._ (2015). Some information on the origin and
production mechanism of the black holes whose binary mergers might contribute
to the diffuse background of GWs may arise from spectral considerations Raidal
_et al._ (2017). Here, we rather focus on the use of the power spectrum of
anisotropies of the stochastic background as a possible tell-tale signature of
the production mechanism of the black holes (for previous studies see e.g.
Bartolo  _et al._ (2020) and references therein).

Previous studies that focused on utilizing the anisotropy of the stochastic
gravitational wave background (SGWB) as a metric to differentiate PBHs from
astrophysical black holes include, for instance, Ref. Wang  _et al._ (2022)
(see also Ref. Bavera  _et al._ (2022) on how to utilize the isotropic SGWB to
constrain the primordial origin of BH formation channel). Here, however, we
take a different perspective, and aim at differentiating one particular aspect
of PBH formation: the portion of the Hubble horizon that collapsed at
formation. Generically, PBHs result from the collapse of patches of the early
universe as large as the Hubble horizon at collapse, but depending on the
formation mechanism the collapsing region can be significantly smaller. We
therefore simply assume, here, in a model-agnostic way, that PBHs result from
the collapse of a fraction $f_{H}<1$ of the Hubble horizon only (we discuss
expectations for $f_{H}$ in different scenarios below), and study the impact
of $f_{H}$ on the anisotropy of the SGWB.

The remainder of this study is as follows: in the next section we discuss PBH
formation from the collapse of curvature perturbations in the standard and in
the presently investigated non-standard scenarios. Sec. III describes the
calculation of the power spectrum of anisotropies of the SGWB, and presents
our results. The final Sec. IV presents our discussion and conclusions.

##  II PBH formation from curvature perturbations

PBHs can form in the early universe as a result of the gravitational collapse
of curvature perturbations $\zeta(x)$. The power spectrum of curvature
perturbations, $P_{\zeta}(k)$, determines the width of the Gaussian
distribution of matter overdensities $\delta\rho(x)/\rho(x)\equiv\delta(x)$.
The distribution of matter overdensities then in turn determines the mass
function of the resulting PBHs. According to the standard scenario, where a
PBH is formed at the horizon re-entry time $t_{H}\sim 1/k_{*}$, with $k_{*}$
the scale of the perturbation that collapses to form the PBH, the mass of the
PBH is directly related, and proportional to, the Hubble mass – the mass
enclosed at that time in a Hubble patch. The later the scale $k_{*}$ reenters
the horizon, the larger the PBH mass. Notably, and generically, however, PBHs
can be also formed at the sub-horizon scales, and the formation time can be
much later than $t_{H}$, modifying the relation between PBH mass and horizon
mass. Here, we explore the implications of sub-horizon-forming PBHs for the
anisotropy of the gravitational wave spectrum induced by the PBH merger.

###  II.1 The standard scenario for PBH formation

The horizon reentry time is defined by, in natural units,

$t_{H}\equiv\frac{1}{k_{*}},$ (1)

where $k_{*}$ is the scale of the over-density that will collapse to form a
PBH, and where we denote the characteristic wavelength of the density
perturbation as $\lambda_{*}\equiv 1/k_{*}$. The size of a Hubble patch at
horizon re-entry is given by

$r_{H}=t_{H}.$ (2)

Thus, we can use Eq. (1) to relate the size of the Hubble patch at the horizon
reentry time and the characteristic wavelength of the density perturbation,

$\lambda_{*}=r_{H}.$ (3)

Statistically-rare perturbations (such that one can neglect hierarchical
random fields) are expected to be approximately spherically symmetric Bardeen
_et al._ (1986). With the assumption that the perturbation is spherically
symmetric, the amplitude of the smoothed density contrast $\bar{\delta}$ is
related to the curvature perturbation, in radiation domination, as

$\bar{\delta}=-\frac{2}{3}r_{m}\zeta^{\prime}(r_{m})(2+r_{m}\zeta^{\prime}(r_{m}))=\delta_{1}-\frac{3}{8}\delta_{1}^{2},$
(4)

where $r_{m}$ is the smoothing scale, the prime is the derivative with respect
to the radial coordinate, and $\delta_{1}$ is the linear component of the
$\bar{\delta}$, given by

$\delta_{1}=-\frac{4}{3}r_{m}\zeta^{\prime}(r_{m}).$ (5)

The distribution of $\delta_{1}$ obeys a Gaussian probability distribution
function (PDF):

${\rm
PDF}(\delta_{1})=\frac{1}{\sqrt{2\pi\sigma^{2}}}\exp\left(-\frac{\delta_{1}^{2}}{2\sigma_{0}^{2}}\right),$
(6)

where the variance $\sigma_{0}^{2}$ characterizes the width of the Gaussian
PDF, and it is obtained from the zeroth moment of the power spectrum of
$\delta_{1}$. The $j$-th moments of the power spectrum $P_{\delta_{1}}$ are
given by

$\sigma_{j}^{2}=\int_{0}^{\infty}\frac{dk}{k}{P}_{\delta_{1}}(k,r_{m})\left(\frac{k}{aH}\right)^{2j}.$
(7)

The variance is

$\sigma^{2}_{0}=\langle\delta_{1}^{2}\rangle=\int_{0}^{\infty}\frac{dk}{k}{P}_{\delta_{1}}(k,r_{m})=\frac{16}{81}\int_{0}^{\infty}\frac{dk}{k}(kr_{m})^{4}\tilde{W}^{2}(k,r_{m})T^{2}(k,r_{m})P_{\zeta}(k),$
(8)

where we choose a Dirac-$\delta$ peak power spectrum of the curvature
perturbation
$P_{\zeta}(k)=A_{s}\delta\left(\ln\frac{k}{k_{*}}\right)+A_{b}~{}(A_{s}\gg
A_{b})$, $A_{s}$ and $A_{b}$ are the amplitudes of the short-wavelength mode
and background mode, respectively, $r_{m}$ is the smoothing scale that
determines how to choose the spatial volume to average the overdensity peak,

$\tilde{W}(k,r_{m})=3\frac{\sin(kr_{m})-kr_{m}\cos(kr_{m})}{(kr_{m})^{3}}$ (9)

is the Fourier transform of the top-hat window function with the time-
dependent smoothing scale $r_{m}$, and

$T(k,r_{m})=3\frac{\sin(kr_{m}/\sqrt{3})-kr_{m}/\sqrt{3}\cos(kr_{m}/\sqrt{3})}{(kr_{m}/\sqrt{3})^{3}}$
(10)

is the linear transfer function. In the literature Young  _et al._ (2019), it
is common to choose the smoothing scale to be the Hubble scale, i.e.,
$r_{m}(z)=r_{H}=(aH)^{-1}$.

When the smoothed density contrast $\bar{\delta}$ is greater than the critical
threshold value $\delta_{c}$, a PBH is formed, with a mass closely related to
the modeling of the gravitational collapse,

$M_{\rm PBH}=\mathcal{K}(\bar{\delta}-\delta_{c})^{\gamma}M_{H},$ (11)

where $M_{H}$ is the horizon mass at the formation time of PBH,
$\mathcal{K}=4$ accounts for most of the shape of the density contract,
$\gamma=0.36$ is the critical exponent dependent on the equation of state of
the universe and here is evaluated in the radiation-dominated era Young  _et
al._ (2019), and $\delta_{c}=0.51$ is chosen for a typical profile shape of
the density perturbation and a real-space top-hat window function at the
horizon reentry time Musco (2019); Young (2019). Given the above PBH formation
model, the distribution of the Gaussian random field $\delta_{1}$ thus
determines the abundance and mass function of the resulting PBHs. By inverting
Eq. (11) and combining with Eq. (4), one can obtain the value of $\delta_{1}$
corresponding to a given PBH mass $M_{\rm PBH}$,

$\delta_{1}(M_{\rm
PBH})=\frac{2}{3}\left(2-\sqrt{4-6\delta_{c}-6\left(\frac{M_{\rm
PBH}}{\mathcal{K}M_{H}}\right)^{1/\gamma}}\right).$ (12)

The initial abundance of PBHs is usually described by the energy density
fraction of PBHs in the universe at a single formation time,

$\beta=\int_{\mu_{c}}^{\mu_{\rm max}}d\mu\frac{M_{\rm PBH}(\mu)}{M_{H}}n_{\rm
peak}(\mu),$ (13)

where $\mu\equiv\delta_{1}/\sigma_{0}$, $\mu\in[\mu_{c},\mu_{\rm max}]$,
$\mu_{c}=\delta_{c,1}/\sigma_{0}$, $\mu_{\rm max}=\delta_{1,{\rm
max}}/\sigma_{0}=4/(3\sigma_{0})$, makes the gravitational collapse happen,
$\delta_{c,1}=\frac{4}{3}\left(1-\sqrt{1-\frac{3}{2}\delta_{c}}\right)$ is the
critical value for gravitational collapse of $\delta_{1}$, $n_{\rm
peak}=\frac{1}{4\pi^{2}}\left(\frac{\sigma_{1}}{\sigma_{0}}\right)^{3}\mu^{3}\exp\left(-\frac{\mu^{2}}{2}\right)$
is the number density of over-density peaks in a comoving volume Bardeen  _et
al._ (1986).

The mass function is defined as

$\Psi(M_{\rm PBH})=\frac{1}{f_{\rm PBH}}\frac{df_{\rm PBH}}{dM_{\rm PBH}},$
(14)

such that $\int\Psi(M_{\rm PBH})dM_{\rm PBH}=1$, where $f_{\rm
PBH}\equiv\Omega_{\rm PBH}/\Omega_{\rm DM}$ is the fraction of the present
dark matter (DM) energy density stored in PBHs. The fraction $f_{\rm PBH}$ can
be calculated from $\beta$ by integrating the contributions from PBHs formed
across the history of the universe until the time of matter-radiation
equality,

$f_{\rm PBH}=\frac{1}{\Omega_{\rm DM}}\int_{M_{H,{\rm min}}}^{M_{H,{\rm
max}}}d(\ln M_{H})\left(\frac{M_{H,{\rm
eq}}}{M_{H}}\right)^{1/2}\beta(M_{H}),$ (15)

where $M_{H,{\rm eq}}=\frac{4\pi}{3}2\rho_{\rm eq}H_{\rm eq}^{-3}\approx
2.8\times 10^{17}M_{\odot}$ is the approximate horizon mass at the matter-
radiation equality. Assuming a Dirac-$\delta$ form for the power spectrum of
the curvature perturbation and assuming as a result that all PBHs form at a
single time, the mass function is given by

$\Psi(M_{\rm PBH})=\frac{1}{\beta}\frac{d\beta}{dM_{\rm PBH}}.$ (16)

Here is the derivation. First, we relate $\beta$, the abundance of PBHs at a
single formation time, to $\Omega_{\rm PBH}$, the energy density parameter of
PBHs observed today:

$\beta=\frac{\rho_{\rm PBH}}{\rho_{r}}=\frac{\rho_{\rm PBH}}{\rho_{\rm
c,0}}\frac{\rho_{\rm c,0}}{\rho_{r}}=\frac{\rho_{\rm
PBH,0}\left(\frac{a_{0}}{a}\right)^{3}}{\rho_{\rm c,0}}\frac{\rho_{\rm
c,0}}{\rho_{r,0}\left(\frac{a_{0}}{a}\right)^{4}\left(\frac{g_{0}}{g}\right)^{1/3}}=\frac{\Omega_{\rm
PBH}}{\Omega_{r}}\left(\frac{a}{a_{0}}\right)\left(\frac{g}{g_{0}}\right)^{1/3},$
(17)

where $\Omega_{r}$ is the energy density parameter of radiation observed
today, $a$ is the scalar factor, $a_{0}=1$ is the scalar factor today, $g$ and
$g_{0}$ are the degree of freedom of the relativistic species at the time of
$a$ and $a_{0}$ respectively. Therefore, the energy density fraction of PBHs
in terms of DM is given by a simplified expression instead of the integral in
Eq. (15),

$f_{\rm PBH}=\frac{\Omega_{\rm PBH}}{\Omega_{\rm
DM}}=\frac{\Omega_{r}}{\Omega_{\rm
DM}}\frac{1}{a}\left(\frac{g_{0}}{g}\right)^{1/3}\beta.$ (18)

Let $\psi(x)=\Psi\times M_{H}$ be a dimensionless mass function:

$\psi(x)=\frac{1}{4\pi^{2}\beta\sigma_{0}}x\left(\frac{\sigma_{1}}{\sigma_{0}}\right)^{3}\left(\frac{\delta_{1}(x)}{\sigma_{0}}\right)^{3}\exp\left(-\frac{\delta_{1}^{2}(x)}{2\sigma_{0}^{2}}\right)\frac{\sqrt{2}(x/\mathcal{K})^{1/\gamma-1}}{\gamma\mathcal{K}\sqrt{2-3\delta_{c}-3(x/\mathcal{K})^{1/\gamma}}},$
(19)

where $x\equiv\frac{M_{\rm PBH}}{M_{H}}$ is the ratio between PBH mass and
horizon mass. The normalization of the dimensionless mass function becomes
$\int\psi(x)dx=1$, and the PBH-mass-to-horizon-mass ratio $x$ ranges from 0 to
2.05, where the upper limit is coming from the that Eq. (4) has a maximum
value $\bar{\delta}=\frac{2}{3}$ such that

$M_{\rm
PBH,max}=\mathcal{K}\left(\frac{2}{3}-\delta_{c}\right)^{\gamma}M_{H}=2.05M_{H},$
(20)

where the second identity is obtained by substituting the parameters of the
PBH formation model used in this paper.

###  II.2 Sub-horizon PBH formation

Generically, a PBH can form from the collapse of a density perturbation of
$\lambda_{*}$ smaller than the size of the horizon by a factor $f_{H}$,

$\lambda_{*}=f_{H}\times r_{H}^{\prime},$ (21)

where we use the symbol $r_{H}^{\prime}$ to indicate the Hubble patch where
the PBH is formed. For example, Ref. Nakama  _et al._ (2017) considers
$f_{H}=1/4$. Several example formation mechanisms with $f_{H}\lesssim 1$
and/or with $f_{H}\ll 1$ have been proposed and discussed in detail in the
literature. For instance, Ref. Cotner  _et al._ (2019) considers PBH
production from the collapse of scalar field lumps (Q-balls or oscillons) that
can potentially temporarily dominate the energy density of the universe
leading to efficient PBH formation; the typical size of the lumps (formed,
here, not because of gravitational instability, but rather due to self-
interactions) is predicted, from numerical simulations, to be a fraction
$f_{Q}\sim 10^{-2}-10^{-1}$ of the horizon size. Bubble collision upon
spontaneous symmetry breaking and a first-order phase transition also
generically seed sub-horizon density perturbations leading to collapse into
black holes with mass much smaller than the horizon mass (see e.g. Hawking
_et al._ (1982)). Similarly, cosmic string loops lead to sub-horizon mass PBH,
depending on the string tension Garriga and Vilenkin (1993); Caldwell and
Casper (1996). Ref. Meissner and Nicolai (2020) entertains a more speculative
possibility of multigravitino Bogomol’nyi-Prasad-Sommerfield (BPS) states that
if sufficiently large could collapse into black holes whose mass is entirely
unrelated to the horizon mass. A somewhat less speculative possibility is the
collapse of Fermi balls consisting of massive fermions trapped into regions of
false vacuum in a first-order phase transition Kawana and Xie (2022). Other
similar possibilities include the formation of PBH from long-range scalar
forces and scalar radiative cooling Flores and Kusenko (2021), and a late-
phase transition in a strongly interacting fermion-scalar fluid Chakraborty
_et al._ (2022). Finally, the collapse of neutron stars in the very late
universe, possibly triggered by the accumulation of dark matter in the neutron
star core, can also lead to the formation of black holes with masses much
smaller than the horizon mass Acevedo  _et al._ (2021); Giffin  _et al._
(2022).

For $f_{H}<1$, one has $r_{H}^{\prime}=r_{H}/f_{H}>r_{H}$, i.e., the Hubble
patch where the PBH is formed is greater than the one in the standard
scenario, and thus the formation time is later than the one in the standard
scenario, $t_{H}^{\prime}>t_{H}$. Assuming PBHs formed in the early universe,
during radiation domination, the formation redshift can be related to the one
in the standard scenario by the following equation:

$\frac{z^{\prime}_{\text{form}}}{z_{\text{form}}}=f_{H}.$ (22)

Figure 1:  The variance $\sigma^{2}_{0}$ of the over-density distribution as a
function of $1/f_{H}$, for different amplitudes of the short-wavelength
perturbation mode $A_{s}=0.8\,(\text{blue})$, $0.08\,(\text{orange})$,
$0.008\,(\text{green})$, at a fixed value for the amplitude of the background
mode $A_{b}=8\times 10^{-4}$.

In the sub-horizon formation scenario, $r_{m}(z)$ in Eq. (8) is evaluated at
time $z^{\prime}_{\text{form}}$ instead of $z_{\text{form}}$, such that
$k_{*}r_{m}(=f_{H}^{-1})\neq 1$. There is one subtlety for the smoothing scale
of the top-hat window function. In the standard case, there is only one
density contrast peak in one Hubble patch, so choosing $r_{m}=r_{H}$ is
reasonable. However, in the sub-horizon formation case, there are $\lfloor
f_{H}^{-1}\rfloor$ density contrast peaks, but we still choose $r_{m}=r_{H}$
as the smoothing scale. We comment later on the impact of choosing different
smoothing radii and how our results depend on the choice of $r_{m}$.

Figure 1 shows $\sigma^{2}_{0}$ as a function of the inverse of the horizon
size fraction $f_{H}^{-1}$ for different values of $A_{s}$ with a fixed value
of $A_{b}=8\times 10^{-4}$. For a given $A_{s}$, $\sigma_{0}^{2}$ increases
monotonically before reaching its maximum at $f_{H}^{-1}\sim 2.5$, and it
rapidly drops until it reaches its minimum at $f_{H}^{-1}\sim 4.5$. The
behavior of $\sigma_{0}^{2}$ is quasi-periodic with respect to $f_{H}^{-1}$
because the window function and transfer function in Eq. (8), containing
$\sin(k_{*}r_{m})=\sin(f_{H}^{-1})$ and $\cos(k_{*}r_{m})=\cos(f_{H}^{-1})$,
are periodic but with a suppressed amplitude at larger $f_{H}^{-1}$. Figure 1
additionally illustrates how the value of $\sigma_{0}^{2}$ decreases globally
for a smaller value of $A_{s}$, because $\sigma_{0}^{2}$ characterizes the
auto-correlation of the overdensity whose amplitude is determined by the
amplitude of the curvature perturbation. We stress that the non-linear
dependence of the value of $\sigma_{0}^{2}$ on $A_{s}$ at some values of
$f_{H}$ stems from the fact that the ratio $A_{s}/A_{b}$ is not a constant. If
$A_{s}/A_{b}$ were a constant, $\sigma_{0}^{2}$ would linearly depend on
$A_{s}$, which can be directly derived from Eq. (8). Starting from the Sec.
III, for simplicity, we choose $A_{b}=0$, and the $j$-th order moments are
simplified to be

$\sigma_{j}^{2}=432A_{s}f_{H}^{8-2j}\left(\sin\left(f_{H}^{-1}\right)-f_{H}^{-1}\cos\left(f_{H}^{-1}\right)\right)^{2}\left(\sin\left(f_{H}^{-1}/\sqrt{3}\right)-f_{H}^{-1}/\sqrt{3}\cos\left(f_{H}^{-1}/\sqrt{3}\right)\right)^{2}.$
(23)

We note that the above dependence on $A_{s}$ and $f_{H}$ brings a degeneracy
of $A_{s}$ when fixing $f_{\rm PBH}$ and $f_{H}$. This is because $f_{\rm
PBH}$ is proportional to $\beta$, which, in turn, is only determined by
$\sigma_{0}^{2}$ and $\sigma_{1}^{2}$ (see Eq. (13)); specifically, $\beta$
has an exponential dependence on $\sigma_{0}^{2}$, and a power-law dependence
on $\sigma_{1}^{2}$. This results in $f_{\rm PBH}$ being approximately fixed
by a given choice of $\sigma_{0}^{2}$. However, Eq. (23) can have multiple
roots for $A_{s}$ when fixing $\sigma_{0}^{2}$ and $f_{H}$, which can also be
appreciated by imagining a horizontal line in Figure 1 and noting how it would
intersect with curves of different colors from a large $\sigma_{0}^{2}$ value
to a smaller value. Therefore, $A_{s}$ has multiple values for given $f_{\rm
PBH}$ and $f_{H}$. Also, due to the quasi-periodic behavior of
$\sigma_{0}^{2}$ versus $f_{H}$, one has the freedom to pick $A_{s}$ such that
$f_{\rm PBH}$ is sufficiently small at a given $f_{H}$. In the calculation of
Sec. III, we fix $f_{\rm PBH}=0.01$ and choose the smallest value of $A_{s}$
for each $f_{H}$. Eq. (23) also manifests choosing different smoothing scales
$r_{m}$ will change the value of $\sigma_{0}^{2}$ and $\sigma_{1}^{2}$. For
the choice of $r_{m}=r_{H}$, $k_{*}r_{m}=f_{H}^{-1}$, while in general, one
can choose $f_{H}r_{H}\leq r_{m}\leq r_{H}$, i.e., $1\leq k_{*}r_{m}\leq
f_{H}^{-1}$.

Figure 2:  The mass function of PBH for different values of $f_{H}$, for
$A_{s}=0.8$ (top panels), 0.08 (middle panels) and 0.008 (bottom panels) with
$A_{b}=8\times 10^{-4}$ in all cases. Left column: linear-linear scale. Right
column: log-log scale.

Figure 2 shows (on a linear (left) and logarithmic (right) scale) the mass
function of Eq. (19) for different values of the variance $\sigma^{2}$ induced
by $A_{s}=0.8$, $0.08$, and $0.008$, again at fixed $A_{b}=8\times 10^{-4}$.
The log-log plots (Figure 2 right panels) illustrate how the mass function
exhibits a power-law behavior when the normalized PBH mass $x<x_{0}$, where
$x_{0}$ is the point at which the mass function reaches its local maximum.
However, we note that for $A_{s}=0.8$, $f_{H}=1/2~{}\rm{and}~{}1/3$, the mass
function does not have a local maximum, but, rather, it monotonically
increases as a function of the PBH mass, while for the other it is
exponentially suppressed when $x>x_{0}$. This behavior can be explained as
follows: the exponential function ${\rm
exp}\left(-\frac{\delta_{1}^{2}(x)}{2\sigma_{0}^{2}}\right)$ in Eq. (19)
decreases monotonically as $x$ increases for a given $\sigma_{0}^{2}$ but with
a different decreasing rate which depends sensitively on the value of
$\sigma_{0}^{2}$. For $f_{H}^{-1}=2$ or 3 when $A_{s}=0.8$, the value of
$\sigma_{0}^{2}\sim 0.8$, which is almost an order of magnitude larger than
$\sigma_{0}^{2}$ at other $f_{H}^{-1}$, such that the decrease in ${\rm
exp}\left(-\frac{\delta_{1}^{2}(x)}{2\sigma_{0}^{2}}\right)$ is counter-
balanced by the increase in the remaining $x$-dependent part in Eq. (19). For
other $f_{H}^{-1}$ values, since the corresponding value of $\sigma_{0}^{2}$
is small ($\sim 0.1$), the behavior of the mass function at large $x$ is
instead mainly determined by the exponential suppression originated from ${\rm
exp}\left(-\frac{\delta_{1}^{2}(x)}{2\sigma_{0}^{2}}\right)$.

Figure 3:  The mass function (left) and the one normalized by the mass
function of $f_{H}=1$ (right) for $f_{H}\leq 0.1$, $A_{s}=0.8$, and
$A_{b}=8\times 10^{-4}$.

We also investigate the mass function at a very small $f_{H}$, i.e., very
large $f_{H}^{-1}$, ranging from $f_{H}^{-1}\in[10,5000]$ in Figure 3 which
shows the mass function at $f_{H}\ll 0.1$ in the left panel and the ratio of
mass functions normalized by the “standard” mass function corresponding to
$f_{H}=1$ in the right panel. We note that for $f_{H}\leq 1/30$, the mass
functions for different values of $f_{H}$ are self-similar. The reason is that
the mass function $\psi(x)$ only depends on the value of
$\sigma_{0}^{2}=\sigma_{0}^{2}(f_{H})$ and when $f_{H}\leq 1/30$,
$\sigma_{0}^{2}\approx 0.00085\sim A_{b}$ does not significantly change.

##  III Anisotropic SGWB energy density power spectrum

If two PBHs, separated by a distance $x$, have an energy density $4M/(3\pi
x^{3})$ larger than the cosmic radiation energy density $\rho_{r}$ (since we
suppose PBHs are formed deep in the radiation era, radiation is the dominant
component of the universe), they decouple from the Hubble flow and become a
gravitationally bound binary system. PBH binaries then radiate GWs during the
inspiral, merger, and ring-down phases Abadie  _et al._ (2011). The incoherent
superposition of the unresolved GWs emitted by such PBH binaries contributes
to SGWBs across several orders of magnitude in frequency, which advanced LIGO
or LISA can detect in the near future Moore  _et al._ (2015). The SGWB
produced by PBH binaries has both an isotropic and an anisotropic component:
while only spectral information can be extracted from the former, we focus on
the latter, the anisotropic component in seeking a way to ascertain whether
the PBH formation history corresponds to a standard horizon-size collapse, or
to the aforementioned sub-horizon-collapse scenarios. As such, the scope of
the present study is to seek for distinguishable features of PBH binary
mergers in the standard versus sub-horizon formation scenarios by studying the
anisotropic component.

As customary Allen and Romano (1999), we define as the isotropic SGWB energy
density spectrum the dimensionless quantity

$\bar{\Omega}_{\rm gw}(\nu)\equiv\frac{\nu}{\rho_{c}}\frac{d\rho_{\rm
gw}}{d\nu},$ (24)

where $\rho_{c}$ is the critical energy density of the universe and $\rho_{\rm
gw}$ is the SGWB energy density. In a more general situation where the SGWB
has directional dependence, the definition is supplemented by

$\Omega_{\rm gw}(\nu,{\bf n})\equiv\frac{1}{\rho_{c}}\frac{d^{3}\rho_{\rm
gw}(\nu,{\bf n})}{d\ln\nu d^{2}{\bf n}}=\frac{\bar{\Omega}_{\rm
gw}}{4\pi}(1+\delta(\nu,{\bf n})),$ (25)

where $\rho(\nu,{\bf n})$ is the SGWB energy density at frequency $\nu$ along
the line-of-sight direction $\bf n$, $\delta(\nu,{\bf n})$ denotes the
anisotropic fluctuations.

Figure 4:  The isotropic component of the SGWB energy spectrum at
$\nu=10^{-4}{\rm Hz}$ with different $M_{\rm
PBH,max}=M_{H}=10^{-6}M_{\odot}/f_{H}^{2}$, $10^{-4}M_{\odot}/f_{H}^{2}$,
$10^{-2}M_{\odot}/f_{H}^{2}$, and $1M_{\odot}/f_{H}^{2}$ (left), and the
stadard case at $M_{\rm PBH,max}=1M_{\odot}/f_{H}^{2}$, for different observed
frequencies $\nu=10^{-4}{\rm Hz}$, $10^{-2}{\rm Hz}$, $1{\rm Hz}$, and
$100{\rm Hz}$ (right). We fix $f_{\rm PBH}=0.05$.

The isotropic component $\bar{\Omega}_{\rm gw}(\nu)$ produced by PBH binary
mergers is given by the following integral:

$\bar{\Omega}_{\rm
gw}(\nu)=\frac{\nu}{\rho_{c}}\int_{0}^{z_{*}}\frac{dz}{(1+z)H(z)}\int
d{\Theta}_{s}\mathcal{R}(\Theta_{s},z)\frac{dE}{d\nu}(\nu_{s},\Theta_{s}),$
(26)

where $z$ is redshift, $z_{*}$ corresponds to the formation time of PBHs,
$H(z)=H_{0}\sqrt{\Omega_{r}(1+z)^{4}+\Omega_{m}(1+z)^{3}+\Omega_{\Lambda}}$ is
the Hubble parameter, $\mathcal{R}$ is the merger rate of PBH binaries,
$dE/d\nu(\nu_{s},\Theta)$ is the energy spectrum at source frequency
$\nu_{s}=\nu\times(1+z)$ ,with source parameters $\Theta_{s}$ describing
masses, eccentricity, and so on. The energy spectrum used in our calculation
is derived from the inspiral-merger-ringdown waveforms for non-spinning PBH
binaries given e.g. in Ref. Ajith  _et al._ (2008, 2011). Figure 4 shows the
isotropic SGWB energy density $\bar{\Omega}_{\rm gw}(\nu;f_{H})$ as a function
of $f_{H}^{-1}$ at different observed frequencies and different $M_{\rm
PBH,max}$, in which we use the merger rate $\mathcal{R}$ (see Eq. (31)) that
is to be introduced later to take into account the effect of the sub-horizon
formation.

The cross-correlation of the anisotropies in SGWB between two different
directions ${\bf n}$ and ${\bf n^{\prime}}$, $\langle\delta(\nu,{\bf
n})\delta(\nu,{\bf n^{\prime}})\rangle$, describes the amplitude of the
statistical fluctuation at a certain angular scale ${\bf n\cdot n^{\prime}}$.
By decomposing it in a spherical harmonic basis, we have

$\langle\delta(\nu,{\bf n})\delta(\nu,{\bf
n^{\prime}})\rangle=\sum_{l}\frac{2l+1}{2\pi}C_{l}P_{l}({\bf n}\cdot{\bf
n^{\prime}}),$ (27)

and the dimensionless angular power spectrum is given by

$C_{l}(\nu)\equiv\frac{2}{\pi}\int d\ln kk^{3}|\delta_{l}(k,\nu)|^{2},$ (28)

where $\delta_{l}(k,\nu)$ is the Fourier transforms of the spatial
anisotropies $\delta(\nu,{\bf n})$, and are given by Cusin  _et al._ (2017),

$\displaystyle\delta_{l}(\nu,k)$ $\displaystyle=$
$\displaystyle\frac{\nu}{\bar{\Omega}_{\rm
gw}(\nu)\rho_{c}}\int_{\eta_{*}}^{\eta_{0}}d\eta\mathcal{L}(\eta,\nu_{s})$
(29)
$\displaystyle\times\Bigg{\\{}\left[4\Psi_{k}(\eta)+4\Pi_{k}(\eta)+b\delta_{m,k}(\eta)\right]j_{l}(k\Delta\eta)-2kv_{k}(\eta)j_{l}^{\prime}(k\Delta\eta)-6\int_{\eta_{0}}^{\eta}d\eta^{\prime}\dot{\Psi}_{k}(\eta^{\prime})j_{l}(k\Delta\eta^{\prime})$
$\displaystyle+\frac{\partial\ln(dE/d\nu_{s})}{\partial\ln\nu_{s}}\left[-\Psi_{k}(\eta)j_{l}(k\Delta\eta)-\Pi_{k}(\eta)j_{l}(k\Delta\eta)+kv_{k}(\eta)j_{l}^{\prime}(k\Delta\eta)+2\int_{\eta_{0}}^{\eta}d\eta^{\prime}{\Psi}_{k}(\eta^{\prime})j_{l}(k\Delta\eta^{\prime})\right]\Bigg{\\}},$

where $\mathcal{L}(\eta,\nu_{s})=a(\eta)\int
d\Theta_{s}\mathcal{R}(\Theta_{s},z(\eta))\frac{dE}{d\nu}(\nu_{s},\Theta_{s})$
is the GW luminosity of the PBH binary merger, the galaxy bias $b=1$ for PBH,
$\Psi_{k}\approx\phi_{k}$ and $\Pi_{k}=0$ are two scalar gravitational
potentials, $\delta_{m,k}$ is the matter density contrast, $v_{k}$ is the
velocity of matter, the dot represents the derivative with respect to
conformal time $\eta$. The corresponding transfer functions are calculated
numerically by solving the linearized cosmological perturbation equations in
CLASS, and the primordial scalar fluctuation power spectrum is
$P_{s}(k)=\mathcal{A}_{s}^{2}k^{-3}(k/k_{0})^{n_{s}-1}$, where
$\mathcal{A}_{s}=2.215\times 10^{-9}$, $n_{s}=0.96$, and $k_{0}=0.05{\rm
Mpc}^{-1}$ are chosen from the Planck 2018 results Aghanim  _et al._ (2020),
and $j_{l}$ is the spherical Bessel function.

By comparing the numerical results of $\delta_{l}$ in Eq. (29) with or without
including the $v_{k}$’s and $\phi_{k}$’s terms, we find that only the term
with the matter density contrast $\delta_{m}$ is dominant in the calculation,
so it can be approximated by

$\delta_{l}(\nu,k)=\frac{\nu}{\bar{\Omega}_{\rm
gw}(\nu)\rho_{c}}\int_{\eta_{*}}^{\eta_{0}}d\eta\mathcal{L}(\eta,\nu_{s})\delta_{m,k}(\eta)j_{l}(k\Delta\eta).$
(30)

We assume that the merger rate of PBH is given by the functional form Raidal
_et al._ (2019),

$\mathcal{R}(\Theta_{s},\tau)=1.6\times
10^{6}\text{Gpc}^{-3}\text{yr}^{-1}~{}f_{\rm
PBH}^{53/37}\eta(m_{1},m_{2})^{-34/37}\left(\frac{m_{1}+m_{2}}{M_{\odot}}\right)^{-32/37}\left(\frac{\tau}{t_{0}}\right)^{-34/37}\Psi(m_{1})\Psi(m_{2})S(f_{\rm
PBH}),$ (31)

where $f_{\rm PBH}$ (not to be confused with $f_{H}$) is the fraction of PBH
at the present time $t_{0}$, $\eta(m_{1},m_{2})=m_{1}m_{2}/(m_{1}+m_{2})^{2}$
is the symmetric factor for the binary, $\Psi=\Psi(m;f_{H})$ is the mass
function which is defined in the previous section, and $S(f_{\rm
PBH})=\left(\frac{5f_{\rm
PBH}^{2}}{6\sigma_{M}^{2}}\right)^{21/74}U\left(\frac{21}{74},\frac{1}{2},\frac{5f_{\rm
PBH}^{2}}{6\sigma_{M}^{2}}\right)$ is the suppression factor with $U$ the
confluent hypergeometric function. Notice that we assume that the form of the
merger rate is unchanged in the sub-horizon formation scenario, i.e.,
$\Psi(m;f_{H})$ can be evaluated when $f_{H}<1$, which is strongly physically
motivated.

Figure 5:  The angular power spectrum at multiple frequencies at a fixed PBH
mass $M_{\rm PBH}=1M_{\odot}$ (left). The relative amplitude of angular power
spectrum at $l=1$ as a function of frequency $\nu$ for three different PBH
masses $M_{\rm PBH}=0.1,1,10~{}M_{\odot}$ (right).

Figure 5, left, shows the typical angular power spectrum of the SGWB
anisotropic component for PBH binary mergers in the standard scenario. We
choose a Dirac-$\delta$ mass function in the merger rate such that $M_{\rm
PBH}=M_{H}=1M_{\odot}$, and thus PBHs are formed deep in the radiation-
dominated era with $z_{\rm form}=1.62\times 10^{12}$. The right panel shows
the frequency dependence of the $l=1$ dipole as a function of frequency for
various choices of the PBH mass $M_{\rm PBH}=0.1,1,10~{}M_{\odot}$,
illustrating how larger chirp masses, $M_{c}=2^{-1/5}M_{\rm PBH}$ for a binary
PBH with the same mass, induce a lower dipole amplitude at low frequency and a
larger amplitude at higher frequencies, as expected.

Figure 6:  The angular power spectrum $C_{l}$ at 100 Hz with PBH masses given
by a mass function ($f_{H}=1$ and $\sigma_{0}^{2}=0.0101$ such that $f_{\rm
PBH}=0.01$) and by a fixed PBH mass $M_{\rm PBH}=1M_{\odot}$ (left). The ratio
of $C_{l}$ between the mass function cases and the single mass case (right).

In Figure 6 we compare the angular power spectrum obtained before to that
corresponding to a mass function arising from a value of $\sigma_{0}^{2}$ such
that $f_{\rm PBH}=0.01$ at a frequency $\nu=100{\rm Hz}$ (the right plot shows
the ratio of the curves at a given $\sigma_{0}^{2}$ to the standard case).
There are no large deviations between the angular power spectra besides an
overall amplitude shift due to the different weighting on the PBH mass
function. The single mass case has a Dirac-$\delta$ peak mass function, while
the other has a mass function determined by $\sigma_{0}^{2}$.

Figure 7:  The angular power spectrum at 100 Hz for different $f_{H}$. $f_{H}^{-1}$ | 1 | 1.2 | 1.4 | 1.6 | 1.8 | 2 | 3 | 4 | 5 | 6 | 7
---|---|---|---|---|---|---|---|---|---|---|---
$\bar{\Omega}_{\rm gw}/10^{-9}$ | 1.26 | 1.61 | 1.98 | 2.36 | 2.73 | 3.11 | 4.93 | 6.49 | 7.77 | 8.63 | 9.04
$\sigma_{0}^{2}$ | 0.0101 | 0.00988 | 0.00972 | 0.00959 | 0.00948 | 0.00938 | 0.00903 | 0.00882 | 0.00866 | 0.00854 | 0.00844
$\sigma_{1}^{2}$ | 0.0101 | 0.0142 | 0.0191 | 0.0245 | 0.0307 | 0.0375 | 0.0813 | 0.141 | 0.217 | 0.307 | 0.413
$A_{s}$ | 0.0668 | 0.0357 | 0.0219 | 0.0150 | 0.0112 | 0.00913 | 0.00885 | 0.0736 | 0.149 | 0.103 | 1.76
Table 1: Values of the isotropic SGWB energy density $\bar{\Omega}_{\rm gw}$
at 100 Hz, the zeroth order and first order width, and the amplitude of
Dirac-$\delta$ curvature power spectrum for a fixed $f_{\rm PBH}$ and
different $f_{H}^{-1}$.

To study the effects of $f_{H}\neq 1$, we study the SGWB from an
$f_{H}$-dependent mass function $\Psi(m;f_{H})$ to investigate the effect of
sub-horizon formation on the angular power spectrum. Figure 7 shows the
angular power spectrum at 100 Hz for $f_{H}=m^{-1}$, $1\leq m\leq 7$ and
$m\in\mathbb{N}$. The amplitude of the angular power spectrum increases as
$f_{H}^{-1}$ increases when $f_{H}^{-1}\in[1,6]$, but slightly decreases when
$f_{H}^{-1}\in(6,7]$. Table 1 summarizes the isotropic SGWB energy density
values, $\sigma_{0}^{2}$, $\sigma_{1}^{2}$, and $A_{s}$ required to match a
given $f_{\rm PBH}=0.01$ for different $f_{H}$ (including some fractional
values of $m$).

Figure 8:  The angular power spectra of $l=1,~{}251,~{}501,~{}751,~{}{\rm
and}~{}1001$ as a function of $f_{H}^{-1}$, normalized by $C_{l}(f_{H}=1)$ for
the corresponding $l$, for {1$M_{\odot}$,100Hz} (left) and
{$10^{-6}M_{\odot}$,$10^{-4}$Hz} (right).

We are also interested in the effect of sub-horizon formation on the angular
power spectrum $C_{l}$ at different multipole $l^{\prime}$s. We show the
normalized $C_{l}$ (i.e. the ratio $C_{l}(f_{H})/C_{l}(f_{H}=1)$) in Figure 8,
for different parameter sets $\\{M_{H}f_{H}^{2},\nu\\}$, in which the first
parameter sets the horizon mass when PBHs form for different $f_{H}$ and the
second parameter is the observed frequency of the angular power spectrum. The
left panel employs {1$M_{\odot}$,100Hz} while the right panel
{$10^{-6}M_{\odot}$,$10^{-4}$Hz}. We also sampled {$10^{-2}M_{\odot}$,1Hz} and
{$10^{-4}M_{\odot}$,$10^{-2}$Hz}, and found almost identical results as for
{$10^{-6}M_{\odot}$,$10^{-4}$Hz}.

The choice of the parameter sets makes the observed frequency below the merger
frequency of the IMR waveform at a redshift $z=0$. We do not find a
significant difference for $C_{l}$ at different multipole $l$ when choosing
different $f_{H}$. The behavior of the $C_{l}$ ratio at 100Hz is different
from the other cases we studied ($\nu=10^{-2}$Hz, $10^{-4}$Hz, and
$10^{-6}$Hz), which is, in turn, similar to the behavior of $C_{l}$ shown in
Figure 5. There exists a turn-around frequency $\nu_{\rm ta}$ for a given
$f_{H}$ and $M_{\rm PBH}$ or effectively $M_{c}$ such that $C_{l}$ decreases
as $\nu$ increases when $\nu<\nu_{\rm ta}$, while $C_{l}$ increases when
$\nu\geq\nu_{\rm ta}$. Similarly, for a fixed observed frequency $\nu$ and
multiple $f_{H}$ and $M_{H}$ or equivalently $M_{\rm PBH,max}=2.05M_{H}$, we
can define a minimal turn-around frequency $\nu_{\rm ta,min}=\min\\{\nu_{\rm
ta}(M_{H,i}f_{H,i}^{2})\\}$ and a maximal turn-around frequency $\nu_{\rm
ta,max}=\max\\{\nu_{\rm ta}(M_{H,i}f_{H,i}^{2})\\}$. If $\nu<\nu_{\rm
ta,min}$, $C_{l}$ will generally decrease as $M_{H}$ increases. If
$\nu>\nu_{\rm ta,max}$, $C_{l}$ will generally increase as $M_{H}$ increases.

Figure 9:  The angular power spectra for $l=1$ as a function of $f_{H}^{-1}$
normalized by $C_{l}(\nu;f_{H}=1)$ for different frequencies. The left panel
is using $M_{H}=1M_{\odot}/f_{H}^{2}$, while the right panel is using
$M_{H}=10^{-6}M_{\odot}/f_{H}^{2}$.

The ratio $C_{l}(f_{H})/C_{l}(f_{H}=1)$ is an $f_{H}$-dependent and frequency-
dependent object, as shown in Figure 9. To analytically understand this
behavior, we first rewrite $j_{l}(x)$ using the Limber approximation,

$j_{l}(x)\approx\sqrt{\frac{\pi}{2l+1}}\delta(x-(l+0.5)),$ (32)

where $\delta(x)$ is the Dirac $\delta$ function. Applying the Limber
approximation to Eq. 30 by replacing $x=k\Delta\eta=k(\eta-\eta_{0})$, we can
integrate the conformal time integral by using the Dirac $\delta$-function,
which gives

$\delta_{l}(\nu,k)=\sqrt{\frac{\pi}{2l+1}}\frac{\nu}{\bar{\Omega}_{\rm
gw}(\nu)\rho_{c}}\frac{\mathcal{L}(\eta_{*},\nu_{s})}{1+z_{*}}\frac{\delta_{m,k}(\eta_{*})}{k},$
(33)

where $\eta_{*}=\eta_{0}-\frac{l+0.5}{k}$ is the conformal time determined by
different sets of $\\{l,k\\}$ such that the argument of the Dirac
$\delta$-function in Eq. (32) is zero, and $z_{*}$ is the redshift evaluated
at $\eta_{*}$. The angular power spectrum is thus simplified in the following,

$C_{l}(\nu)=\frac{2}{2l+1}\left(\frac{\nu}{\bar{\Omega}_{\rm
gw}(\nu)\rho_{c}}\right)^{2}\int
dk\frac{\mathcal{L}^{2}(\eta_{*},\nu_{s})}{(1+z_{*})^{2}}\delta_{m,k}^{2}(\eta_{*}).$
(34)

For a fixed $l$ and a fixed observed frequency $\nu$ but different $f_{H}$,
the ratio of $C_{l}$ is given by,

$\frac{C_{l}(\nu;f_{H})}{C_{l}(\nu;f_{H}=1)}=\left(\frac{\bar{\Omega}_{\rm
gw}(\nu;f_{H}=1)}{\bar{\Omega}_{\rm
gw}(\nu;f_{H})}\right)^{2}\left(\frac{I(\nu;f_{H})}{I(\nu;f_{H}=1)}\right)^{2},$
(35)

where $I(\nu;f_{H})=\int
dk\frac{\mathcal{L}^{2}(\eta_{*},\nu_{s};f_{H})}{(1+z_{*})^{2}}\delta_{m,k}^{2}(\eta_{*})$.

In the present scenario under consideration, the horizon mass is increased by
a factor of $f_{H}^{-2}$ when PBHs are formed in the radiation-dominated era,
and thus the maximal PBH mass is also increased by the same factor. Note that
in early matter domination, the horizon mass would also be increased, but by
the smaller factor $f_{H}^{-3/2}$; in a general cosmology where the equation
of state of the dominant species is $P=w\rho$, the horizon mass varies as
$f_{H}^{-3(w+1)/2}$. The mass function also depends on $f_{H}$ since the width
of the distribution PDF$(\delta_{1})$ depends on $f_{H}$, i.e.,
$\sigma_{0}=\sigma_{0}(f_{H})$ if $A_{s}$ is given. Since PDF$(\delta_{1})$ is
a Gaussian distribution, the shape of the distribution is sensitive to the
change in $\sigma_{0}$ Therefore, the merger rate and the GW energy spectrum
depend on $f_{H}$ non-trivially in general: For example, the merger rate
depends on the total mass of the binary system by $M_{\rm tot}^{-32/37}$, in
which we suppose the two BHs have the same mass, i.e.,
$\eta(m_{1},m_{2})=0.25=\eta_{\rm max}$, but for different $f_{H}$ the mass
function is different, and it can peak at different $x=M_{\rm PBH}/M_{H}$ (see
Figure 2), so only using $f_{H}$ cannot analytically predict the behavior of
the merger rate. However, by requiring a fixed $f_{\rm PBH}$, the
$\sigma_{0}^{2}$ is almost determined, and varying $f_{H}$ only affects the
horizon mass when PBHs form and PBH formation time. For the GW energy
spectrum, varying $f_{H}$ will directly vary the spectrum amplitude – the
spectrum amplitude is proportional to $M_{\rm tot}^{5/3}$ if $m_{1}=m_{2}$ for
the inspiral phase where $\nu$ is smaller than the merger frequency, and the
shape of the spectrum, e.g., the cutoff frequency shifts to a lower frequency
for a larger chirp mass $M_{c}$.

Figure 10:  The absolute value of angular power spectra as a function of
$f_{H}^{-1}$ at $10^{-4}$Hz (left) and 100Hz (right) for $l=1$ and different
values of $f_{H}^{2}M_{H}$ normalized by the corresponding
$C_{l=1}\times\bar{\Omega}_{\rm gw}^{2}\big{|}_{f_{H}=1}$, with ${f_{H}}_{\rm
min}=5000^{-1}$ (left) and $f_{H,{\rm min}}=300^{-1}$ (right). The black solid
line in the right panel indicates the value of 1.

Even though $\sigma_{0}^{2}$ is highly suppressed at a large
$f_{H}^{-1}\gtrsim 10$ for a given $A_{s}$, we further study the regime where
$f_{H}^{-1}\gg 1$ for completeness, noting that this regime is also
theoretically motivated by several deeply sub-horizon formation mechanisms
Cotner  _et al._ (2019); Hawking  _et al._ (1982). We note that since we
require $f_{\rm PBH}=0.01$ in the plots, the value of $\sigma_{0}^{2}$ is
almost unchanged even for large $f_{H}^{-1}$ by manually tuning $A_{s}$.
Figure 10 (left) shows the absolute angular power spectra
$C_{l}\times\bar{\Omega}_{\rm gw}$ of $l=1$ and $f_{\rm PBH}=0.05$ as a
function of $f_{H}^{-1}$ at $\nu=10^{-4}{\rm Hz}$ and
$f_{H}^{2}M_{H}\in\\{10^{-6}M_{\odot},10^{-5}M_{\odot},10^{-4}M_{\odot}\\}$,
while Figure 10 (right) shows the similar angular power spectra at
$\nu=100{\rm Hz}$ and
$f_{H}^{2}M_{H}\in\\{0.1M_{\odot},1M_{\odot},10M_{\odot}\\}$. There is a
manifest power-law behavior for the left subplot across the whole range of
$f_{H}^{-1}$ and for the right subplot right before the spectrum drops
significantly due to the fact that the observed frequency is far above the
cutoff frequency of the GW energy spectrum.

The $f_{H}$-dependence or effectively $M_{\rm PBH}$-dependence of the absolute
angular power spectrum can be understood as follows: First, when considering
the observed frequency in the inspiral phase of the IMR spectrum and the two
PBHs in the binary have the same mass, we notice the $f_{H}$-dependence of the
isotropic SGWB energy density :

$\bar{\Omega}_{\rm
gw,inspiral}\propto\mathcal{L}\propto\frac{dE}{d\nu}\bigg{|}_{\rm
inspiral}\times\mathcal{R}\propto f_{H}^{-1.6},$ (36)

where the power law index,
$1.6\approx\left(\frac{5}{3}-\frac{32}{37}\right)\times 2$, is determined by
the binary PBH mass dependence from both the merger rate and GW energy
spectrum, i.e., $\mathcal{R}\propto M_{b}^{-32/37}$ and
$\frac{dE}{d\nu}\bigg{|}_{\rm inspiral}\propto M_{b}^{5/3}$. By using the
above power-law $f_{H}$-dependence, we can further show that the
$C_{l}\times\bar{\Omega}_{\rm gw}^{2}$ ratio (shown in Figure 10) has the
power-law behavior with respect to $f_{H}$:

$C_{l}\times\bar{\Omega}_{\rm
gw}^{2}(f_{H})\propto\mathcal{L}^{2}\propto\left(\frac{dE}{d\nu}\bigg{|}_{\rm
inspiral}\times\mathcal{R}\right)^{2}\propto f_{H}^{-3.2},$ (37)

where we have used the fact that $f_{H}$-dependence on $C_{l}$ is sub-dominant
compared to the one on $\bar{\Omega}_{\rm gw}^{2}$, which can be found by
considering the cancellation of $f_{H}$-dependence between $\bar{\Omega}_{\rm
gw}$ and $\mathcal{L}$ in Eq. (35) and can be cross-checked by comparing
Figure 8 and Figure 10 for $\nu=100$ Hz and $10^{-4}$Hz.

Figure 11:  The angular power spectra as a function of mutipole $l$ for
$f_{H}^{-1}=100,700,5000$ with $\\{\nu=10^{-4}{\rm
Hz},~{}M_{H}=10^{-6}M_{\odot}\\}$ (left) and for $f_{H}^{-1}=10,30,100$ with
$\\{\nu=10{\rm Hz},~{}M_{H}=1M_{\odot}\\}$ (right).

Finally, we show in Figure 11 the angular power spectra as a function of $l$
given large values of $f_{H}^{-1}=100$, $700$, and $5000$ at
$\\{\nu=10^{-4}{\rm Hz},~{}M_{H}=10^{-6}M_{\odot}\\}$ (left panel), and the
corresponding ones for $f_{H}^{-1}=10$, $30$, and $100$ at $\\{\nu=10{\rm
Hz},~{}M_{H}=1M_{\odot}\\}$ in Figure 11 (right panel). There is no unexpected
behavior in the $C_{l}$ ratio for different $l$ except the case for
$f_{H}^{-1}=100$ in the right subplot. This is, again, because $M_{\rm PBH}$
increases with $f_{H}^{-2}$, and as a result the cutoff frequency of the IMR
waveform shifts to a lower frequency which is far below the observed
frequency.

##  IV Discussion and Conclusion

We performed, for the first time, a comprehensive study of the behavior of the
isotropic SGWB energy density spectrum and the angular power spectrum of PBH
binary mergers for PBH forming from the collapse of sub-horizon patches. We
introduced a model-independent parameter $f_{H}$ characterizing the fraction
of the wavelength of the perturbation mode in units of the horizon radius
where the patch gravitationally collapses into a BH. The standard scenario is
recovered for $f_{H}=1$, whereas values $f_{H}<1$ correspond to sub-horizon
collapsing regions. We illustrated our findings using a simple Dirac-$\delta$
peak curvature spectrum, but the spectrum shape can be modified easily
depending on the specific sub-horizon PBH formation mechanism. One of the
limitations in our results is the assumption of a spherically symmetric
distribution of overdensities, which in general can be non-spherically
distributed when entering the regime of very small $f_{H}$.

We found that the sub-horizon PBH formation in general enhances the isotropic
SGWB energy density $\bar{\Omega}_{\rm gw}$ and the absolute angular power
spectrum $C_{l}\times\bar{\Omega}^{2}_{\rm gw}$. However, the almost monotonic
increases in both $\bar{\Omega}_{\rm gw}$ and
$C_{l}\times\bar{\Omega}^{2}_{\rm gw}$, as $f_{H}$ decreases, cease when the
chirp mass of binary PBHs reaches a mass threshold determined by a given
observed frequency, such that the isotropic SGWB energy density spectrum
significantly drops above that specific cutoff frequency. The important effect
of the sub-horizon formation is changing the PBH mass (distribution) and also
the formation time or redshift of PBHs, which in turn affects the GW
observables, including both $\bar{\Omega}_{\rm gw}$ and
$C_{l}\times\bar{\Omega}^{2}_{\rm gw}$.

We investigated the isotropic SGWB energy density spectrum and the angular
power spectrum at various frequencies, PBH masses, and horizon size fractions
$f_{H}$ during PBH formation. When the observed frequency sits at the
frequency range of the inspiral phase of the IMR waveform, the isotropic SGWB
energy density spectrum and the absolute angular power spectrum at that
frequency have $f_{H}$-dependent power-law behaviors, i.e., $\bar{\Omega}_{\rm
gw}\propto f_{H}^{-1.6}$ and $C_{l}\times\bar{\Omega}^{2}_{\rm gw}\propto
f_{H}^{-3.2}$.

By introducing a sub-horizon formation scenario in the calculation of
$\bar{\Omega}_{\rm gw}$ and $C_{l}\times\bar{\Omega}^{2}_{\rm gw}$, one can
study the rich phenomenology of PBHs formed by non-standard mechanisms across
the universe’s history and provide a way to potentially test PBH formation
mechanisms upon the hopefully forthcoming detection of an SGWB signal.

## Acknowledgments

This work is supported in part by the U.S. Department of Energy grant number
de-sc0010107 (SP). F.W.Y. is supported by the U.S. Department of Energy under
Award No. DESC0009959.

## References

  * Hawking (1971) S. Hawking, Mon. Not. Roy. Astron. Soc. 152, 75 (1971).
  * Carr and Kuhnel (2020) B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Smyth  _et al._ (2022) N. Smyth, L. Santos-Olmsted, and S. Profumo, JCAP 03, 013 (2022), arXiv:2110.14660 [hep-ph] .
  * Morrison  _et al._ (2019) L. Morrison, S. Profumo, and Y. Yu, JCAP 05, 005 (2019), arXiv:1812.10606 [astro-ph.CO] .
  * Bean and Magueijo (2002) R. Bean and J. Magueijo, Phys. Rev. D 66, 063505 (2002), arXiv:astro-ph/0204486 .
  * Bird  _et al._ (2016) S. Bird, I. Cholis, J. B. Muñoz, Y. Ali-Haïmoud, M. Kamionkowski, E. D. Kovetz, A. Raccanelli, and A. G. Riess, Phys. Rev. Lett. 116, 201301 (2016), arXiv:1603.00464 [astro-ph.CO] .
  * Belczynski  _et al._ (2010) K. Belczynski, T. Bulik, C. L. Fryer, A. Ruiter, J. S. Vink, and J. R. Hurley, Astrophys. J. 714, 1217 (2010), arXiv:0904.2784 [astro-ph.SR] .
  * Giffin  _et al._ (2022) P. Giffin, J. Lloyd, S. D. McDermott, and S. Profumo, Phys. Rev. D 105, 123030 (2022), arXiv:2105.06504 [hep-ph] .
  * Fernandez and Profumo (2019) N. Fernandez and S. Profumo, JCAP 08, 022 (2019), arXiv:1905.13019 [astro-ph.HE] .
  * Arzoumanian  _et al._ (2020) Z. Arzoumanian _et al._ (NANOGrav), Astrophys. J. Lett. 905, L34 (2020), arXiv:2009.04496 [astro-ph.HE] .
  * Shannon  _et al._ (2015) R. M. Shannon _et al._ , Science 349, 1522 (2015), arXiv:1509.07320 [astro-ph.CO] .
  * Raidal  _et al._ (2017) M. Raidal, V. Vaskonen, and H. Veermäe, JCAP 09, 037 (2017), arXiv:1707.01480 [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, JCAP 02, 028 (2020), arXiv:1909.12619 [astro-ph.CO] .
  * Wang  _et al._ (2022) S. Wang, V. Vardanyan, and K. Kohri, Phys. Rev. D 106, 123511 (2022), arXiv:2107.01935 [gr-qc] .
  * Bavera  _et al._ (2022) S. S. Bavera, G. Franciolini, G. Cusin, A. Riotto, M. Zevin, and T. Fragos, Astron. Astrophys. 660, A26 (2022), arXiv:2109.05836 [astro-ph.CO] .
  * Bardeen  _et al._ (1986) J. M. Bardeen, J. R. Bond, N. Kaiser, and A. S. Szalay, Astrophys. J. 304, 15 (1986).
  * Young  _et al._ (2019) S. Young, I. Musco, and C. T. Byrnes, JCAP 11, 012 (2019), arXiv:1904.00984 [astro-ph.CO] .
  * Musco (2019) I. Musco, Phys. Rev. D 100, 123524 (2019), arXiv:1809.02127 [gr-qc] .
  * Young (2019) S. Young, Int. J. Mod. Phys. D 29, 2030002 (2019), arXiv:1905.01230 [astro-ph.CO] .
  * Nakama  _et al._ (2017) T. Nakama, J. Silk, and M. Kamionkowski, Phys. Rev. D 95, 043511 (2017), arXiv:1612.06264 [astro-ph.CO] .
  * Cotner  _et al._ (2019) E. Cotner, A. Kusenko, M. Sasaki, and V. Takhistov, JCAP 10, 077 (2019), arXiv:1907.10613 [astro-ph.CO] .
  * Hawking  _et al._ (1982) S. W. Hawking, I. G. Moss, and J. M. Stewart, Phys. Rev. D 26, 2681 (1982).
  * Garriga and Vilenkin (1993) J. Garriga and A. Vilenkin, Phys. Rev. D 47, 3265 (1993).
  * Caldwell and Casper (1996) R. R. Caldwell and P. Casper, Phys. Rev. D 53, 3002 (1996).
  * Meissner and Nicolai (2020) K. A. Meissner and H. Nicolai, Phys. Rev. D 102, 103008 (2020), arXiv:2007.11889 [astro-ph.HE] .
  * Kawana and Xie (2022) K. Kawana and K.-P. Xie, Phys. Lett. B 824, 136791 (2022), arXiv:2106.00111 [astro-ph.CO] .
  * Flores and Kusenko (2021) M. M. Flores and A. Kusenko, Phys. Rev. Lett. 126, 041101 (2021), arXiv:2008.12456 [astro-ph.CO] .
  * Chakraborty  _et al._ (2022) A. Chakraborty, P. K. Chanda, K. L. Pandey, and S. Das, Astrophys. J. 932, 119 (2022), arXiv:2204.09628 [astro-ph.CO] .
  * Acevedo  _et al._ (2021) J. F. Acevedo, J. Bramante, A. Goodman, J. Kopp, and T. Opferkuch, JCAP 04, 026 (2021), arXiv:2012.09176 [hep-ph] .
  * Abadie  _et al._ (2011) J. Abadie _et al._ (LIGO Scientific, VIRGO), Phys. Rev. D 83, 122005 (2011), [Erratum: Phys.Rev.D 86, 069903 (2012)], arXiv:1102.3781 [gr-qc] .
  * Moore  _et al._ (2015) C. J. Moore, R. H. Cole, and C. P. L. Berry, Class. Quant. Grav. 32, 015014 (2015), arXiv:1408.0740 [gr-qc] .
  * Allen and Romano (1999) B. Allen and J. D. Romano, Phys. Rev. D 59, 102001 (1999), arXiv:gr-qc/9710117 .
  * Ajith  _et al._ (2008) P. Ajith _et al._ , Phys. Rev. D 77, 104017 (2008), [Erratum: Phys.Rev.D 79, 129901 (2009)], arXiv:0710.2335 [gr-qc] .
  * Ajith  _et al._ (2011) P. Ajith _et al._ , Phys. Rev. Lett. 106, 241101 (2011), arXiv:0909.2867 [gr-qc] .
  * Cusin  _et al._ (2017) G. Cusin, C. Pitrou, and J.-P. Uzan, Phys. Rev. D 96, 103019 (2017).
  * Aghanim  _et al._ (2020) N. Aghanim _et al._ (Planck), Astron. Astrophys. 641, A6 (2020), [Erratum: Astron.Astrophys. 652, C4 (2021)], arXiv:1807.06209 [astro-ph.CO] .
  * Raidal  _et al._ (2019) M. Raidal, C. Spethmann, V. Vaskonen, and H. Veermäe, JCAP 02, 018 (2019), arXiv:1812.01930 [astro-ph.CO] .