Gravitational wave probes of dark matter: challenges and opportunities

In this white paper, we discuss the prospects for characterizing and identifying dark matter using gravitational waves, covering a wide range of dark matter candidate types and signals. We argue that present and upcoming gravitational wave probes offer unprecedented opportunities for unraveling the nature of dark matter and we identify the most urgent challenges and open problems with the aim of encouraging a strong community effort at the interface between these two exciting fields of research.


Introduction
The direct detection of gravitational waves (GWs) has opened a new era and window into the study of the Universe [1][2][3][4], allowing us to probe in exquisite detail the nature of compact astrophysical objects and their environments. In the coming decades, a range of current and future experiments will continue to study GWs, including LIGO/Virgo [5], LISA [6,7], the Einstein Telescope [8], Pulsar Timing Arrays [9,10], and others. These experiments promise to address, and possibly solve, a variety of longstanding problems in astrophysics, cosmology, and particle physics (see e.g. Ref. [11] and references therein).
Here, we present a survey of the most promising prospects for characterizing dark matter (DM) using gravitational wave detectors. Exploiting synergies between these two very active fields of research, we review the status of six different areas of research. We show that tremendous advances in our understanding are likely using this combined approach. Our paper is organized around the following topics, which are also summarized in Fig. 1  • Primordial black holes. In Sec. 2, we discuss the connections between DM and primordial black holes (PBHs). We briefly review the currently proposed PBH formation scenarios, the strategies to discriminate this population from astrophysical black holes using GWs and other probes, and the implications that their discovery would have for particle DM.
• Environmental effects. In Sec. 3, we review the so-called "environmental" effects induced by DM around compact objects -such as the de-phasing of gravitational waveforms and the spin-down of black holes due to super-radiance -and the prospects for probing them with GW experiments.
• Exotic compact objects. In Sec. 4, we discuss the possibility that DM forms compact objects, or that it accretes onto and modifies astrophysical objects, and explore the ensuing gravitational wave signatures.
• Direct detection with GW experiments. In Sec. 5, we review ideas on how to repurpose GW experiments (inteferometers and Pulsar Timing Arrays) to search for the signature of DM not via gravitational wave signals, but "directly," e.g. through the effects induced by DM particle physics couplings, or time-varying gravitational potentials.
• Non-perturbative DM dynamics. In Sec. 6, we consider scenarios in which DM particles are produced non-perturbatively and discuss the detectability of the stochastic gravitational wave background that they would produce.
• Phase transitions. Finally, in Sec. 7, we review the prospects for detecting a stochastic GW background induced by cosmological phase transitions, and the implications that this detection would have for DM.
For each of these areas, we identify the most urgent challenges, open problems and the prospects for resolution for these outstanding questions. We argue that present and upcoming gravitational wave probes offer unprecedented opportunities for DM studies, and we encourage a strong community effort at the interface between these two exciting fields of research.

Primordial Black Holes
Primordial black holes (PBHs) have been of longstanding theoretical interest [12,13], particularly as a possible component of DM [14]. Such interest has been reinvigorated with the dawn of GW astrophysics, as PBHs have the potential to produce signals for current and future GW experiments [15][16][17]. The existence of non-baryonic matter (i.e. matter that does not appreciably interact with the electroweak or strong sectors of the Standard Model) can be seen as early as Big Bang Nucleosynthesis [18], indicating that in order for PBHs to be a viable DM candidate they must have formed less than ∼1 second after the Big Bang. PBHs with masses below ∼ 10 −16 M cannot constitute a substantial fraction of the DM, because they are not stable to Hawking evaporation on cosmological timescales and this can drastically change cosmological observables [19]. Meanwhile, the existence of DM dominated ∼ 10 6 M dwarf galaxies like Segue 2 [20] indicate that PBHs cannot constitute a large fraction of the DM above this mass scale. Between these limiting masses, there are many constraints on the fraction of DM that can be comprised of PBHs, f PBH : these limits on f PBH arise from constraints on extragalactic gamma rays [19]; gravitational microlensing [21][22][23] and lensing of type Ia supernovae [24]; the survival of star clusters in dwarf galaxies [25] and of wide stellar binaries [26]; and the cosmic microwave background [27,28]. For the remainder of this section, we focus on the possibility of detecting PBHs via their GW signatures and the implications of such a discovery for particle DM candidates.

Formation
PBHs have a number of formation scenarios, each of which could leave slightly different GW signatures. Large curvature perturbations can be generated in hybrid [29,30] and axioncurvaton [31] models of inflation, which can lead to the formation of PBHs. Phase transitions in the early universe can lead to topological defects such as cosmic strings, loops of which can collapse into PBHs [32,33] (and the dynamics of the loops may be a potential source of gravitational radiation). PBHs may also be formed by the fragmentation of extended solitonic and pseudo-solitonic field configurations [34][35][36], or as a byproduct of the metastability of the electroweak vacuum [37,38]. These scenarios can possibly be distinguished by the shape of the inferred PBH mass spectrum, which impacts any potential GW signal from PBH merger events. Further, PBHs formed in a matter-dominated era (e.g. [34][35][36]) are expected to possess large spins [39], unlike PBHs formed in the radiation-era as with standard inflationary perturbations; the LIGO Collaboration has placed some limits on the spin distribution of binary black holes through the merger events detected in Observing Runs 1 and 2 [40]. Additionally, some formation scenarios are accompanied by a tensor mode [41], which impacts the gravitational wave background and can thus be probed with pulsar timing and interferometric searches [42][43][44][45][46][47].

Detection
GWs have the potential to provide evidence for PBHs in the near future, as there are at least two detection scenarios that appear incompatible with black holes of astrophysical origin. The first scenario is the detection of nearby mergers of sub-solar mass black holes (BHs) with LIGO/Virgo: according to the standard theory of stellar evolution, black holes do not form below the Chandrasekhar mass of ∼ 1.4 M , and therefore the detection of the merger of subsolar mass BHs using LIGO and Virgo would strongly suggest a primordial origin (though see Refs. [48,49] for alternative exotic formation mechanisms). The second detection scenario is the observation of BH-BH mergers at high redshift with next generation interferometers such as Einstein Telescope (ET) [50] and Cosmic Explorer (CE) [8]: above a redshift of z ∼ 40, the merger rate of astrophysical BHs should be negligible [51]. Very high-redshift mergers would therefore point towards a primordial BH population which formed much earlier [52].
Other strategies to discriminate between primordial and astrophysical BHs include: the analysis of their eccentricities [66], mass function [67], as well as angular momentum distributions [68,69], which are already being constrained via BH-BH merger waveforms [40]; the study of the spatial distribution and mass function of radio and X-ray sources powered by the accretion of gas onto Galactic BHs [70][71][72]; and the search for a population of compact objects with a large optical depth to GW lensing [73].

Implications for particle DM
If PBHs do not make up all of the DM, they are generically expected to accrete a dense mini-halo of DM. This begins before matter-radiation equality, as the sphere of gravitational influence of each PBH grows [74,75]. Analytic considerations suggest that such mini-halos would have a very steep central density profile ρ(r) ∝ r −9/4 [76], which has recently been confirmed in numerical simulations [77,78]. The growth of the halo is expected to continue at least until PBHs are subsumed into bound structures, at which point the DM halo may be 100 times more massive than the PBH itself [74].
Such large over-densities would have profound implications for particle DM. Enhanced DM annihilation within mini-halos would potentially make solar-mass PBHs bright γ-ray sources [79][80][81], as well as contributing to the diffuse γ-ray background [82][83][84]. By comparing with γ-ray observations, the abundance of PBHs can be constrained below the level of f PBH 10 −8 , assuming that the DM is a Weakly Interacting Massive Particle (WIMP) produced as a thermal relic [78,[85][86][87]. This suggests that PBHs and WIMPs are fundamentally incompatible. More recently, it has been pointed out that if PBHs are detected in near-future GW searches (as discussed above), this would point towards an abundance greater than f PBH 10 −5 [88]. This in turn would place stringent constraints on models of WIMP DM, ruling out thermal WIMP DM with GeV-scale masses and above. Such a detection would also strongly constrain large regions of the parameter space for models such as the minimal supersymmetric standard model, even when they predict only a sub-dominant population of WIMPs.
The existence of PBHs may also have more indirect implications. In scenarios where PBHs are formed from enhanced primordial perturbations (Sec. 2.1), the same rare, large density fluctuations which produce PBHs should also lead to the formation of gravitationally bound ultra-compact mini-halos (UCMHs) of particle DM [77,[89][90][91][92]. This in turn would lead to enhanced lensing and annihilation signatures due to these UCMHs (see e.g. [93,94]). Finally, the detection of PBHs would shed light on a number of models that suggest a common origin for PBHs and particle DM [95,96].

Environmental Effects
Future GW experiments should allow for precise reconstructions of the properties of inspiraling and merging compact objects, including mass measurements at the sub-percent level [97]. With such exquisite precision, these experiments would be sensitive to tiny deviations in the gravitational waveforms which may be induced by matter in the violent environments of the merging BHs and NSs [98,99]. These "environmental effects," once detected, could perhaps signal the presence of DM, although similar effects due to baryons must also be taken into account.

Cold DM
As discussed above, large overdensities of cold DM may form around primordial black holes. Numerical and analytical studies of the mergers of "dressed" black holes show that the distribution of DM around them dramatically affects the dynamical evolution of the binaries [62].
Cold DM overdensities are also expected to inevitably form around astrophysical black holes, although the slope and normalization of the DM density profile depends strongly on the formation mechanism. This process has been explored in the context of supermassive black holes lying at the centers of galaxies [100][101][102][103][104][105][106], as well as for intermediate-mass BHs [107][108][109].
The dynamical friction induced by DM particles is expected to modify the dynamics of the merger, possibly leaving an imprint on the gravitational waveform, in the form of a change in phase relative to the inspiral without DM [11,98,99,[110][111][112]. Much remains however to be understood about the evolution of these systems, as only approximate solutions have been obtained thus far. The full problem of evolving the BH-BH pair -including post-Newtonian corrections as well as the gravitational feedback induced on the DM distribution -has yet to be solved. We also stress that the dephasing induced by a distribution of DM might substantially alter the waveform, to the point that dedicated templates would be needed in order to extract the predicted signal from the noise.
If this dephasing effect is observed, it would potentially allow for the detection of DM around black holes, as well as a measurement of the DM density. It has been suggested that these "dark dresses" are incompatible with light bosonic and fermionic DM and with selfannihilating DM [113] and their observation may therefore hint at the nature of the DM particle.

Ultralight bosons
Another attractive class of DM models consists of ultralight bosonic particles such as the axion, axion-like particles (ALPs), and dark photons. For DM candidates with masses below ∼ keV, the occupation numbers are high enough that one can describe the DM as a classical field. These light bosonic fields can form gravitational bound states around spinning black holes, and subsequently extract energy and angular momentum from the BH through a process known as "rotational superradiance" (see [114] for a review). Notably, this process depends only on the gravitational interaction between the bosonic fields and the black hole. 1 Super-radiance drives an exponential growth of bosonic clouds around the BH, potentially reaching masses of up to ∼ 10% of the BH, which causes the BH to spin-down [117][118][119][120]. Therefore, the existence of these bosonic fields can be probed indirectly by measuring the masses and spins of BHs. Since this process is most efficient when the Compton wavelength of the bosonic field is of comparable size to the BH, observations of stellar mass black holes (supermassive black holes) set limits on the scalar masses that apply in the range ∼ 10 −13 − 10 −11 eV (10 −18 − 10 −16 eV) [121,122]. The superradiance process is more effective for vector particles, and in that case the limits shift to ∼ 10 −14 − 10 −11 eV (10 −20 − 10 −17 eV), for stellar mass (supermassive) black holes, respectively [123].
In addition, the oscillations of the bosonic clouds can source GWs that are detectable with LIGO, LISA, and other gravitational wave detectors [121][122][123]. There are three types of signals that can arise from the bosonic cloud: graviton emission from level transitions, boson annihilations into gravitons, and a bosenova collapse of the boson cloud [117]. The first two processes result in monochromatic GWs. The frequency of the signal is determined by the masses of the boson and black hole. For stellar mass black holes, this corresponds to the frequencies probed by Advanced LIGO while supermassive black holes correspond to the frequencies of LISA. The third process is a consequence of the self-interactions; if the attractive self-interactions are stronger than the gravitational binding energy, the bosonic cloud collapses, resulting in a burst of GWs. The amplitude of these waves is smaller, but may be observable for supermassive black holes. The gravitational wave signatures are also altered if the black hole is a part of a binary system. Most notably, the existence of a companion black hole induces resonant mixing between the growing and decaying modes of the bosonic cloud [124,125]. Some consequences on the gravitational wave signal are a Doppler modulation of the frequency, as well as modifications to the waveforms from the cloud's multipole moments and tidal deformations caused by the companion.
DM in the form of ALPs may also give rise to electromagnetic signatures around compact objects, creating a new opportunity for multimessenger astronomy. For example, dense clouds of ALPs (grown by superradiance around rotating BHs) may lead to stimulated photon emission. This gives rise to a periodic radio signal, relevant for axion masses above ∼10 −8 eV and BH masses below ∼0.01 M [126,127]. Another possibility is that ALPs may be converted into radio photons via the Primakoff effect [128]. This process requires strong magnetic fields, making the magnetospheres of neutrons stars (NSs) a promising environment for this axionphoton conversion [129][130][131]. Around intermediate mass black holes (IMBHs), the DM density is expected to be significantly enhanced (see Sec. 2.3 and Sec. 3.1) leading to an even larger predicted signal in NS-IMBH binaries. GW observations of such a system could provide a measurement of the DM density (through the dephasing effect), fixing the normalization of the expected radio signal [132]. Joint observations of these systems using LISA and the Square Kilometer Array [133] would thus probe the natural parameter space of the QCD axion in the mass range m a ∈ [10 −7 , 10 −5 ] eV.

Baryons
Baryons can also play a major role in driving the coalescence of Massive Black Hole Binaries (MBHBs), whether in the form of stars or gas. The dominant mechanism for driving a change in the orbital dynamics of MBHBs has yet to be fully understood, but several channels are possible. Stars that are ejected from the binary loss cone via gravitational slingshot can lead to hardening of the binary, as the ejected star carries away energy [134,135]. This process will stall if the loss cone is depleted, but the loss cone can be efficiently replenished in sufficiently non-spherical galaxies [136][137][138]. Gas disks can also play a role via tidal torquing of the binary wave emission is relevant in this case.
if the mass of the disk is similar to the mass of the secondary BH, although this is complicated by star formation in the disk [139,140]. If the binary accretes infalling molecular clouds, this can also efficiently drive MBHBs into the regime of gravitational wave domination [141].
Baryons can have a residual effect on the MBHB even after it has entered the GW-dominated phase. For instance, while the emission of GWs tends to circularize elliptical orbits, interactions with a circumbinary gaseous disk [142] or stars [143,144] can increase the orbital eccentricity, particularly for binaries with a large mass discrepancy. Such effects may be a fingerprint for baryonic influence on the binaries, but may also make it harder to identify the subtle environmental effects expected from Dark Matter.

Exotic Binary Mergers
If dark particles coalesce into exotic compact objects (ECOs) of astrophysical size, they may form new binary systems. The gravitational radiation emitted when such binary systems merge may be observed in gravitational wave experiments. Binary black hole mergers are often separated into three sequential phases: the inspiral, the merger, and the ringdown. Gravitational wave emission builds up in the inspiral phase and peaks during the merger. The end of the inspiral phase is characterized by the innermost stable circular orbit (ISCO), which can be defined for ECOs in analogy with black holes, [145] where M 1,2 are the masses of the stars in the binary and C * = G N M * /R * is the typical compactness of an ECO of mass M * and radius R * , limited from above by black holes (C * = 1/2). ECOs with C * > 4/9 are in violation of Buchdahl's theorem [146,147] 2 ; we will describe some examples in Section 4.3. The best detection prospect for an ECO merger is for f ISCO within the frequency window of an experiment. This implies that ground-based interferometers are most sensitive to 10 < f ISCO /Hz < 10 3 , which includes BNS, BBH, and solar mass-sized ECOs, whereas Pulsar Timing Arrays (PTAs) are sensitive to supermassive BHs and other compact objects with M * ∼ 10 6 M . The experimental prospects also rely on the formation history and abundance of the ECOs. ECOs may form at redshifts long before first-star formation, and depending on the fraction of DM that they constitute, may not follow a Navarro-Frenk-White profile [148]. ECOs with masses close to solar masses are also constrainable via the non-observation of microlensing events.

Exotic stars
Real and complex scalar fields provide excellent DM candidates; these fields can support ultracompact, coherent solitonic configurations held together by gravity (oscillatons & boson stars) [149,150], self-interactions (oscillons & Q-balls) [151][152][153][154], or both. When gravity is the dominant force holding the solitons together, their mass is given by M 0.6m 2 pl /m φ [155] where m φ is the mass of the scalar. Their compactness C * can be comparable to that of neutron stars (C * ∼ 0.1). For sufficiently compact solitons, and a small mass m φ (M ∼ m 2 pl /m φ ), such objects can source detectable GWs for LIGO/LISA via their mergers with one another, or with other compact objects (see for example, [156][157][158][159][160]).
Boson stars stabilized by a repulsive self-interaction can form compact, astronomically sized Bose-Einstein condensates with mass M * = O(10 −2 ) λ φ M 3 p /m 2 φ [161]. The maximum compactness of such spherically symmetric stable stars is C * 0.16 [162]. In this case, it is seen that the binary merger signal peaks at f peak ∼ 10 −14 (m φ /eV) 2 / λ φ Hz, which falls within the LISA window for m φ ∼ λ The differences in tidal deformability of scalar stars [164][165][166] and the possibility of energy loss through scalar radiation provide avenues for distinguishing their gravitational wave signatures from those of neutron stars and black holes. In certain cases, the gravitational wave output from the mergers of such scalar stars can exceed that of their corresponding mass black-hole counterparts [158,159].
Finally, asymmetric DM comprised of fermions can also form stable compact objects [167,168]. Typically, self interactions are required to allow for efficient formation, see for example Ref. [169] for a recent discussion. The resulting fermion stars, due to Fermi repulsion, are in general less compact compared to boson stars in the absence of additional strong attractive interactions. As for the case of boson stars these scenarios can be distinguished through differences in tidal deformability [170].

Captured DM
Small (sub-Chandrasekhar-mass) PBHs constituting a significant fraction of DM could be efficiently captured in DM-rich environments by neutron stars (or white dwarfs). A captured black hole, growing within the star, will eventually consume the host. This can lead to novel multi-messenger coincidence signatures with GWs (e.g. a kilonova without a binary merger GW counterpart) as well as the appearance of binaries with solar-mass black holes that typically do not arise in standard astrophysics [171][172][173].
Alternatively, neutron star properties can be modified in the presence of particle DM. This can lead to observable deviations in the BNS mergers probed by both current and future GW observatories. Effects considered thus far in the literature include modifications to the tidal deformability through the presence of DM in the NS core [174,175] as well as by extended DM clouds around coalescing NSs [176], axionic induced fifth-forces [177,178] and long-range dark forces affecting the inspiral phase [179][180][181][182][183][184]. In addition, many of these effects will also be present in extreme-mass-ratio inspirals (EMRIs), where the duration of the waveform will be observable for significantly longer periods of time in comparison to LIGO/VIRGO NS-NS binaries.

Black hole mimickers
Several proposed horizonless compact objects mimic the space-time of black holes at large distances as well as in the vicinity of the photon ring. Such objects are referred to as clean-photon sphere objects (ClePhOs). ClePhOs violate the Buchdahl limit by violating one or several of its axioms [147], such as in alternative gravitational theories [185]. Abundant ClePhOs may form a sub-fraction of DM, and will be subject to the same constraints as PBHs; we will describe a few examples here. Gravastars [186] are supported by negative pressure from radiative QFT effects in curved space-times, and do not necessarily require new physics. Stable gravastars may exist for a range in compactness [187]. Wormholes [188] connect different regions of space-time. Solutions with different geometries exist, and can have any mass or compactness. Their stability and formation depends on the theory of (modified) gravity, though generically they are unstable [189]. Anisotropic stars are ECOs subject to large anisotropic stresses. Although covariant studies of anisotropic stars are challenging, it is believed they can exist for a wide range in mass and compactness [190].
ClePhOs may partake in binary mergers, and can be distinguished from black holes through their nonzero tidal Love numbers [191], which leads to higher order post-Newtonian corrections (tidal heating at 2.5PN, tidal deformability at 5PN) [192]. EMRIs hold a unique potential to probe black hole mimickers, as any nonzero tidal Love number of the central object leads to large post-Newtonian deviations of the gravitational waveform in the long inspiral phase, creating a clear signature of a non-standard BH-BH inspiral scenario [193].

Searches with interferometers
Cold DM particles with sub-eV masses feature in general large occupation numbers of lowmomentum states. This is a consequence of the high number densities required to yield the observed DM energy density. Sub-eV cold dark matter (CDM) is hence typically described in terms of classical fields rather than distinct particles (as already discussed in Sec. 3.2). For DM masses around 10 −13 to 10 −12 eV, the DM field oscillation frequencies match the best sensitivity of LIGO (∼ 100 Hz). The DM fields are expected to be in coherent oscillation over length scales of coh. ∼ 10 6 km in the Milky Way, and might exhibit topological defects depending on the details of the production mechanism [194]. A wide class of sub-eV CDM models exist, and, as we will discuss below, the optical cavities of gravitational wave detectors have a unique potential to provide complementary probes to previously unconstrained parts of the sub-eV DM parameter space. However, there are also situations where significantly heavier DM can be constrained using interferometers. For instance, if DM is composite with a super-Planckian mass (e.g., if DM is made of "dark blobs") and interacts with baryons via a longrange mediator, it can induce an observable acceleration on interferometer elements [195]. Furthermore, Ref. [196] recently proposed optical cavities as sensitive probes to detect the Brownian motion caused by interacting electroweak-scale DM particles.
Dark photons. Ultra-light dark photons, produced via the mis-alignment mechanism, can constitute all of dark matter [197][198][199]. If this dark photon DM (DPDM) is associated with the U(1) B or U(1) B−L symmetries, DPDM couples directly to baryon or neutron number and hence acts as a fifth force. The strongest constraints on such scenarios come from tests of the weak equivalence principle [200,201]. In Ref. [202] it was shown that DPDM with masses around 10 −13 to 10 −12 eV can potentially lead to GW signatures that are similar to monochromatic stochastic GWs. Since the coherence length of DPDM in this mass-range is much larger than the separation between various GW detectors on Earth (or the future LISA detector), the crosscorrelation between measurements from individual GW detectors can be used significantly enhance the signal-to-noise ratio of the otherwise stochastic signal. The first searches for DPDM, using LIGO's first observing run from the detectors in Hanford and Livingston, O1, have been conducted [203], and the constraints obtained already exceed those of fifth force searches for DM masses m ∼ 10 −14 − 10 −13 eV. Both methodological improvements of the analysis technique and additional data have the potential to further strengthen constraints on the coupling parameter ε 2 by more than an order of magnitude. On the other hand, the future LISA detector will probe DM masses around m ∼ 10 −17 eV.
Axions. The presence of axion DM affects the propagation of photons by inducing minuscule changes in the phase velocity of circularly polarized light. The authors of Ref. [204] (see also Refs. [205][206][207]) have proposed a new scheme for exploiting and enhancing this effect in the linear optical cavities of gravitational wave detection experiments. The basic idea is to inject a linearly polarized laser beam into the optical cavities and search for the generation of orthorgonally polarized light, which would be a signature of axion DM. The detection strategy is thought to be resilient against the most common systematics that plague GW detectors, since they would affect both circular polarizations in similar ways. Projected sensitivies indicate that this search strategy can probe axion masses below around 10 −11 eV-up to three orders of magnitude below the current helioscope bound from the CAST experiment [208,209].
General searches. Ref. [210] (see also Ref. [211]) studies the effects of bosonic sub-eV DM fields coupling to the photon kinetic or the fermion mass terms. Coherent oscillations of dark matter fields or the collisional encounters of topological defects induce spatial and temporal variations of physical constants. Recently, various clock-clock comparison experiments have been conducted to search for time varying physical constants [212,213]. These experiments are mainly sensitive to sub-Hz frequencies, with DM masses m 10 −15 eV. Searches with the cryogenic resonant-mass detector AURIGA [214] provide the most sensitive constraints on a narrow frequency around a kHz. In the optical cavities of GW detectors, the variation of physical constants leads (like in the case of DPDM) to additional forces on the freely suspended mirrors, as well as variations of the extent of the beam splitter (through variations of the Bohr radius) and its optical index. Ref. [210] shows that (Fabry-Perot-)Michelson interferometers, like GEO 600 and Advanced LIGO (for DM masses m = 10 −13 − 10 −11 eV) or the future LISA (for m = 10 −16 −10 −18 eV), have the potential to probe scalar DM models in previously unconstrained regions of the parameter space. Minor modifications of the experimental configurations (e.g. changing the thickness of Fabry-Perot mirrors in one of the arms) can lead to significant enhancements of the DM reach.

Searches with PTAs
Millisecond pulsars are excellent clocks, emitting radiation at regular intervals over long periods of time. Given their stability, these objects are ideal for searches of GWs in the nHz frequency range [215,216]. The presence of GWs causes variations in the time of arrival of the electromagnetic pulses, and worldwide efforts are underway to detect the stochastic GW background using pulsar timing arrays (PTAs). Notably, the International Pulsar Timing Array (IPTA) [217] is a consortium of three collaborations: the North American Nanohertz Observatory for gravitational waves [218], the European Pulsar Timing Array [219], and the Parkes Pulsar Timing Array [220,221]. The first IPTA data release contains 49 millisecond pulsars that have been observed for 5-30 years with µs precision [217]. Future experiments such as the Square Kilometre Array (SKA) will achieve a much greater sensitivity over current PTAs by finding many new stable millisecond pulsars with better timing precision [222]. While one of the main goals of PTAs is to detect the stochastic GW background, expected to arise primarily from binaries of supermassive black holes, PTAs are also sensitive to the time-varying gravitational potentials that arise in certain models of DM.
Substructure and compact objects. Massive objects, such as DM subhalos, UCMHs, and PBHs, passing near the Earth-pulsar system can cause a shift in the expected time of arrival of the pulses through two possible mechanisms: the Shapiro time delay [223] and the Doppler effect [224]. The Shapiro delay occurs when the travel time of a pulse is altered due to the presence of a gravitational potential of a transiting object [223,[225][226][227][228][229], while the Doppler effect arises from the object accelerating the Earth or pulsar as it passes by [224,225,229]. These two detection strategies cover a wide range of object masses, from 10 −12 M to above 100 M , complementing lensing searches at low masses and LIGO and LISA searches at high masses.
Fuzzy DM. Another possibility is DM as an ultralight scalar field. For particles of mass m and velocity v, the corresponding de Broglie wavelength is λ dB ≈ 600 pc 10 −23 eV m The wave nature of the DM stabilizes it from collapse on scales of λ dB , smoothing out inhomogeneities on smaller scales and thereby suppressing structure [230]. For DM within the Galaxy, the scalar field behaves classically. Its pressure oscillates with an angular frequency ω ≈ 2m and induces oscillations in the gravitational potential [231]. DM masses of m ∼ 10 −23 eV are particularly interesting, since the oscillation frequency f ∼ 5 × 10 −8 (m/10 −22 eV) Hz is in the sensitivity range of PTAs. Current limits constrains the ultralight DM density to be below 6 GeV/cm 3 for masses m 10 −23 eV [232].

Non-perturbative DM dynamics
Non-perturbative production and dynamics of DM can give rise to an observable stochastic GW background. In this section we consider two main possibilities: GWs sourced by the breakdown of coherent oscillations of a scalar DM candidate, and GWs sourced by the non-perturbative production of DM. In each of these scenarios, further dedicated studies will be necessary to pin down the exact relationship between the GW spectrum and the mass and dynamics of the particular DM candidate in question. Nevertheless, the detection of a signal will be a clear sign of new physics. Homogeneous, oscillating scalar fields provide attractive DM candidates. In the absence of self-interactions, their energy density clusters gravitationally on cosmological time-scales, essentially behaving as CDM on scales larger than the de Broglie scale [230]. However, when such fields have significant attractive self-interactions, their homogeneous oscillations are unstable to spatial perturbations, leading to self-resonance (parametric and tachyonic), rapid fragmentation, and spatial clustering in the condensate (including formation of solitonic configurations) [233][234][235][236]. Such rapid fragmentation can source a stochastic background of GWs [237][238][239].
Bosonic DM can also be efficiently produced as a daughter field through non-perturbative mechanisms akin to preheating in the early universe (for a review, see [240]). For example, such mechanisms have been proposed as an efficient means of producing cold vector DM as light as 10 −20 eV [199,229,241,242]. Such exponential particle production and ensuing nonlinear dynamics can again source a stochastic gravitational wave background due to time dependent anisotropies in the energy-momentum tensor [243][244][245]. A Chern-Simons coupling between an axion and a dark-photon field strength tensor of the form (a/ f a )X µνX µν , where f a is the axion decay constant, causes a tachyonic instability in one of the dark-photon helicities. This leads to a spectrum of chiral GWs where the amplitude is controlled by f a /m pl , while the frequency is determined by the axion mass. In particular, the QCD axion can be probed by PTAs if the Chern-Simons coupling is large enough [246].
Nonlinear dynamics of energy densities resulting from the fragmentation of a condensate and production of daughter fields can lead to a detectable stochastic background only if (i) a substantial fraction of the total energy density participates in the production of GWs; (ii) it does so at sufficiently late times (low energies); and (iii) the fragmentation scale is not too small compared to the size of the horizon at the time of production [238,240,247,248]. For fields that eventually constitute all of the DM in the late Universe, however, it is non-trivial to get their energy fraction to be significant in the early Universe.

Phase Transitions
Phase transitions (PTs) occur when the vacuum state of a theory changes, for example, when a symmetry breaks spontaneously. Phase transitions that feature a discontinuity in the first derivative of the free energy are first-order and inhomogeneous. Bubbles of the new vacuum nucleate in a background of the old vacuum, and as the new vacuum is energetically favored, they expand. Gravitational wave radiation is associated with the collisions of bubbles, as well as the acoustic waves and turbulence in the plasma coupled to the bubble wall.
The GW spectrum from a first order phase transition (FOPT) is expected to follow a broken power-law, which peaks at a frequency roughly set by the inverse size of the bubbles at collision redshifted to the present time, f peak ∼ R −1 * (a * /a 0 ) [249][250][251][252][253][254][255]. For transitions which occur during radiation-domination, such as weakly-first-order phase transitions, the peak frequency is predominantly set by the nucleation temperature. PTAs are sensitive to T N ∼ 10 −6 − 10 −4 GeV, space-based interferometers such as LISA are sensitive to PTs around the EW-scale (T N ∼ 10 −1 − 10 3 GeV), and ground-based interferometers are sensitive to PTs at higher scales T N ∼ 10 5 GeV [256].
Because the Standard Model (SM) does not feature a first-order phase transition [257,258], an observed GW background of this kind would point uniquely to new physics (for examples see Refs. [259,260]). An observable background also implies that at least an O(10 −2 − 10 −1 ) fraction of the energy density of the Universe was coupled to the order parameter. In this section we review scenarios in which dark matter or a hidden sector generates a FOPT.

The electroweak phase transition
A strongly first-order electroweak phase transition is a necessary ingredient for electroweak baryogenesis. As such, there are a variety of well-studied mechanisms for generating such a phase transition. For example, it is well known that the addition of a scalar singlet to the SM can promote the SM electroweak phase transition from second-order to first-order by providing an additional cubic term to the effective potential (see e.g. [261][262][263]). If this singlet has an additional Z 2 -symmetry, it can also serve as a DM candidate [264], although this condition necessitates that the singlet is a sub-dominant component of the DM relic abundance [265]. A singlet scalar can generate Majorana neutrino mass and produce sterile neutrinos with the correct DM relic abundance [266][267][268].
Models in which new scalars are charged under the SM gauge groups can also impact the electroweak phase transition; the simplest realization of this type of scenario is the addition of a second scalar doublet to the SM Higgs sector, known as the two-Higgs-doublet model (2HDM). There are many variations of 2HDMs, several of which contain DM candidates [269] and can lead to a sizable gravitational wave signature [270]. Thermal loops of bosons from beyond-Standard-Model (BSM) theories are also a source for a large cubic term in the finite temperature effective potential, and occur in the minimal supersymmetric extension of the SM (MSSM) through the stop squark loop [271], provided that the mass of the lightest stop is below that of the top quark. However, this minimal scenario is in tension with LHC data. Instead, one can consider singlet extensions of the MSSM, for which there are still viable regions of parameter space in which strong electroweak phase transitions are allowed [272]. Finally, the electroweak phase transition can be modified to feature a nearly-conformal potential, inducing a very strong first-order electroweak phase transition. Such possibilities include which also include a dark matter state are composite or Randall-Sundrum models [273][274][275] and extended gauge sector models [276,277]. The key point is that for typical electroweak phase transitions, the relevant dimensionful parameters are near the electroweak scale, and thus could potentially be probed by space-based interferometers like LISA.

Phase transitions in hidden sectors
There are a variety of dynamics associated with a FOPT in a hidden sector that can directly affect the production of the observed DM relic abundance. The first, most widely explored opportunity, is a dark Higgs generating masses in the hidden sector. In this context, the phase transition is not necessarily first-order, nor does it provide direct information into the DM micro-physics. However, models featuring a large number of either gauge bosons or scalars thermally induce sizable potential barriers yielding strong first-order phase transitions [278,279]. The gauge structure of such models, as well as any additional field content of the theory, all lead to observable deviations in the GW signal [279]. Note that producing an observable GW signal requires that the hidden sector is at least partially thermalized, so BBN and ∆N eff constraints place strong lower bounds on the scale at which the phase transition can occur [280,281].
Another distinct opportunity relies on phase transitions that occur at a much lower temperature than the critical temperature-a phenomenon termed supercooling. Given sufficient supercooling, the phase transition drives a period of inflation which dilutes the DM relic density, followed by insufficient reheating to re-thermalize the DM [282][283][284]. For example, the supercooling of the electroweak phase transition down to QCD temperatures leads to sizable GW signatures and affects the abundance of QCD axion DM [283,285,286]. Alternatively, the collision of sufficiently energetic bubble walls can lead to non-thermal production of DM [287][288][289][290]. In the latter case, DM significantly heavier than the scale of the electroweak phase transition (M DM 10 8 GeV) can be produced.

Hidden sector confinement
If the dark sector features a gauge coupling that grows large in the IR, i.e. dark QCD, confinement of dark quarks or gauge bosons will occur. An analytic argument by Pisarski and Wilczek in 1983 [291] determines the confining phase transition to be first order for N F ≥ 3 light quark flavors at the time of confinement. 3 Hidden sector confinement has been studied recently in the context of dark quark nuggets [293], as well as solutions to the strong CP problem [294,295]. Models of confinement invariably predict states around the confinement scale, as well as lighter pions. Therefore, a GW spectrum from a confining phase transition could motivate a collider search for states charged under dark QCD.
Another potential avenue for studies of DM utilizes the accumulation of DM in front of the bubble wall. As the bubble walls collide, the DM becomes trapped, leading to the formation of large bound states, which can have implications for a number of search strategies such as microlensing and gravitational wave detection [293,296,297].

Discussion and conclusions
We summarise current and future avenues for constraining dark matter using gravitational waves in Fig. 2. The range of possible constraints is impressive, as it spans a wide range of detection techniques and DM candidates and almost 90 orders of magnitude in DM candidate mass.
We also remark that GWs allow us to constrain alternative theories of gravity, some of which have been proposed as alternatives to particle DM. The coincident observation of electromagnetic radiation and GWs from the merger of two neutron stars [2] has for instance already provided stringent constraints on theories of modified gravity in which photons travel on different geodesics with respect to GWs [298][299][300].
The connection between DM and GWs remains, however, largely unexplored, and much
remains to be done to improve existing constraints and to assess in a more complete and robust fashion the prospects for identifying DM using future experiments. We highlight here a list of challenges and open questions for each of the research areas discussed above: • Primordial black holes. What is the initial distribution of PBHs and how does it affect the GW signal and other probes? Can observations help us discriminate between different PBH production scenarios? How can we probe sub-lunar mass PBHs constituting DM? What type of radiation is emitted by accreting PBHs? Are PBHs the seeds of supermassive BHs? How do constraints on PBHs change in the presence of other DM candidates?
• Environmental effects. What is the effect of an inspiraling object on the DM halo surrounding an intermediate-mass or supermassive black hole? Is the halo disrupted, or destroyed? How common and robust is such a DM overdensity? And would its presence be detectable through a dephasing in the gravitational waveform?
• Exotic Compact Objects. What is their formation history and their distribution in the Galaxy? What is the impact on the number of predicted events and on the amplitude of the stochastic gravitational wave background they induce?
• Direct detection with GW experiments. What modifications of current and upcoming interferometers would maximize the sensitivity to DM models without jeopardizing gravitational wave searches? What would we learn about the nature of DM from the detection of subhalos with PTAs?
• Non-perturbative DM dynamics. Thorough analyses have been performed mostly in the context of the very early Universe, where the typical GW frequencies are much higher than any current experiment. Is it possible to get a measurable GW signal in the case of fragmentation of a coherent field or non-perturbative particle production in the relatively late Universe?
• Phase transitions. A broken power law in the stochastic background would imply new physics comprising a large fraction of the energy density at the time of the transition. Can we break the degeneracy among models predicting this feature, e.g. by searching for new states with masses at scales similar to that of the nucleation temperature of the phase transition?
In conclusion, we believe that the time is ripe for investigating the many connections between GWs and DM, and we encourage a strong community effort to further explore the interface between these two exciting fields of research.