# Exotic Compact Objects: The Dark White Dwarf

Michael Ryan<EMAIL_ADDRESS>Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA  Department of
Physics, The Pennsylvania State University, University Park, PA, 16802, USA
David Radice<EMAIL_ADDRESS>Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA  Department of
Physics, The Pennsylvania State University, University Park, PA, 16802, USA
Department of Astronomy & Astrophysics, The Pennsylvania State University,
University Park, PA, 16802, USA

###### Abstract

Several dark matter models allow for the intriguing possibility of exotic
compact object formation. These objects might have unique characteristics that
set them apart from their baryonic counterparts. Furthermore, gravitational
wave observations of their mergers may provide the only direct window on a
potentially entirely hidden sector. Here we discuss dark white dwarfs,
starting with an overview of the microphysical model and analytic scaling
relations of macroscopic properties derived from the non-relativistic limit.
We use the full relativistic formalism to confirm these scaling relations and
demonstrate that dark white dwarfs, if they exist, would have masses and tidal
deformabilities that are very different from those of baryonic compact
objects. Further, and most importantly, we demonstrate that dark white dwarf
mergers would be detectable by current or planned gravitational observatories
across several orders of magnitude in the particle-mass parameter space.
Lastly, we find universal relations analogous to the compactness-Love and
binary Love relations in neutron star literature. Using these results, we show
that gravitational wave observations would constrain the properties of the
dark matter particles constituting these objects.

cosmology: theory – dark matter – compact objects

##  I Introduction

Current dark matter search techniques focus on two primary channels: large-
scale structure constraints (e.g. [1, 2]) and direct and indirect detection
experiments (e.g. [3, 4, 5, 6, 7]). While these have constrained several of
the bulk properties of dark matter, i.e. that dark matter is cold,
particulate, and effectively collisionless on large scales, the current lack
of dark matter detection or production has not helped to narrow the space of
models. Likewise, the field of astro-particle indirect detection (e.g. [8, 9,
10]), while showing possible signals of interest [11], has also not yet
produced results. On the other hand, the advent of gravitational wave
observations has opened a new window on the universe that could illuminate the
dark sector in a completely novel manner.

Several promising alternative dark matter models have a “complex” (two or more
massive particles) particle zoo, and dissipative interactions create the
potential for gravitationally-bound macroscopic structures. Many of these
models even form exotic compact objects, with both “dark” black holes (a
normal black hole formed from dark matter instead of baryonic matter) [12, 13,
14, 15, 16, 17, 18, 19] and dark (neutron) stars [20, 21, 22, 23] having been
proposed. The merger of these exotic compact objects with each other, or with
astrophysical compact objects could be revealed by gravitational wave
detectors, such as LIGO. Alternatively, dark matter capture by ordinary
compact objects could result in the formation of compact, dark matter cores in
their interior[24, 25, 26, 27, 28, 29], creating hybrid dark/baryonic objects
ranging from planetary to stellar masses.

While dark black holes and neutron stars are the obvious counterparts to
ordinary black holes and neutron stars, little attention has been paid to the
dark equivalent of the third member of the ordinary compact object trio: the
(dark) white dwarf (DWD). In the simplest sense, a white dwarf is a compact
object predominantly composed of degenerate light electrons and massive
nuclei. In white dwarfs, a balance is struck between gravitational and
(electron-dominated) fermion degeneracy pressure forces. This balance sets
many of their macroscopic properties, like the radius and compactness. Above a
certain mass (the well-known Chandrasekhar limit), this balance is broken and
no static configurations exist. Likewise, above a certain density (the onset
of neutronization) these objects again become dynamically unstable to
collapse. Here, we consider the dark matter analogs to white dwarfs. These
would be compact objects, predominantly composed of two or more fermion
species of dark matter, where the pressure support is primarily provided by
fermion degeneracy pressure. Importantly, in hidden, multi-particle dark
matter models that lack “nuclear” reactions, DWDs may be the only sub-
Chandrasekhar-mass compact objects that can form. Also note that, with this
definition, in the limit where the particle species contribute equally to the
pressure and mass (the single particle limit), these objects may be
macroscopically indistinguishable from dark neutron stars in models with weak
or no nuclear interactions (e.g. [30, 22]).

The literature on objects fitting the definition above is sparse, with few
articles discussing objects that fit all three criteria. For example, Narain
et al.[30] describe a general, particle-mass-independent framework for dark,
fermionic compact objects, assuming single species composition, even studying
the effects of a generic interaction term. Gross et al.[31] studied DWD-like
objects as a possible final state of collapsed dark quark pockets (assuming
multiple species, but identical masses) in an analogue of primordial black
holes, computing their mass, radii, and long-term stability with that method.
Using a less general approach, we extend these analyses to a two-particle,
fermionic gas with potentially varying masses, and include a brief examination
of additional binary properties, like the tidal deformability and potential
universal relations. Hippert et al.[21] mentions the high plausibility of dark
white dwarf formation in the Twin Higgs mirror model, but focuses on neutron
stars instead. Brandt et al. [32] consider the stability of a multi-species
model for what they refer to as a pion star, focusing on using lattice QCD
methods to obtain a precise equation of state. Meanwhile, the discussion on
dark planets and similar low-mass, multi-component objects([24, 25, 27])
requires the mixing of ordinary matter with dark matter. Consequently, the DWD
space has remained unexplored.

We hasten to mention that, like most dark, non-primordial black hole and
neutron star models, these objects do not necessarily constitute the entirety
of dark matter. While the formation of these objects is outside the scope of
this work, we note that regardless of the pathway, the population must obey
the constraints imposed by microlensing and similar primordial black hole and
massive, compact halo object searches(e.g. [33, 34]). These searches impose
constraints on the fraction of dark matter contained in compact objects versus
all dark matter. For the masses considered here, this corresponds to less than
$\mathcal{O}(0.01)$ for subsolar mass DWDs, decreasing to
$\mathcal{O}(10^{-4})$ for objects above $10\text{\,}\mathrm{M_{\odot}}$.
Given that ordinary white dwarfs in the Milky Way correspond to a similar
$\mathcal{O}(0.01)$ fraction of the ordinary matter [35, 36], this constraint
seems reasonable for models that predict a dark astrophysical formation like
[21].

In the following sections we examine the properties of the most basic DWD
model following our definition above: two particle species forming a compact
object, with fermion degeneracy pressure providing the dominant support
against gravitational collapse. We start with a discussion of the basic
properties of DWD that can be inferred analytically in Section II, including
the calculation of an equation of state and scaling relations for the mass,
radius, and compactness in the non-relativistic limit. We then discuss the
results of the fully relativistic hydrostatic-equilibrium calculations across
the particle-mass parameter space, highlighting four example parameter cases
in Section III and examining several of the macroscopic attributes, potential
universal relations, and implications for DWD merger observations.
Importantly, we demonstrate that DWD mergers should be detectable by current
and planned gravitational wave observatories across much of the dark parameter
space and observations can be used to constrain the dark microphysics. Lastly,
we conclude in Section IV.

##  II Analytic Scaling Relations

###  II.1 Equation of State

We consider a simplified, cold compact object comprised of a cloud of
degenerate, fundamental fermionic particles, $L$, and $H$. The particles have
masses $m_{H}$ and $m_{L}$, defined such that $m_{H}\geq m_{L}$, analogous to
the Standard Model proton and electron. We will use the notation

$\displaystyle\mathcal{r}_{H}=\frac{m_{H}}{m_{p}},\;\;\text{ and
}\;\;\mathcal{r}_{L}=\frac{m_{L}}{m_{e}},$ (1)

throughout this and following sections, where $m_{p}$ and $m_{e}$ are the
proton and electron masses. We will assume approximately neutral “charge” in
bulk, i.e. an equal numbers of $L$ and $H$ particles.

The basic thermodynamic properties of such a cloud are well known ([37, 38]),
with the number density, pressure, and energy density of a single fermionic
particle $f$ given by:

$n_{f}=\frac{8\pi}{3h^{3}}p_{\rm
fermi}^{3}=\frac{x^{3}}{3\pi^{2}\lambda_{f}^{3}}$ (2) $\displaystyle P_{f}$
$\displaystyle=\frac{8\pi
m_{f}^{4}c^{5}}{3h^{3}}\int_{0}^{x}\frac{y^{4}}{(1+y^{2})^{1/2}}\text{d}y$
$\displaystyle=\frac{m_{f}c^{2}}{\lambda_{f}^{3}}\Bigl{[}\frac{1}{8\pi^{2}}\bigl{(}x(1+x^{2})^{1/2}\left(\frac{2}{3}x^{2}-1\right)$
$\displaystyle\qquad+\ln\left(x+(1+x^{2})^{1/2}\right)\bigr{)}\Bigr{]}$
$\displaystyle=\frac{m_{f}c^{2}}{\lambda_{f}^{3}}\phi(x)$ (3)
$\displaystyle\varepsilon_{f}$ $\displaystyle=4\pi
c^{5}m_{f}^{4}\int_{0}^{x}y^{2}\sqrt{y^{2}+1}\text{d}y$
$\displaystyle=\frac{m_{f}c^{2}}{\lambda_{f}^{3}}\Bigl{[}\frac{1}{8\pi^{2}}\bigl{(}x(1+x^{2})^{1/2}\left(1+2x^{2}\right)$
$\displaystyle\qquad-\ln\left(x+(1+x^{2})^{1/2}\right)\bigr{)}\Bigl{]}$
$\displaystyle=\frac{m_{f}c^{2}}{\lambda_{f}^{3}}\chi(x).$ (4)

Here, $\lambda_{f}=\hbar/(m_{f}c)$ is the Compton wavelength, $x$ is the
dimensionless Fermi momentum, defined using Equation (2), and
$y=(pc)/(m_{f}c^{2})$. The total pressure in a two-component gas is then just
$P_{\rm tot}=P_{H}+P_{L}$, while the total energy density is
$\epsilon_{tot}=\epsilon_{H}+\epsilon_{L}$.

We will assume that interparticle interactions contribute at most a small
correction to the energy density and will not include their effects here. This
follows the standard white dwarf model, where the electrostatic interaction
contributes the dominant term, at least at high densities, with a correction
to the pressure on the order of the fine structure constant, $(P+P_{\rm
correction})/P\approx$0.4\text{\,}\mathrm{\char 37\relax}$$ for a hydrogen
white dwarf [38]. We have chosen this path to preserve the generalizability of
these results, since the type of interaction changes between dark matter
models.

###  II.2 Polytropic approximation

Since the full EoS is complicated, especially when considered across the
$m_{L}-m_{H}$ parameter space, it is convenient to expand the EoS as a power
series in $x$, keeping only the dominant term. Then the EoS can be
approximated as a polytropic function, $P(x)=Kx^{\Gamma/3}$, where $K$ and
$\Gamma$ are the polytropic constants. Commonly, this is rewritten in terms of
the rest mass density ($\rho_{0}=\sum m_{f}n_{f}\approx m_{H}n_{H}\approx
m_{H}n_{L}$) and the polytropic index $n$, as $P(\rho_{0})=K\rho_{0}^{1+1/n}$.

The polytropic approximation then falls into one of four limiting cases,
depending on whether the particles are highly relativistic ($x_{f}\gg 1$) or
non-relativistic ($x_{f}\ll 1$), and whether the particle masses are
substantially different ($m_{L}\ll m_{H}$) or similar ($m_{L}\approx m_{H}$).
Generally, the heavy particles are non-relativistic except in the similar-
mass, relativistic electron limit and the similar mass limit can also be
thought of as the single particle limit, approximately obtainable with the
replacement $L\rightarrow H$. From Equation (2), the relativity condition can
be written as a condition on $\rho_{0}$, with the non-relativistic (light
particle) limit as
$\rho_{0}\ll${10}^{6}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\mathcal{r}_{L}^{3}\,\mathcal{r}_{H}$
and the highly relativistic limit
$\rho_{0}\gg${10}^{6}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\mathcal{r}_{L}^{3}\,\mathcal{r}_{H}$.

Using the notation from above and defining $K_{\rm WD}(n)$ as the
$n$-dependent polytropic constant for the ordinary white dwarf, the four cases
can be written as

$P_{\rm
tot}\approx\begin{cases}\mathcal{r}_{L}^{-1}\,\mathcal{r}_{H}^{-5/3}K_{\rm
WD}(\frac{3}{2})\,\rho_{0}^{5/3}&(a)\\\ \mathcal{r}_{H}^{-4/3}K_{\rm
WD}(3)\,\rho_{0}^{4/3}&(b)\\\ 2\,\mathcal{r}_{H}^{-8/3}K_{\rm
WD}(\frac{3}{2})\,\rho_{0}^{5/3}&(c)\\\ 2\,\mathcal{r}_{H}^{-4/3}K_{\rm
WD}(3)\,\rho_{0}^{4/3}&(d)\end{cases},$ (5)

with $(a)$ being
$\rho_{0}\ll${10}^{6}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\mathcal{r}_{L}^{3}\,\mathcal{r}_{H}$,
$m_{L}\ll m_{H}$, $(b)$ being
$\rho_{0}\gg${10}^{6}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\mathcal{r}_{L}^{3}\,\mathcal{r}_{H}$,
$m_{L}\ll m_{H}$, $(c)$ being
$\rho_{0}\ll${10}^{6}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\mathcal{r}_{H}^{4}$,
$m_{L}\approx m_{H}$, and $(d)$ being
$\rho_{0}\gg${10}^{6}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$\mathcal{r}_{H}^{4}$,
$m_{L}\approx m_{H}$. Of note, the $(a)$ and $(b)$ cases correspond to the
ordinary white dwarf when $\mathcal{r}_{H}=\mathcal{r}_{L}=1$.

###  II.3 Newtonian Hydrostatic Approximation

Next, we examine the parametric dependencies of the DWD mass, radius, and
compactness. This can be accomplished by solving the Newtonian hydrostatic
equilibrium equations. Defining $m(r)$ as the total mass contained within
radius $r$, $p(r)$ as the net outward pressure, and gravitational constant
$G$, we have

$\displaystyle\frac{dp}{dr}$ $\displaystyle=-\frac{Gm}{r^{2}}\rho(r),$ (6a)
$\displaystyle\frac{dm}{dr}$ $\displaystyle=4\pi r^{2}\rho(r).$ (6b)

With the inclusion of a polytropic EoS and the definitions,
$\rho_{c}=\rho(r=0)$, $\rho=\rho_{c}\theta^{n}$, $r=a\xi$, and
$a=[(n+1)K\rho_{c}^{1/n-1}/(4\pi G)]^{1/2}$, Equations (6a) and (6b) can be
combined into the Lane-Emden equation,

$\frac{1}{\xi^{2}}\frac{d}{d\xi}\xi^{2}\frac{d\theta}{d\xi}=-\theta^{n}.$ (7)

Note that the only remaining polytropic parameter is the index; this equation
is otherwise independent of the EoS and thus mass dependencies. Numeric
integration of Equation (7) with the boundary conditions $\theta(0)=1$,
$\theta^{\prime}(0)=0$, gives the point $\theta(\xi_{1})=0$, which corresponds
to the surface of the object. Undoing the previous transformations provides
solutions for the final radius ($R$) and mass ($M$) of the DWD,

$\displaystyle R_{\rm DWD}$ $\displaystyle=\left(\frac{(n+1)K}{4\pi
G}\right)^{1/2}\rho_{c}^{(1-n)/2n}\xi_{1}$ $\displaystyle M_{\rm DWD}$
$\displaystyle=4\pi\left(\frac{(n+1)K}{4\pi
G}\right)^{3/2}\rho_{c}^{(3-n)/2n}\xi_{1}^{2}|\theta^{\prime}(\xi_{1})|$ . (8)

These equations can be rewritten in terms of the ordinary white dwarf mass and
radius using the density scaling $\rho_{\rm
DWD}=\mathcal{r}_{L}^{3}\,\mathcal{r}_{H}\rho_{WD}$ and polytropic constant
scaling

$\mathcal{r}_{K}=\frac{K_{\rm DWD}(n)}{K_{\rm WD}(n)},$ (9)

with $K_{\rm DWD}$ from Equation (5), giving

$\displaystyle R_{\rm DWD}$
$\displaystyle=\left(\mathcal{r}_{L}^{3}\mathcal{r}_{H}\right)^{(1-n)/2n}\mathcal{r}_{K}^{1/2}R_{\rm
WD}(\rho_{c})$ $\displaystyle=\mathcal{r}_{L}^{-1}\,\mathcal{r}_{H}^{-1}R_{\rm
WD}(\rho_{c})$ (10) $\displaystyle M_{\rm DWD}$
$\displaystyle=\left(\mathcal{r}_{L}^{3}\mathcal{r}_{H}\right)^{(3-n)/2n}\mathcal{r}_{K}^{3/2}M_{\rm
WD}(\rho_{c})$ $\displaystyle=\mathcal{r}_{H}^{-2}M_{\rm WD}(\rho_{c}).$ (11)

Unsurprisingly, we recover the classic Chandrasekhar mass limit scaling, Eq.
((11)), commonly seen in the literature on exotic compact objects (e.g. [20,
17, 18]) as either a lower mass bound, when discussing black holes, or an
upper mass bound when discussing other types of objects. Lastly, the
compactness is given by

$\displaystyle C_{\rm
DWD}(\rho_{c})=\left.\frac{M(\rho_{c})}{R(\rho_{c})}\right|_{\rm DWD}$
$\displaystyle=\frac{\mathcal{r}_{L}}{\mathcal{r}_{H}}C_{\rm WD}(\rho_{c}).$
(12)

As expected, when $m_{L}\to m_{H}$ we recover the single-particle limit,

$\displaystyle R_{\rm DWD}$ $\displaystyle\propto\mathcal{r}_{H}^{-2}R_{\rm
SP}(\rho_{c})$ (13a) $\displaystyle M_{\rm DWD}$
$\displaystyle\propto\mathcal{r}_{H}^{-2}M_{\rm SP}(\rho_{c})$ (13b)
$\displaystyle C_{\rm DWD}$ $\displaystyle\propto C_{\rm SP}(\rho_{c})$ (13c)

with the $m_{L}^{-2}=m_{H}^{-2}$ scaling seen in the literature [39]. Note
that while Equations 13a–13c display the scaling for dark neutron stars, whose
Fermi pressure and energy density are dominated by the neutron terms, they
would only be useful for order of magnitude estimation, because the properties
of dark neutron stars, like ordinary neutron stars, are heavily influenced by
inter-particle interactions that are lacking in this model(see e.g. [15, 22,
16, 21]).

##  III Numerical Results

Now that we have established approximate scaling relations for DWDs, we
proceed to compute their properties using a fully relativistic treatment. In
particular, the Tolman-Oppenheimer-Volkoff (TOV) equations[38, 39],

$\displaystyle\frac{dp}{dr}$
$\displaystyle=-\frac{Gm}{r^{2}}\left[\frac{\epsilon}{c^{2}}\right]\left[1+\frac{p}{\epsilon}\right]\times$
$\displaystyle\qquad\left[1+\frac{4\pi
r^{3}p}{mc^{2}}\right]\left[1-\frac{2Gm}{c^{2}r}\right]^{-1},$ (14a)
$\displaystyle\frac{dm}{dr}$ $\displaystyle=4\pi
r^{2}\left[\frac{\epsilon}{c^{\textbf{{\color[rgb]{0,0,1}2}}}}\right],$ (14b)

update the Newtonian hydrostatic equilibrium equations (Equations 6a–6b) to
add corrections due to general relativity (square brackets), necessary to
describe the gravitational field of compact objects. We have suppressed the
$r$-dependence of $\epsilon$ and $p$ for clarity. Like the Newtonian
approximation, the TOV equation can be non-dimensionalized and solved
numerically to obtain the approximate scaling relations from Section II.3,
given the corresponding polytrope[40]. As the matter density in an actual DWD
could fall anywhere between the non- and highly-relativistic limiting cases,
we need to use the full EoS and so a numerical TOV solver. To this aim we use
a modified version of the TOVL solver developed by [41] and [42].

Solving either the Newtonian approximation or the TOV equation across a range
of central densities (i.e. $\rho_{c}=\rho(r=0)$, as before) for a given
$m_{L}$ and $m_{H}$ generates a relationship in the $M-R$ space known as a
mass-radius relation. In Fig. 1, we plot some of these mass-radius relations
in the slice of the $m_{L}-m_{H}$ parameter space specified by
$m_{L}=$4.1\text{\,}\mathrm{MeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$ and
$m_{H}=$4.1\text{\times}{10}^{-3}94\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$
on logarithmic $R-M$ axes. Clearly visible are the wide ranges in $M$ and $R$,
even over a small range in the $m_{H}$ parameter space. Conversely, the
overall shape of the $M-R$ curve remains similar over that same range. There
is a noticeable plateau that appears in the $m_{H}\gg m_{L}$ regime and
disappears as $m_{H}\to m_{L}$. The plateau is due to the fact that the light
particle becomes ultrarelativistic in the core of these DWDs and the equation
of state becomes a polytrope with $\Gamma=4/3$ (see Section II.2). As
$m_{H}\to m_{L}$, the maximum mass is achieved before the particles enter the
highly relativistic regime. For example, the maximum mass for
$\\{m_{L},m_{H}\\}=\\{4.8,5.6\\}$\mathrm{MeV}$$ (solid blue line) occurs at
$x\approx 1.1$, well within the transition regime. The transition to the
single particle limit can be seen in the behavior of the radius scaling in
Figure 1b. As $m_{H}\to m_{L}$, the radii begins changing by a factor
approaching $2^{2/3}$.

(a)

(b)

Figure 1:  Example mass-radius relations for varying parameter values. Each
line in Figure 1a corresponds to the mass-radius relation for a single value
of $m_{H}$ and the fixed value
$m_{L}=$4.8\text{\,}\mathrm{MeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$
resulting from solving the Tolman-Oppenheimer-Volkoff equation over a range of
central densities. The blue, solid, markerless line is
$m_{H}=$5.6\text{\times}{10}^{-3}\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$,
with increasing $m_{H}$ shifting the relation’s maximum mass down and to the
left. Note that the maximum mass clearly scales as described in Eq. ((11)),
and the $M-R$ curves exhibit the classic white dwarf shape, even when
approaching the single-particle limit,
$m_{H}=$5.6\text{\times}{10}^{-3}\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$
(solid blue, no marker). Above the densities corresponding to the maximum mass
(to the left on the plot), the DWD is gravitationally unstable. The cutoffs on
the right simply correspond to the minimum density plotted. Figure 1b is
identical to Figure 1a, but with radius rescaled as
$R\to\mathcal{r}_{L}\,\mathcal{r}_{H}R$ and mass rescaled as
$M\to\mathcal{r}_{H}^{2}M$.

The dimensionless, gravito-electric quadrupolar tidal deformability
($\Lambda_{2}$;[43]) is also of interest for comparison with gravitational
wave observations. The calculation of $\Lambda_{2}$ requires the numeric
solution of the reduced, relativistic, quadrupole gravitational potential
($y$) differential equation from [44, 41, 45]

$\displaystyle\frac{dy}{dr}$ $\displaystyle=-\frac{y^{2}}{r}-\frac{r+4\pi
r^{3}(p-\epsilon)}{r(r-2m)}y+\frac{4(m+4\pi r^{3}p)^{2}}{r(r-2m)^{2}}$
$\displaystyle\qquad+\frac{6}{r-2m}-\frac{4\pi
r^{2}}{r-2m}\left[5\epsilon+9p+\frac{(\epsilon+p)}{(dp/d\epsilon)^{2}}\right]$
(15)

This is solved in parallel with the TOV equation to find the value of $y$ at
the surface of the object, $Y=y(R)$.

The quantity $\Lambda_{2}$ is then defined by

$\Lambda_{2}=\frac{2}{3}\frac{k_{2}}{C^{5}}$ (16)

where $k_{2}$ is the tidal apsidal constant,

$\displaystyle k_{2}$
$\displaystyle=\frac{8C^{5}}{5}\left[(1-2C)^{2}(2C(Y-1)-Y+2)\right]\times$
$\displaystyle\qquad\Bigl{[}6C(2-Y+C(5Y-8))+4C^{3}[13-11Y+C(3Y-2)+2C^{2}(1+Y)]$
$\displaystyle\qquad\qquad+3(1-2C)^{2}[2-Y+2C(Y-1)]\log(1-2C)\Bigr{]}^{-1}.$
(17)

The $\Lambda_{2}$ parameter specifies how much the object deforms in the tidal
field of a companion star and is directly related to the compactness. Black
holes, for example, have $\Lambda_{2}=0$, LIGO-detectable neutron stars in
binaries are in the range $1{10}^{4}$, and white dwarfs are $>{10}^{10}$[43,
46]. The usage of $\Lambda_{2}$ in regards to observations is explained in
section III.3.

Note that we have included two modifications to the TOVL solver from [41, 42].
First, the TOVL solver, as written, computes the solution to the second-order
Regge-Wheeler equation, instead of Equation 15. As the numerical solver had
difficulty converging on a solution in parts of the parameter space, and for
improved numerical efficiency, we use the equivalent, first-order Equation 15
instead[44, 45]. Second, Equation 17 runs into numerical precision
difficulties when the compactness is small, with, for instance, terms in the
denominator (not) canceling as they should, leading to negative values of
$k_{2}$. We replaced the analytic form of Equation 17 with a series expansion
around both $C=0$ and $Y=1$ out to fifth order. This introduces an error of
$<1\%$ into the calculations across the entire parameter space defined below.

###  III.1 Parameter Space

In Figure 2, we plot the dark white dwarf mass, compactness, and tidal
deformability computed by solving the TOV and $\Lambda_{2}$ equations across
the values $m_{L}=$
$0.511\text{\,}\mathrm{keV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ to
$5.11\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ and $m_{H}=$
$93.8\text{\,}\mathrm{keV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ to
$93.8\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ (corresponding
to $\mathcal{r}_{L}=$ ${10}^{-3}{10}^{4}$, $\mathcal{r}_{H}=$
${10}^{-4}{10}^{2}$), with the restriction $m_{L}\leq m_{H}$ (the behavior is
symmetric across the $m_{L}=m_{H}$ line). At each sampled $m_{H}-m_{L}$ point
in the parameter space, the TOV and $\Lambda_{2}$ equations are solved for
$\rho_{c}$ ranging from
${10}^{-5}{10}^{25}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$. Figure 2a
shows the mass-radius relations for three parameter cases, somewhat
representative of the three parameter space corners: light-light (yellow),
heavy-light (blue), and heavy-heavy (purple), and the Standard Model (SM)
relation (red), where light (heavy) corresponds to significantly below (above)
the Standard Model value. In Fig. 2b, we plot the maximum mass obtained for
each mass-radius relation, while in Fig.s 2c-2d, we plot $C$ and $\Lambda_{2}$
at the value of $\rho_{c}$ corresponding to said maximum mass.

Of note, the maximum mass scales predominantly with $m_{H}$, as seen by
comparing the light-light and SM cases and SM and heavy-heavy cases, and
ranges from
$3\text{\times}{10}^{-4}5\text{\times}{10}^{8}\text{\,}\mathrm{M_{\odot}}$.
The compactness also scales approximately with the ratio $m_{L}/m_{H}$, as
shown by comparing the light-light and heavy-heavy cases and predicted in the
Newtonian approximation, ranging from ${10}^{-6}0.09$. The scaling with the
ratio generally holds for $\Lambda_{2}$ as well, due to the strong dependence
on $C$, leading to a minimum of $\sim$500$$ near the $m_{L}\approx m_{H}$ line
and a maximum of $4\text{\times}{10}^{26}$ in the upper right corner. The
maximum mass configuration is attained at a density that depends on the values
of $m_{H}$ and $m_{L}$, corresponding to a central Fermi momentum of $x\sim
m_{H}/m_{L}$, rather than a fixed value. This additional factor corresponds to
a change in the central density scaling used in Equations 10–12 from
$\rho_{c}\propto r_{L}^{3}r_{H}$ to $\rho_{c}\propto r_{L}^{2}r_{H}^{2}$.
Substitution gives $C_{\rm
DWD}\propto\left(\mathcal{r}_{L}/\mathcal{r}_{H}\right)^{2/3}$, which is what
we see in Figure 2c for $m_{L}\ll m_{H}$. For $m_{L}\to m_{H}$, $C$ approaches
the single particle limit, $C=0.114$[30].

(a) Example Cases

(b) $M_{\rm max}$ in $\mathrm{M_{\odot}}$

(c) Compactness at $M_{\rm max}$

(d) $\Lambda_{2}$ at $M_{\rm max}$

Figure 2: Tollman-Oppenheimer-Volkoff and electric quadrupolar tidal
deformability ($\Lambda_{2}$) solution results for the parameter range
$m_{L}=$ $0.511\text{\,}\mathrm{keV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ to
$5.11\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ and $m_{H}=$
$93.8\text{\,}\mathrm{keV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ to
$93.8\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$, with values
for $m_{L}>m_{H}$ ignored. In Fig. 2a, we display the mass-radius relations
near the three “corners” of the $m_{H}-m_{L}$ parameter space (with
corresponding points in the parameter space subfigure) as well as
(approximately) the Standard Model white dwarf (thick, solid/largest).
Notably, the masses and radii of these objects span many orders of magnitude.
Panels (a)-(d) show the maximum mass, compactness and tidal deformability
found. The $C$ and $\Lambda_{2}$ values were plotted at the central density
corresponding to the maximum mass achieved at that value of $(m_{H},m_{L})$,
and are the maximum (minimum) possible compactness (tidal deformability) for
that parameter set.

###  III.2 Implications for Gravitational Wave Observations

In many of the dark matter models the dark sector is mostly or entirely
hidden, only observable through gravitational interactions. Thus, DWD
observations may be limited to purely gravitational techniques, like, for
example, the detection of gravitational waves from the merger of a DWD and
some other compact object. As a first step, and since the DWDs span a large
range in both mass and compactness, it is worthwhile to determine the
detectability across the microphysical parameter space. A gravitational wave
observation is detected if the signal-to-noise ratio (SNR), defined as[47]

$\langle\text{SNR}^{2}\rangle=4\int_{0}^{\infty}\frac{|\tilde{h}(f)|^{2}}{S_{n}(f)}\text{d}f,$
(18)

for a signal with strain $h(t)$ and Fourier transform $\tilde{h}(f)$ observed
by detector with sensitivity curve $S_{n}(f)$, achieves a specified detection
threshold. As the choice of threshold is somewhat arbitrary, $\text{SNR}\geq
8$ is used here to match recent LIGO usage[48, 49].

With the assumption that the strain can be approximated as originated from a
quadrupole source and truncated to Newtonian order, its Fourier transform can
be written as

$\tilde{h}(f)\approx\frac{\sqrt{5/24}}{\pi^{2/3}D_{L}}M_{C}^{5/6}f^{-7/6},$
(19)

where $M_{C}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$ is the chirp mass and
$D_{L}$ is the luminosity distance of the merger. We will restrict the
analysis to the in-spiral phase of the merger; postmerger components,
especially for high mass objects, will require numerical relativity
simulations. This reduces the integral bounds to $0<f<f_{\rm contact}$, where
$f_{\rm contact}$ is the contact frequency, i.e. the binary orbital frequency
at the termination of the in-spiral period. Using $f_{\rm contact}$ will
possibly overestimate the final SNR, since it does not, for example, account
for tidal effects (see, e.g., [50] for other choices), but it does provide a
reasonable, simple estimate[51]. In Fig. 3, we consider two identical,
maximum-mass DWDs, and plot the contact frequency given by[41, 42]

$\displaystyle f_{\rm contact}$
$\displaystyle=\sqrt{\frac{G}{4\pi^{2}}\frac{2M_{\rm max}}{(2R_{M_{\rm
max}})^{3}}}.$ (20)

Using identical, maximum-mass DWDs to calculate $f_{\rm contact}$ provides an
optimistic estimate of the maximum frequency emitted during the merger. As
real DWDs are not likely to be at the maximum mass, actual contact frequencies
will be lower. From Figure 3, we can see that the contact frequency in the
optimistic case ranges from $8\text{\times}{10}^{-8}\text{\,}\mathrm{Hz}$, for
the most massive DWDs in the top right corner, to $600\text{\,}\mathrm{kHz}$,
for the least massive DWDs in the bottom left.

Figure 3: Frequency at merger contact in $\mathrm{Hz}$ for two identical,
maximum-mass dark white dwarfs as a function of $m_{H}$ and $m_{L}$.

Figure 4 demonstrates the application of Equations 18, 19 and 20 using the
current sensitivity curves for Advanced LIGO, and the design sensitivity
curves for LISA and DECIGO [52, 53, 54]. We compute the SNR at a luminosity
distance of $\\{100,250,450\\}$\mathrm{Mpc}$$, multiplied by a factor $4/25$
to include an averaging over sky position, inclination, and polarization and
assuming a source dominated by quadrupolar radiation [55], and shade the
region satisfying the condition $\text{SNR}>8$. The contact frequency was
computed in the same manner as in Fig. 3, i.e. assuming two, identical,
maximum-mass DWDs. As hinted at in Fig. 3, since the different gravitational
wave observatories are sensitive over different frequency ranges, the
different parameter cases will be visible by different observatories. While
LIGO can observe only the $m_{L}\in$0.011\text{\,}\mathrm{GeV}$$,
$m_{H}\in$0.13\text{\,}\mathrm{GeV}$$ region, corresponding to ordinary-
neutron-star-like DWDs and the heavy-heavy case, LISA and especially DECIGO
should be able to explore a much wider range of parameter space. LISA can
probe the $m_{L}\sim${10}^{-6}0.01\text{\,}\mathrm{GeV}$$,
$m_{H}\sim${10}^{-4}1\text{\,}\mathrm{GeV}$$ regime, encompassing the light-
light and SM cases, and nicely complimenting the LIGO region. Finally DECIGO
would be able to explore below
${10}^{-6}\text{\,}\mathrm{GeV}$(${10}^{-4}\text{\,}\mathrm{GeV}$) and up to
$1\text{\,}\mathrm{GeV}$($30\text{\,}\mathrm{GeV}$) in $m_{L}$($m_{H}$) space,
verifying LIGO/LISA results and even including the light-heavy cases.

Figure 4: Gravitational wave detectability (SNR) of identical, maximum-mass
DWD mergers. Shaded regions correspond to $\text{SNR}\geq 8$. SNR contours are
derived from Eqs. ((18)) and ((19)) with a frequency cutoff of $f_{\rm
contact}$ (Eq. (20)) and plotted at three luminosity distances,
$D_{L}=\\{100,250,450\\}$\mathrm{Mpc}$$ (lighter to darker) using design
sensitivity curves for LIGO, LISA and DECIGO [52, 53, 54]. Clearly, the
different detectors will probe different regions of the parameter space in a
complementary fashion and combined, LIGO and LISA will probe much of the
$m_{L}\sim${10}^{-3}1\text{\,}\mathrm{GeV}$$,
$m_{H}\sim${10}^{-4}3\text{\,}\mathrm{GeV}$$ parameter space.

###  III.3 Universal Relations

Potential universal relations, especially those of the electric quadrupolar
tidal deformability, are of further interest for potential gravitational wave
observations. A universal relation is a relation between two or more
macrophysical properties that is generally independent of the equation of
state, and, more importantly in our case, allows the breaking of observational
degeneracies. Consider an example DWD merger, with DWD 1 having macroscopic
parameters $\\{M_{1},C_{1},\Lambda_{2,1}\\}$ and DWD 2 having
$\\{M_{2},C_{2},\Lambda_{2,2}\\}$ (assume $M_{2}<M_{1}$). The corresponding
gravitational wave detection would observe the chirp mass,
$\mathcal{M}=(M_{1}M_{2})^{3/5}/(M_{1}+M_{2})^{1/5}$, mass ratio,
$q=M_{2}/M_{1}$, and reduced tidal deformability [56],

$\displaystyle\tilde{\Lambda}$
$\displaystyle=\frac{16}{13}\frac{1}{(1+q)^{5}}\Bigl{[}\bigl{(}q^{4}(q+12)-(1+12q)\bigr{)}\Lambda_{A}$
$\displaystyle\qquad+\bigl{(}q^{4}(12+q)+(1+12q)\bigr{)}\Lambda_{S}\Bigr{]},$
(21)

where $\Lambda_{S}$ and $\Lambda_{A}$ are the symmetric and anti-symmetric
tidal deformabilities,

$\Lambda_{S,A}=\frac{1}{2}(\Lambda_{2,2}\pm\Lambda_{2,1}).$ (22)

While the mass ratio can be measured directly from the gravitational wave
signal, the tidal deformability enters at leading order in the phasing only
through $\tilde{\Lambda}$. Further, the radii and compactness of the two
objects do not directly enter into the phasing or magnitude of the signal.
This is where universal relations are useful: breaking the
$\Lambda_{A},\Lambda_{S}$ degeneracy in Equation 21 and calculating $C$ and
$R$ from the $\Lambda_{2}$ of the individual DWDs. First, the binary love
relation, discovered by Yagi and Yunes in 2015 while examining binary neutron
star properties [57], with the form

$\displaystyle\Lambda_{A}(q,\Lambda_{S})$
$\displaystyle=F_{n}(q)\frac{1+\sum_{i=1}^{3}\sum_{j=1}^{2}b_{ij}q^{j}\Lambda_{S}^{i/5}}{1+\sum_{i=1}^{3}\sum_{j=1}^{2}c_{ij}q^{j}\Lambda_{S}^{i/5}}\Lambda_{S}$
(23)

where $b_{ij}$ and $c_{ij}$ are numerical fitting coefficients and
$F_{n}(q)=(1-q^{10/(3-n)})/(1+q^{10/(3-n)})$ is a polytropic-index-dependent
controlling factor, lets Equation (21) be rewritten as a function of $q$ and
$\Lambda_{S}$ only. From this, one can solve for $\Lambda_{S}$ and then
$\Lambda_{A}$. Then the individual $\Lambda_{2,1},\Lambda_{2,2}$ can be
computed using the definitions of $\Lambda_{S},\Lambda_{A}$.

Second, in 2013, Yagi and Yunes [58] and Maselli et al.[59] demonstrated that
the relation between $C$ and $\Lambda_{2}$ was also universal. This relation,
which follows mostly from Equation (16) provides an estimate of $C_{1}$ and
$C_{2}$ from $\Lambda_{2,1}$ and $\Lambda_{2,2}$. The radii then follow
directly from the definition of $C$.

In Fig. 5a, we demonstrate that the $\Lambda_{A}(q,\Lambda_{S})$ function from
Equation 23 also applies to DWDs. Here, we considered 2025 $m_{L}{-}m_{H}$
pairs in the ranges
$m_{L}\in\\{$5.11\text{\times}{10}^{-7}$,$5.11$\\}$\mathrm{GeV}$$ and
$m_{H}\in\\{$9.38\text{\times}{10}^{-3}$,93.8\\}$\mathrm{GeV}$$. For each
$m_{L}{-}m_{H}$ pair, we picked 20 random pairs of central densities and
computed
$\\{M_{1},M_{2},\Lambda_{2,1},\Lambda_{2,2},q,\Lambda_{S},\Lambda_{A}\\}$, as
well as $\Lambda_{A,{\rm fit}}$. The fit is computed using the functional form
from Equation 23 with new $b_{ij}$ and $c_{ij}$ coefficients given in Table 1.
Note that for our basic model the $F_{n}(q)$ controlling factor is $\approx 1$
and has been dropped.

| $\square_{1,1}$ | $\square_{1,2}$ | $\square_{2,1}$ | $\square_{2,2}$ | $\square_{3,1}$ | $\square_{3,2}$
---|---|---|---|---|---|---
$b_{ij}$ | $1.73\text{\times}{10}^{00}$ | $-1.57\text{\times}{10}^{00}$ | $5.48\text{\times}{10}^{-02}$ | $-5.10\text{\times}{10}^{-02}$ | $1.27\text{\times}{10}^{-06}$ | $-7.15\text{\times}{10}^{-07}$
$c_{ij}$ | $1.68\text{\times}{10}^{00}$ | $-1.42\text{\times}{10}^{00}$ | $5.47\text{\times}{10}^{-02}$ | $-5.01\text{\times}{10}^{-02}$ | $1.27\text{\times}{10}^{-06}$ | $-6.80\text{\times}{10}^{-07}$
Table 1: Fitting coefficients for the binary love relation given in Equation
(23).

(a) $\Lambda_{A}$ vs $\Lambda_{A,{\rm fit}}$

(b) $\Lambda$ vs C

Figure 5: DWD universal relations. The first panel shows that the binary Love
relation, $\Lambda_{A}=\Lambda_{A}(q,\Lambda_{S})$, is reasonably approximated
by a functional form from the neutron star literature [57, 60, 46], using the
new coefficients from Table 1, whereas the second panel demonstrates that for
DWDs $\Lambda_{2}$ can be well approximated by the fit given in Equation (24),
effectively $\Lambda_{2}\propto C^{5.1}$.

In Fig. 5b we plot $\Lambda_{2}$ versus $C$ for the 2025 $m_{L}{-}m_{H}$ pairs
mentioned previously, at logrithmically spaced central densities from
${10}^{-5}\text{\,}\mathrm{g}\text{\,}{\mathrm{cm}}^{-3}$ to the density
corresponding to the maximum mass (approximately $56\,000$ DWDs). Comparing
with the simple linear log fit,

$\log_{10}C=-0.1958\log_{10}(\Lambda_{2})-0.3931,$ (24)

which provides an excellent fit over the entire range, we see that
$\Lambda_{2}\propto C^{-5.1}$ to good approximation. This should not be
surprising; even though Equation (17) appears to show $k_{2}$ has a strong
($C^{5}$) dependence on $C$, in reality, the dependence is actually relatively
weak,111The lower order terms in the denominator cancel, leaving $C^{5}$ as
the lowest surviving term. This then cancels with the $C^{5}$ in the numerator
and $k_{2}$ approaches a nonzero constant as $C\rightarrow 0$. functionally
leaving $\Lambda\sim C^{-5}$. It is important to point out this relationship
is effectively independent of the dark parameters for this simple model and
thus also “universal”.

With these relations, we can substitute some numbers and see what sort of
constraints we might be able to place on the dark parameters from an
observation. For example, let us consider an hypothetical scenario in which a
binary DWD is detected with $q=$0.705\pm 0.04$$, $\mathcal{M}=$19\pm
4\text{\,}\mathrm{M_{\odot}}$$, and $\tilde{\Lambda}=$(9.09\pm
0.9)\text{\times}{10}^{4}$$. Using Equations (21) and (23), and simply
propagating the bounds, we would obtain $\Lambda_{S}=$(1.6\pm
0.3)\text{\times}{10}^{5}$$ and $\Lambda_{A}=1.5_{-0.3}^{+0.4}\times 10^{5}$.
By definition, we would then have $\Lambda_{2,1}=9.52_{-10}^{+6}\times 10^{3}$
and $\Lambda_{2,2}=3.1_{-0.6}^{+0.8}\times 10^{5}$, which, using our
$C-\Lambda_{2}$ relation would give $C_{1}=$0.0887\pm 0.02$$ and
$C_{2}=0.0443_{-0.009}^{+0.01}$. Using $m_{1}=q^{-3/5}(1+q)^{1/5}\mathcal{M}$
and $m_{2}=q^{2/5}(1+q)^{1/5}\mathcal{M}$, the chirp mass and mass ratio would
give $m_{1}=$26\pm 6\text{\,}\mathrm{M_{\odot}}$$ and $m_{2}=$18.3\pm
4\text{\,}\mathrm{M_{\odot}}$$. From this, the maximum mass would be at least
$20\text{\,}\mathrm{M_{\odot}}$, and we could constrain the heavy particle
mass, $m_{H}<$0.45\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$.
Further, while the conservative compactness of $C_{max}<0.1087$ would provide
a minimal constraint on the heavy to light particle mass ratio
($m_{H}/m_{L}>2$), the other end would suggest a not dissimilar constraint of
$m_{H}/m_{L}>5$, so $m_{L}$ would likely be
$90200\text{\,}\mathrm{MeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$ (assuming
$m_{H}=$0.45\text{\,}\mathrm{GeV}\text{\,}{\mathrm{\text{$c$}}}^{-2}$$).

##  IV Conclusion

We present a first look at the DWD, using a basic model of two, different-mass
fundamental fermions to explore some of the possibilities of this exotic
compact object. We determine analytic scaling relations for the mass, radius,
and compactness of the DWD as function of the Standard Model white dwarf and
the fermion masses (Equations (2)-(4) and (11)-(12)) in the non-relativistic
limit. We accomplish this by solving the Newtonian hydrostatic-equilibrium
approximation using the well-known equation of state of fermionic ideal
gasses. As expected, we recover the scaling relations found in the literature
upon approaching the single-particle limit.

We then solve the Tollman-Oppenheimer-Volkoff and tidal-deformability
differential equations numerically to obtain fully relativistic versions of
the Newtonian approximations. Using the relativistic formalism confirms the
approximate Newtonian scaling as well as highlights the large span of the
macrophysical and even binary attributes (Figures 1-3).

We further find universal relations between macroscopic properties of DWDs
that are analogous to those found for neutron stars. In particular, we
investigate the $C$ vs $\Lambda_{2}$ and binary love universal relations, with
the net result that the $\Lambda_{2}$-$C$ relationship can be well
approximated by a simple power law and that binary love relation can be well
approximated by fits from the neutron star literature (Figure 5). These
relations could be used to determine the radii of DWD from gravitational wave
observations of their mergers, thus directly constraining the masses of the
dark particles.

Lastly, we discuss the detectability of DWD binary mergers across the fermion
mass parameter space. We show that not only are DWD mergers detectable but
that, assuming design sensitivity, different gravitational wave observatories
would probe different regions of the space (Figure 4). For example, LIGO
should be able to detect mergers of high-compactness, lower-mass DWDs
corresponding to a dark light particle in the mass range
$0.011\text{\,}\mathrm{GeV}$ and dark heavy particle in the range
$0.13\text{\,}\mathrm{GeV}$, while LISA could detect both higher-mass, high-
compactness and lower-mass, lower-compactness DWDs, corresponding to light
particles that are ${10}^{-6}0.01\text{\,}\mathrm{GeV}$ and heavy particles
that are ${10}^{-4}1\text{\,}\mathrm{GeV}$. Later-generation space-based
detectors like DECIGO may be able to detect mergers across an even larger part
of the parameter space.

We have left four significant topics to future work, though two of those are
interrelated. First, is the effect of inter-particle interactions, both those
that do and do not change particle type, on the equation of state. While it is
reasonable to ignore particle-conversion interactions, as many dark matter
models do not contain them, interactions such as the dark electromagnetism of
atomic dark matter[61] or the Yukawa interactions of the model in [62] should
also be studied. The lack of such an interaction term is not fatal; after all,
the model presented here works quite well for estimating ordinary white dwarf
properties and should also work well in cases with weak dark interactions.
Conversely, dissipative interactions are the entire point of dissipative dark
matter models, necessitating their inclusion in follow-up work. Doing so in a
general fashion is non-trivial, but including additional polynomial terms into
the total pressure and energy equations as in [30], similar to the
electrostatic correction in ordinary white dwarfs[38] or the Yukawa term in
Kouvaris and Nielsen[22] would likely be a good first step.

Second, we have restricted our gravitational wave signal analysis to the
detectability of the in-spiral portion only. As demonstrated in Figure 2d, the
tidal deformability of these DWDs can be significantly larger than that of
ordinary neutron stars and black holes. As such, general template bank
searches based on ordinary binary black hole or binary neutron star mergers
may not find these objects. Additionally, the post-merger portion of the
signal may contain various features that could be strongly dependent on the
equation of state and the dark microphysics. Computing the full merger and
post-merger signal using numerical relativity simulations at several points in
the parameter space would help resolve both of these issues, providing data
for both a more targeted search and demonstrating any potential microphysical
dependence.

The remaining two issues concern the formation and populations of these
objects and using gravitational wave observations to constrain their
properties. Just as there are a number of dark matter models that have the
particle types to create DWDs, so are there a number of possible formation
mechanisms, ranging from primordial direct collapse, as in [31], to
astrophysical direct collapse, analogous to the dark black hole formation or
asymmetric stars in e.g. [12, 17, 63], to the astrophysical remnant of a dark
star, suggested in [21]. This makes estimating the DWD population highly
model-dependent. The possibility that such binaries might not be able to form
naturally at all also cannot be excluded. On the other hand, determining the
current constraints on the merger rates from LIGO observations should be
tractable, assuming the current state of nondetection holds. Likewise,
identifying a particular merger as a possible DWD merger (as opposed to an
ordinary object merger) should not be difficult, given the significant
discrepancies between DWD and ordinary compact object characteristics across
the majority of the parameter space. Distinguishing between a DWD and a dark
neutron star, or determining the dark composition pose a much higher level of
difficulty, however, given the potential overlap in macroscopic traits, and
may require population analysis, circling back to the formation problem.

While there are several major questions left to be resolved, the potential for
DWDs and their mergers to shine a light on the dark sector strongly motivates
the development of targeted search strategies in gravitational wave detectors
data.

###### Acknowledgements.

Funding for this work was provided by the Charles E. Kaufman Foundation of the
Pittsburgh Foundation. The authors also thank Rahul Kashyap and Daniel
Godzieba for their input on the TOVL and $\Lambda_{2}$ calculations. D.R.
acknowledges funding from the U.S. Department of Energy, Office of Science,
Division of Nuclear Physics under Award Number(s) DE-SC0021177 and from the
National Science Foundation under Grants No. PHY-2011725, PHY-2020275,
PHY-2116686, and AST-2108467.

## References

  * Aghanim  _et al._ [2020] N. Aghanim _et al._ (Planck Collaboration), Planck2018 results, Astron. Astrophys. 641, A1 (2020).
  * Markevitch  _et al._ [2004] M. Markevitch, A. H. Gonzalez, D. Clowe, A. Vikhlinin, L. David, W. Forman, C. Jones, S. Murray, and W. Tucker, Direct constraints on the dark matter self-interaction cross-section from the merging galaxy cluster 1E0657-56, Astrophys. J. 606, 819 (2004), arXiv:astro-ph/0309303 [astro-ph] .
  * Clark  _et al._ [2020] M. Clark, A. Depoian, B. Elshimy, A. Kopec, R. F. Lang, and J. Qin, Direct detection limits on heavy dark matter, Phys. Rev. D 102, 123026 (2020), arXiv:2009.07909 [hep-ph] .
  * Akerib  _et al._ [2020] D. S. Akerib _et al._ (LUX Collaboration), First direct detection constraint on mirror dark matter kinetic mixing using lux 2013 data, Phys. Rev. D 101, 012003 (2020), arXiv:1908.03479 [hep-ex] .
  * Berlin and Kling [2019] A. Berlin and F. Kling, Inelastic dark matter at the LHC lifetime frontier: ATLAS, CMS, LHCb, CODEX-b, FASER, and MATHUSLA, Phys. Rev. D 99, 015021 (2019), arXiv:1810.01879 [hep-ph] .
  * Undagoitia and Rauch [2015] T. M. Undagoitia and L. Rauch, Dark matter direct-detection experiments, J. Phys. G Nucl. Partic. 43, 013001 (2015), arXiv:1509.08767 [physics.ins-det] .
  * Mitsou [2015] V. A. Mitsou, Overview of searches for dark matter at the LHC, J. Phys. Conf. Ser. 651, 012023 (2015), arXiv:1402.3673 [hep-ex] .
  * Guépin  _et al._ [2021] C. Guépin, R. Aloisio, A. Cummings, L. A. Anchordoqui, J. F. Krizmanic, A. V. Olinto, M. H. Reno, and T. M. Venters, Indirect dark matter searches at ultrahigh energy neutrino detectors, Phys. Rev. D 104, 083002 (2021), arXiv:2106.04446 [hep-ph] .
  * Génolini  _et al._ [2021] Y. Génolini, M. Boudaud, M. Cirelli, L. Derome, J. Lavalle, D. Maurin, P. Salati, and N. Weinrich, New minimal, median, and maximal propagation models for dark matter searches with galactic cosmic rays, Phys. Rev. D 104, 083005 (2021), arXiv:2103.04108 [astro-ph.HE] .
  * Maity and Queiroz [2021] T. N. Maity and F. S. Queiroz, Detecting bosonic dark matter with neutron stars, Phys. Rev. D 104, 083019 (2021), arXiv:2104.02700 [hep-ph] .
  * Leane [2020] R. K. Leane, Indirect detection of dark matter in the galaxy, arXiv:2006.00513 [hep-ph] (2020).
  * Buckley and DiFranzo [2018] M. R. Buckley and A. DiFranzo, Collapsed dark matter structures, Phys. Rev. Lett. 120, 10.1103/physrevlett.120.051102 (2018).
  * D’Amico  _et al._ [2017] G. D’Amico, P. Panci, A. Lupi, S. Bovino, and J. Silk, Massive black holes from dissipative dark matter, Mon. Not. R. Astron. Soc. 473, 328 (2017), arXiv:1707.03419 [astro-ph.CO] .
  * de Lavallaz and Fairbairn [2010] A. de Lavallaz and M. Fairbairn, Neutron stars as dark matter probes, Phys. Rev. D 81, 123521 (2010), arXiv:1004.0629 [astro-ph.GA] .
  * Kouvaris and Tinyakov [2011] C. Kouvaris and P. Tinyakov, Constraining asymmetric dark matter through observations of compact stars, Phys. Rev. D 83, 083512 (2011), arXiv:1012.2039 [astro-ph.HE] .
  * Kouvaris  _et al._ [2018] C. Kouvaris, P. Tinyakov, and M. H. G. Tytgat, Nonprimordial solar mass black holes, Phys. Rev. Lett. 121, 221102 (2018), arXiv:1804.06740 [astro-ph.HE] .
  * Shandera  _et al._ [2018] S. Shandera, D. Jeong, and H. S. G. Gebhardt, Gravitational waves from binary mergers of subsolar mass dark black holes, Phys. Rev. Lett. 120, 241102 (2018), arXiv:1802.08206 [astro-ph.CO] .
  * Singh  _et al._ [2021] D. Singh, M. Ryan, R. Magee, T. Akhter, S. Shandera, D. Jeong, and C. Hanna, Gravitational-wave limit on the chandrasekhar mass of dark matter, Phys. Rev. D 104, arXiv:2009.05209 (2021), arXiv:2009.05209 [astro-ph.CO] .
  * Choquette  _et al._ [2019] J. Choquette, J. M. Cline, and J. M. Cornell, Early formation of supermassive black holes via dark matter self-interactions, JCAP 2019, 036 (2019).
  * Giudice  _et al._ [2016] G. F. Giudice, M. McCullough, and A. Urbano, Hunting for dark particles with gravitational waves, JCAP 2016, 001 (2016), arXiv:1605.01209 [hep-ph] .
  * Hippert  _et al._ [2021] M. Hippert, J. Setford, H. Tan, D. Curtin, J. Noronha-Hostler, and N. Yunes, Mirror neutron stars, arXiv:2103.01965 [astro-ph.HE] (2021).
  * Kouvaris and Nielsen [2015] C. Kouvaris and N. G. Nielsen, Asymmetric dark matter stars, Phys. Rev. D 92, 063526 (2015), arXiv:1507.00959 [hep-ph] .
  * Maselli  _et al._ [2017] A. Maselli, P. Pnigouras, N. G. Nielsen, C. Kouvaris, and K. D. Kokkotas, Dark stars: gravitational and electromagnetic observables, Phys. Rev. D 96, 023005 (2017), arXiv:1704.07286 [astro-ph.HE] .
  * Leung  _et al._ [2013] S.-C. Leung, M.-C. Chu, L.-M. Lin, and K.-W. Wong, Dark-matter admixed white dwarfs, Phys. Rev. D 87, 123506 (2013), arXiv:1305.6142 [astro-ph.CO] .
  * Tolos and Schaffner-Bielich [2015] L. Tolos and J. Schaffner-Bielich, Dark compact planets, Phys. Rev. D 92, 123002 (2015), arXiv:1507.08197 [astro-ph.HE] .
  * Bauswein  _et al._ [2020] A. Bauswein, G. Guo, J.-H. Lien, Y.-H. Lin, and M.-R. Wu, Compact dark objects in neutron star mergers, arXiv:2012.11908 [astro-ph.HE] (2020).
  * Dengler  _et al._ [2021] Y. Dengler, J. Schaffner-Bielich, and L. Tolos, The second love number of dark compact planets and neutron stars with dark matter, arXiv:2111.06197 [astro-ph.HE] (2021).
  * Gleason  _et al._ [2022] T. Gleason, B. Brown, and B. Kain, Dynamical evolution of dark matter admixed neutron stars, Phys. Rev. D 105, 023010 (2022), arXiv:2201.02274 [gr-qc] .
  * Dasgupta  _et al._ [2021] B. Dasgupta, R. Laha, and A. Ray, Low mass black holes from dark core collapse, Phys. Rev. Lett. 126, 141105 (2021), arXiv:2009.01825 [astro-ph.HE] .
  * Narain  _et al._ [2006] G. Narain, J. Schaffner-Bielich, and I. N. Mishustin, Compact stars made of fermionic dark matter, Phys. Rev. D 74, 063003 (2006), arXiv:astro-ph/0605724 [astro-ph] .
  * Gross  _et al._ [2021] C. Gross, G. Landini, A. Strumia, and D. Teresi, Dark matter as dark dwarfs and other macroscopic objects: multiverse relics?, arXiv:2105.02840 [hep-ph] (2021).
  * Brandt  _et al._ [2018] B. Brandt, G. Endrődi, E. Fraga, M. Hippert, J. Schaffner-Bielich, and S. Schmalzbauer, New class of compact stars: Pion stars, Phys. Rev. D 98, 094510 (2018), arXiv:1802.06685 [hep-ph] .
  * Carr  _et al._ [2020] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, Constraints on primordial black holes, arXiv:2002.12778 [astro-ph.CO] (2020).
  * Katz  _et al._ [2020] A. Katz, J. Kopp, S. Sibiryakov, and W. Xue, Looking for MACHOs in the spectra of fast radio bursts, Mon. Not. R. Astron. Soc. 496, 564 (2020), arXiv:1912.07620 [astro-ph.CO] .
  * Napiwotzki [2009] R. Napiwotzki, The galactic population of white dwarfs, J. Phys. Conf. Ser. 172, 012004 (2009), arXiv:0903.2159 [astro-ph.SR] .
  * McMillan [2016] P. J. McMillan, The mass distribution and gravitational potential of the milky way, Mon. Not. R. Astron. Soc. 465, 76 (2016), arXiv:1608.00971 [astro-ph.GA] .
  * Chandrasekhar [2010] S. Chandrasekhar,  _An Introduction to the Study of Stellar Structure_ (Dover Publications, Inc., 2010).
  * Shapiro and Teukolsky [1983] S. L. Shapiro and S. A. Teukolsky,  _Black Holes, White Dwarfs, and Neutron Stars_ (Wiley, 1983).
  * Reisenegger and Zepeda [2016] A. Reisenegger and F. S. Zepeda, Order-of-magnitude physics of neutron stars, Eur. Phys. J. A. 52, 10.1140/epja/i2016-16052-y (2016), arXiv:1511.08813 [astro-ph.SR] .
  * Tooper [1965] R. F. Tooper, Adiabatic fluid spheres in general relativity., Astrophys. J. 142, 1541 (1965).
  * Damour and Nagar [2009] T. Damour and A. Nagar, Relativistic tidal properties of neutron stars, Phys. Rev. D80, 084035 (2009), arXiv:0906.0096 [gr-qc] .
  * Bernuzzi and Nagar [2008] S. Bernuzzi and A. Nagar, Gravitational waves from pulsations of neutron stars described by realistic Equations of State, Phys. Rev. D78, 024024 (2008), arXiv:0803.3804 [gr-qc] .
  * Binnington and Poisson [2009] T. Binnington and E. Poisson, Relativistic theory of tidal love numbers, Phys. Rev. D 80, 084018 (2009), arXiv:0906.1366 [gr-qc] .
  * Hinderer [2008] T. Hinderer, Tidal love numbers of neutron stars, Astrophys. J. 677, 1216 (2008), arXiv:0711.2420 [astro-ph] .
  * Lindblom and Indik [2014] L. Lindblom and N. M. Indik, Spectral approach to the relativistic inverse stellar structure problem II, Phys. Rev. D 89, 064003 (2014).
  * Godzieba  _et al._ [2021] D. A. Godzieba, R. Gamba, D. Radice, and S. Bernuzzi, Updated universal relations for tidal deformabilities of neutron stars from phenomenological equations of state, Phys. Rev. D 103, 063036 (2021), arXiv:2012.12151 [astro-ph.HE] .
  * Jaranowski and Królak [2012] P. Jaranowski and A. Królak, Gravitational-wave data analysis. formalism and sample applications: The gaussian case, Living Rev. Relativ. 15, 10.12942/lrr-2012-4 (2012), arXiv:0711.1115 [gr-qc] .
  * Abbott  _et al._ [2020a] R. Abbott _et al._ , GWTC-2: Compact binary coalescences observed by LIGO and virgo during the first half of the third observing run, Phys. Rev. X 11, 021053 (2020a), arXiv:2010.14527 [gr-qc] .
  * Abbott  _et al._ [2020b] R. Abbott _et al._ (LIGO Scientific, Virgo), Population properties of compact objects from the second ligo-virgo gravitational-wave transient catalog, Phys. Rev. X 11913, L7 (2020b), arXiv:2010.14533 [astro-ph.HE] .
  * De Luca and Pani [2021] V. De Luca and P. Pani, Tidal deformability of dressed black holes and tests of ultralight bosons in extended mass ranges, JCAP 2021, 032 (2021), arXiv:2106.14428 [gr-qc] .
  * Bernuzzi  _et al._ [2012] S. Bernuzzi, A. Nagar, M. Thierfelder, and B. Bruegmann, Tidal effects in binary neutron star coalescence, Phys. Rev. D 86, 044030 (2012), arXiv:1205.3403 [gr-qc] .
  * Barsotti  _et al._ [2018] L. Barsotti, S. Gras, M. Evans, and P. Fritschel,  _The updated Advanced LIGO design curve_ , Tech. Rep. T1800044-v5 (The LIGO Project, 2018).
  * Robson  _et al._ [2019] T. Robson, N. J. Cornish, and C. Liu, The construction and use of LISA sensitivity curves, Classical Quant. Grav. 36, 105011 (2019), arXiv:1803.01944 [astro-ph.HE] .
  * Kawamura  _et al._ [2008] S. Kawamura, M. Ando, T. Nakamura, K. Tsubono, T. Tanaka, _et al._ , The japanese space gravitational wave antenna - DECIGO, J. Phys. Conf. Ser. 122, 012006 (2008).
  * Dominik  _et al._ [2015] M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. Holz, T. Bulik, and F. Pannarale, Double compact objects iii: Gravitational wave detection rates, Astrophys. J. 806, 263 (2015), arXiv:1405.7016 [astro-ph.HE] .
  * Favata [2014] M. Favata, Systematic parameter errors in inspiraling neutron star binaries, Phys. Rev. Lett. 112, 101101 (2014), arXiv:1310.8288 [gr-qc] .
  * Yagi and Yunes [2015] K. Yagi and N. Yunes, Binary love relations, Classical Quant. Grav. 33, 13LT01 (2015), arXiv:1512.02639 [gr-qc] .
  * Yagi and Yunes [2013] K. Yagi and N. Yunes, I-love-q relations in neutron stars and their applications to astrophysics, gravitational waves, and fundamental physics, Phys. Rev. D 88, 023009 (2013), arXiv:1303.1528 [gr-qc] .
  * Maselli  _et al._ [2013] A. Maselli, V. Cardoso, V. Ferrari, L. Gualtieri, and P. Pani, Equation-of-state-independent relations in neutron stars, Phys. Rev. D 88, 023007 (2013), arXiv:1304.2052 [gr-qc] .
  * Yagi and Yunes [2016] K. Yagi and N. Yunes, Approximate universal relations for neutron stars and quark stars, Phys. Rep. 681, 1 (2016), arXiv:1608.02582 [gr-qc] .
  * Kaplan  _et al._ [2010] D. E. Kaplan, G. Z. Krnjaic, K. R. Rehermann, and C. M. Wells, Atomic dark matter, JCAP 2010, 021 (2010).
  * Kaplan  _et al._ [2009] D. E. Kaplan, M. A. Luty, and K. M. Zurek, Asymmetric dark matter, Phys. Rev. D 79, 10.1103/physrevd.79.115016 (2009).
  * Chang  _et al._ [2019] J. H. Chang, D. Egana-Ugrinovic, R. Essig, and C. Kouvaris, Structure formation and exotic compact objects in a dissipative dark sector, JCAP 2019, 036 (2019).