††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS>Corresponding
author<EMAIL_ADDRESS>
# Primordial Black Hole Mass Function with Mass Gap

Xiao-Ming Bi1 Lu Chen2 Ke Wang1,3,4 1School of Physical Science and
Technology, Lanzhou University, Lanzhou 730000, China  2School of Physics and
Electronics, Shandong Normal University, Jinan 250014, China  3Lanzhou Center
for Theoretical Physics, Key Laboratory of Theoretical Physics of Gansu
Province, Lanzhou University, Lanzhou 730000, China  4Institute of Theoretical
Physics & Research Center of Gravitation, Lanzhou University, Lanzhou 730000,
China

###### Abstract

In this paper, we investigate the primordial black hole (PBH) mass function
with mass gap. Firstly, to obtain a data-supported PBH mass function with mass
gap for subsolar masses PBHs, we fine-tune the coefficients of a model-
independent power spectrum of primordial curvature perturbations. Then we take
this unique PBH mass function into consideration and calculate the energy
density spectrum of the stochastic gravitational wave background from PBH
mergers. We find the location of its first peak almost has no relationship
with the mass gap and is only determined by the probability distribution of
frequencies at which PBH binaries merge. Apart from the first peak, there must
be an accompanying smaller trough at higher frequency resulting from the mass
gap. Therefore, the detection of this smaller trough will provide more
information about inflation and PBH formation.

###### pacs:

???

##  I Introduction

Since LIGO detected the first binary black hole(BH) merger
LIGOScientific:2016aoc , it has become possible to directly constrain
primordial black holes (PBHs) using gravitational wave (GW) observations
Wang:2016ana ; LIGOScientific:2018glc ; LIGOScientific:2019kan . In contrast
to astrophysical BHs, the formation of PBHs is different and is not subject to
mass constraints from the stellar evolution model Rhoades:1974fn . More
precisely, PBHs whose masses range from Planck mass to masses of supermassive
BHs are possible to form through the gravitational collapse of primordial
overdensities in the early universe Hawking:1971ei . Therefore, PBHs can sweep
across the entire frequency band of most current and future GW detectors, such
as LIGO-Virgo network KAGRA:2013rdx ; LIGOScientific:2014pky ; VIRGO:2014yos
and LISA Robson:2018ifk . Also, there are many other talented methods to
constrain PBHs. For reviews of constraints on PBHs, see Ref. Carr:2020gox .
More importantly, as one of the most promising candidates for cold dark matter
(CDM), the constraints on the abundance of PBHs may help us end the debate
that which one is the main composition of CDM: WIMPs PandaX-II:2016vec ;
LUX:2015abn ; ATLASCMS ; AMS:2014bun ; Fermi-LAT:2011baq , Axion Marsh:2015xka
; Khmelnitsky:2013lxt ; Aoki:2016kwl ; Blas:2016ddr ; Boskovic:2018rub  or
PBHs Carr:2016drx .

Instead of a monochromatic mass distribution, PBHs probably span an extended
range of masses. There are three common types of mass function: the lognormal
mass function Dolgov:1992pu , the power-law mass function Carr:1975qj  and the
critical collapse mass function Yokoyama:1998xd ; Niemeyer:1999ak ;
Musco:2012au ; Carr:2016hva . All of these heuristic mass functions are
relatively simple. They are too smooth to directly describe a mass gap in mass
function, for example. In fact, mass gap is a common feature for astrophysical
BH mass distribution: on the one hand, stellar evolution models predict the
existence of a gap in the BH mass spectrum from about $55M_{\odot}$ to
$120M_{\odot}$ due to pair-instability supernovae Edelman:2021fik ;
Woosley:2019hmi ; Spera:2017fyx ; on the other hand, astronomical measurements
and independent sophisticated statistical analyses have found that there are
no observed BHs in the mass range about $2–5M_{\odot}$ Bailyn:1997xt ;
Ozel:2010su ; Farr:2010tu ; Belczynski:2011bn . Therefore, it is acceptable
that there may be also mass gaps in the PBH mass function resulting from some
independent PBH production mechanisms.

Consequently, PBH mass function with mass gap would leave some unique
footprints on the observations, such as stochastic gravitational wave
background (SGWB). Although there are a large number of sources of SGWB, most
of them just produce a smooth SGWB energy-density spectrum within a certain
frequency interval. Usually, these smooth SGWB energy-density spectra are
degenerate with each other. So it is hard for us to distinguish one source of
SGWB from the others. Of course, there are also some sources which can predict
an unsmooth SGWB energy-density spectrum with more than one peak Saito:2008jc
; Liu:2017hua ; Cai:2019jah . However, their unsmooth features in the SGWB
energy-density spectrum are different from the ones naturally resulting from
PBH mass function with mass gap, as shown in our paper.

In this paper, we first propose a model-independent power spectrum of
primordial curvature perturbations on PBH scales. Then we follow the procedure
of Wang:2019kaf  to relate the power spectrum to PBH mass function. As a
result, we can obtain any type of PBH mass function, especially the one with
mass gap, by fine-tuning the coefficients of power spectrum. Lastly, we use
subsolar masses PBHs as an example to investigate the effect of PBH mass
function with mass gap. More precisely, we calculate the SGWB energy-density
spectrum by integrating the contribution from all possible PBH binaries, as
did in Braglia:2021wwa . We find that mass gap serves as a natural suppressor
to suppress the contribution from certain PBH binaries and makes the SGWB
energy-density spectrum unsmooth.

This paper is organized as follows. In section II, a data-supported PBH mass
function with mass gap is given. In section III, the SGWB energy-density
spectrum from PBH mergers is calculated. Finally, a brief summary and
discussions are included in section IV. We adopt natural units $c=\hbar=1$.

##  II Mass Function with Mass Gap

On the cosmic microwave background (CMB) scales, the power spectrum of
primordial curvature perturbations is a quasi scale-invariant spectrum with
the amplitude $10^{9}A_{s}=2.1$ and the spectral index $n_{s}=0.96$
Planck:2018vyg . In order for PBHs to form, however, the amplitude of power
spectrum on PBH scales should be orders of magnitude larger. As for the shape
of power spectrum on PBH scales, it depends on the inflation model. Here, we
propose a model-independent power spectrum on PBH scales

$P_{\zeta}(k)=\sum_{i=1}^{N}A_{i}\delta\left(\ln k-\ln k_{i}\right),$ (1)

where $\delta(\ln k)$ is a delta function of $\ln k$, $A_{i}$ is the
dimensionless amplitude at the given wavenumber $\ln k_{i}$. Therefore, we can
mimic the power spectrum on PBH scales from any inflation model by fine-tuning
$A_{i}$ and $\ln k_{i}$.

Next, we can relate the variance of the primordial density perturbations in
the early Universe to the power spectrum of primordial curvature perturbations
as

$\displaystyle\sigma^{2}(k)$ $\displaystyle=$
$\displaystyle\int_{-\infty}^{+\infty}d\ln
q~{}w^{2}\left(\frac{4}{9}\right)^{2}\left(\frac{q}{k}\right)^{4}T^{2}P_{\zeta}(q),$
(2) $\displaystyle w$ $\displaystyle=$
$\displaystyle\exp\left(-\frac{q^{2}}{2k^{2}}\right),$ $\displaystyle T$
$\displaystyle=$ $\displaystyle 3\frac{\sin y-y\cos y}{y^{3}},$ $\displaystyle
y$ $\displaystyle=$ $\displaystyle\frac{q}{\sqrt{3}k},$

where is $w$ a Gaussian window function and $T$ is a transfer function between
the primordial density perturbations and the primordial curvature
perturbations Ando:2018qdb . If the primordial density perturbations with
wavenumber $k$ have a Gaussian distribution, then the probability distribution
of the smoothed density contrast is given by

$\mathcal{P}_{M_{H}}(\delta(M))=\frac{1}{\sqrt{2\pi\sigma^{2}\left(k\left(M_{H}\right)\right)}}\exp\left(-\frac{\delta^{2}(M)}{2\sigma^{2}\left(k\left(M_{H}\right)\right)}\right).$
(3)

The wavenumber $k$ is a function of the horizon mass $M_{H}$ Wang:2019kaf

$\displaystyle\frac{k}{k_{*}}=$ $\displaystyle 7.49\times
10^{7}\left(\frac{M_{\odot}}{M_{H}}\right)^{1/2}\left(\frac{g_{*,\rho}\left(T\left(M_{H}\right)\right)}{106.75}\right)^{1/4}\times$
$\displaystyle\left(\frac{g_{*,s}\left(T\left(M_{H}\right)\right)}{106.75}\right)^{-1/3},$
(4)

where $k_{*}=0.05{\rm Mpc}^{-1}$, the temperature $T$ is also a function of
$M_{H}$ Wang:2019kaf , $g_{*\rho}$ and $g_{*s}$ are the effective degrees of
freedom for the energy density and for the entropy density in the Standard
Model respectively. We assume that the PBHs are formed in the early Universe
through the critical collapse. Then the density contrast $\delta$ and the
horizon mass $M_{H}$ can determine the mass of PBHs $M$ as

$\displaystyle M=KM_{H}\left(\delta-\delta_{c}\right)^{\gamma},$ (5)

where $K=3.3$, $\gamma=0.36$ and $\delta_{c}=0.45$ are set by numerical
simulations. In the Press-Schechter formalism Press:1973iz , the probability
of the PBH production is related to the probability distribution of the
density contrast as

$\displaystyle\beta_{M_{H}}$
$\displaystyle=\int_{\delta_{c}}^{\infty}d\delta(M)\frac{M}{M_{H}}\mathcal{P}_{M_{H}}(\delta(M))$
$\displaystyle=\int_{-\infty}^{\infty}d\ln M\frac{d\delta(M)}{d\ln
M}\frac{M}{M_{H}}\mathcal{P}_{M_{H}}(\delta(M))$
$\displaystyle\equiv\int_{-\infty}^{\infty}d\ln
M\widetilde{\beta}_{M_{H}}(M),$ (6)

where $\widetilde{\beta}_{M_{H}}(M)$ has following explicit form

$\displaystyle\widetilde{\beta}_{M_{H}}(M)=$
$\displaystyle\frac{K}{\sqrt{2\pi}\gamma\sigma\left(k\left(M_{H}\right)\right)}\left(\frac{M}{KM_{H}}\right)^{1+\frac{1}{\gamma}}\times$
$\displaystyle\exp\left(-\frac{1}{2\sigma^{2}\left(k\left(M_{H}\right)\right)}\left(\delta_{c}+\left(\frac{M}{KM_{H}}\right)^{\frac{1}{\gamma}}\right)^{2}\right).$
(7)

Given the definition of mass function of PBHs $f(M)\equiv\frac{1}{\Omega_{\rm
CDM}}\frac{d\Omega_{\rm PBH}}{d\ln M}$, we can obtain the mass function of
PBHs as

$\displaystyle f(M)=$
$\displaystyle\frac{\Omega_{\mathrm{m}}}{\Omega_{\mathrm{CDM}}}\int_{-\infty}^{\infty}d\ln
M_{H}\widetilde{\beta}_{M_{H}}(M)\times$
$\displaystyle\left(\frac{g_{*,\rho}\left(T\left(M_{H}\right)\right)}{g_{*,\rho}\left(T_{\mathrm{eq}}\right)}\frac{g_{*,s}\left(T_{\mathrm{eq}}\right)}{g_{*,s}\left(T\left(M_{H}\right)\right)}\frac{T\left(M_{H}\right)}{T_{\mathrm{eq}}}\right),$
(8)

where $T_{\rm eq}$ is the temperature of the epoch of matter-radiation
equality. Then we can obtain the abundance of PBHs in CDM as $f_{\rm PBH}=\int
f(M)d\ln(M/M_{\odot})$. For subsolar masses PBHs, there are constraints on
$f_{\rm PBH}$ for monochromatic mass function from the microlensing
observations of Subaru/HSC Niikura:2017zjd , OGLE Niikura:2019kqi , EROS-2
EROS-2:2006ryy , MACHO Macho:2000nvd  and the caustic crossing Oguri:2017ock
as shown by the red dashed line in Fig. 1. Since the power spectrum in Eq. (1)
is model-independent, we can obtain a just data-supported $f(M)$ (blue solid)
by seting $N=18$ and $k_{i}(M_{H})=k(10^{(2-i)/2}\times M_{\odot}$) and fine-
tuning $A_{i}$. To obtain a mass function with mass gap, we only need to set
certain $A_{i}$s equal to $0$: $A_{8}=0$ leads to a little gap (red solid);
$A_{8}=A_{9}=0$ leads to a middle gap (orange solid);
$A_{8}=A_{9}=A_{10}=A_{11}=0$ leads to a large gap (purple solid).

Figure 1: Mass functions $f(M)$ derived from the model-independent power
spectrum $P_{\zeta}(k)$, where $N=18$ and $k_{i}(M_{H})=k(10^{(2-i)/2}\times
M_{\odot}$). $A_{i}$ is fine-tuned so that $f(M)$ is similar to the
constraints (blue solid). $A_{8}=0$ leads to a little gap (red solid).
$A_{8}=A_{9}=0$ leads to a middle gap (orange solid).
$A_{8}=A_{9}=A_{10}=A_{11}=0$ leads to a large gap (purple solid). For all
cases. the abundance of PBHs in CDM $f_{\rm PBH}=\int
f(M)d\ln(M/M_{\odot})<1$. The constraints on $f_{\rm PBH}$ for monochromatic
mass function from the microlensing observations of Subaru/HSC Niikura:2017zjd
, OGLE Niikura:2019kqi , EROS-2 EROS-2:2006ryy , MACHO Macho:2000nvd  and the
caustic crossing Oguri:2017ock  are also shown (red dashed).

##  III Stochastic Gravitational Wave Background from PBH Mergers

The present day SGWB energy-density spectrum is given by following integral
Rosado:2011kv

$\displaystyle\Omega_{\mathrm{GW}}(f)=$
$\displaystyle\frac{f}{\rho_{c}}\iint\mathrm{d}\lg m_{1}\mathrm{d}\lg
m_{2}\int_{z_{\min}}^{z_{\max}}\mathrm{d}z^{\prime}\frac{1}{\left(1+z^{\prime}\right)H\left(z^{\prime}\right)}\times$
$\displaystyle\frac{\mathrm{d}^{2}\tau_{\operatorname{merg}}\left(z^{\prime},m_{1},m_{2}\right)}{\mathrm{d}\lg
m_{1}\mathrm{d}\lg
m_{2}}\frac{\mathrm{d}E_{\mathrm{GW}}\left(f_{s}\right)}{\mathrm{d}f_{s}},$
(9)

where $f_{s}=f(1+z)$ is frequency in the source frame,
$\rho_{c}=3H_{0}^{2}/8\pi G$ is the critical density of the Universe and the
Hubble parameter $H(z)$ is calculated by $\Lambda$CDM model Planck:2018vyg .
The merger rate Braglia:2021wwa  is given by

$\displaystyle\frac{\mathrm{d}^{2}\tau\left(z,m_{1},m_{2}\right)}{\mathrm{d}\lg
m_{1}\mathrm{~{}d}\lg m_{2}}=$ $\displaystyle
R_{\mathrm{clust}}\frac{\left(m_{1}+m_{2}\right)^{10/7}}{\left(m_{1}m_{2}\right)^{5/7}}\frac{\left(1+z\right)^{\alpha_{z}}}{\left(1+M_{\mathrm{{tot}}}/M_{*}\right)^{\alpha_{c}}}\times$
$\displaystyle f\left(m_{1}\right)f\left(m_{2}\right),$ (10)

where $M_{\rm tot}=m_{1}+m_{2}$ is the total mass of a PBH binary, $R_{\rm
clust}$ is equal to $3.3\times 10^{4}{\rm yr}^{-1}{\rm Gpc}^{-3}$ for subsoler
masses PBHs or $2.1\times 10^{4}{\rm yr}^{-1}{\rm Gpc}^{-3}$ for mass bin
$[5,100]M_{\odot}$, $M_{*}=10^{5}M_{\odot}$ is a cutoff to suppress the merger
rate when $M_{\rm tot}>M_{*}$, $\alpha_{z}$ and $\alpha_{c}$ are positive
indexes to weight the corresponding contribution. There are also other methods
to give the merger rate Ali-Haimoud:2017rtz ; Chen:2018czv ; Raidal:2018bbj .
The energy spectrum of a single PBH binary Ajith:2009bn  is given by

$\displaystyle\frac{\mathrm{d}E_{\mathrm{GW}}(f)}{\mathrm{d}f}=$
$\displaystyle\frac{\pi^{2/3}}{3}\left(G\mathcal{M}_{c}\right)^{5/3}f^{-1/3}\times$
$\displaystyle\left\\{\begin{array}[]{ll}\left(1+\alpha_{2}u^{2}\right)^{2}&\text{
for }f<f_{1},\\\
w_{1}f\left(1+\epsilon_{1}u+\epsilon_{2}u^{2}\right)^{2}&\text{ for }f_{1}\leq
f<f_{2},\\\
w_{2}f^{7/3}\frac{f_{4}^{4}}{\left(4\left(f-f_{2}\right)^{2}+f_{4}^{2}\right)^{2}}&\text{
for }f_{2}\leq f<f_{3},\\\ 0&\text{ for }f_{3}\leq f,\end{array}\right.$ (15)

where $\mathcal{M}_{c}=(m_{1}m_{2})^{3/5}/(m_{1}+m_{2})^{1/5}$,
$u_{\left(i\right)}\equiv\left(\pi
M_{\mathrm{tot}}Gf_{\left(i\right)}\right)^{1/3},\eta=m_{1}m_{2}/M_{\mathrm{tot}}^{2},\alpha_{2}=-323/224+451/168\eta,\epsilon_{1}=-1.8897,\epsilon_{2}=1.6557$,

$\displaystyle w_{1}$ $\displaystyle=$ $\displaystyle
f_{1}^{-1}\dfrac{\left[1+\alpha_{2}u_{1}^{2}\right]^{2}}{\left[1+\epsilon_{1}u_{1}+\epsilon_{2}u_{1}^{2}\right]^{2}},$
$\displaystyle w_{2}$ $\displaystyle=$ $\displaystyle
w_{1}f_{2}^{-4/3}\left[1+\epsilon_{1}u_{2}+\epsilon_{2}u_{2}^{2}\right]^{2},$
(16)

and

$\displaystyle u_{1}^{3}$ $\displaystyle=$ $\displaystyle
0.066+0.6437\eta-0.05822\eta^{2}-7.092\eta^{3},$ $\displaystyle u_{2}^{3}$
$\displaystyle=$ $\displaystyle
0.37/2+0.1469\eta-0.0249\eta^{2}+2.325\eta^{3},$ $\displaystyle u_{3}^{3}$
$\displaystyle=$ $\displaystyle
0.3236-0.1331\eta-0.2714\eta^{2}+4.922\eta^{3},$ $\displaystyle u_{4}^{3}$
$\displaystyle=$
$\displaystyle(1-0.63)/4-0.4098\eta+1.829\eta^{2}-2.87\eta^{3}.$ (17)

In Fig. 2, we show the probability density of $f_{i}$ for a given mass span.
According to Eq. (III), we know that only these PBH binaries whose $f_{i}$ are
larger than $f$ may contribute to $\Omega_{\mathrm{GW}}(f)$. Therefore, for
smaller $f$, all the PBH binaries may contribute to $\Omega_{\mathrm{GW}}(f)$.
As $f$ increases, $f$ is larger than $f_{i}$ of enough PBH binaries, which
gives the location of the first peak of $\Omega_{\mathrm{GW}}(f)$. As shown in
Fig. 2 and Fig. 3, the peak of the probability density of $f_{i}$ coincides
with the first peak of $\Omega_{\mathrm{GW}}(f)$. As $f$ continue to increase,
the higher $f$, the more PBH binaries whose $f_{i}$ are smaller than $f$, as
shown in Fig. 2. That is to say, there are more PBH binaries that will not
contribute to $\Omega_{\mathrm{GW}}(f)$. So, after the first peak,
$\Omega_{\mathrm{GW}}(f)$ decreases with $f$, as shown by the blue solid line
in Fig. 3. When we take the mass gap into consideration, there are some PBH
binaries that meet not only the condition of $f_{i}>f$ but also the condition
of $f(m_{1})f(m_{2})\sim 0$. Therefore, these PBH binaries also don’t
contribute to $\Omega_{\mathrm{GW}}(f)$ and there is a trough in
$\Omega_{\mathrm{GW}}(f)$. In Fig. 3, we plot $\Omega_{\mathrm{GW}}(f)$ for
PBH mass function with no gap (blue solid), little gap (red solid), middle gap
(orange solid) and large gap (purple solid) individually by setting $R_{\rm
clust}=3.3\times 10^{4}{\rm yr}^{-1}{\rm Gpc}^{-3}$, $M_{*}=10^{5}M_{\odot}$,
$\alpha_{z}=0$ and $\alpha_{c}=0$. We also show the sensitivity curves for
DECIGO/BBO Sato:2017dkf ; Harry:2006fi  (green dashed), ET+CE Punturo:2010zz ;
Reitze:2019iox  (cyan dashed) and LISA Robson:2018ifk  (purple dashed)
individually. We find that these detectors just can probe the behavior of
$f^{2/3}$ of $\Omega_{\mathrm{GW}}(f)$. And the first peak of
$\Omega_{\mathrm{GW}}(f)$ is out of the reach of these detectors, let alone
the footprints of mass gap.

Figure 2: The probability density of $f_{i}$ for a given mass span. Only these
PBH binaries whose $f_{i}$ are larger than $f$ may contribute to
$\Omega_{\mathrm{GW}}(f)$. Figure 3: The present day SGWB energy-density
spectra for mass function with no gap (blue solid), little gap (red solid),
middle gap (orange solid) and large gap (purple solid), where $R_{\rm
clust}=3.3\times 10^{4}{\rm yr}^{-1}{\rm Gpc}^{-3}$, $M_{*}=10^{5}M_{\odot}$,
$\alpha_{z}=0$ and $\alpha_{c}=0$. We also show the sensitivity curves for
DECIGO/BBO Sato:2017dkf ; Harry:2006fi  (green dashed), ET+CE Punturo:2010zz ;
Reitze:2019iox  (cyan dashed) and LISA Robson:2018ifk  (purple dashed)
individually.

##  IV Summary and Discussions

In this paper, we first construct a model-independent power spectrum of
primordial curvature perturbations on PBH scales through the superposition of
a series of delta functions. This power spectrum can not only mimic the other
power spectra from certain inflation models but also be easily related to the
variance of the primordial density perturbations. Then we follow the procedure
of Wang:2019kaf  to utilize the Press-Schechter formalism Press:1973iz  and
obtain the PBH mass function. By fine-tuning the amplitude and location of
$18$ delta functions, we can not only make the mass function satisfied with
the constraints on subsolar masses PBHs, but also give the mass function
arbitrary mass gap. Lastly, we look for the footprints of mass gap by
comparing the present day SGWB energy-density spectra derived from mass
functions with no gap or different gap. We find that there is a trough in the
SGWB energy-density spectra due to mass gap, which contains the information
about the power spectrum of primordial curvature perturbations. Unfortunately,
such unique feature is out of the reach of detectors.

For mass span $[0.1,10^{9}]M_{\odot}$, there are four stronger constraints on
$f_{\rm PBH}$ from CMB Serpico:2020ehh , X-ray Inoue:2017csr , dynamical
friction Carr:1997cn  and GW LIGOScientific:2018glc ; LIGOScientific:2019kan .
Compared with the other three constraints, the constraints from CMB are much
stronger, reaching the level $f_{\rm PBH}<3\times 10^{-9}$ around
$10^{4}M_{\odot}$. We assume a mass function shares the same shape with the
combination of these four constraints. So the mass function has an approximate
mass gap around $10^{4}M_{\odot}$, as shown in Fig. 10 and Fig. 15 of
Carr:2020gox . Using this mass function, we calculate the present day SGWB
energy-density spectrum, as shown in Fig. 4. We find both of the first peak
and the trough can be detected. Of course, here we have ignored the
contributions of astrophysical BHs.

Figure 4: The present day SGWB energy-density spectrum for PBH mass function
summarized from the constraints on $f_{\rm PBH}$ for mass span
$[0.1,10^{9}]M_{\odot}$ Serpico:2020ehh ; Inoue:2017csr ; Carr:1997cn ;
LIGOScientific:2018glc ; LIGOScientific:2019kan , where $R_{\rm
clust}=2.1\times 10^{4}{\rm yr}^{-1}{\rm Gpc}^{-3}$, $M_{*}=10^{5}M_{\odot}$,
$\alpha_{z}=0$ and $\alpha_{c}=1$. We also show the sensitivity curves for
DECIGO/BBO Sato:2017dkf ; Harry:2006fi  (green dashed), ET+CE Punturo:2010zz ;
Reitze:2019iox  (cyan dashed) and LISA Robson:2018ifk  (purple dashed)
individually.

###### Acknowledgements.

Ke Wang is supported by grants from National Key Research and Development
Program of China Grant No.2020YFC2201502 and grants from NSFC (grant No.
12005084 and grant No. 12047501). Lu Chen is supported by grants from NSFC
(grant No. 12105164). This work has also received funding from project
ZR2021QA021 supported by Shandong Provincial Natural Science Foundation.

## References

  * (1) B. P. Abbott et al. [LIGO Scientific and Virgo], “Observation of Gravitational Waves from a Binary Black Hole Merger,” Phys. Rev. Lett. 116, no.6, 061102 (2016) [arXiv:1602.03837 [gr-qc]].
  * (2) S. Wang, Y. F. Wang, Q. G. Huang and T. G. F. Li, “Constraints on the Primordial Black Hole Abundance from the First Advanced LIGO Observation Run Using the Stochastic Gravitational-Wave Background,” Phys. Rev. Lett. 120, no.19, 191102 (2018) [arXiv:1610.08725 [astro-ph.CO]].
  * (3) B. P. Abbott et al. [LIGO Scientific and Virgo], “Search for Subsolar-Mass Ultracompact Binaries in Advanced LIGO’s First Observing Run,” Phys. Rev. Lett. 121, no.23, 231103 (2018) [arXiv:1808.04771 [astro-ph.CO]].
  * (4) B. P. Abbott et al. [LIGO Scientific and Virgo], “Search for Subsolar Mass Ultracompact Binaries in Advanced LIGO’s Second Observing Run,” Phys. Rev. Lett. 123, no.16, 161102 (2019) [arXiv:1904.08976 [astro-ph.CO]].
  * (5) D. Prialnik, “An Introduction to the Theory of Stellar Structure and Evolution”, Cambridge University Press, ISBN 0-521-65065-8 (2000)
  * (6) S. Hawking, “Gravitationally collapsed objects of very low mass,” Mon. Not. Roy. Astron. Soc. 152, 75 (1971)
  * (7) B. P. Abbott et al. [KAGRA, LIGO Scientific, Virgo and VIRGO], “Prospects for observing and localizing gravitational-wave transients with Advanced LIGO, Advanced Virgo and KAGRA,” Living Rev. Rel. 21, no.1, 3 (2018) [arXiv:1304.0670 [gr-qc]].
  * (8) J. Aasi et al. [LIGO Scientific], “Advanced LIGO,” Class. Quant. Grav. 32, 074001 (2015) [arXiv:1411.4547 [gr-qc]].
  * (9) F. Acernese et al. [VIRGO], “Advanced Virgo: a second-generation interferometric gravitational wave detector,” Class. Quant. Grav. 32, no.2, 024001 (2015) [arXiv:1408.3978 [gr-qc]].
  * (10) T. Robson, N. J. Cornish and C. Liu, “The construction and use of LISA sensitivity curves,” Class. Quant. Grav. 36, no.10, 105011 (2019) [arXiv:1803.01944 [astro-ph.HE]].
  * (11) B. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, “Constraints on primordial black holes,” Rept. Prog. Phys. 84, no.11, 116902 (2021) [arXiv:2002.12778 [astro-ph.CO]].
  * (12) A. Tan et al. [PandaX-II], “Dark Matter Results from First 98.7 Days of Data from the PandaX-II Experiment,” Phys. Rev. Lett. 117, no.12, 121303 (2016) [arXiv:1607.07400 [hep-ex]].
  * (13) D. S. Akerib et al. [LUX], “Improved Limits on Scattering of Weakly Interacting Massive Particles from Reanalysis of 2013 LUX Data,” Phys. Rev. Lett. 116, no.16, 161301 (2016) [arXiv:1512.03506 [astro-ph.CO]].
  * (14) “ATLAS Collaboration and CMS Collaboration (ATLAS, CMS) (2016).”
  * (15) L. Accardo et al. [AMS], “High Statistics Measurement of the Positron Fraction in Primary Cosmic Rays of 0.5–500 GeV with the Alpha Magnetic Spectrometer on the International Space Station,” Phys. Rev. Lett. 113, 121101 (2014)
  * (16) M. Ackermann et al. [Fermi-LAT], “Measurement of separate cosmic-ray electron and positron spectra with the Fermi Large Area Telescope,” Phys. Rev. Lett. 108, 011103 (2012) [arXiv:1109.0521 [astro-ph.HE]].
  * (17) D. J. E. Marsh, “Axion Cosmology,” Phys. Rept. 643, 1-79 (2016) [arXiv:1510.07633 [astro-ph.CO]].
  * (18) A. Khmelnitsky and V. Rubakov, “Pulsar timing signal from ultralight scalar dark matter,” JCAP 02, 019 (2014) [arXiv:1309.5888 [astro-ph.CO]].
  * (19) A. Aoki and J. Soda, “Detecting ultralight axion dark matter wind with laser interferometers,” Int. J. Mod. Phys. D 26, no.07, 1750063 (2016) [arXiv:1608.05933 [astro-ph.CO]].
  * (20) D. Blas, D. L. Nacir and S. Sibiryakov, “Ultralight Dark Matter Resonates with Binary Pulsars,” Phys. Rev. Lett. 118, no.26, 261102 (2017) [arXiv:1612.06789 [hep-ph]].
  * (21) M. Bošković, F. Duque, M. C. Ferreira, F. S. Miguel and V. Cardoso, “Motion in time-periodic backgrounds with applications to ultralight dark matter haloes at galactic centers,” Phys. Rev. D 98, 024037 (2018) [arXiv:1806.07331 [gr-qc]].
  * (22) B. Carr, F. Kuhnel and M. Sandstad, “Primordial Black Holes as Dark Matter,” Phys. Rev. D 94, no.8, 083504 (2016) [arXiv:1607.06077 [astro-ph.CO]].
  * (23) A. Dolgov and J. Silk, “Baryon isocurvature fluctuations at small scales and baryonic dark matter,” Phys. Rev. D 47, 4244-4255 (1993)
  * (24) B. J. Carr, “The Primordial black hole mass spectrum,” Astrophys. J. 201, 1-19 (1975)
  * (25) J. Yokoyama, “Cosmological constraints on primordial black holes produced in the near critical gravitational collapse,” Phys. Rev. D 58, 107502 (1998) [arXiv:gr-qc/9804041 [gr-qc]].
  * (26) J. C. Niemeyer and K. Jedamzik, “Dynamics of primordial black hole formation,” Phys. Rev. D 59, 124013 (1999) [arXiv:astro-ph/9901292 [astro-ph]].
  * (27) I. Musco and J. C. Miller, “Primordial black hole formation in the early universe: critical behaviour and self-similarity,” Class. Quant. Grav. 30, 145009 (2013) [arXiv:1201.2379 [gr-qc]].
  * (28) B. J. Carr, K. Kohri, Y. Sendouda and J. Yokoyama, “Constraints on primordial black holes from the Galactic gamma-ray background,” Phys. Rev. D 94, no.4, 044029 (2016) [arXiv:1604.05349 [astro-ph.CO]].
  * (29) B. Edelman, Z. Doctor and B. Farr, “Poking Holes: Looking for Gaps in LIGO/Virgo’s Black Hole Population,” Astrophys. J. Lett. 913, no.2, L23 (2021) [arXiv:2104.07783 [astro-ph.HE]].
  * (30) S. E. Woosley, “The Evolution of Massive Helium Stars, Including Mass Loss,” Astrophys. J. 878, no.1, 49 (2019)
  * (31) M. Spera and M. Mapelli, “Very massive stars, pair-instability supernovae and intermediate-mass black holes with the SEVN code,” Mon. Not. Roy. Astron. Soc. 470, no.4, 4739-4749 (2017) [arXiv:1706.06109 [astro-ph.SR]].
  * (32) C. D. Bailyn, R. K. Jain, P. Coppi and J. A. Orosz, “The Mass distribution of stellar black holes,” Astrophys. J. 499, 367 (1998) [arXiv:astro-ph/9708032 [astro-ph]].
  * (33) F. Ozel, D. Psaltis, R. Narayan and J. E. McClintock, “The Black Hole Mass Distribution in the Galaxy,” Astrophys. J. 725, 1918-1927 (2010) [arXiv:1006.2834 [astro-ph.GA]].
  * (34) W. M. Farr, N. Sravan, A. Cantrell, L. Kreidberg, C. D. Bailyn, I. Mandel and V. Kalogera, “The Mass Distribution of Stellar-Mass Black Holes,” Astrophys. J. 741, 103 (2011) [arXiv:1011.1459 [astro-ph.GA]].
  * (35) K. Belczynski, G. Wiktorowicz, C. Fryer, D. Holz and V. Kalogera, “Missing Black Holes Unveil The Supernova Explosion Mechanism,” Astrophys. J. 757, 91 (2012) [arXiv:1110.1635 [astro-ph.GA]].
  * (36) R. Saito and J. 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]].
  * (37) J. Liu, Z. K. Guo, R. G. Cai and G. Shiu, “Gravitational Waves from Oscillons with Cuspy Potentials,” Phys. Rev. Lett. 120, no.3, 031301 (2018) [arXiv:1707.09841 [astro-ph.CO]].
  * (38) Y. F. Cai, C. Chen, X. Tong, D. G. Wang and S. F. Yan, “When Primordial Black Holes from Sound Speed Resonance Meet a Stochastic Background of Gravitational Waves,” Phys. Rev. D 100, no.4, 043518 (2019) [arXiv:1902.08187 [astro-ph.CO]].
  * (39) S. Wang, T. Terada and K. 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, no.10, 103531 (2019) [erratum: Phys. Rev. D 101, no.6, 069901 (2020)] [arXiv:1903.05924 [astro-ph.CO]].
  * (40) M. Braglia, J. Garcia-Bellido and S. Kuroyanagi, “Testing Primordial Black Holes with multi-band observations of the stochastic gravitational wave background,” JCAP 12, no.12, 012 (2021) [arXiv:2110.07488 [astro-ph.CO]].
  * (41) N. Aghanim et al. [Planck], “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641, A6 (2020) [erratum: Astron. Astrophys. 652, C4 (2021)] [arXiv:1807.06209 [astro-ph.CO]].
  * (42) K. Ando, K. Inomata and M. Kawasaki, “Primordial black holes and uncertainties in the choice of the window function,” Phys. Rev. D 97, no.10, 103528 (2018) [arXiv:1802.06393 [astro-ph.CO]].
  * (43) W. H. Press and P. Schechter, “Formation of galaxies and clusters of galaxies by selfsimilar gravitational condensation,” Astrophys. J. 187, 425-438 (1974)
  * (44) H. Niikura, M. Takada, N. Yasuda, R. H. Lupton, T. Sumi, S. More, T. Kurita, S. Sugiyama, A. More and M. Oguri, et al. “Microlensing constraints on primordial black holes with Subaru/HSC Andromeda observations,” Nature Astron. 3, no.6, 524-534 (2019) [arXiv:1701.02151 [astro-ph.CO]].
  * (45) H. Niikura, M. Takada, S. Yokoyama, T. Sumi and S. Masaki, “Constraints on Earth-mass primordial black holes from OGLE 5-year microlensing events,” Phys. Rev. D 99, no.8, 083503 (2019) [arXiv:1901.07120 [astro-ph.CO]].
  * (46) 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 [astro-ph]].
  * (47) R. A. Allsman et al. [Macho], “MACHO project limits on black hole dark matter in the 1-30 solar mass range,” Astrophys. J. Lett. 550, L169 (2001) [arXiv:astro-ph/0011506 [astro-ph]].
  * (48) M. Oguri, J. M. Diego, N. Kaiser, P. L. Kelly and T. Broadhurst, “Understanding caustic crossings in giant arcs: characteristic scales, event rates, and constraints on compact dark matter,” Phys. Rev. D 97, no.2, 023518 (2018) [arXiv:1710.00148 [astro-ph.CO]].
  * (49) P. A. Rosado, “Gravitational wave background from binary systems,” Phys. Rev. D 84, 084004 (2011) [arXiv:1106.5795 [gr-qc]].
  * (50) Y. Ali-Haïmoud, E. D. Kovetz and M. Kamionkowski, “Merger rate of primordial black-hole binaries,” Phys. Rev. D 96, no.12, 123523 (2017) [arXiv:1709.06576 [astro-ph.CO]].
  * (51) Z. C. Chen and Q. G. Huang, “Merger Rate Distribution of Primordial-Black-Hole Binaries,” Astrophys. J. 864, no.1, 61 (2018) [arXiv:1801.10327 [astro-ph.CO]].
  * (52) M. Raidal, C. Spethmann, V. Vaskonen and H. Veermäe, “Formation and Evolution of Primordial Black Hole Binaries in the Early Universe,” JCAP 02, 018 (2019) [arXiv:1812.01930 [astro-ph.CO]].
  * (53) P. Ajith, M. Hannam, S. Husa, Y. Chen, B. Bruegmann, N. Dorband, D. Muller, F. Ohme, D. Pollney and C. Reisswig, et al. “Inspiral-merger-ringdown waveforms for black-hole binaries with non-precessing spins,” Phys. Rev. Lett. 106, 241101 (2011) [arXiv:0909.2867 [gr-qc]].
  * (54) S. Sato, S. Kawamura, M. Ando, T. Nakamura, K. Tsubono, A. Araya, I. Funaki, K. Ioka, N. Kanda and S. Moriwaki, et al. “The status of DECIGO,” J. Phys. Conf. Ser. 840, no.1, 012010 (2017)
  * (55) G. M. Harry, P. Fritschel, D. A. Shaddock, W. Folkner and E. S. Phinney, “Laser interferometry for the big bang observer,” Class. Quant. Grav. 23, 4887-4894 (2006) [erratum: Class. Quant. Grav. 23, 7361 (2006)]
  * (56) M. Punturo, M. Abernathy, F. Acernese, B. Allen, N. Andersson, K. Arun, F. Barone, B. Barr, M. Barsuglia and M. Beker, et al. “The Einstein Telescope: A third-generation gravitational wave observatory,” Class. Quant. Grav. 27, 194002 (2010)
  * (57) D. Reitze, R. X. Adhikari, S. Ballmer, B. Barish, L. Barsotti, G. Billingsley, D. A. Brown, Y. Chen, D. Coyne and R. Eisenstein, et al. “Cosmic Explorer: The U.S. Contribution to Gravitational-Wave Astronomy beyond LIGO,” Bull. Am. Astron. Soc. 51, no.7, 035 (2019)
  * (58) P. D. Serpico, V. Poulin, D. Inman and K. Kohri, “Cosmic microwave background bounds on primordial black holes including dark matter halo accretion,” Phys. Rev. Res. 2, no.2, 023204 (2020) [arXiv:2002.10771 [astro-ph.CO]].
  * (59) Y. Inoue and A. Kusenko, “New X-ray bound on density of primordial black holes,” JCAP 10, 034 (2017) [arXiv:1705.00791 [astro-ph.CO]].
  * (60) B. J. Carr and M. Sakellariadou, “Dynamical constraints on dark compact objects,” Astrophys. J. 516, 195-220 (1999)