Unique Multimessenger Signal of QCD Axion Dark Matter

We propose a multimessenger probe of QCD axion dark matter based on observations of black hole-neutron star binary inspirals. It is suggested that a dense dark matter spike may grow around intermediate mass black holes ( 10 3 – 10 5 M ⊙ ). The presence of such a spike produces two unique effects: a distinct phase shift in the gravitational wave strain during the inspiral and an enhancement of the radio emission due to the resonant axion-photon conversion occurring in the neutron star magnetosphere throughout the inspiral and merger. Remarkably, the observation of the gravitational wave signal can be used to infer the dark matter density and, consequently, to predict the radio emission. We study the projected reach of the LISA interferometer and next-generation radio telescopes such as the Square Kilometre Array. Given a sufficiently nearby system, such observations will potentially allow for the detection of QCD axion dark matter in the mass range 10 − 7 eV to 10 − 5 eV. DOI: 10.1103/PhysRevLett.124.161101

Introduction.-Theparticle nature of dark matter (DM) remains a mystery [1,2] despite efforts to detect it through astrophysical and lab-based observations [3][4][5][6].Another indication of new physics comes from the strong charge parity (CP) problem of quantum chromodynamics (QCD) [7].The nonobservation of the neutron electric dipole moment [8] constrains the CP-violating θ parameter in the QCD sector to be surprisingly small, jθj ≲ 10 −10 , while it could generically be Oð1Þ.A popular solution to this finetuning issue is the Peccei-Quinn mechanism, which predicts the existence of the axion [9][10][11][12].Axionlike particles are also predicted in several extensions of the standard model (SM), as well as in string theory [13].However, in the case of the QCD axion, there is a tight relation between its mass and its couplings with ordinary matter [14][15][16][17].
Gravitational waves (GWs) have provided a new observational portal into extreme astrophysical environments [51].The detection of the binary NS merger GW170817 and electromagnetic (EM) counterparts further revolutionized astrophysics and multimessenger astronomy [52,53].
Distortions to GWs from binary black holes, caused by finite size effects of superradiant clouds, have recently been shown to provide a new probe of beyond the standard model (BSM) physics [54,55].References [56,57] demonstrated that a DM minispike around an intermediate mass black hole (IMBH) can dramatically affect the GW waveform through dynamical friction, providing yet another direct probe of BSM physics.DM environmental effects on GW signals were studied more generally in Refs.[58][59][60].
In this Letter, we explore the possibility of probing QCD axion DM with multimessenger astronomy.We study the combined signal of GWs and radio emission from a NS inspiraling toward an IMBH surrounded by a dense spike of axion DM as sketched in Fig. 1.We show that by measuring the spike profile from the GW signal using the planned space-based observatory LISA [61] (Fig. 2), we can predict the gradual evolution of the radio signal during the yearslong inspiral phase.Most importantly, the increased density from the minispike amplifies the radio signal, allowing for the potential detection of QCD axions with photon couplings expected in the most commonly studied models [14][15][16][17].This is illustrated for radio observations with the Square Kilometre Array (SKA) [62,63] in Fig. 3.
These IMBHs may exist in DM halos [79][80][81][82].It has been shown that for a BH undergoing adiabatic growth [83] at the centre of such a halo, the surrounding DM would form a dense spike whose profile ρ DM ðrÞ is approximately a power law with index α [84][85][86][87][88][89] ρ DM ðrÞ ¼ The Navarro-Frenk-White (NFW) parameters ρ s and r s [90] are related to the cosmological history and mass of the halo, for which we follow Ref. [57], assuming a formation redshift z f ¼ 20 and total halo mass 10 6 M ⊙ .The radius of the BH's innermost stable circular orbit (ISCO) is denoted r ISCO .To solve for the spike parameters, we set where r sp ∼ 0.2r h [57].The spike profile varies with the initial DM profile.For an initially NFW-like profile, the spike slope is α ¼ 7=3 (our benchmark scenario).
These systems are speculative from both the perspective of IMBH formation as well as the presence and properties (such as α) of the spike [80,[91][92][93].For instance, for the spike to be preserved, the BH must not have undergone any mergers in its recent past [92], nor should it be in a dense baryonic environment [91].So if these systems form, their most likely location is in globular clusters [94].
In addition to the IMBH with a DM spike, we consider an inspiraling NS (on a circular orbit, for concreteness).NSs can have extremely high magnetic fields (10 9 -10 15 G), allowing for efficient axion-photon conversion close to the NS surface.NSs readily form in globular clusters [95] FIG. 2. Error on the DM density from GW measurements of α. Green bands show 1σ uncertainties on the reconstructed DM density from analyzing the GW waveform (for a system at d ¼ 0.01 Gpc, representing a signal-to-noise ratio for LISA of ∼92) over 10 bins in radius (measured from the position of the 10 4 M ⊙ IMBH).The fiducial density profile is shown as a blue dashed line.Along the top axis, we also label the approximate time-to-merger in the vacuum case.FIG. 3. Projected reach in axion-photon coupling from radio observations.SKA2 sensitivity (100 hours) to the axion-photon coupling for different orbital separations r and IMBH-DM-NS distances d, assuming θ ¼ 90°.The QCD axion parameter space is represented by the blue band, while the vertical and horizontal gray bands show the ADMX [30,31] and CAST [32] limits, respectively.and, therefore, are plausible candidates for mergers with IMBHs.We refer to the total system as IMBH-DM-NS.
Reference [80] argues that there are many IMBHs within our own Galactic halo.For an IMBH-DM-NS system to form, the IMBH must capture a NS.This process is very uncertain, relying on tracing formation models from the early Universe to today [94,96].Reference [94] suggests a detection rate density in LISA of approximately R ∼ 3-10 Gpc −3 yr −1 .Therefore, we consider two scenarios, one in which the IMBH-DM-NS system is close, at 0.01 Gpc, and one in which the system is further away, at 1 Gpc.The former is an optimistic scenario in terms of the strength of the radio signal, whereas many of the farther systems are likely to be observed over a ten year observing period.(Note that 1 Gpc corresponds to z ≈ 0.25 and a signal-to-noise ratio of ∼1 in LISA.)Importantly, these events would be dominated by IMBHs with masses 10 3 -10 4 M ⊙ .For concreteness, we consider an IMBH of 10 4 M ⊙ since the additional gravitational potential of the BH preserves the structure of the spike for longer times [97].
The properties of NSs in globular clusters are uncertain.However, they are thought to be much older than normal pulsars in galactic disks, and it has been suggested that most are formed from electron-capture supernova processes due to their minimal kick velocities [50,95].For the inspiraling NS, we take the magnetic field strength B 0 ¼ 10 12 G and spin period P ¼ 10 s [49,95]; NSs with similar properties have been found in observed globular clusters [98][99][100][101][102].We assume M NS ¼ 1.4 M ⊙ and r NS ¼ 10 km as benchmark values for the NS mass and radius.
Gravitational wave signal.-Invacuum, intermediate mass-ratio inspiral produces sub-Hertz gravitational waves.In the IMBH-DM-NS system, the dominant effect causing a deviation from the vacuum inspiral is the gravitational interaction between the DM halo and the NS passing through it, known as dynamical friction (DF) [103][104][105].Dynamical friction exerts a drag force on the NS where v NS is the velocity of the NS, and we take ln Λ ∼ 3 for the Coulomb logarithm.This force causes a loss of orbital energy, dE DF =dt ¼ v NS f DF , changing the accumulated phase of the GW signal and eventually reducing the inspiral time before merger with respect to the vacuum waveform.We see, from Eq. ( 2), that this force grows as the NS inspirals, (The NS orbital velocity grows roughly as r −1=2 , so that the DF force scales roughly as r −αþ1 .)although so, too, does the radiation reaction force due to GW emission.
In the Newtonian regime, the waveform of the IMBH-DM-NS system is computed by solving the energy balance equation, taking into account the effect of both DF and GW emission on the orbital energy E orbit of the system [57] − For circular orbits in the Newtonian regime, the energy loss due to GW emission is where ω s is the orbital frequency and r is the orbital radius.
The resulting phase difference with respect to the vacuum inspiral signal depends on the chirp mass M c ≃ M 3=5 NS M 2=5 BH , on the individual masses M BH and M NS , and on the density of DM. (Note that we only consider low redshift systems, z ≪ 1, therefore, we ignore any difference between lab and system frame.) Figure 2 shows the constraints on α as a function of radius from the IMBH recast as an error on the DM density.To calculate the error, we take ten log-spaced radial bins and integrate the noise-weighted inner product between the associated lower and upper frequencies, f i l and f i u , respectively.The error on α [106] (using the Fisher information) is then given by where S n ðfÞ is the LISA noise spectral density taken from Ref. [57] and h is the GW strain.We assume a five-year observation with LISA, beginning at a frequency of 0.04 Hz at r ≈ 1.5 × 10 −8 pc and ending at 0.44 Hz at the ISCO.We neglect any errors from the correlation between different parameters, which are expected to be small for α [57].We assume that all quantities (for example, spins and masses) can be measured precisely and do not contribute significantly to the error on the DM density.Note that higher order post-Newtonian effects on the inspiral will be important in breaking the degeneracy between M BH , M NS , and M c , as well as deducing the spins of the NS and IMBH.This degeneracy breaking has been demonstrated for current and future experiments in Refs.[107][108][109], although the impact on our projections should be tested in future work.
Figure 2 shows the 1σ uncertainty on the density reconstruction.At radii larger than r ≳ 6 × 10 −9 pc, the DM density can be constrained to better than 10%, but as the separation of the binary decreases, the uncertainty on the DM density increases.This is due to three effects; first, as the system approaches merger, GW emission (and not DF) begins to dominate the phase evolution of the waveform.Second, time spent at a given radius is not evenly distributed, as can be seen in the upper y axis of Fig. 2. Finally, the LISA sensitivity decreases at higher frequencies, weakening the constraining power at small r.
PHYSICAL REVIEW LETTERS 124, 161101 (2020) 161101-3 Therefore, the phase evolution of the waveform is very sensitive to dynamical friction (and thus, the DM density) predominantly when r is large.
While the formation, properties, and survival of DM spikes is not guaranteed [80,[91][92][93], GW observations can be used to confirm (or disfavor) the presence of a spike in a given system.DM profiles with α ¼ 1.5 will still induce a detectable phase shift [57], though these shallower slopes will significantly degrade constraints on α (see Supplemental Material I [110]).In any case, the resulting density constraints can be fed directly into the EM signal calculation, predicting the expected radio emission.
Radio signal.-Theradio signal arises from resonant axion-photon conversion occurring when the axion mass m a matches the frequency ω p of the plasma surrounding the NS with n c the number density of charged particles with mass m c .As a concrete example, we consider the Goldreich-Julian model for the NS plasma [111] which, in the nonrelativistic limit, provides where Ω ¼ 2π=P is the angular velocity and B is the magnetic field, which we consider to be in a dipole configuration with the axis aligned to the rotation axis.
Resonant axion-photon conversion then occurs at a specific radial distance from the NS center, which is given by [49] where θ is the polar angle with respect to the rotation axis.Equation ( 8) is obtained by setting ω p ¼ m a =2π and considering electrons or positions plasma (m c ¼ m e ).Following Ref. [49], the radiated power is given by where ρ DM ðr c Þ and v c are the DM density and velocity at the conversion radius.The energy transfer function p aγ is obtained using the WKB and stationary phase approximations to give where g aγγ is the strength of the coupling that leads to axion-photon conversion through the interaction L ¼ −g aγγ aE • B=4.
We use Eddington's formula to calculate the phase-space distribution of the DM in the BH frame [112,113], assuming isotropy and spherical symmetry (see Supplemental Material II [110] for more details).This distribution fðEÞ depends on the relative energy E ¼ ΨðrÞ − 1 2 v 2 and the gravitational potential ΨðrÞ ¼ Φ 0 − ΦðrÞ, relative to the potential at the minispike boundary, Φ 0 .For r ≲ 10 −8 pc (the point at which the GW signal would become observable) the enclosed mass is dominated by the BH, and therefore, we neglect the contribution of the minispike to the relative potential: Nearby DM particles are accelerated under gravity as they infall toward the NS.Particles with initial velocity v reach velocity ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi v 2 þ 2Ψ NS p at the conversion radius, where the NS potential is Ψ NS ¼ G N M NS =r c .Applying Liouville's theorem [114], we find the DM density at the conversion radius as We assume that the amplitude of the radiated power is dominated by the peak of the velocity distribution (We do not consider the boost to the NS frame since the NS orbital velocity is subdominant with respect to the DM peak velocity.) Finally, the flux density of the radio signal is given by where d is the distance to the system and B is the signal bandwidth (calculated as the 90% containment region of the DM velocity distribution far from the NS).Given the central frequency f of the radio signal, we find B=f to be 0.06 and 0.12 at an orbit of r ¼ 6 × 10 −9 pc and r ¼ 3 × 10 −9 pc, respectively.Figure 3 shows the projected reach of the future SKA2 telescope in the axion parameter space, obtained by considering the minimum detectable flux density S which provides a signal-to-noise ratio (SNR) equal to one.(An SNR of one was chosen to directly compare with Ref. [49].Increasing the required SNR for detection leads to a corresponding increase in the couplings we can probe as g aγγ ∝ ffiffiffiffiffiffiffiffiffi ffi SNR p .)In particular, for a radio telescope where n pol ¼ 2 is the number of polarizations, SEFD ¼ 0.098 is the SKA2 system-equivalent flux density as estimated in Ref. [50], and we assume an observation time of Δt obs ¼ 100 hours.This is roughly the time spent by the system from the closest orbit we consider (r ¼ 3 × 10 −9 pc) until merger.We have assumed that the sensitivity is limited by the detector's thermal noise.Since IMBHs have not yet been conclusively detected, it is thought that they do not have an appreciable accretion disk, meaning that there is unlikely to be any background radio emission during the inspiral.The sensitivity curves are valid in the mass range 10 −7 eV ≤ m a ≤ 1.4 × 10 −5 eV.The lower limit is set by the lower cutoff frequency potentially probed by SKA, (We note that, for m a > 10 −7 eV, the corresponding axion Compton wavelength is smaller than 0.01 km, 3 orders of magnitude smaller than the size of the NS.This allows us to neglect quantum effects when computing both signatures.) while the upper limit comes from the requirement that axion-photon conversion occurs outside the NS, r c ≥ r NS .
Figure 3 shows that a crucial parameter is the distance of the IMBH-DM-NS system, since the flux density depends on its inverse square.On the other hand, the sensitivity does not strongly depend on the BH-NS separation r.Radio observations taken when r ∼ 3 × 10 −9 pc (solid lines) yield sensitivities to g aγγ which are roughly a factor of 2 stronger than for r ∼ 6 × 10 −9 pc (dashed lines).In Fig. 3, we have fixed ρ DM to the fiducial density profile.However, as we saw in Fig. 2, the DM density is likely to be more poorly constrained at smaller radii, making the radio sensitivity at large r substantially more robust (though not substantially weaker).
Discussion.-With a sufficiently nearby detection of an IMBH-DM-NS system, it will be possible to probe the parameter space of QCD axion DM.We find roughly a 0.05% probability of a detection closer than d ¼ 0.01 Gpc over 10 years, using predicted LISA detection rates for such systems [94] (though these typically come with large uncertainties).Instead, out to d ¼ 1 Gpc, we expect a few tens of detections per year.
We emphasize that setting upper limits on g aγγ from the nondetection of a radio signal is hampered by uncertainties in the individual NS properties and magnetosphere modeling (see Supplemental Material III [110]).Detecting and studying a larger population of such systems would perhaps allow for robust limit setting, through modeling of the expected properties of the NS population.Nevertheless, a joint GW þ EM detection is within reach of upcoming experiments and would be a striking confirmation of axion Dark Matter.GW observations can provide constraints on the DM density around BHs, as in Fig. 2, with the better estimation of the density at larger separations reducing uncertainties on the expected radio signal.
Above around m a ∼ 10 −6 eV, these broadband sensitivities would be complementary to current and proposed axion haloscope experiments [30,31,[34][35][36][37] (some of which are plotted in Fig. 3).These are sensitive to the density of DM local to Earth, which carries its own uncertainties [115].Such uncertainties could be mitigated in our scenario by combining information from GW and radio emission.Therefore, multimessenger observations of black hole-dark matter-neutron star systems have the potential to detect QCD axion dark matter for masses between 10 −7 eV and 10 −5 eV.
We thank David Nichols, Tanja Hinderer, Mikael Leroy, and Gianfranco Bertone for fruitful discussions.Finally, we thank the python scientific computing packages numpy [116] and scipy [117].

FIG. 1 .
FIG. 1. Illustration of the IMBH-DM-NS system.A DM halo of axions a around the intermediate mass black hole (IMBH) produces a phase shift in the GW signal and radio emission due to its conversion into photons γ in the neutron star (NS) magnetosphere.
This research is funded by NWO through the VIDI research program "Probing the Genesis of Dark Matter" [Grant No. 680-47-532; (T.E., C. W.)]. S. N. is grateful for support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) through the VIDI and Projectruimte grants (PI Nissanke)