# Constraining properties of asymmetric dark matter candidates from
gravitational-wave observations

Divya Singh<EMAIL_ADDRESS>Institute for Gravitation and the Cosmos,
Department of Physics, Pennsylvania State University, University Park, PA,
16802, USA  Anuradha Gupta  Department of Physics and Astronomy, The
University of Mississippi, University MS 38677, USA  Emanuele Berti
Department of Physics and Astronomy, Johns Hopkins University, 3400 N. Charles
Street, Baltimore, Maryland, 21218, USA  Sanjay Reddy  Institute for Nuclear
Theory, University of Washington, Seattle, WA USA  B. S. Sathyaprakash
Institute for Gravitation and the Cosmos, Department of Physics, Pennsylvania
State University, University Park, PA, 16802, USA  Department of Astronomy &
Astrophysics, Pennsylvania State University, University Park, PA, 16802, USA
School of Physics and Astronomy, Cardiff University, Cardiff, UK, CF24 3AA

###### Abstract

The accumulation of certain types of dark matter particles in neutron star
cores due to accretion over long timescales can lead to the formation of a
mini black hole. In this scenario, the neutron star is destabilized and
implodes to form a black hole without significantly increasing its mass. When
this process occurs in neutron stars in coalescing binaries, one or both stars
might be converted to a black hole before they merge. Thus, in the mass range
of $\sim\mbox{1--2}\,M_{\odot},$ the Universe might contain three distinct
populations of compact binaries: one containing only neutron stars, the second
population of only black holes, and a third, mixed population consisting of a
neutron star and a black hole. However, it is unlikely to have a mixed
population as the various timescales allow for both neutron stars to remain or
collapse within a short timescale. In this paper, we explore the capability of
future gravitational-wave detector networks, including upgrades of Advanced
LIGO and Virgo, and new facilities such as the Cosmic Explorer and Einstein
Telescope (XG network), to discriminate between different populations by
measuring the effective tidal deformability of the binary, which is zero for
binary black holes but nonzero for binary neutron stars. Furthermore, we show
that observing the relative abundances of the different populations can be
used to infer the timescale for neutron stars to implode into black holes, and
in turn, provide constraints on the particle nature of dark matter. The XG
network will infer the implosion timescale to within an accuracy of 0.01 Gyr
at 90% credible interval and determine the dark matter mass and interaction
cross section to within a factor of 2 GeV and 10 cm-2, respectively.

Gravitational Waves; Astrophysics, Black Holes, Star Clusters, Dark Matter

##  I Introduction and background

The origin and properties of dark matter (DM) have been long-standing problems
in fundamental physics and cosmology. Astronomical observations have
increasingly provided evidence for a non-baryonic component of matter that
either does not interact electromagnetically with baryons or has a negligibly
small interaction cross-section. Consequently, the presence of DM is inferred
due to its gravitational effect on baryonic matter. Laboratory experiments to
detect DM particles from their weak interaction with baryons have so far
produced null results, as have the observations of decay products that would
result from the annihilation of certain types of DM particles. Although there
are a few plausible DM candidates in the Standard Model, theoretical insight
into what they might be in theories beyond the Standard Model is plentiful and
not very constraining. Currently, there is an effort to look for DM over sixty
orders of magnitude in mass, with candidates ranging from wave-like [1] and
particle DM [2] to macroscopic objects such as primordial black holes (BHs)
[3].

Figure 1: Two plausible scenarios for the formation of BHs by imploding NSs
without significantly changing their mass. In the first scenario, DM particles
accumulate, thermalize and form a self-gravitating object, which collapses to
a BH if the number of DM particles exceeds the Chandrasekhar limit. In the
second scenario, when a sufficiently large number of DM particles accumulates,
a Bose-Einstein condensate (BEC) could form under favorable conditions and
then collapse to a BH. Once a BH is assembled at the core, it can grow by
accretion of NS matter, eventually leading to the implosion of the NS. The BH
forms over a shorter timescale through BEC formation for DM particles of mass
$m_{\chi}\leq 2\times 10^{4}$ GeV. Therefore, we show the timescales with
$m_{\chi}=2\times 10^{4}$ GeV for the channel where the BH forms without a BEC
(blue), and with $m_{\chi}=1$ GeV for the BEC channel (green), assuming a
scattering cross-section $\sigma_{\chi}=2\times 10^{-45}\rm cm^{-2}$ and
ambient dark matter density $\rho_{\chi}=1\ {\rm GeV/cm^{3}}$.

Observations of gravitational waves (GWs) by the Laser Interferometer
Gravitational-Wave Observatory (LIGO) and Virgo over the past seven years [4,
5, 6] have opened up a new avenue for exploring DM. On the one hand, detecting
BHs of unusually large masses (compared to astrophysical BHs observed until
then) could hint at their primordial origin [7]. This remains a possibility,
although several astrophysical models can account for the broad range of BH
masses detected by LIGO and Virgo (see, e.g., [8, 9]). The search for GWs from
sub-solar mass BHs has so far been unsuccessful, severely constraining the
fraction of total DM content in primordial BHs [10, 11]. If the Universe has
no primordial BHs with masses of ${\cal O}(10\,M_{\odot})$, XG detectors can
set upper limit on their abundance as a fraction of DM energy density as low
as $f_{\rm PBH}\sim{\cal O}(10^{-5})$, about two orders of magnitude lower
than current upper limits in this mass range; if instead $f_{\rm PBH}\gtrsim
10^{-4}$, future GW observations would exclude $f_{\rm PBH}=0$ at the $95\%$
credible interval [12]. The minimum testable abundance as a fraction of DM
energy density depends on the primordial BH mass, and can be as low as $f_{\rm
PBH}\sim{\cal O}(10^{-10})$ (see e.g. Fig. 5 of [13]).

An alternative way to constrain the presence of DM would be to look for its
gravitational drag on the orbits of BHs and neutron stars (NSs) [14, 15, 16,
17, 18, 19, 20, 21, 22, 23, 24]. Additionally, the presence of an axionic
cloud around BHs could extract the rotational energy from BHs, thereby
affecting the spin distribution of BHs or producing continuous GWs from newly
formed BHs [25, 26]. Several authors have explored the prospect of making such
observations [27, 28, 29, 30]. In fact, next-generation ground-based GW
observatories, with the prospect of detecting several binary black hole (BBH)
inspiral events each year with large signal-to-noise-ratios (SNRs), could
observe dozens of post-merger axionic signals [31], confirming or constraining
bosons in the mass range $\sim[7\times 10^{-14},2\times 10^{-11}]$ eV [32].

In this paper, we explore the accumulation of bosonic DM in NS cores that
could eventually form a stable, mini-BH, grow by Bondi-Hoyle accretion and
eventually lead NSs to implode and form BHs, without significantly changing
their mass. Two plausible scenarios are described in Fig. 1. These mechanisms
could be particularly efficient in regions of large DM densities, such as the
central cores of large galaxies. It has been suggested that the lack of a
sizeable population of pulsars in the core of the Milky Way, where the density
of DM is expected to be particularly high, is because most of them have
imploded to form BHs [33]. While this explanation might not be the root cause
of why the Galactic center is deficient of pulsars, future GW observations
could test if the implosion mechanism operates in NSs, as described below.

The timescale $t_{c}(\rho_{\chi},m_{\chi},\sigma_{\chi})$ over which the
accumulation of DM eventually makes NSs implode to form BHs depends on the DM
density $\rho_{\chi}$ in the neighborhood of NSs, its mass $m_{\chi}$, and its
interaction cross-section with hadrons $\sigma_{\chi}.$ NSs that live for a
time longer than $t_{c}$ will get converted to BHs, and those that live for a
shorter duration won’t. Although isolated NSs can last forever, those in a
merging binary would only live for a time $t_{d}$, called the _delay time_ ,
before they inspiral and merge to either form (rarely) supermassive NSs or
(frequently) BHs. The delay time depends on the companion masses, and the
periapsis and eccentricity of the binary at the time when it first forms. The
delay time, therefore, is not the same for all binary neutron stars (BNSs).
Instead, the NS binary population is characterized by a certain delay-time
distribution $P(t_{d})$.

The probability distribution $P(t_{d})$ is not well known, but it is often
assumed to scale like $P(t_{d})\propto 1/t_{d}$ [34, 35, 36, 37, 38]. The
fraction of the population for which $t_{d}>t_{c}$ will be converted to BHs
and the rest will remain as NSs. For the population of mergers detected with a
sufficiently large SNR, GW observations can determine the fraction of the BNS
population that has been converted to BBHs. This fraction, if nonzero, can be
used to infer the implosion timescale
$t_{c}(\rho_{\chi},m_{\chi},\sigma_{\chi})$ and hence constrain the parameter
space of local DM density, DM mass, and interaction cross-section. If the
population does not contain any BBHs, then it will be possible to set limits
on the very same quantities.

It is quite possible, although unlikely, that BHs in the mass range of NSs of
1–3 $M_{\odot}$ are produced by stellar evolution or, alternatively, they
could be primordial in origin [39, 40]. The current consensus is that massive
stars up to $\sim 23\,M_{\odot}$ leave behind NSs of masses in the range
1.2–2.0 $M_{\odot}$ at the end of their lives, while more massive stars are
likely to leave behind a BH of mass greater than about $\sim 5\,M_{\odot}$
[41, 42, 43, 44, 45]. Although the primordial Universe could produce BHs in
the mass range of NSs, they should also produce sub-solar mass BHs. A
detection of sub-solar mass BHs could hint at the early-Universe origin of BHs
with NS masses. Moreover, primordial BHs are expected to have small spin
magnitudes if they do not increase their spin by coherent accretion [46],
while BHs formed from imploding NSs could have nonzero spins (see, e.g., [47,
48]). Consequently, it might be possible to discriminate between the two
populations from their spin distributions [49]. Another proposal to
distinguish primordial BHs from astrophysical compact objects is to use their
mass distribution and the redshift evolution of the merger rates by using LVK
and A+ detections of the stochastic GW background [50], and possibly (in the
future) Cosmic Explorer and Einstein Telescope observations of the stochastic
background produced by sub-solar mass compact objects [51].

The above argument assumes that either both NSs will implode, or neither does.
If the time difference between the formation of the two NSs is large compared
to the delay time $t_{d}$, then it is possible that only one of the NSs gets
converted to a BH, but not its companion. However, this scenario is likely to
be very rare.

Figure 2: The fraction of BBH and NSBH binaries formed from the implosion of
one or both NSs, compared to the fraction of BNS systems as a function of the
implosion time-scale $t_{c}$. As a function of $t_{c}$ the NSBH fraction
remains negligibly small and can be ignored.

The complex evolutionary process leading to NS formation is not completely
understood, but stellar evolution models broadly suggest that NSs form from
progenitors whose mass $M$ lies in the range $8\,M_{\odot}\lesssim M\lesssim
23\,M_{\odot}$ [45]. The lifetime of such progenitor stars in the main-
sequence, which varies as $10(M/M_{\odot})^{-2.5}\,\rm Gyr$ [52], would be in
the range 4–55 Myr. The heavier companion would evolve through the main
sequence first to form a NS, followed by the lighter progenitor after a delay
$\Delta t_{\rm MS}.$ Assuming that the progenitors are drawn from the Salpeter
mass function, i.e., $P(m)\propto m^{-2.3}$ [53], we find that $\Delta t_{\rm
MS}$ has a median value of $\sim 14$ Myr, which is smaller than the smallest
delay time $t_{d}^{\rm min}=20\,\rm Myr$ that we will be using in this work,
and likely much smaller than the implosion timescale $t_{c},$ which could be
as large as billions of years. A specific binary with delay time $t_{d}$ will
be seen as a BBH if $t_{c}<t_{d},$ as a BNS if $t_{c}>t_{d}+\Delta t_{\rm
MS}$, or as a neutron star-black hole (NSBH) system otherwise. In Fig. 2 we
plot the three fractions, and we see that the NSBH population constitutes at
best about 6% of the total, and only over a small range of values of $t_{c}.$
Consequently, we can safely assume that either both NSs in a binary would
implode to form BHs, or neither would.

Observations of GWs could potentially discriminate the BBH population from
that of BNSs. In the final moments before a BNS coalesces, each star
experiences the tidal field of its companion, inducing a time-varying
quadrupole deformation and associated emission of GWs. This is a high-order
post-Newtonian effect – technically a fifth post-Newtonian effect or a
$(v/c)^{10}$ correction in the orbital phase evolution of the binary, where
$v$ is the orbital speed and $c$ is the speed of light [54, 55] – which
becomes important before the two NSs merge with each other, and will be absent
in the case of BBHs [56, 57, 58, 59]. Thus, by measuring the _tidal
polarizability,_ often referred to as the _effective tidal deformability_ , it
will be possible to ascertain the fraction of the two populations. It will not
be possible to measure tidal deformability with sufficiently good accuracy for
the entire observed population, but only for a fraction with sufficiently
large SNR. In this paper we explore the sensitivity of future ground-based GW
detectors to constrain the properties of a class of DM particles that can
accumulate in NS cores and cause implosion, given sufficiently large time.

The rest of the paper is organized as follows. In Sec. II, we recall how tidal
effects are encoded in GWs from BNSs. In Sec. III we introduce the various GW
detector networks and the waveform model used in this study, and compute the
accuracy with which the tidal deformability of a BNS can be measured using the
Fisher matrix formalism. This will be followed by a computation in Sec. IV of
the relative rates of BBHs and BNSs, as a function of the unknown implosion
timescale $t_{c}$ needed to form this novel population of BBHs. It turns out
that the most important hyperparameter needed in the computation of relative
rates is the implosion timescale $t_{c}.$ Thus, from the measured relative
rates we can infer the implosion timescale, as shown in Sec. V. In Sec. VI, we
discuss the physics of implosion from DM accretion and the various timescales
involved in the problem, from accumulation of DM particles to form a mini-BH,
through self-gravitation with or without the formation of a Bose-Einstein
condensate, its growth, and the final implosion of the NS (cf. Fig. 1). In
Sec. VII we derive the constraints that can be placed on DM particles if the
proposed analysis does not find a single BBH event in the mass range of 1–3
$M_{\odot}$ (which would imply that the implosion timescale is larger than the
Hubble time) and obtain the properties of bosonic DM particles assuming a
collapse time of 1 Gyr. In Sec. VIII we briefly summarize our findings, as
well as our plans to apply this technique to the known population of LIGO-
Virgo binary mergers for a number of implosion scenarios and different DM
candidates.

##  II Tidal interaction to distinguish binary neutron stars from binary black
holes

In this Section we will discuss how to distinguish BNS mergers from BBH
mergers using GW observations. An important difference in the GWs from the
coalescence of BBH and BNS systems is that waves from BNS mergers have
imprinted in them the tidal interaction between the two bodies, while BBH
mergers will have no such signature [56, 60]. Additionally, while BBH mergers
leave behind a BH remnant, BNS mergers could either promptly form a BH or
leave behind a long-lived NS remnant with neutron-rich relativistic ejecta and
a thermonuclear fireball. In this work, we will only consider the tidal
interaction between the two bodies during the adiabatic inspiral regime. A
merger accompanied by an electromagnetic afterglow essentially rules out a BBH
merger.

###  II.1 Tidal deformability

A massive body produces a tidal field. The deformation induced by the tidal
field on other bodies can be expressed as a multipole expansion, the
quadrupole being the dominant multipole. Consider a spherically symmetric NS
of radius $R$ in the tidal field ${\cal E}_{ij}$ of its companion NS. The
quadrupole deformation $Q_{ij}$ induced in the star is related to the tidal
field via the tidal deformability as

$Q_{ij}=-\lambda{\cal E}_{ij},$ (1)

where $\lambda\equiv-\frac{2}{3G}k_{2}\,R^{5},$ $R$ is the star’s radius, and
$k_{2}$ is the dimensionless tidal Love number. The tidal Love number, which
measures a body’s rigidity, depends on the equation of state (EOS) of the NS
via its compactness $C\equiv Gm/(c^{2}R),$ where $m$ is the mass of the star
[61]. For NS equations of state considered in this paper, we have $k_{2}\sim
0.1$ [62].

In a binary system of stars orbiting each other the above quadrupole
deformation is a function of time, which generates gravitational radiation,
modifying the emitted signal at the fifth post-Newtonian order, inducing a
$(v/c)^{10}$ correction to the dynamics of the system beyond the dominant
quadrupole radiation reaction. In other words, the tidal interaction
dissipates additional orbital energy into GWs, thus changing the orbital phase
evolution of the waves at $(v/c)^{10}$ order beyond the quadrupole. This
modification is significant in the final few cycles of the inspiral and
coalescence of a binary, and can be detected if the signal is observed with a
high SNR.

The tidal deformability $\lambda$ has dimensions of $\rm kg~{}m^{2}~{}s^{2}$,
but what appears in the post-Newtonian dynamics is the dimensionless tidal
deformability, defined by

$\Lambda\equiv\frac{c^{10}\lambda}{G^{4}m^{5}}=\frac{2}{3}k_{2}{\cal C}^{-5}.$
(2)

As mentioned before, the Love number $k_{2}$ generally decreases with
increasing compactness, thus the tidal deformability falls of steeper than
$C^{-5}.$ For candidate equations of state of NSs $\Lambda$ decreases with the
NS’s mass and varies over the range $\sim[100,4000]$ for NS masses in the
range 1.1-1.5 $M_{\odot}$ considered in this study (see, e.g., Fig. 1 of Ref.
[63]), the smallest values corresponding to largest NS masses and softer
equations of state, and largest values corresponding to smallest masses and
stiffer equations of state. For BHs, $\Lambda=0$ [56, 57, 58]: this is the key
to distinguishing BBH mergers from NSBH and BNS mergers [64, 65, 66].

###  II.2 Tidal signature in neutron-star binary signal

The GWs produced by BNSs are accurately described by post-Newtonian theory. We
assume NSs have negligibly small spins and are on quasi-circular orbits. These
are reasonable assumptions, as companions in Galactic double NS systems have
negligible spins (based on Ref. [67], see also Fig. 2.17 of Ref. [68]) and
gravitational radiation back reaction causes orbital eccentricity to decay
more rapidly compared to the orbital separation [69]. In the Fourier domain,
the strain amplitude $\tilde{h}(f)$ measured by an interferometric GW detector
in response to an incident BNS signal on a quasi-circular orbit is given by:

$f\tilde{h}(f)={\cal A}(f)\,e^{i\psi_{\rm PP}(f)+i\psi_{\rm Tidal}(f)},$ (3)

where $\psi_{\rm PP}$ and $\psi_{\rm Tidal}$ are contributions to the Fourier
phase from the point-particle approximation and tidal effects, respectively,
and a factor of $f$ is included to make the right-hand side dimensionless. The
amplitude ${\cal A}(f)$ and the phase $\psi_{\rm PP}(f)$ are given by [70]:

$\displaystyle{\cal A}(f)$ $\displaystyle=$
$\displaystyle\sqrt{\frac{5\nu}{6}}\frac{M}{D_{\rm eff}}\left(\pi
Mf\right)^{-1/6},$ (4) $\displaystyle\psi_{\rm PP}(f)$ $\displaystyle=$
$\displaystyle\frac{5}{128\nu}\sum_{k=-5}^{2}\left[\alpha_{k}+\alpha_{kl}\log\frac{v}{v_{0}}\right](\pi
Mf)^{k/3},$ (5) $\displaystyle D_{\rm eff}$ $\displaystyle=$
$\displaystyle\frac{4D_{L}}{\sqrt{F_{+}^{2}(1+\cos^{2}\iota)^{2}+4F^{2}_{\times}\cos\iota}}.$
(6)

Here $M\equiv m_{1}+m_{2}$ is the binary’s total mass, $m_{1}$ and $m_{2}$ are
masses of the companion stars, $F_{+}(\theta,\phi,\psi)$ and
$F_{\times}(\theta,\phi,\psi)$ are the detector antenna pattern functions
[see, e.g., [71]], $D_{L}$ is the luminosity distance to the source, $D_{\rm
eff}>D_{L}$ is the effective distance, $(\theta,\phi)$ describe the position
of the source in the sky, $\iota$ is the angle between the line sight to the
binary and the orbital angular momentum, and $\psi$ is the polarization angle.
The post-Newtonian coefficients $\alpha_{k}$ and $\alpha_{kl}$ depend on the
symmetric mass ratio $\nu\equiv m_{1}m_{2}/M^{2}$, except for the log-terms in
the post-Newtonian expansion and the tidal terms, which do depend on the total
mass. The post-Newtonian expansion is carried out in powers of $v/c,$ where
$v=(\pi Mf)^{1/3}$, and $f$ is the GW frequency. Relative to the dominant
quadrupole term, the tidal terms occur at the fifth post-Newtonian order and
higher, i.e., a $(v/c)^{10}$ effect or, equivalently, $(\pi Mf)^{10/3}$ term
in the Fourier phase in Eq. (5) relative to the dominant term [72]. The
dominant tidal contribution to the phase and the first post-Newtonian
corrections are (see e.g. [73]).

$\psi_{\rm
Tidal}(f)=-\frac{39}{2}\tilde{\Lambda}v^{10}+\left(\frac{6595}{364}\delta\tilde{\Lambda}-\frac{3115}{64}\tilde{\Lambda}\right)v^{12},$
(7)

The dominant tidal term at the fifth post-Newtonian order depends on the mass
ratio-weighted sum of the individual tidal deformabilities, defined as

$\tilde{\Lambda}=\frac{1}{26}\left[(1+12q)\Lambda_{1}+\left(1+\frac{12}{q}\right)\Lambda_{2}\right],$
(8a) while the correction at the sixth post-Newtonian order depends on
$\tilde{\Lambda}$ as well as the difference in the tidal deformabilities,
given by $\displaystyle\delta\tilde{\Lambda}$
$\displaystyle=\sqrt{1-4\nu}\left(1-\frac{13272}{1319}\nu+\frac{8944}{1319}\nu^{2}\right)\frac{(\Lambda_{2}+\Lambda_{1})}{2}$
(8b)
$\displaystyle+\left(1-\frac{15910}{1319}\nu+\frac{32850}{1319}\nu^{2}+\frac{3380}{1319}\nu^{3}\right)\frac{(\Lambda_{2}-\Lambda_{1})}{2}.$

Note that for NSs of comparable masses, i.e., $\nu\simeq 1/4,$ we have
$\Lambda_{1}\simeq\Lambda_{2},$ $\tilde{\Lambda}\simeq\Lambda_{1,2}$, and
$\delta\tilde{\Lambda}\simeq 0.$ Thus, the second term is not only a sub-
dominant post-Newtonian effect – the expansion coefficient is itself small for
most BNSs. Therefore, we neglect the sixth post-Newtonian correction in our
computations.

While NSs have non-zero tidal deformability, the tidal deformability for BHs
is zero. Tidal parameters can be inferred from GW signals, but current
measurements have large uncertainties [74, 75]. Using GW observations, one can
compute the value of $\tilde{\Lambda}$ for the binary system and infer whether
the system is a BNS or a BBH [76, 77, 66].

If the system is a BBH, it could have formed through the imploding DM channel.
This allows us to find constraints on DM properties. On the other hand, if
none of the systems are concluded to be BBH or NSBH binaries, we can still get
limits for the DM properties.

##  III Tidal deformability measurement with gravitational-wave detector
networks

We now discuss the accuracy with which the effective tidal deformability
$\tilde{\Lambda}$ can be measured using the GW signals emitted by coalescing
BNSs. We will estimate the accuracy within the Fisher information matrix
formalism [78, 79], as implemented in gwbench [80]. The two ingredients needed
for the measurement of $\tilde{\Lambda}$ are (i) a GW detector network, which
we introduce in Sec. III.1, and (ii) the waveform model used in the Fisher
matrix, which we discuss briefly in Sec. III.2. Section III.3 describes the
accuracy with which $\tilde{\Lambda}$ can be measured with the detector
networks considered in this study.

###  III.1 Detector networks

Several authors have studied the capabilities of Advanced LIGO and Advanced
Virgo in measuring the tidal deformability of NSs: for a review, see e.g.
[61]. During the second observing run, LIGO and Virgo observed the first BNS
inspiral event with a joint SNR of 33 [81] which allowed the measurement of
$\tilde{\Lambda}$ to within an accuracy of $\sigma_{\tilde{\Lambda}}=630$ at
90% credible interval [82], with some authors ruling out the possibility that
this was a BBH [74], especially when combined with optical and infrared
observations [83]. The second BNS merger event, GW190425, was observed with a
far lower SNR of 12.9, and it did not allow to place any meaningful bounds on
the tidal deformability of NSs [84]. The upcoming year-long fourth observing
run of the LIGO, Virgo and KAGRA detectors is expected to detect a handful of
BNS mergers but at the current sensitivity this network won’t measure
$\tilde{\Lambda}$ with an accuracy good enough to conclusively say that tidal
effects are absent. We will therefore restrict ourselves to future upgrades of
LIGO, Virgo, KAGRA and LIGO-Aundh,111LIGO-Aundh is the preferred name for the
new LIGO-India observatory, coming up near the town of Aundh in central India.
as well as next generation (XG) observatories such as Cosmic Explorer [85] and
the Einstein Telescope [86].

More precisely, we consider three ground-based detector networks—A+, Voyager
and XG—to determine the accuracy of measuring the effective tidal
deformability, with particular interest in the XG network, since these
observatories will have the required sensitivities to obtain an informative
estimate of the tidal deformability parameter for BNSs.

  * •

The A+ network comprises five detectors: LIGO-Hanford, LIGO-Livingston, Virgo,
KAGRA and LIGO-Aundh at A+ sensitivity [87].

  * •

The Voyager network consists of LIGO-Hanford, LIGO-Livingston and LIGO-Aundh
at Voyager sensitivity [88], with Virgo and KAGRA at A+ sensitivities.

  * •

Finally, the XG network includes the Einstein Telescope, one Cosmic Explorer
in the US, and another Cosmic Explorer in Australia, as in [80].

Table 1: Parameter space of BNSs used in this study. In addition to the detector-frame companion masses and the tidal deformability of NSs, there are four angles [two describing the orientation of the binary’s orbit relative to the detector frame $(\iota,\,\psi)$ and two for the sky position of the source $(\alpha,\,\delta)$], the luminosity distance (or, equivalently, redshift $z$), a fiducial ‘arrival time’ $t_{C}$ when the strain amplitude of the signal is largest, and the phase of the signal at that time $\phi_{C}$. Component mass, ${m_{1}}$ and ${m_{2}}$ | $\mbox{[1, 2]}\,M_{\odot}$
---|---
Tidal Parameters222These parameters don’t directly enter the waveform, but only the effective tidal deformability $\tilde{\Lambda}.$, $\Lambda_{1}$ and $\Lambda_{2}$ | 2.0
Effective tidal parameter, $\tilde{\Lambda}$ | $(6q^{2}+q+6)\Lambda_{1}/(13q)$
Right ascension, $\alpha$ | [0, 2$\pi)$
Declination, $\delta$ | $[-\pi/2,\,\pi/2]$
Inclination, $\iota$ | [0, $\pi]$
Polarization, $\psi$ | [0, $2\pi]$
Redshift, ${z}$ | [0, 10]
Fiducial time of arrival, $t_{C}$ | 0
Constant phase offset, $\phi_{C}$ | 0
Figure 3: These plots show $\sigma_{\tilde{\Lambda}}^{90\%},$ the 90%
confidence interval in the accuracy of measurement of $\tilde{\Lambda},$ for
A+ (left panel), Voyager (middle panel) and XG (right panel), as a function of
$\tilde{\Lambda}$ and the SNR (color bar). The measurement accuracy is an
order of magnitude better for the XG network compared to the A+ network. In XG
detectors, half of all events within $z=10$ have their tidal deformability
constrained to within $\sigma^{90\%}_{\tilde{\Lambda}}\lesssim 100,$ an order
of magnitude larger for than A+.

###  III.2 Binary NS waveform model

As discussed in Sec. II, the signature of tidal deformation is imprinted in
the GWs emitted by an inspiralling BNS system. The simplest modification of
the point-particle post-Newtonian model in the frequency domain – cf. Eq. (5)
– captures the essence of the tidal effects, but it is not in agreement with
waveforms obtained from numerical simulations of BNS mergers. We adopt the
IMRPhenomD_NRTidalv2 waveform model, in which the Fourier amplitude and phase
are given algebraically in closed form [89, 90, 91]. This choice improves the
accuracy in the calculation of derivatives of the waveform with respect to the
various parameters, which are needed to compute the Fisher matrix.

The tidal terms in this model are valid at frequencies as low as $\sim 50$–100
Hz where Cosmic Explorer and Einstein Telescope have good sensitivity. The
model does not incorporate the post-merger signal that could be present if the
remnant is a long-lived NS, nor dynamical tides, which could also be important
[92, 93, 94, 95, 96, 97, 98, 99, 100, 101], although sub-dominant compared to
the static tides that are included in the model.

The waveform model is described by ten parameters (see Table 1): the companion
masses $(m_{1},m_{2}),$ two angles describing the position of the source in
the sky $(\alpha,\delta)$, two angles describing the inclination of the
binary’s orbit and the wave polarization $(\iota,\psi),$ the luminosity
distance of the source $D_{L}$ (or, equivalently, its redshift $z$), the
effective tidal parameter $\tilde{\Lambda}$, the coalescence time $t_{C}$ and
the coalescence phase $\phi_{C}$.

In the Fisher matrix approximation, which involves derivatives with respect to
the parameters of the waveform, the absolute error in the measurement of
$\tilde{\Lambda}$ at leading post-Newtonian order is the same for all values
of $\tilde{\Lambda}.$ This is because the effective tidal parameter appears
linearly in the waveform model and we have neglected the sub-dominant tidal
terms. Thus, the error $\sigma_{\tilde{\Lambda}}$ is determined by the
correlation of $\tilde{\Lambda}$ with the other waveform parameters and the
loudness of the signal and not any particular value of $\tilde{\Lambda}.$ To
compute the Fisher matrix, without loss of generality we use a tidal
deformability of $\Lambda_{\rm 1,2}=2$, because a nonzero tidal deformability
value is required to use IMRPhenomD_NRTidalv2.

The NS mass distribution is not known very well, but accurate measurement of
NSs in radio binary pulsars seems to suggest that they are normally
distributed with a standard deviation that’s small compared to the mean (see
[102, 103]). The two BNS mergers discovered by LIGO and Virgo, in particular
GW190425 [84], already indicate that NS masses in merging binaries could be
different from those of galactic BNSs. There is currently no concrete
distribution that we could use from GW measurements (see, however, [104]) and
hence we draw the companion masses from a Gaussian distribution with a mean
value of $1.3\ M_{\odot}$ and a standard deviation of $0.09\,M_{\odot}$ [103],
but with the constraint that $1\leq m_{1},m_{2}\leq 2M_{\odot}.$

If $\Lambda_{1}=\Lambda_{2}=\Lambda,$ then Eq. (8a) leads to

$\tilde{\Lambda}=\frac{6q^{2}+q+6}{13q}\Lambda.$

This implies that for $0.5\leq q\leq 1,$ which is the range of $q=m_{2}/m_{1}$
allowed by the Gaussian distribution and the hard upper- and lower-cutoff of
NS masses, $2\leq\tilde{\Lambda}\lesssim 2.5.$ The range of $\tilde{\Lambda}$
is a proxy for the mass ratio in our sample, and the variation in the
measurement error in $\tilde{\Lambda}$ is determined largely by its
correlation with other parameters and by the SNR.

###  III.3 Measurement accuracy of $\tilde{\Lambda}$

Bayesian inference is the preferred method to estimate the error in the
measurement of parameters, but the long Markov chains needed for the
convergence of the posterior distribution are expensive and time consuming.
For the exploratory work carried out in this paper it suffices to employ the
faster Fisher matrix approach. In this approach, one first computes the Fisher
information matrix $\Gamma_{mn}$ defined by

$\Gamma_{mn}\equiv\left<\partial_{m}\tilde{h},\,\partial_{n}\tilde{h}\right>,$
(9)

where $\partial_{m}$ denotes the derivative of the waveform with respect to
the parameter $\lambda_{m}$ of the waveform, $\tilde{h}(f)$ is the Fourier
transform of the detector response $h(t)\equiv
F_{+}h_{+}+F_{\times}h_{\times}$, $(h_{+},\,h_{\times})$ are the two
polarization strain amplitudes, and $\left<a,\,b\right>$ is the scalar product
of waveforms $a$ and $b$, defined as

$\left<a,\,b\right>\equiv
4\Re\int_{f_{L}}^{f_{H}}a(f)\,b^{*}(f)\frac{df}{S_{h}(f)}.$ (10)

Here $S_{h}(f)$ is the noise power spectral density of the detector in
question, and $f_{L}$ and $f_{H}$ are suitably chosen lower and upper
frequency cutoffs. We choose $f_{L}=20$ Hz for the A+ and Voyager networks,
and $f_{L}=5$ Hz in the case of the XG network. The upper frequency cutoff is
chosen to be the Nyquist frequency with a sampling rate of 4096 Hz.

The information matrix of a detector network is just the sum of the individual
information matrices: $\Gamma_{mn}\equiv\sum_{A}\Gamma^{A}_{mn}.$ The
covariance matrix is the inverse of $\Gamma:$ $C_{mn}=(\Gamma^{-1})_{mn}.$ By
definition the Fisher matrix is symmetric, and so is the covariance matrix.
Its diagonal elements $C_{mm}$ are the variances in the inference of
parameters $\lambda_{m}$, and the off-diagonal elements $C_{mn}$ are the
covariances in parameters $\lambda_{m}$ and $\lambda_{n}.$

We perform the analysis for a population of BBH sources with component masses
$m_{1},m_{2}\in$ [$1\,M_{\odot}$, $2\,M_{\odot}$], distributed in redshift as
described in Sec. IV, up to a maximum redshift of $z=10.$ Table 1 lists the
parameters of the binary population considered in this study. The errors on
the effective tidal deformability $\sigma_{\tilde{\Lambda}}$ for this
population were computed using the gwbench toolkit [80], which performs Fisher
analysis to provide measurement errors on GW parameters given a network of GW
detectors. This provides the distribution of inferred values of
$\tilde{\Lambda}$ for the cosmic population of sources. We use
$\sigma_{\tilde{\Lambda}}^{90\%}$ as the criterion to differentiate between
the population of BBHs and BNSs in this mass-range because $\tilde{\Lambda}=0$
for BHs. Therefore, the confidence with which we can classify a binary as a
BBH or BNS is inversely proportional to the measurement error, given some
$\tilde{\Lambda}$.

The 90% confidence interval in the measurement of the effective tidal
deformability for the full population, $\sigma_{\tilde{\Lambda}}^{90\%}$, is
shown in Fig. 3 for the three networks considered in this study, with the
color representing the SNR of the events. In Fig. 4 we plot the cumulative
distribution of $\sigma_{\tilde{\Lambda}}^{90\%}$ for the three networks. The
A+ and Voyager networks can determine the tidal deformability to within
$\sigma_{\tilde{\Lambda}}^{90\%}\simeq 100$ for 0.05% and 0.5% of the
population. The XG network, on the other hand, can determine $\tilde{\Lambda}$
to the same accuracy for 30% of the full population. This is a good enough
accuracy to distinguish BNSs from BBHs if the preferred EOS is stiff, such as
ALF2, producing larger NS radii and greater tidal deformabilities (e.g.,
$\tilde{\Lambda}\sim{\color[rgb]{1,0,0}300}$ for a $1.4\,M_{\odot}$ NS, which
is true for most NS masses considered in our study). Smaller errors
$\sigma^{90\%}_{\tilde{\Lambda}}\lesssim{\color[rgb]{1,0,0}20}$ would be
required if NSs are described by a softer EOS, such as APR4, with tidal
deformabilities $\tilde{\Lambda}\sim{\color[rgb]{1,0,0}50}.$ A greater
measurement accuracy (i.e., smaller values of
$\sigma^{90\%}_{\tilde{\Lambda}}$) requires louder signals, which means fewer
systems can be classified as belonging to one of the two classes. For example,
only the XG network can measure $\tilde{\Lambda}$ to better than
$\sigma_{\tilde{\Lambda}}^{90\%}<20.$ Thus, the measurement accuracy of
effective tidal deformability directly impacts how well we can determine the
collapse time discussed in Sec. V.

Figure 4: Cumulative distribution of the error in the measurement of effective
tidal deformability in the three detector networks considered in this study.
Notice that the fraction of events detected with an error
$\sigma_{\tilde{\Lambda}}^{90\%}\simeq 100$ is 0.05%, 0.5% and 30% in the A+,
Voyager and XG networks, respectively.

##  IV Merger rates

In this section we discuss how to deduce constraints on DM mass and
interaction cross-section based on the observed merger rate of BNS and BBH
systems in the NS mass range. The constraints follow by comparing the time-
scale for conversion of NSs to BHs by DM accumulation, $t_{c},$ to the time-
scale for coalescence of NSs by gravitational radiation backreaction, $t_{d}$.
The timescales in the problem suggest that either both or neither of the NSs
in a binary will be converted to BHs before they inspiral and merge, so we
will not consider NSBH binaries, but we will derive an equation that relates
an upper limit on the merger rate of BBHs to the properties of DM particles.

Figure 5: Left: Merger rate density of BBHs formed by implosion as a function
of redshift for several values of the collapse time $t_{c}$ from 0.02 Gyr to
10 Gyr, as given in the legend (in Gyr). Right: The number of BNS and BBH
mergers per year $N=\dot{N}T,$ $T=1\,\rm yr,$ in Earth’s frame as a function
of the collapse time $t_{c}.$ Not all of these would be observable by a
detector network, but only those that above a certain SNR threshold. The
shaded regions in both panels correspond to the 90% credible interval for the
local BNS merger rate for the multi-source model in [105].

###  IV.1 Collapse timescales and merger time delays

As we shall discuss in Sec. VI, the time-scale $t_{c}$ for NSs to implode to
form BHs largely depends on the properties of the DM particles: (i) the
interaction cross section $\sigma_{\chi},$ (ii) the DM density at the site of
BNSs and their dispersion velocity (which are both determined by the location
of the binary within a galaxy, being larger at the galactic core and smaller
in the halo), and (iii) the mass of the DM particles $m_{\chi}.$ On the other
hand, the time-scale $t_{d}$ for NSs to coalesce depends on (i) the
eccentricity and semi-major axis when the BNS first forms, and (ii) the masses
of the two NSs [69].

The two NSs in a binary do not form from their stellar progenitors at the same
time. The delay in the formation of the second NS with respect to the first
could be substantial if the masses of the parent stars are very different.
However, the timescale arguments in Sec. I imply that binaries with component
masses $1\,M_{\odot}\leq m\leq 2\,M_{\odot}$ will either be BNS or BBH
binaries, and not mixed (NSBH) binaries, so we will ignore mixed binaries from
now on.

###  IV.2 Binary neutron star merger rate

LIGO and Virgo have so far observed two BNS mergers: GW170817 at a distance of
$40^{+8}_{-14}$ Mpc [81] and GW190425 at a distance of $159^{+69}_{-72}$ Mpc
[84], both at 90% credible interval (CI). Since they are both at very low
redshift, the merger rate determined from them is essentially the local (i.e.,
$z=0$) rate. The local rate for BNSs inferred from the third Gravitational
Wave Transient Catalog-3 GWTC-3 [6] is $R_{0}=660_{-530}^{+1040}{\rm
Gpc^{-3}\,yr^{-1}}$ at 90% CI [105], under the assumption of a multi-source
model including BNS, NSBH, and BBH sub-populations.

The merger rate evolves with redshift because (i) the star formation rate
varies as a function of redshift; (ii) BNSs that form at a certain redshift
don’t merge immediately but only after a delay time $t_{d}$, and hence at a
different redshift; and (iii) the metallicity evolves with redshift, affecting
the mass function and formation rate of compact binaries. In this study, we
will ignore the effect of metallicity, as it plays a greater role in the case
of BBHs and is less likely to affect the merger rate of BNSs [34, 106, 107,
108]. We will assume the star formation rate $\psi(z)$ (SFR) given by Ref.
[109]:

$\psi(z|\alpha_{\rm{F}},\beta_{\rm{F}},C_{\rm{F}})\propto\frac{(1+z)^{\alpha_{\rm{F}}}}{1+\left(\frac{1+z}{C_{\rm{F}}}\right)^{\beta_{\rm{F}}}}\,,$
(11)

with $(\alpha_{\rm{F}},\beta_{\rm{F}},C_{\rm{F}})=(2.7,5.6,2.9)$. We assume
that the merger rate is the same as the SFR except that binaries that form at
redshift $z_{f}$ merge at redshift $z$ after a delay $t_{d},$ with a
corresponding redshift interval $\Delta z=z_{f}-z.$ Given a redshift $z$ at
which merger rate is required and the time delay $t_{d},$ the redshift $z_{f}$
at which the binary forms can be found by solving

$t_{d}=\frac{1}{H_{0}}\int_{z}^{z_{f}}\frac{dz^{\prime}}{(1+z^{\prime})E(z^{\prime})},\quad
E(z)=\Omega_{\Lambda}+\Omega_{M}(1+z)^{3},$ (12)

where $\Omega_{M}$ and $\Omega_{\Lambda}$ are the DM and dark energy
densities, respectively, and we have assumed a flat Universe in which dark
energy is interpreted as a cosmological constant [110]. Now, the merger rate
density $\dot{n}(z)$ as a function of redshift can be computed by integrating
the SFR over all delay times, the delay-time probability function $P(t_{d})$
serving as a weighting factor:

$\displaystyle\dot{n}(z)=A\int_{t_{d}^{\rm min}}^{t_{d}^{\rm
max}}\psi(z_{f}(z,t_{d}))\mathcal{P}(t_{d})\,dt_{d},$ (13)

where $z_{f}$ is found using Eq. (12). The delay-time distribution is not very
well known, and we take it to be the Jeffrey’s prior
$\mathcal{P}(t_{d})\propto t_{d}^{-1}$ extending over the range $[t_{d}^{\rm
min},\,t_{d}^{\rm max}],$ with the normalization

$\int_{t_{d}^{\rm min}}^{t_{d}^{\rm max}}P(t_{d})\,dt_{d}=1.$

The constant $A$ in Eq. (13) is chosen such that $\dot{n}(0)=R_{0}:$

$A=\frac{R_{0}}{\int_{t_{d}^{\rm min}}^{t_{d}^{\rm
max}}\psi(z_{f}(0,t_{d}))P(t_{d})\,dt_{d}}.$ (14)

Figure 6: The number of BBH mergers as a function of collapse time $t_{c}$ for
several values of the threshold $\sigma_{\tilde{\Lambda}_{T}}$ on the
estimated error in $\tilde{\Lambda},$ for 5 years of observing time.

If NSs in a binary implode to form BHs within a collapse time $t_{c}$ that
lies between the minimum and maximum delay times, i.e.,
$t_{d}^{\rm{min}}<t_{c}<t_{d}^{\rm{max}}$, then the total merger rate density
in Eq. (13) is the sum of the BNS and BBH merger rates, i.e., $\dot{n}_{\rm
total}=\dot{n}_{\rm BNS}+\dot{n}_{\rm BBH}$, where

$\displaystyle\dot{n}(z)_{\rm BNS}$ $\displaystyle=$ $\displaystyle
A\int_{t_{d}^{\rm
min}}^{t_{c}}\psi(z_{f}(z,t_{d}))\mathcal{P}(t_{d})\,dt_{d},$ (15)
$\displaystyle\dot{n}(z)_{\rm BBH}$ $\displaystyle=$ $\displaystyle
A\int_{t_{c}}^{t_{d}^{\rm
max}}\psi(z_{f}(z,t_{d}))\mathcal{P}(t_{d})\,dt_{d}.$ (16)

The left panel of Fig. 5 shows the merger rate of BBHs formed by implosion as
a function of redshift for different choices of the collapse time $t_{c}$ (in
Gyr). The shaded region for each choice of $t_{c}$ represents the uncertainty
in the local merger rate of BNS found using the latest GW catalog, which is in
the range $130\leq R_{0}\leq 1700\,{\rm Gpc^{-3}\,yr^{-1}}$ [105]. If the
collapse time is less than the smallest delay time, then all BNSs are
converted to BBHs, and the entire BNS population will be observed as BBHs. If,
on the other hand, the collapse time is larger than the largest delay time, no
BNS is converted to BBH. The observed population will be a mixture of BNS and
BBH for values of the collapse time in between these extremes.

The merger rate $d\dot{N}$ in a cosmological volume $dV_{c}$ (redshift range
$dz$), as measured by an observer at $z=0,$ is
$\dot{N}=\dot{n}(z)\,dV_{c}/(1+z),$ where a factor $(1+z)$ accounts for the
time dilation between observers at redshift $z$ and redshift $0.$ Thus, the
merger rate within some redshift $z$ as measured by an observer at $z=0$ is
given by

$\dot{N}=\int_{0}^{z}\frac{\dot{n}(z^{\prime})}{1+z^{\prime}}\frac{dV_{c}}{dz^{\prime}}\,dz^{\prime}.$
(17)

In our simulations we distribute sources as a function of redshift using the
above equation. The right panel of Fig. 5 shows the number of BNS and BBH
mergers in the Universe per year $N=\dot{N}T,$ with $T=1\,\rm yr,$ in the
Earth’s frame, as a function of the collapse time $t_{c}.$ If
$t_{c}<t_{d}^{\rm min}$ every BNS will be converted to a BBH, while
$t_{c}>t_{d}^{\rm max}$ would imply none will be converted. At intermediate
values of $t_{c}$ we expect to observe, among all mergers, some that are BBHs.

A detector network only observes those mergers that stand above the detector
noise at an acceptably low false alarm rate (say, one false event per year).
Thus, the merger rate _observable_ by a detector network, $\dot{N}_{\rm obs}$,
is

$\dot{N}_{\rm
obs}=\int_{0}^{z}\frac{\dot{n}(z^{\prime})}{1+z^{\prime}}\frac{dV_{c}}{dz^{\prime}}\epsilon(z^{\prime})dz^{\prime}.$
(18)

Here $\epsilon(z)$ is the detection efficiency of a network as a function of
redshift, to be discussed below. It measures the fraction of events detected
by a given detector network subject to one or more conditions. At the outset,
we require an event’s SNR to be larger than a preset threshold $\rho_{T}$ to
make a high-confidence detection. For the current network of LIGO/Virgo/KAGRA
detectors, a SNR threshold of $\rho_{T}=10$ assures that the false alarm rate
is no more than a few per year.

In this work we not only wish to make high-confidence detections, but
additionally select only those events for which the effective tidal
deformability can be measured with a good accuracy—good enough to distinguish
BBH mergers (for which $\tilde{\Lambda}=0$) from BNS mergers (for which
$\tilde{\Lambda}>\tilde{\Lambda}^{\rm NS}_{\rm min}\gg 0$). To this end, we
require the 90% credible interval $\sigma_{\tilde{\Lambda}}^{90\%}$ in the
measurement of $\tilde{\Lambda}$ to be less than a preset threshold
$\sigma_{\tilde{\Lambda}_{T}}$, i.e.,
$\sigma_{\tilde{\Lambda}}^{90\%}<\sigma_{\tilde{\Lambda}_{T}}$, and choose the
threshold so that it is (significantly) smaller than the smallest value of the
effective tidal deformability of BNSs:

$\sigma_{\tilde{\Lambda}_{T}}\ll\tilde{\Lambda}^{\rm NS}_{\rm
min}\equiv\min_{m_{1},m_{2}}\tilde{\Lambda}(m_{1},m_{2}),$

where $m_{1}$ and $m_{2}$ are the companion masses. This condition guarantees
that a merger is correctly classified as a BNS or a BBH merger with high
confidence. Typically, NSs with the largest masses and softest equations of
state have the smallest tidal deformability. For equations of state that are
still viable, the smallest value of $\tilde{\Lambda}$ varies over the range
$\tilde{\Lambda}\in[10,200]$ [63]. We present our results for several values
of the threshold $\sigma_{\tilde{\Lambda}_{T}}$ restricting the value to be
less than 100.

With the conditions on SNR and $\tilde{\Lambda}$ imposed, the efficiency of a
detector network is given by:

$\epsilon(z)=\frac{1}{N}\sum^{N}_{i=1}\Pi\left(\frac{\rho}{\rho_{T}}-1|z\right)\Pi\left(\frac{\sigma_{\tilde{\Lambda}_{T}}}{\sigma_{\tilde{\Lambda}}}-1|z\right),$
(19)

where the sum is over the full population of BNSs and $\Pi(x)$ is the step
function: $\Pi(x)=0$ if $x\leq 0$, and $\Pi(x)=1$ if $x>0$.

Figure 7: Left: Gaussian distribution for the observed number of BBHs for
$\sigma_{{\tilde{\Lambda}_{T}}}^{90\%}=100$ and $t_{c}=1~{}{\rm Gyr}$. Right:
Inferred distribution of $t_{c}$ from $N_{\rm BBH}$ for the three detector
networks under consideration and an observing period of 5 years. The dashed
lines show the 90% credible intervals for the distributions of $t_{c}$ from
$N_{\rm BBH}$.

The number of BBH mergers observed up to some redshift, $z$, in an observing
period $T,$ is given by

$N_{\rm BBH}=\dot{N}_{\rm BBH}T=T\int_{0}^{z}\frac{\dot{n}(z^{\prime})_{\rm
BBH}}{1+z^{\prime}}\,\frac{dV_{c}}{dz^{\prime}}\,\epsilon(z^{\prime})dz^{\prime},$
(20)

and similarly for BNSs. Therefore, the number of BBHs we expect to observe in
this mass spectrum not only depends on the sensitivity of GW detectors and
searches, but also on how well we can measure the tidal effects from these
observations. A single detection with some assumed measurement efficiency for
$\tilde{\Lambda}$ provides a limit on $t_{c}$ which can be used to infer DM
properties using Eqs. (33) and (39).

Figure 6 shows the number of BBH mergers _observable_ by a detector network
over a five-year duration, $N_{\rm BBH}=\dot{N}_{\rm BBH}T$ with $T=5\,\rm
yr,$ as a function of the collapse time $t_{c}$ for different choices of the
threshold $\sigma_{\tilde{\Lambda}_{T}}$ (from 10 to 100) and for the A+,
Voyager and XG detector networks described in Sec III. The imminent upgrade of
LIGO and Virgo detectors could identify tens to hundreds of of BBHs if the
collapse time is in the range of 100 Myr to 1 Gyr over an observational period
of 5 years but this network is more likely to acquire about a year’s worth of
data. The Voyager network increases these numbers by an order of magnitude
while the XG network of Einstein Telescope and Cosmic Explorer will observe
several orders of magnitude more BH binaries compared to the other networks.
As we shall see later, the larger fraction of systems that can be clearly
identified as BBH helps in placing a tighter constraints on DM mass and
interaction cross section.

##  V Inference of collapse time

In this section we elucidate how to infer the average time it takes for a NS
to collapse to a BH due to accretion of DM. Let us first note a caveat in our
argument: what GWs _can_ infer is the relative abundance of BBHs versus BNSs.
If BHs with companion masses in the 1–2 $M_{\odot}$ range form by unknown
astrophysical processes or in the primordial Universe they will be part of the
BH population in this mass range, and will be indistinguishable from the
population that formed from the implosion of NSs due to accumulation of DM.
First, we assume that stellar evolution cannot produce BHs in this mass range.
Second, if primordial BHs exist then there is no fundamental reason they
should only appear in this mass range. In particular, we expect sub-solar mass
BBHs to exist as well [111, 112]. In the absence of such a population we can
be fairly confident (although not certain) that stellar mass BHs in the 1–2
$M_{\odot}$ range are not of primordial origin. At present, there is no
preference for any of these scenarios, so we will proceed with the assumption
that any detections of BBHs with BNS masses formed by the implosion of NSs.

Under this assumption, it is straightforward to deduce the collapse time from
the observed population of BBHs and BNSs. For illustration, we assume that the
true collapse time is $t_{c}=1\,\rm Gyr.$ Given $t_{c},$ Eq. (16) gives the
local merger rate density, which could then be used in Eq. (20) to compute the
number of BBHs expected to be detected over an observational period $T$ in a
detector network whose detection efficiency is $\epsilon(z)$, as defined in
Eq. (19). Since the rate is Poisson distributed, the number of detections
$N_{\rm BBH}$ would be uncertain by $\sqrt{N_{\rm BBH}}$, and the relative
error scales like $1/\sqrt{N_{\rm BBH}}.$

The left panel of Fig. 7 shows the expected number of BBHs in different
detector networks, together with the rate uncertainty in our simulation. The
relative uncertainty will, obviously, be larger for less sensitive detectors,
and this impacts how well the collapse time can be deduced. The right panel of
Fig. 7 shows the collapse time deduced from the rate posterior plotted on the
left. In reality, we would determine the number of BBH mergers given the
collapse time and network efficiency [i.e., $N_{\rm BBH}(t_{c},\epsilon(z))]$
and interpolate this function to determine $t_{c}$ from $N_{\rm BBH}$ for a
given network efficiency $\epsilon(z).$ It is clear that XG detectors will be
able to constrain the collapse time far better than the A+ network.

##  VI Implosion of neutron stars by accumulation of dark matter

The inferred collapse time of NSs into BHs from GW observations is model-
agnostic. These limits can constrain the particle properties of DM in
scenarios where DM particles get captured in NS cores. In this section, we
illustrate how we can constrain the asymmetric DM scenario under certain
simplified assumptions using previously derived limits on the collapse time.

In the asymmetric DM scenario, DM particles do not self-annihilate due to the
assumed asymmetry between the number density of particles and anti-particles.
Therefore, the capture and accumulation of DM particles in the core of NSs
could lead to the formation of a BH at the NS core, leading to the potential
implosion of the host star. Here we examine this scenario in the case of non-
interacting bosonic DM, because the bosonic Chandrasekhar limit is much
greater than that for fermions, implying that bosonic DM undergoes
gravitational collapse sooner than fermionic DM. Indeed, the number of
fermionic and bosonic particles corresponding to the Chandrasekhar mass for a
DM particle of mass $m_{\chi}$ are given, respectively, by:

$\displaystyle N^{\rm fermionic}_{\rm Chandra}\simeq 1.8\times
10^{57}\left(\frac{{\rm GeV}}{m_{\chi}}\right)^{3},$ (21) $\displaystyle
N^{\rm bosonic}_{\rm Chandra}\simeq 1.5\times 10^{38}\left(\frac{{\rm
GeV}}{m_{\chi}}\right)^{2}.$ (22)

We also note that due to the dispersion velocities of DM, the NS gravity alone
cannot capture DM particles, and some dissipative mechanism involving the
interaction of DM with hadrons would be needed. For the species of DM
considered in this paper, DM particles are assumed to interact with hadrons
through the weak interaction. However, the accumulation of DM particles over
the lifetime of the NS for the accretion rate considered in this paper will
not significantly increase its mass.

###  VI.1 Dark matter capture by neutron stars

The ambient DM attracted by the NS’s gravity is captured if its trajectory
intersects the star, and it loses energy through its interactions with baryons
and leptons. The capture rate of DM by gravitating bodies such as the earth
and the sun was first computed in [113] and a general analytic theory was
developed in [114]. These calculations have been refined to include the
effects of matter degeneracy in estimating the capture rates of NSs [115,
116]. Most recently, the calculations have been further improved to include
general relativistic corrections and provide a consistent treatment of several
different operator structures that define the interaction between nucleons and
DM [117].

The capture rate, in general, depends on the ambient energy density of DM,
$\rho_{\chi}$; the DM scattering cross-section off targets in the NSs,
$\sigma$; the number density of targets, $n_{\rm t}(r)$; the escape velocity
in the NS at radius $r$, $v(r)$; the DM velocity dispersion far away from the
star, $v_{\chi}$; and the velocity of the NS, $v_{*}$. In what follows, we
adopt the result derived in [117] to obtain the total capture rate for DM in
the mass range $1-10^{6}$ GeV and for a constant cross-section $\sigma$. In
this mass range, Pauli blocking suppression of scattering is unimportant and
the total capture rate is

$C_{\rm NS}=C_{\rm geom}~{}{\rm Min}\left[\frac{\sigma}{\sigma_{\rm
th}},1\right]\,,\ $ (23)

where

$C_{\rm geom}=\pi R_{*}^{2}\frac{\rho_{\chi}}{m_{\chi}}v_{*}~{}{\rm
Erf}\left[\sqrt{\frac{3}{2}}\frac{v_{*}}{v_{\chi}}\right]\frac{(v^{2}_{\rm
esc}/v_{*}^{2})}{1-(v^{2}_{\rm esc}/c^{2})}$ (24)

is the geometric capture rate that includes the effect of gravitational
focusing in general relativity, and

$\sigma_{\rm th}=\frac{\pi R_{*}^{2}}{N_{\rm t}~{}\xi}$ (25)

is the threshold value of the cross-section required to ensure that the DM
particles traversing the NS have at least one collision. Here, $N_{\rm t}$ is
the total number of target particles in the NS, and

$\xi=\frac{4\pi}{N_{\rm t}}\int_{0}^{R_{*}}r^{2}dr\sqrt{g_{rr}(t)}n_{\rm
t}(r)\frac{1-g_{tt}(r)}{1-g_{tt}(R_{*})}\frac{g_{tt}(R_{*})}{g_{tt}(r)}\,,$
(26)

where $g_{tt}$ and $g_{rr}$ are temporal and spatial components of the
Schwarzschild metric. The dimensionless number $\xi\approx{\cal O}(1)$, and
its precise value depends on the structure of the NS and the density profile
of targets in the NS interior.

For a typical NS with mass $M_{*}\simeq 1.4~{}M_{\odot}$, radius $R_{*}=12$
km, and total baryon number $N_{t}\simeq 2\times 10^{57}$, assuming that
$\xi\approx 1$, the threshold cross-section $\sigma_{\rm th}\approx 2\times
10^{-45}~{}{\rm cm}^{2}$. For this case, the accumulation of DM in a NS at
rest, i.e., $v_{*}=0$, is given by

$N_{\chi}(t)=\sqrt{6\pi}R_{*}^{2}\frac{\rho_{\chi}}{m_{\chi}}\frac{(v^{2}_{\rm
esc}/v_{\chi})}{1-(v^{2}_{\rm esc}/c^{2})}{\rm
Min}\left[\frac{\sigma}{\sigma_{\rm th}},1\right]~{}t\,.$ (27)

Using $v_{\rm esc}=\sqrt{2GM_{*}/R_{*}}=0.6c$, and a DM velocity dispersion
$v_{\chi}=220$ km/s, we obtain

$N_{\chi}(t)=2.9\times 10^{43}\frac{\rho_{\chi}}{\rho_{0}}\left(\frac{{\rm
GeV}}{m_{\chi}}\right){\rm Min}\left[\frac{\sigma}{\sigma_{\rm
th}},1\right]~{}\frac{t}{10^{10}~{}{\rm yr}}\,,$ (28)

where $\rho_{0}=1$ GeV/cm3 is our choice for the fiducial DM density.

###  VI.2 Thermalization and formation of a
mini-black hole

The captured DM particles continue to scatter and eventually thermalize with
NS matter. A general analysis of thermalization is challenging because the DM
energy decreases over several orders of magnitude during this process. At low
energy, matter degeneracy and correlations due to strong interactions between
baryons can substantially alter the scattering rate [118, 119]. For scattering
off neutrons with a constant cross-section, the thermalization time for the DM
masses of interest is

$t_{\rm th}\simeq 3750\
\text{years}\frac{\gamma}{(1+\gamma)^{2}}\left(\frac{2\times
10^{-45}\rm{cm^{2}}}{\sigma_{\rm\chi}}\right)\left(\frac{10^{5}\rm K}{\rm
T}\right)\,,$ (29)

where $\gamma=m_{\chi}/m_{\rm n}$, and $m_{\rm n}$ is the neutron mass [118].

Upon thermalization, the radius of the DM sphere is determined by the
temperature and the gravitational potential of the NS [118]

$r_{\rm{th}}\approx 2.2\ {\rm m}\left(\frac{{\rm T}}{10^{5}\ {\rm
K}}\right)^{1/2}\left(\frac{{1~{}\rm GeV}}{m_{\chi}}\right)^{1/2}\,.$ (30)

When the mass density of bosonic DM exceeds that of baryons, the DM becomes
self-gravitating. This occurs when the total number of DM particles exceeds

$N_{\rm{self}}\simeq 4.8\times
10^{46}\left(\frac{\rm{GeV}}{m_{\chi}}\right)^{5/2}\left(\frac{\rm{T}}{10^{5}\
\rm{K}}\right)^{3/2}\,.$ (31)

The self-gravitating DM sphere collapses to form a BH when its mass exceeds
the Chandrasekhar limit for bosonic matter, i.e., if $N_{\rm{self}}>N_{\rm
Chandra}^{\rm bosonic}$. Using Eq. (22), we can deduce that $N_{\rm
self}>N_{\rm Chandra}^{\rm bosonic}$ for $m_{\chi}\leq 10^{17}~{}\rm
GeV(T/10^{5}~{}\rm K)^{3}$. Thus, for the mass range $m_{\chi}\in[1,\,10^{6}]$
GeV considered here, a mini-BH of mass

$M_{\rm BH}=m_{\chi}N_{\rm{self}}\simeq 4.8\times
10^{46}\left(\frac{\rm{GeV}}{m_{\chi}}\right)^{3/2}\left(\frac{\rm{T}}{10^{5}\
\rm{K}}\right)^{3/2}{\rm GeV}$ (32)

forms at time

$\begin{split}t_{\rm BH}&=\frac{N_{\rm{self}}}{C_{\rm geom}}{\rm
Max}\left[\frac{\sigma_{\rm th}}{\sigma},1\right]\\\ &\approx
10^{13}~{}\left(\frac{\rho_{0}}{\rho_{\chi}}\right)\left(\frac{\rm
GeV}{m_{\chi}}\right)^{3/2}{\rm Max}\left[\frac{\sigma_{\rm
th}}{\sigma},1\right]~{}{\rm yr}\end{split}$ (33)

when $N_{\chi}(t_{\rm BH})=N_{\rm self}$. For $\rho_{\chi}\simeq\rho_{0}$,
collapse occurs at $t<10^{10}$ yrs for $m_{\chi}\gtrsim 100$ GeV.

###  VI.3 Bose-Einstein Condensation

In the preceding discussion, we neglected the formation of a Bose-Einstein
Condensate (BEC) of DM. Earlier work has shown that a BEC can accelerate BH
formation for bosonic DM with mass $m_{\chi}<m_{\rm BEC}$, where $m_{\rm
BEC}\simeq 2\times 10^{4}$ GeV for $T=10^{5}$ K and $m_{\rm BEC}\simeq 5\times
10^{3}$ GeV for $T=10^{6}$ K [120, 115]. We briefly summarize their findings
below.

When the critical temperature for BEC of DM

$T_{C}=\frac{2\pi}{m_{\chi}}\left(\frac{3N_{\chi}(t)}{4\pi r_{\rm
th}^{3}~{}\zeta(3/2)}\right)^{2/3}$ (34)

exceeds the NS temperature, a condensate forms at the core within a radius

$r_{\rm BEC}=1.5\times 10^{-4}{\rm cm}\left(\frac{1{\rm
GeV}}{m_{\chi}}\right)^{1/2}$ (35)

which is much smaller than $r_{\rm th}$. This occurs when the total number of
accumulated DM,

$N_{\chi}(t)=N_{\rm C}=10^{36}~{}\left(\frac{T}{10^{5}~{}{\rm
K}}\right)^{3}\,,$ (36)

and the number of DM particles in the condensate,

$N_{\rm
BEC}(t)=N_{\chi}(t)-N_{\chi}(t)\left(\frac{T}{T_{C}}\right)^{3/2}=N_{\chi}(t)-N_{C}\,,$
(37)

increase at the accretion rate on timescales that are large compared to the
thermalization time. The BEC becomes self-gravitating when $N_{\rm BEC}>4\pi
r_{\rm BEC}^{3}\rho_{B}/m_{\chi}\approx 10^{28}({\rm GeV}/m_{\chi})$. Since
$N_{\rm BEC}\ll N^{\rm bosonic}_{\rm Chandra}$ for the mass range of interest,
the self-gravitating BEC is stable, and the time required for it to collapse
is obtained as a solution to $N_{\chi}(t^{\rm BEC}_{\rm BH})=N^{\rm
bosonic}_{\rm Chandra}+N_{C}$. We find that a BH of mass

$M^{\rm BEC}_{\rm BH}=m_{\chi}N^{\rm bosonic}_{\rm Chandra}\simeq 1.5\times
10^{38}\left(\frac{\rm GeV}{m_{\chi}}\right)\rm GeV$ (38)

forms at a time

$\begin{split}t^{\rm BEC}_{\rm BH}&=\frac{N^{\rm bosonic}_{\rm
Chandra}+N_{C}}{C_{\rm geom}}{\rm Max}\left[\frac{\sigma_{\rm
th}}{\sigma},1\right]\\\ &\approx 5.2\times
10^{4}\left(\frac{\rho_{0}}{\rho_{\chi}}\right)\left(\frac{\rm
GeV}{m_{\chi}}\right){\rm Max}\left[\frac{\sigma_{\rm th}}{\sigma},1\right]\\\
&\times\left[1+\frac{2}{3}\left(\frac{m_{\chi}}{10~{}\rm
GeV}\right)^{2}\left(\frac{T}{10^{5}~{}{\rm K}}\right)^{3}\right]{\rm
yr}\,.\end{split}$ (39)

At $T=10^{5}$ K, this equation implies $t^{\rm BEC}_{\rm BH}<t_{\rm BH}$ for
$m_{\chi}<2\times 10^{4}$ GeV, and BEC greatly reduces the time for BH
formation in this mass range. For example, when $m_{\chi}=1$ GeV, the BEC
reduces the collapse time by $\approx 10^{8}$ yrs. Consequently, mini-BH
formation is possible within the Universe’s lifetime for the range of DM
considered in this study ($m_{\chi}=1-10^{6}$ GeV).

###  VI.4 Growth of the black hole due to accretion of baryons

The mini-BH continues to grow if the Bondi-Hoyle accretion rate,
$\dot{M}_{\rm{acc}}=4\pi\lambda_{s}(GM_{\rm{BH}}/v_{\rm{s}}^{2})^{2}\rho_{\rm{C}}v_{\rm{s}}$,
exceeds the Hawking evaporation rate, $\dot{M}_{\rm{evap}}=(15360\pi
G^{2}M^{2}_{\rm{BH}})^{-1}$, where $v_{s}$ is the sound speed, $\rho_{\rm C}$
is the central baryon density of the NS surrounding the BH, $M_{\rm BH}$ is
the mass of the BH at the star’s core, and $\lambda_{s}$ is a dimensionless
constant that depends on the EOS of matter in the NS core. The growth is
possible only if the BH reaches a critical mass,

$M_{\rm{BH}}^{\rm{crit}}\simeq 1.3\times
10^{37}\left(\frac{v_{s}}{0.3~{}c}\right)^{3/4}\left(\frac{10^{15}{\rm
g/cm}^{3}}{\lambda_{s}\rho_{B}}\right)^{1/4}{\rm GeV}\,.$ (40)

The analysis of BH growth rates presented in Refs. [115, 120] use
$\lambda_{s}=1/4$, corresponding to a polytropic index $\Gamma=5/3$. Since the
EOS of the NS core is expected to be stiff, with $\Gamma\simeq 2$ and
$v_{s}\simeq c/\sqrt{3}$, there has been much recent work on understanding
Bondi-Hoyle accretion under these conditions [121, 122, 123, 124, 125, 126].
These studies suggest $\lambda_{s}\simeq 1.3$ for realistic NS EOSs.

Figure 8: Distribution for the mass $m_{\chi}\in[10^{-1}$, $2.5\times 10^{5}]\
{\rm GeV}$ and scattering cross-section $\sigma_{\chi}$ inferred from the
collapse time derived from the observed number of BBHs by the detector
networks considered in Fig. 7. Here we consider the scenario where the NS
implodes without the formation of a BEC state for $m_{\chi}\geq 2\times
10^{4}\ {\rm GeV}$, shown in grey. For $m_{\chi}<2\times 10^{4}\ {\rm GeV}$,
NS implosion occurs at shorter timescales through the formation of a BEC. We
use $\rho_{\chi}=[1,10,100]\ {\rm GeV/cm^{3}}$.

For the scenario in which a BH forms without a BEC intermediate state, i.e.,
for $m_{\chi}\gtrsim m_{\rm BEC}$, from Eq. (32) and Eq. (40) we find that
$M_{\rm BH}>M_{\rm{BH}}^{\rm{crit}}$ for the mass range of interest. In the
scenario that involves a BEC intermediate state, i.e., for $m_{\chi}<m_{\rm
BEC}$, comparing Eq. (38) and Eq. (40) we find that $M_{\rm
BH}>M_{\rm{BH}}^{\rm{crit}}$ is only satisfied for $m_{\chi}\lesssim 10$ GeV.
For larger $m_{\chi}$ the BH evaporates, because the accretion of baryons
cannot keep pace with mass loss from Hawking radiation.

Figure 9: Comparison between the constraints obtained for
$m_{\chi}\in[1,10^{4}]$ GeV with the lowest ambient DM density considered in
this work ($\rho_{\chi}=1\ {\rm Gev/cm^{3}}$) with the latest constraints from
the direct detection experiment, LZ [127].

###  VI.5 Growth due to dark matter accretion

If the thermalization time is short compared to the BH Hawking evaporation
time $t_{\rm Haw}=15360\pi G^{2}M^{3}_{BH}/3\simeq 5\times 10^{4}(100~{}{\rm
GeV}/m_{\chi})^{3}$, DM accreting onto the NS reaches the BEC efficiently. The
maximal impact parameter for DM capture is larger than the radius of the BEC,
and DM particles reaching the BEC can feed the growth of the BH at a rate
$\dot{M}_{\rm DM}=C_{\rm NS}t$ if Hawking radiation does not disrupt the BEC
[115]. Direct heating of DM will be absent if the dark sector does not contain
mediator particles with mass less than the Hawking temperature $T_{\rm
haw}=1/(8\pi GM_{\rm BH})$. Further, even when DM is directly heated by
Hawking radiation, the change in its temperature would be negligible if it can
thermalize with baryons quickly [115]. In this case, $\dot{M}_{\rm DM}=C_{\rm
NS}t>\dot{M}_{\rm{evap}}$, and the accretion of DM prevents BH evaporation for
heavier masses [115]. Comparing the timescales for thermalization and
evaporation, we find that DM accretion can prevent BH evaporation for
$m_{\chi}\lesssim 4\times 10^{3}~{}(\sigma_{\chi}/\sigma_{\rm
th})(T/10^{5}~{}{\rm K})$ GeV.

###  VI.6 Neutron star implosion timescale

The timescale for a mini-BH to devour the entire NS by Bondi-Hoyle accretion
has been computed recently for the stiff EOS expected in NS cores [122]. For a
realistic NS model, they estimate this timescale to be given by

$t_{\rm dev}\simeq 0.3\left(\frac{10^{46}~{}{\rm GeV}}{M_{\rm BH}}\right){\rm
yr}\,,$ (41)

where $M_{\rm BH}$ is the mass of the initial mini-BH. For large DM mass,
i.e., $m_{\chi}>m_{\rm BEC}$ where $m_{\rm BEC}\simeq 2\times 10^{4}$ GeV at
$T=10^{5}$ K, the NS is destroyed on a timescale $t\approx 6\times
10^{-2}(m_{\chi}/{\rm GeV})^{3/2}(10^{5}~{}{\rm K}/T)^{3/2}~{}{\rm yr}$. For
$m_{\chi}<m_{\rm BEC}$ (which involves a BEC intermediate state), the NS is
destroyed on a timescale $t\approx 3.2\times 10^{7}(m_{\chi}/{\rm GeV})~{}{\rm
yr}$.

##  VII Constraining dark matter properties from the implosion timescale

As discussed before, GW observations can determine the relative abundance of
BBHs and BNSs, which allows us to infer the implosion timescale $t_{c}.$ The
posterior distribution of $t_{c}$, or the lower bound on $t_{c}$ if GW
observations can not conclusively measure a nonzero BBH rate, can be used to
constrain a region in the space of DM mass, interaction cross section and
dispersion velocity in the vicinity of merger. To illustrate the sensitivity
of this method, we assume that the time it takes for DM particle capture to
form a BH of critical mass, $M^{\rm{crit}}_{\rm{BH}}$, is approximately equal
to the total collapse time $t_{c}$ of the NS to form a BH through this
channel. Furthermore, for illustration, we assume that $t_{c}=10\,\rm Gyr.$

Figure 8 shows the constraint on the DM particle mass, $m_{\chi}$, and the DM-
baryon interaction cross-section, $\sigma_{\chi}$, for the inferred collapse
times shown in Fig. 7 and three values of the DM density: $\rho_{\chi}=1\,{\rm
GeV/m^{3}},$ ${\rm 10\,GeV/m^{3}}$, and $100\,{\rm GeV/m^{3}}.$ In all cases
the dispersion velocity is assumed to be 200 km/s. The precision of the
inference of DM parameters is directly related to how precisely we can deduce
the collapse time from the observed number of BBH mergers and the measured
$\tilde{\Lambda}$ from the GW signal. Compared to the A+ network, the XG
network of Cosmic Explorer and Einstein Telescope can improve the width of the
distribution by almost a factor of 10. From Fig. 9, these constraints are
competitive with those of direct detection DM experiments, such as the LUX-
ZEPLIN (LZ) experiment [127], over the entire range of $m_{\chi}$.

##  VIII Conclusions and Future prospects

In the last decade, direct detection experiments for DM have proven to be
successful in constraining the parameter space for WIMPs as well as axionic DM
[127, 128, 129, 130, 131]. Microlensing surveys [132, 133], calculations from
dwarf-galaxy dynamics [134, 135], as well as searches for sub-solar mass
compact binary mergers in GW data [136, 10, 11, 112] have probed the DM
compact object parameter space, deriving limits on the abundance of DM in
these objects and on their mass spectrum.

In this work, we present a method to combine measurements from GW observations
with the particle properties of DM, especially in the WIMP mass range, through
the observation (or lack thereof) of a novel population of BBHs in the mass
range $1M_{\odot}\leq m\leq 2M_{\odot}$. With XG ground-based GW detectors,
our ability to measure the effective tidal deformability will improve
tremendously, as shown in Fig. 4. With this refinement, we expect to identify
a significant number of BBH mergers in this mass range, if this population
exists in the Universe. The observed number of mergers can potentially
constrain the implosion time of NSs if these BHs form from the collapse of old
NSs due to the presence of a mini-BH at their cores. If the mini-BH forms
through the accumulation of DM in the cores of NSs over their lifetime, the
collapse time can then inform us on DM particle properties.

We illustrate how this mechanism works. We report limits on the interaction
cross-section of DM particles with baryons and particle masses for a
simplistic scenario where asymmetric DM interacts weakly with hadrons and gets
captured through scattering within NS cores to form a self-gravitating mass
that forms a BH. As expected, the constraints are more precise in the case of
XG GW detectors, which is a direct consequence of how precisely the binary’s
effective tidal deformability, and hence the collapse time, is inferred from
the observed number of BBH mergers.

One potential issue is the degeneracy with primordial BHs, which could also
populate this spectrum of masses. One possible way to remove the degeneracy
involves accurate spin measurements to distinguish between the two formation
channels. Furthermore, if primordial BHs contribute to the super-solar mass
range, we would also expect to observe some sub-solar mass primordial BHs. The
absence of observed old NSs also complements the observation of BHs formed
from imploding NSs. In this work, we show how effectively GW detections can
not only constrain DM interactions, but also rule out models that allow for
NSs to implode through DM accumulation if no such population is observed in
the future.

## Acknowledgments

We thank K. Belczynski for useful discussions on BNS delay times, as well as
M. Baryakhtar and T. Slatyer for discussions on the DM scenarios considered in
this paper. D.S. and B.S.S. were supported in part by NSF grant No.
PHY-1836779, PHY-2012083, AST-2006384 and PHY-2207638. A.G. is supported by
NSF grant No. AST-2205920. E.B. is supported by NSF Grants No. AST-2006538,
PHY-2207502, PHY-090003 and PHY20043, and NASA Grants No. 19-ATP19-0051,
20-LPS20-0011 and 21-ATP21-0010. This research project was conducted using
computational resources at the Maryland Advanced Research Computing Center
(MARCC). S.R. is supported by the U.S. Department of Energy under Grant No.
DE-FG02- 00ER41132 and National Science Foundation’s Physics Frontier Center:
The Network for Neutrinos, Nuclear Astrophysics, and Symmetries. This paper
has the INT preprint number INT-PUB-22-028. Part of E.B.’s and B.S.S.’s work
was performed at the Aspen Center for Physics, which is supported by National
Science Foundation grant PHY-1607611. This research was also supported in part
by the National Science Foundation under Grant No. NSF PHY-1748958.

## References

  * Marsh and Hoof [2021] D. J. E. Marsh and S. Hoof, (2021), arXiv:2106.08797 [hep-ph] .
  * Bertone  _et al._ [2005] G. Bertone, D. Hooper, and J. Silk, Phys. Rept. 405, 279 (2005), arXiv:hep-ph/0404175 .
  * Carr and Kuhnel [2020] B. Carr and F. Kuhnel, Ann. Rev. Nucl. Part. Sci. 70, 355 (2020), arXiv:2006.02838 [astro-ph.CO] .
  * Abbott  _et al._ [2019a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 031040 (2019a), arXiv:1811.12907 [astro-ph.HE] .
  * Abbott  _et al._ [2021a] R. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 11, 021053 (2021a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2021b] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021b), arXiv:2111.03606 [gr-qc] .
  * 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] .
  * Wong  _et al._ [2021] K. W. K. Wong, G. Franciolini, V. De Luca, V. Baibhav, E. Berti, P. Pani, and A. Riotto, Phys. Rev. D 103, 023026 (2021), arXiv:2011.01865 [gr-qc] .
  * Franciolini  _et al._ [2021] G. Franciolini, V. Baibhav, V. De Luca, K. K. Y. Ng, K. W. K. Wong, E. Berti, P. Pani, A. Riotto, and S. Vitale, (2021), arXiv:2105.03349 [gr-qc] .
  * Abbott  _et al._ [2018a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 231103 (2018a), arXiv:1808.04771 [astro-ph.CO] .
  * Abbott  _et al._ [2019b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 123, 161102 (2019b), arXiv:1904.08976 [astro-ph.CO] .
  * Ng  _et al._ [2022] K. K. Y. Ng, G. Franciolini, E. Berti, P. Pani, A. Riotto, and S. Vitale, Astrophys. J. Lett. 933, L41 (2022), arXiv:2204.11864 [astro-ph.CO] .
  * De Luca  _et al._ [2021] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, JCAP 11, 039 (2021), arXiv:2106.13769 [astro-ph.CO] .
  * Eda  _et al._ [2013] K. Eda, Y. Itoh, S. Kuroyanagi, and J. Silk, Phys. Rev. Lett. 110, 221101 (2013), arXiv:1301.5971 [gr-qc] .
  * Eda  _et al._ [2015] K. Eda, Y. Itoh, S. Kuroyanagi, and J. Silk, Phys. Rev. D 91, 044045 (2015), arXiv:1408.3534 [gr-qc] .
  * Hannuksela  _et al._ [2019] O. A. Hannuksela, K. W. K. Wong, R. Brito, E. Berti, and T. G. F. Li, Nature Astron. 3, 447 (2019), arXiv:1804.09659 [astro-ph.HE] .
  * Hannuksela  _et al._ [2020] O. A. Hannuksela, K. C. Y. Ng, and T. G. F. Li, Phys. Rev. D 102, 103022 (2020), arXiv:1906.11845 [astro-ph.CO] .
  * Kavanagh  _et al._ [2020] B. J. Kavanagh, D. A. Nichols, G. Bertone, and D. Gaggero, Phys. Rev. D 102, 083006 (2020), arXiv:2002.12811 [gr-qc] .
  * Annulli  _et al._ [2020] L. Annulli, V. Cardoso, and R. Vicente, Phys. Rev. D 102, 063022 (2020), arXiv:2009.00012 [gr-qc] .
  * Traykova  _et al._ [2021] D. Traykova, K. Clough, T. Helfer, E. Berti, P. G. Ferreira, and L. Hui, Phys. Rev. D 104, 103014 (2021), arXiv:2106.08280 [gr-qc] .
  * Coogan  _et al._ [2022] A. Coogan, G. Bertone, D. Gaggero, B. J. Kavanagh, and D. A. Nichols, Phys. Rev. D 105, 043009 (2022), arXiv:2108.04154 [gr-qc] .
  * Vicente and Cardoso [2022] R. Vicente and V. Cardoso, Phys. Rev. D 105, 083008 (2022), arXiv:2201.08854 [gr-qc] .
  * Speeney  _et al._ [2022] N. Speeney, A. Antonelli, V. Baibhav, and E. Berti, Phys. Rev. D 106, 044027 (2022), arXiv:2204.12508 [gr-qc] .
  * Cole  _et al._ [2022] P. S. Cole, A. Coogan, B. J. Kavanagh, and G. Bertone, (2022), arXiv:2207.07576 [astro-ph.CO] .
  * Arvanitaki and Dubovsky [2011] A. Arvanitaki and S. Dubovsky, Phys. Rev. D 83, 044026 (2011), arXiv:1004.3558 [hep-th] .
  * Brito  _et al._ [2015] R. Brito, V. Cardoso, and P. Pani, Lect. Notes Phys. 906, pp.1 (2015), arXiv:1501.06570 [gr-qc] .
  * Arvanitaki  _et al._ [2015] A. Arvanitaki, M. Baryakhtar, and X. Huang, Phys. Rev. D 91, 084011 (2015), arXiv:1411.2263 [hep-ph] .
  * Brito  _et al._ [2017a] R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. Lett. 119, 131101 (2017a), arXiv:1706.05097 [gr-qc] .
  * Brito  _et al._ [2017b] R. Brito, S. Ghosh, E. Barausse, E. Berti, V. Cardoso, I. Dvorkin, A. Klein, and P. Pani, Phys. Rev. D 96, 064050 (2017b), arXiv:1706.06311 [gr-qc] .
  * Ng  _et al._ [2021] K. K. Y. Ng, S. Vitale, O. A. Hannuksela, and T. G. F. Li, Phys. Rev. Lett. 126, 151102 (2021), arXiv:2011.06010 [gr-qc] .
  * Ghosh  _et al._ [2019] S. Ghosh, E. Berti, R. Brito, and M. Richartz, Phys. Rev. D 99, 104030 (2019), arXiv:1812.01620 [gr-qc] .
  * Yuan  _et al._ [2021] C. Yuan, R. Brito, and V. Cardoso, (2021), arXiv:2106.00021 [gr-qc] .
  * Bramante and Linden [2014] J. Bramante and T. Linden, Phys. Rev. Lett. 113, 191301 (2014), arXiv:1405.1031 [astro-ph.HE] .
  * Dominik  _et al._ [2012] M. Dominik, K. Belczynski, C. Fryer, D. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Astrophys. J. 759, 52 (2012), arXiv:1202.4901 [astro-ph.HE] .
  * Safarzadeh and Berger [2019] M. Safarzadeh and E. Berger, Astrophys. J. Lett. 878, L12 (2019), arXiv:1904.08436 [astro-ph.GA] .
  * Safarzadeh  _et al._ [2019] M. Safarzadeh, E. Berger, K. K. Y. Ng, H.-Y. Chen, S. Vitale, C. Whittle, and E. Scannapieco, Astrophys. J. Lett. 878, L13 (2019), arXiv:1904.10976 [astro-ph.HE] .
  * McCarthy  _et al._ [2020] K. S. McCarthy, Z. Zheng, and E. Ramirez-Ruiz, Mon. Not. Roy. Astron. Soc. 499, 5220 (2020), arXiv:2007.15024 [astro-ph.GA] .
  * Greggio  _et al._ [2020] L. Greggio, P. Simonetti, and F. Matteucci, Mon. Not. Roy. Astron. Soc. 500, 1755 (2020), arXiv:2009.13138 [astro-ph.HE] .
  * Takhistov  _et al._ [2021] V. Takhistov, G. M. Fuller, and A. Kusenko, Phys. Rev. Lett. 126, 071101 (2021), arXiv:2008.12780 [astro-ph.HE] .
  * Dasgupta  _et al._ [2021] B. Dasgupta, R. Laha, and A. Ray, Phys. Rev. Lett. 126, 141105 (2021), arXiv:2009.01825 [astro-ph.HE] .
  * Bailyn  _et al._ [1998] C. D. Bailyn, R. K. Jain, P. Coppi, and J. A. Orosz, Astrophys. J. 499, 367 (1998), arXiv:astro-ph/9708032 .
  * Özel  _et al._ [2010] F. Özel, D. Psaltis, R. Narayan, and J. E. McClintock, ApJ 725, 1918 (2010), arXiv:1006.2834 [astro-ph.GA] .
  * Farr  _et al._ [2011] W. M. Farr, N. Sravan, A. Cantrell, L. Kreidberg, C. D. Bailyn, I. Mandel, and V. Kalogera, ApJ 741, 103 (2011), arXiv:1011.1459 [astro-ph.GA] .
  * Belczynski  _et al._ [2012] K. Belczynski, G. Wiktorowicz, C. L. Fryer, D. E. Holz, and V. Kalogera, ApJ 757, 91 (2012), arXiv:1110.1635 [astro-ph.GA] .
  * Fryer  _et al._ [2022] C. L. Fryer, A. Olejak, and K. Belczynski, Astrophys. J. 931, 94 (2022), arXiv:2204.13025 [astro-ph.HE] .
  * De Luca  _et al._ [2020] V. De Luca, G. Franciolini, P. Pani, and A. Riotto, Phys. Rev. D 102, 043505 (2020), arXiv:2003.12589 [astro-ph.CO] .
  * Gerosa  _et al._ [2018] D. Gerosa, E. Berti, R. O’Shaughnessy, K. Belczynski, M. Kesden, D. Wysocki, and W. Gladysz, Phys. Rev. D 98, 084036 (2018), arXiv:1808.02491 [astro-ph.HE] .
  * Belczynski  _et al._ [2020] K. Belczynski _et al._ , Astron. Astrophys. 636, A104 (2020), arXiv:1706.07053 [astro-ph.HE] .
  * Franciolini  _et al._ [2022] G. Franciolini, R. Cotesta, N. Loutrel, E. Berti, P. Pani, and A. Riotto, Phys. Rev. D 105, 063510 (2022), arXiv:2112.10660 [astro-ph.CO] .
  * Mukherjee and Silk [2021] S. Mukherjee and J. Silk, Mon. Not. Roy. Astron. Soc. 506, 3977 (2021), arXiv:2105.11139 [gr-qc] .
  * Mukherjee  _et al._ [2022] S. Mukherjee, M. S. P. Meinema, and J. Silk, Mon. Not. Roy. Astron. Soc. 510, 6218 (2022), arXiv:2107.02181 [astro-ph.CO] .
  * Hansen and Kawaler [1994] C. J. Hansen and S. D. Kawaler,  _Stellar Interiors: Physical Principles, Structure, and Evolution_ (Birkhäuser, 1994).
  * Chabrier [2005] G. Chabrier, Astrophys. Space Sci. Libr. 327, 41 (2005), arXiv:astro-ph/0409465 .
  * Damour [1986] T. Damour, in  _300 Years of Gravity: A Conference to Mark the 300th Anniversary of the Publication of Newton’s Principia_ (1986).
  * Mora and Will [2004] T. Mora and C. M. Will, Phys. Rev. D 69, 104021 (2004), [Erratum: Phys.Rev.D 71, 129901 (2005)], arXiv:gr-qc/0312082 .
  * Damour and Nagar [2009] T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009), arXiv:0906.0096 [gr-qc] .
  * Landry and Poisson [2015] P. Landry and E. Poisson, Phys. Rev. D 91, 104018 (2015), arXiv:1503.07366 [gr-qc] .
  * Chia [2021] H. S. Chia, Phys. Rev. D 104, 024013 (2021), arXiv:2010.07300 [gr-qc] .
  * Poisson [2021] E. Poisson, Phys. Rev. D 104, 104062 (2021), arXiv:2108.07328 [gr-qc] .
  * Binnington and Poisson [2009] T. Binnington and E. Poisson, Phys. Rev. D 80, 084018 (2009), arXiv:0906.1366 [gr-qc] .
  * Chatziioannou [2020] K. Chatziioannou, Gen. Rel. Grav. 52, 109 (2020), arXiv:2006.03168 [gr-qc] .
  * Hinderer  _et al._ [2010] T. Hinderer, B. D. Lackey, R. N. Lang, and J. S. Read, Phys. Rev. D 81, 123016 (2010), arXiv:0911.3535 [astro-ph.HE] .
  * Abbott  _et al._ [2020a] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Class. Quant. Grav. 37, 045006 (2020a), arXiv:1908.01012 [gr-qc] .
  * Yang  _et al._ [2018] H. Yang, W. E. East, and L. Lehner, Astrophys. J. 856, 110 (2018), [Erratum: Astrophys.J. 870, 139 (2019)], arXiv:1710.05891 [gr-qc] .
  * Chen and Chatziioannou [2020] H.-Y. Chen and K. Chatziioannou, Astrophys. J. Lett. 893, L41 (2020), arXiv:1903.11197 [astro-ph.HE] .
  * Fasano  _et al._ [2020] M. Fasano, K. W. K. Wong, A. Maselli, E. Berti, V. Ferrari, and B. S. Sathyaprakash, Phys. Rev. D 102, 023025 (2020), arXiv:2005.01726 [astro-ph.HE] .
  * Manchester  _et al._ [2005] R. N. Manchester, G. B. Hobbs, A. Teoh, and M. Hobbs, Astron. J. 129, 1993 (2005), arXiv:astro-ph/0412641 .
  * Breton [2008] R. P. Breton,  _Radio Pulsars in Binary Systems_ , Ph.D. thesis, McGill U. (2008), arXiv:0907.2623 [astro-ph.SR] .
  * Peters and Mathews [1963] P. C. Peters and J. Mathews, Phys. Rev. 131, 435 (1963).
  * Sathyaprakash and Dhurandhar [1991] B. S. Sathyaprakash and S. V. Dhurandhar, Phys. Rev. D 44, 3819 (1991).
  * Sathyaprakash and Schutz [2009] B. S. Sathyaprakash and B. F. Schutz, Living Rev. Rel. 12, 2 (2009), arXiv:0903.0338 [gr-qc] .
  * Flanagan and Hinderer [2008] E. E. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008), arXiv:0709.1915 [astro-ph] .
  * Favata [2014] M. Favata, Phys. Rev. Lett. 112, 101101 (2014), arXiv:1310.8288 [gr-qc] .
  * De  _et al._ [2018] S. De, D. Finstad, J. M. Lattimer, D. A. Brown, E. Berger, and C. M. Biwer, Phys. Rev. Lett. 121, 091102 (2018), [Erratum: Phys.Rev.Lett. 121, 259902 (2018)], arXiv:1804.08583 [astro-ph.HE] .
  * Abbott  _et al._ [2018b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018b), arXiv:1805.11581 [gr-qc] .
  * Johnson-Mcdaniel  _et al._ [2020] N. K. Johnson-Mcdaniel, A. Mukherjee, R. Kashyap, P. Ajith, W. Del Pozzo, and S. Vitale, Phys. Rev. D 102, 123010 (2020), arXiv:1804.08026 [gr-qc] .
  * Chen  _et al._ [2020] A. Chen, N. K. Johnson-McDaniel, T. Dietrich, and R. Dudi, Phys. Rev. D 101, 103008 (2020), arXiv:2001.11470 [astro-ph.HE] .
  * Finn [1992] L. S. Finn, Phys. Rev. D 46, 5236 (1992), arXiv:gr-qc/9209010 .
  * Vallisneri [2008] M. Vallisneri, Phys. Rev. D 77, 042001 (2008), arXiv:gr-qc/0703086 .
  * Borhanian [2020] S. Borhanian, (2020), arXiv:2010.15202 [gr-qc] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Abbott  _et al._ [2019c] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. X 9, 011001 (2019c), arXiv:1805.11579 [gr-qc] .
  * Radice  _et al._ [2018] D. Radice, A. Perego, F. Zappa, and S. Bernuzzi, Astrophys. J. Lett. 852, L29 (2018), arXiv:1711.03647 [astro-ph.HE] .
  * Abbott  _et al._ [2020b] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Astrophys. J. Lett. 892, L3 (2020b), arXiv:2001.01761 [astro-ph.HE] .
  * Reitze  _et al._ [2019] D. Reitze _et al._ , Bull. Am. Astron. Soc. 51, 035 (2019), arXiv:1907.04833 [astro-ph.IM] .
  * Punturo  _et al._ [2010] M. Punturo _et al._ , Class. Quant. Grav. 27, 084007 (2010).
  * Barsotti  _et al._ [2018] L. Barsotti, L. McCuller, M.Evans, and P. Fritschel, “The A+ Design Curve,”  (2018).
  * Adhikari  _et al._ [2019] R. X. Adhikari _et al._ , Class. Quant. Grav. 36, 245010 (2019), arXiv:1905.02842 [astro-ph.HE] .
  * Dietrich  _et al._ [2017] T. Dietrich, S. Bernuzzi, and W. Tichy, Phys. Rev. D 96, 121501 (2017), arXiv:1706.02969 [gr-qc] .
  * Dietrich  _et al._ [2019a] T. Dietrich _et al._ , Phys. Rev. D 99, 024029 (2019a), arXiv:1804.02235 [gr-qc] .
  * Dietrich  _et al._ [2019b] T. Dietrich, A. Samajdar, S. Khan, N. K. Johnson-McDaniel, R. Dudi, and W. Tichy, Phys. Rev. D 100, 044003 (2019b), arXiv:1905.06011 [gr-qc] .
  * Kokkotas and Schaefer [1995] K. D. Kokkotas and G. Schaefer, Mon. Not. Roy. Astron. Soc. 275, 301 (1995), arXiv:gr-qc/9502034 .
  * Gualtieri  _et al._ [2001] L. Gualtieri, E. Berti, J. A. Pons, G. Miniutti, and V. Ferrari, Phys. Rev. D 64, 104007 (2001), arXiv:gr-qc/0107046 .
  * Pons  _et al._ [2002] J. A. Pons, E. Berti, L. Gualtieri, G. Miniutti, and V. Ferrari, Phys. Rev. D 65, 104021 (2002), arXiv:gr-qc/0111104 .
  * Berti  _et al._ [2002] E. Berti, J. A. Pons, G. Miniutti, L. Gualtieri, and V. Ferrari, Phys. Rev. D 66, 064013 (2002), arXiv:gr-qc/0208011 .
  * Steinhoff  _et al._ [2016] J. Steinhoff, T. Hinderer, A. Buonanno, and A. Taracchini, Phys. Rev. D 94, 104028 (2016), arXiv:1608.01907 [gr-qc] .
  * Andersson and Pnigouras [2021] N. Andersson and P. Pnigouras, Mon. Not. Roy. Astron. Soc. 503, 533 (2021), arXiv:1905.00012 [gr-qc] .
  * Schmidt and Hinderer [2019] P. Schmidt and T. Hinderer, Phys. Rev. D 100, 021501 (2019), arXiv:1905.00818 [gr-qc] .
  * Pratten  _et al._ [2022] G. Pratten, P. Schmidt, and N. Williams, Phys. Rev. Lett. 129, 081102 (2022), arXiv:2109.07566 [astro-ph.HE] .
  * Gamba and Bernuzzi [2022] R. Gamba and S. Bernuzzi, (2022), arXiv:2207.13106 [gr-qc] .
  * Williams  _et al._ [2022] N. Williams, G. Pratten, and P. Schmidt, Phys. Rev. D 105, 123032 (2022), arXiv:2203.00623 [astro-ph.HE] .
  * Kiziltan  _et al._ [2013] B. Kiziltan, A. Kottas, M. De Yoreo, and S. E. Thorsett, Astrophys. J. 778, 66 (2013), arXiv:1309.6635 [astro-ph.SR] .
  * Farrow  _et al._ [2019] N. Farrow, X.-J. Zhu, and E. Thrane, Astrophys. J. 876, 18 (2019), arXiv:1902.03300 [astro-ph.HE] .
  * Landry and Read [2021] P. Landry and J. S. Read, Astrophys. J. Lett. 921, L25 (2021), arXiv:2107.04559 [astro-ph.HE] .
  * Abbott  _et al._ [2021c] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), (2021c), arXiv:2111.03634 [astro-ph.HE] .
  * Dominik  _et al._ [2013] M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Astrophys. J. 779, 72 (2013), arXiv:1308.1546 [astro-ph.HE] .
  * Dominik  _et al._ [2015] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik, and F. Pannarale, Astrophys. J. 806, 263 (2015), arXiv:1405.7016 [astro-ph.HE] .
  * Santoliquido  _et al._ [2021] F. Santoliquido, M. Mapelli, N. Giacobbo, Y. Bouffanais, and M. C. Artale, Mon. Not. Roy. Astron. Soc. 502, 4877 (2021), arXiv:2009.03911 [astro-ph.HE] .
  * Madau and Dickinson [2014] P. Madau and M. Dickinson, Ann. Rev. Astron. Astrophys. 52, 415 (2014), arXiv:1403.0007 [astro-ph.CO] .
  * Sahni and Starobinsky [2000] V. Sahni and A. A. Starobinsky, Int. J. Mod. Phys. D 9, 373 (2000), arXiv:astro-ph/9904398 .
  * Byrnes  _et al._ [2018] C. T. Byrnes, M. Hindmarsh, S. Young, and M. R. S. Hawkins, JCAP 08, 041 (2018), arXiv:1801.06138 [astro-ph.CO] .
  * Abbott  _et al._ [2022] R. Abbott _et al._ (LIGO Scientific, VIRGO, KAGRA), Phys. Rev. Lett. 129, 061104 (2022), arXiv:2109.12197 [astro-ph.CO] .
  * Press and Spergel [1985] W. H. Press and D. N. Spergel, Astrophys. J. 296, 679 (1985).
  * Gould [1987] A. Gould, Astrophys. J. 321, 571 (1987).
  * McDermott  _et al._ [2012] S. D. McDermott, H.-B. Yu, and K. M. Zurek, Phys. Rev. D 85, 023519 (2012), arXiv:1103.5472 [hep-ph] .
  * Garani  _et al._ [2019] R. Garani, Y. Genolini, and T. Hambye, JCAP 05, 035 (2019), arXiv:1812.08773 [hep-ph] .
  * Bell  _et al._ [2020] N. F. Bell, G. Busoni, S. Robles, and M. Virgato, JCAP 09, 028 (2020), arXiv:2004.14888 [hep-ph] .
  * Bertoni  _et al._ [2013] B. Bertoni, A. E. Nelson, and S. Reddy, Phys. Rev. D 88, 123505 (2013), arXiv:1309.1721 [hep-ph] .
  * Garani  _et al._ [2021] R. Garani, A. Gupta, and N. Raj, Phys. Rev. D 103, 043019 (2021), arXiv:2009.10728 [hep-ph] .
  * Kouvaris and Tinyakov [2011] C. Kouvaris and P. Tinyakov, Phys. Rev. Lett. 107, 091301 (2011), arXiv:1104.0382 [astro-ph.CO] .
  * East and Lehner [2019] W. E. East and L. Lehner, Phys. Rev. D 100, 124026 (2019), arXiv:1909.07968 [gr-qc] .
  * Baumgarte and Shapiro [2021] T. W. Baumgarte and S. L. Shapiro, Phys. Rev. D 103, L081303 (2021), arXiv:2101.12220 [astro-ph.HE] .
  * Richards  _et al._ [2021a] C. B. Richards, T. W. Baumgarte, and S. L. Shapiro, Phys. Rev. D 103, 104009 (2021a), arXiv:2102.09574 [astro-ph.HE] .
  * Richards  _et al._ [2021b] C. B. Richards, T. W. Baumgarte, and S. L. Shapiro, Mon. Not. Roy. Astron. Soc. 502, 3003 (2021b), arXiv:2101.08797 [astro-ph.HE] .
  * Aguayo-Ortiz  _et al._ [2021] A. Aguayo-Ortiz, E. Tejeda, O. Sarbach, and D. López-Cámara,  (2021), 10.1093/mnras/stab1127, arXiv:2102.12529 [astro-ph.HE] .
  * Giffin  _et al._ [2021] P. Giffin, J. Lloyd, S. D. McDermott, and S. Profumo, (2021), arXiv:2105.06504 [hep-ph] .
  * Aalbers  _et al._ [2022] J. Aalbers _et al._ (LZ), (2022), arXiv:2207.03764 [hep-ex] .
  * Agnes  _et al._ [2022] P. Agnes _et al._ (DarkSide-50), (2022), arXiv:2207.11966 [hep-ex] .
  * Aprile  _et al._ [2022] E. Aprile _et al._ ((XENON Collaboration)††, XENON), Phys. Rev. Lett. 129, 161805 (2022), arXiv:2207.11330 [hep-ex] .
  * Barak  _et al._ [2020] L. Barak _et al._ (SENSEI), Phys. Rev. Lett. 125, 171802 (2020), arXiv:2004.11378 [astro-ph.CO] .
  * Aprile  _et al._ [2018] E. Aprile _et al._ (XENON), Phys. Rev. Lett. 121, 111302 (2018), arXiv:1805.12562 [astro-ph.CO] .
  * Allsman  _et al._ [2001] R. A. Allsman _et al._ (Macho), Astrophys. J. Lett. 550, L169 (2001), arXiv:astro-ph/0011506 .
  * Tisserand  _et al._ [2007] P. Tisserand _et al._ (EROS-2), Astron. Astrophys. 469, 387 (2007), arXiv:astro-ph/0607207 .
  * Wyrzykowski  _et al._ [2011] L. Wyrzykowski, J. Skowron, S. Kozłowski, A. Udalski, M. K. Szymański, M. Kubiak, G. Pietrzyński, I. Soszyński, O. Szewczyk, K. Ulaczyk, R. Poleski, and P. Tisserand, Monthly Notices of the Royal Astronomical Society 416, 2949 (2011), arXiv:1106.2925 [astro-ph.GA] .
  * Koushiappas and Loeb [2017] S. M. Koushiappas and A. Loeb, Phys. Rev. Lett. 119, 041102 (2017), arXiv:1704.01668 [astro-ph.GA] .
  * Abbott  _et al._ [2005] B. Abbott _et al._ (LIGO Scientific), Phys. Rev. D 72, 082002 (2005), arXiv:gr-qc/0505042 .