A Unique Multi-Messenger Signal of QCD Axion Dark Matter

TL;DR

This paper proposes a multi-messenger approach combining gravitational wave and radio observations to detect QCD axion dark matter around black holes, offering a new method to identify dark matter properties.

Contribution

It introduces a novel multi-messenger method using gravitational waves and radio signals from black hole-neutron star systems to detect QCD axion dark matter.

Findings

Potential detection of QCD axion dark matter in the mass range $10^{-7}$ to $10^{-5}$ eV.

Projected sensitivity of LISA and SKA for observing these signals.

Distinct phase shifts and radio emission enhancements as signatures of dark matter presence.

Abstract

We propose a multi-messenger 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} \mathrm{\,M_{\odot}}$). 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}\,\mathrm{eV}$ to $10^{-5}\,\mathrm{eV}$.

Equations

$${\rho }_{\text{DM}}(r)=\{\begin{array}{ll}{\rho }_{\text{sp}}{\left(\frac{r_{\text{sp}}}{r}\right)}^{\alpha },& r_{\text{ISCO}}<r\le r_{\text{sp}}\\ \frac{{\rho }_{\text{s}}}{(r/r_{\text{s}})(1+r/r_{\text{s}}{)}^2},& r>r_{\text{sp}}.\end{array}$$

$$f_{\text{DF}}=4\pi G_N^2{M}_{\mathrm{NS}}^2\frac{{\rho }_{\text{DM}}(r)}{v_{\text{NS}}^2(r)}\ln \Lambda ,$$

$$-\frac{\mathrm{d}{E}_{\text{orbit }}}{\mathrm{d}t}=\frac{\mathrm{d}{E}_{\mathrm{DF}}}{\mathrm{d}t}+\frac{\mathrm{d}{E}_{\mathrm{GW}}}{\mathrm{d}t}.$$

$$\frac{\mathrm{d}{E}_{\mathrm{GW}}}{\mathrm{d}t}=\frac{32}{5}\frac{G_N{M}_{\mathrm{NS}}^2}{c^5}{r}^4{\omega }_s^6,$$

$$\frac{\Delta \alpha }{\alpha }={\left[4\mathrm{Re}\left({\int }_{f_l^i}^{f_u^i}\frac{\frac{\partial h}{\partial \ln \alpha }\frac{\partial h^{\ast }}{\partial \ln \alpha }}{S_n(f)}\mathrm{d}f\right)\right]}^{-\frac{1}{2}},$$

$${\omega }_p\approx \sqrt{\frac{4\pi {\alpha }_{\mathrm{EM}}{n}_c}{m_c}},$$

$$n_c=\frac{2\mathbf{\Omega }\cdot \mathbf{B}}{e},$$

$$r_c$$

$$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}\Omega }\sim 2\times p_{a\gamma }{\rho }_{\text{DM}}(r_c)v_c{r}_c^2,$$

$$p_{a\gamma }\sim \frac{\pi }{12}\frac{g_{a\gamma \gamma }^2{B}_0^2{r}_{\mathrm{NS}}}{m_a}{\left(\frac{r_{\mathrm{NS}}}{r_c}\right)}^5\left(3{\cos }^2\theta +1\right),$$

$$\begin{array}{rl}{\rho }_{\text{DM}}(r_c)& =\sqrt{\frac{2}{\pi }}\frac{{\rho }_{\mathrm{sp}}{r}_{\mathrm{sp}}{}^{\alpha }}{(G_N{M}_{\mathrm{BH}}{)}^{\alpha }}\frac{\alpha (\alpha -1)\Gamma (\alpha -1)}{\Gamma (\alpha -\frac{1}{2})}\times \\ & {\int }_{v_{\min }}^{v_{\max }}{\left[{\Psi }_{\mathrm{BH}}+{\Psi }_{\mathrm{NS}}-\frac{v^2}{2}\right]}^{\alpha -\frac{3}{2}}{v}^2\mathrm{d}v,\end{array}$$

$$v_c^2\sim \frac{2G_N{M}_{\mathrm{BH}}}{r}{\left[\alpha -\frac{1}{2}\right]}^{-1}+\frac{2G_N{M}_{\mathrm{NS}}}{r_c}.$$

$$S=\frac{1}{\mathcal{B}{d}^2}\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}\Omega },$$

$$\mathrm{SNR}=\frac{S}{\mathrm{SEFD}}\sqrt{n_{\mathrm{pol}}\mathcal{B}\Delta t_{\mathrm{obs}}},$$

$$\Delta \psi =\tilde{\varphi }-\varphi ,$$

$$f(\mathcal{E})=\frac{1}{\sqrt{8}{\pi }^2}{\int }_0^{\mathcal{E}}\frac{\mathrm{d}\Psi }{\sqrt{\mathcal{E}-\Psi }}\frac{{\mathrm{d}}^2\rho }{\mathrm{d}{\Psi }^2}.$$

$$\rho (\Psi )$$

$$\begin{array}{rl}f(\mathcal{E})& =\frac{\alpha (\alpha -1)}{\sqrt{8}{\pi }^2}{\rho }_{\mathrm{sp}}{\left(\frac{r_{\mathrm{sp}}}{G_N{M}_{\mathrm{BH}}}\right)}^{\alpha }{\int }_0^{\mathcal{E}}{\Psi }^{\alpha -2}\frac{\mathrm{d}\Psi }{\sqrt{\mathcal{E}-\Psi }}\\ & =\frac{\alpha (\alpha -1)}{(2\pi {)}^{3/2}}{\rho }_{\mathrm{sp}}{\left(\frac{r_{\mathrm{sp}}}{G_N{M}_{\mathrm{BH}}}\right)}^{\alpha }\frac{\Gamma (\alpha -1)}{\Gamma (\alpha -\frac{1}{2})}{\mathcal{E}}^{\alpha -3/2}.\end{array}$$

$$f(v|r)$$

$${\int }_0^{v_{\max }(r)}f(v|r)\mathrm{d}v=1.$$

$$v_{\mathrm{peak}}^2=v_{\max }^2{\left[\alpha -\frac{1}{2}\right]}^{-1}=\frac{2G_N{M}_{\mathrm{BH}}}{r}{\left[\alpha -\frac{1}{2}\right]}^{-1}.$$

$$\frac{\mathrm{d}\mathcal{P}}{\mathrm{d}\Omega }\propto B_{0}P\left(\frac{3{\cos }^2\theta +1}{\left|3{\cos }^2\theta -1\right|}\right)\left[g_{a\gamma \gamma }^2{m}_a{\rho }_{\mathrm{DM}}(r_c)v_c\right].$$

Full text

Introduction — The particle nature of Dark Matter (DM) remains a mystery Bertone et al. (2005a); Bertone and Tait (2018) despite efforts to detect it through astrophysical and lab-based observations Duffy and van Bibber (2009); Marrodán Undagoitia and Rauch (2016); Gaskins (2016); Kahlhoefer (2017). Another indication of New Physics comes from the Strong charge parity (CP) problem of quantum chromodynamics (QCD) Dine (2000). The non-observation of the neutron electric dipole moment Pendlebury et al. (2015) constrains the CP-violating $\theta$-parameter in the QCD sector to be surprisingly small, $\left|\theta\right|\lesssim 10^{-10}$, while it could generically be $\mathcal{O}(1)$. A popular solution to this fine-tuning issue is the Peccei-Quinn mechanism, which predicts the existence of the axion Peccei and Quinn (1977a, b); Weinberg (1978); Wilczek (1978). Axion-like particles are also predicted in several extensions of the Standard Model, as well as in string theory Arvanitaki et al. (2010). However, in the case of the QCD axion there a tight relation between its mass and its couplings with ordinary matter Kim (1979); Shifman et al. (1980); Zhitnitsky (1980); Dine et al. (1981).

These two fundamental issues can be addressed simultaneously by treating the QCD axion as DM Preskill et al. (1983); Abbott and Sikivie (1983); Dine and Fischler (1983) (see Ref. Marsh (2016) for a review). Axions may be produced with the correct relic abundance through the misalignment mechanism Wantz and Shellard (2010), through the decay of topological defects such as strings and domain walls Davis (1986); Harari and Sikivie (1987); Hagmann et al. (2001); Hiramatsu et al. (2012); Kawasaki et al. (2015), or via thermal production Turner (1987); Salvio et al. (2014). So far, only a small part of the QCD axion parameter space has been explored  Asztalos et al. (2010); Du et al. (2018); Anastassopoulos et al. (2017); Hoof et al. (2019), though new experimental techniques have been recently proposed Brubaker et al. (2017); Caldwell et al. (2017); Zhong et al. (2018); Brun et al. (2019); Lawson et al. (2019); Kahn et al. (2016); McAllister et al. (2017); Shokair et al. (2014); Al Kenany et al. (2017); Alesini et al. (2017); Caputo et al. (2018, 2019) (see Ref. Irastorza and Redondo (2018) for a comprehensive review). Furthermore, it has been noted that the Primakov effect can efficiently convert axions to photons in the magnetic fields of Neutron Stars (NSs) Pshirkov and Popov (2009). These photons are potentially observable with current and future radio telescopes  Huang et al. (2018); Hook et al. (2018); Safdi et al. (2018).

Gravitational Waves (GWs) has provided a new observational portal into extreme astrophysical environments Abbott et al. (2016a). The detection of the binary NS merger GW170817 and electromagnetic counterparts further revolutionised astrophysics and multi-messenger astronomy Abbott et al. (2017a, b). 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  Baumann et al. (2019); Hannuksela et al. (2019a). References Eda et al. (2013, 2015) demonstrated that a DM mini-spike 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. Macedo et al. (2013); Barausse et al. (2014); Yue and Han (2018).

In this Letter, we explore the possibility of probing QCD axion DM with multi-messenger astronomy. We study the combined signal of GWs and radio emission from a NS inspiraling towards 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 Amaro-Seoane et al. (2017) (Fig. 2), we can predict the gradual evolution of the radio signal during the years-long inspiral phase. Most importantly, the increased density from the mini-spike amplifies the radio signal, allowing for the potential detection of QCD axions with photon couplings expected in the most commonly studied models Kim (1979); Shifman et al. (1980); Zhitnitsky (1980); Dine et al. (1981). This is illustrated for radio observations with the Square Kilometre Array (SKA) Bull et al. (2018); Dewdney in Fig. 3.

Astrophysical System — IMBHs have masses $M_{\text{BH}}=10^{3}-10^{5}M_{\odot}$. Thought to reside in the centres of smaller spiral galaxies ($\lesssim 10^{9}\,M_{\odot}$ Greene (2012)), as well as in dense stellar environments such as globular clusters Miller and Hamilton (2002), a growing number of observations point toward the existence of IMBHs in nature Webb et al. (2012); Ballone et al. (2018); Takekawa et al. (2019); Woo et al. (2019). There are multiple proposed formation mechanisms, including runaway growth through the mergers of stellar mass objects Taniguchi et al. (2000); Portegies Zwart and McMillan (2002); Portegies Zwart et al. (2004); the direct collapse of gas clouds at high redshift Begelman et al. (2006); Agarwal et al. (2012); or primordial formation from large density perturbations Carr (1975); Kawaguchi et al. (2008); Carr and Kuhnel (2019); Carr et al. (2019).

These IMBHs may exist in DM halos Islam et al. (2003); Zhao and Silk (2005); Bertone et al. (2005b); Rashkov and Madau (2014). It has been shown that for a BH undergoing adiabatic growth Sigurdsson (2003) at the centre of such a halo, the surrounding DM would form a dense spike whose profile $\rho_{\text{DM}}(r)$ is approximately a power law with index $\alpha$ Blumenthal et al. (1986); Quinlan et al. (1995); Gondolo and Silk (1999); Bertone et al. (2002, 2009); Sadeghian et al. (2013):

[TABLE]

The NFW parameters $\rho_{\text{s}}$ and $r_{\text{s}}$ Navarro et al. (1996) are related to the cosmological history and mass of the halo, for which we follow Ref. Eda et al. (2015), assuming a formation redshift $z_{f}=20$ and total halo mass $10^{6}\,M_{\odot}$. The radius of the BH’s inner-most stable circular orbit (ISCO) is denoted $r_{\mathrm{ISCO}}$. To solve for the spike parameters we set $M(<r_{\mathrm{h}})=4\pi\int_{r_{\mathrm{ISCO}}}^{r_{\mathrm{h}}}\rho_{\mathrm{DM}}r^{2}\,\mathrm{d}r=2M_{\mathrm{BH}}$, where $r_{\mathrm{sp}}\sim 0.2\,r_{\mathrm{h}}$ Eda et al. (2015). The spike profile varies with the initial DM profile. For an initially NFW-like profile, the spike slope is $\alpha=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 $\alpha$) of the spike Ullio et al. (2001); Merritt et al. (2002); Zhao and Silk (2005); Hannuksela et al. (2019b). For instance, for the spike to be preserved, the BH must not have undergone any mergers in its recent past Merritt et al. (2002), nor should it be in a dense baryonic environment Ullio et al. (2001). So if these systems form, their most likely location is in globular clusters Fragione et al. (2018).

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}\mathrm{\,G}$), allowing for efficient axion-photon conversion close to the NS surface. NSs readily form in globular clusters Freire (2013) and are therefore plausible candidates for mergers with IMBHs. We refer to the total system as IMBH-DM-NS.

Reference Zhao and Silk (2005) 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 Mandel et al. (2008); Fragione et al. (2018). Reference Fragione et al. (2018) suggests a detection rate density in LISA of approximately $\mathcal{R}\sim 3-10\,\mathrm{Gpc^{-3}\,yr^{-1}}$. We therefore 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.111Note that 1 Gpc corresponds to $z\approx 0.25$ and a signal-to-noise ratio of $\sim 1$ in LISA. Importantly, these events would be dominated by IMBHs with masses $10^{3}-10^{4}\,\mathrm{M}_{\odot}$. For concreteness we consider an IMBH of $10^{4}\,\mathrm{M}_{\odot}$ since the additional gravitational potential of the BH preserves the structure of the spike for longer times Kavanagh et al. (2020).

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 Freire (2013); Safdi et al. (2018). For the inspiraling NS, we take the magnetic field strength $B_{0}=10^{12}\mathrm{\,G}$ and spin period $P=10\mathrm{\,s}$ Freire (2013); Hook et al. (2018); NSs with similar properties have been found in observed globular clusters Lyne et al. (1996); Biggs et al. (1994); Lyne et al. (1993); Ivanova et al. (2008); Boyles et al. (2011). We assume $M_{\rm NS}=1.4\,\mathrm{M}_{\odot}$ and $r_{\rm NS}=10$ km as benchmark values for the NS mass and radius.

Gravitational Wave signal — In vacuum, 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) Chandrasekhar (1943a, b, c). Dynamical friction exerts a drag force on the NS:

[TABLE]

where $v_{\text{NS}}$ is the velocity of the NS, and we take $\ln\Lambda\sim 3$ for the Coulomb logarithm. This force causes a loss of orbital energy, $\mathrm{d}E_{\mathrm{DF}}/\mathrm{d}t=v_{\mathrm{NS}}f_{\mathrm{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,222The NS orbital velocity grows roughly as $r^{-1/2}$, so that the DF force scales roughly as $r^{-\alpha+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_{\text{orbit }}$ of the system Eda et al. (2015):

[TABLE]

For circular orbits in the Newtonian regime, the energy loss due to GW emission is

[TABLE]

where $\omega_{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 $\mathcal{M}_{c}\simeq M_{\mathrm{NS}}^{3/5}\,M_{\mathrm{BH}}^{2/5}$, on the individual masses $M_{\mathrm{BH}}$ and $M_{\mathrm{NS}}$, and on the density of DM.333Note that we only consider low redshift systems, $z\ll 1$, therefore we ignore any difference between lab and system frame.

Figure 2 shows the constraints on $\alpha$ 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 $\alpha$ Cutler and Flanagan (1994) (using the Fisher information) is then given by:

[TABLE]

where $S_{n}(f)$ is the LISA noise spectral density taken from Ref. Eda et al. (2015) and $h$ is the GW strain. We assume a 5-year observation with LISA, beginning at a frequency of 0.04 Hz at $r\approx 1.5\times 10^{-8}\,\mathrm{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 $\alpha$ Eda et al. (2015). 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_{\mathrm{BH}}$, $M_{\mathrm{NS}}$ and $\mathcal{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. Abbott et al. (2016b); Barack and Cutler (2004); Abbott et al. (2016c), although the impact on our projections should be tested in future work.

Figure 2 shows the $1\sigma$ uncertainty on the density reconstruction. At radii larger than $r\gtrsim 6\times 10^{-9}\,\mathrm{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; firstly, as the system approaches merger, GW emission (and not DF) begins to dominate the phase evolution of the waveform. Secondly, 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$. The phase evolution of the waveform is therefore 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 Ullio et al. (2001); Merritt et al. (2002); Zhao and Silk (2005); Hannuksela et al. (2019b), GW observations can be used to confirm (or disfavor) the presence of a spike in a given system. DM profiles with $\alpha=1.5$ will still induce a detectable phase shift Eda et al. (2015), though these shallower slopes will significantly degrade constraints on $\alpha$ (see Supplementary Material I). In any case, the resulting density constraints can be fed directly into the EM signal calculation, predicting the expected radio emission.

Radio Signal — The radio signal arises from resonant axion-photon conversion occurring when the axion mass $m_{a}$ matches the frequency $\omega_{p}$ of the plasma surrounding the NS:

[TABLE]

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 Goldreich and Julian (1969), which in the non-relativistic limit provides

[TABLE]

where $\Omega=2\pi/P$ is the angular velocity and $\mathbf{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 Hook et al. (2018)

[TABLE]

where $\theta$ is the polar angle with respect to the rotation axis. Equation (8) is obtained by setting $\omega_{p}=m_{a}/2\pi$ and considering electrons/positions plasma ($m_{c}=m_{e}$).

Following Ref. Hook et al. (2018), the radiated power is given by

[TABLE]

where $\rho_{\rm DM}(r_{c})$ and $v_{c}$ are the DM density and velocity at the conversion radius. The energy transfer function $p_{a\gamma}$ is obtained using the WKB and stationary phase approximations to give

[TABLE]

where $g_{a\gamma\gamma}$ is the strength of the coupling that leads to axion-photon conversion through the interaction $\mathcal{L}=-g_{a\gamma\gamma}\,a\,\mathbf{E}\cdot\mathbf{B}/4$.

We use Eddington’s formula to calculate the phase-space distribution of the DM in the BH frame Binney and Tremaine (2008); Catena and Ullio (2012), assuming isotropy and spherical symmetry (see Supplementary Material II for more details). This distribution $f(\mathcal{E})$ depends on the relative energy $\mathcal{E}=\Psi(r)-\frac{1}{2}v^{2}$ and the gravitational potential $\Psi(r)=\Phi_{0}-\Phi(r)$, relative to the potential at the mini-spike boundary, $\Phi_{0}$. For $r\lesssim 10^{-8}$ pc (the point at which the GW signal would become observable) the enclosed mass is dominated by the BH and we therefore neglect the contribution of the mini-spike to the relative potential: $\Psi=\Psi_{\mathrm{BH}}=G_{N}\,M_{\mathrm{BH}}/r$. In this case, we find $f(\mathcal{E})\propto\mathcal{E}^{\alpha-3/2}$ (for $\mathcal{E}>0$).

Nearby DM particles are accelerated under gravity as they infall toward the NS. Particles with initial velocity $v$ reach velocity $\sqrt{v^{2}+2\Psi_{\mathrm{NS}}}$ at the conversion radius, where the NS potential is $\Psi_{\mathrm{NS}}=G_{N}M_{\mathrm{NS}}/r_{c}$. Applying Liouville’s theorem Liouville (1838), we find the DM density at the conversion radius as,

[TABLE]

where $v_{\mathrm{min}}=\sqrt{2\Psi_{\mathrm{NS}}}$ and $v_{\mathrm{max}}=\sqrt{2\left(\Psi_{\mathrm{BH}}+\Psi_{\mathrm{NS}}\right)}$. We assume that the amplitude of the radiated power is dominated by the peak of the velocity distribution:444We do not consider the boost to the NS frame since the NS orbital velocity is subdominant with respect to the DM peak velocity.

[TABLE]

Finally, the flux density of the radio signal is given by

[TABLE]

where $d$ is the distance to the system and $\mathcal{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 $\mathcal{B}/f$ to be 0.06 and 0.12 at an orbit of $r=6\times 10^{-9}~{}\mathrm{pc}$ and $r=3\times 10^{-9}~{}\mathrm{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.555An SNR of one was chosen to directly compare with Ref. Hook et al. (2018). Increasing the required SNR for detection leads to a corresponding increase in the couplings we can probe as $g_{a\gamma\gamma}\propto\sqrt{\mathrm{SNR}}$. In particular, for a radio telescope,

[TABLE]

where $n_{\rm pol}=2$ is the number of polarizations, $\rm SEFD=0.098$ is the SKA2 system-equivalent flux density as estimated in Ref. Safdi et al. (2018) and we assume an observation time of $\Delta t_{\rm obs}=100$ hours. This is roughly the time spent by the system from the closest orbit we consider ($r=3\times 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}\mathrm{\,eV}\leq m_{a}\leq 1.4\times 10^{-5}\mathrm{\,eV}$. The lower limit is set by the lower cut-off frequency potentially probed by SKA,666We note that for $m_{a}>10^{-7}\mathrm{\,eV}$ the corresponding axion Compton wavelength is smaller than $0.01\mathrm{\,km}$, three 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}\geq r_{\rm 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\sim 3\times 10^{-9}\,\mathrm{pc}$ (solid lines) yield sensitivities to $g_{a\gamma\gamma}$ which are roughly a factor of 2 stronger than for $r\sim 6\times 10^{-9}\,\mathrm{pc}$ (dashed lines). In Fig. 3, we have fixed $\rho_{\mathrm{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\,\mathrm{Gpc}$ over 10 years, using predicted LISA detection rates for such systems Fragione et al. (2018) (though these typically come with large uncertainties). Instead, out to $d=1\,\mathrm{Gpc}$, we expect a few tens of detections per year.

We emphasise that setting upper limits on $g_{a\gamma\gamma}$ from the non-detection of a radio signal is hampered by uncertainties in the individual NS properties and magnetosphere modelling (see Supplementary Material III). Detecting and studying a larger population of such systems would perhaps allow for robust limit-setting, through modelling 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}\sim 10^{-6}\,\mathrm{eV}$, these broadband sensitivities would be complementary to current and proposed axion haloscope experiments Asztalos et al. (2010); Du et al. (2018); Brubaker et al. (2017); Caldwell et al. (2017); Zhong et al. (2018); Brun et al. (2019) (some of which are plotted in Fig. 3). These are sensitive to the density of DM local to Earth, which carries its own uncertainties Read (2014). Such uncertainties could be mitigated in our scenario by combining information from GW and radio emission. Multi-messenger observations of Black Hole - Dark Matter - Neutron Star systems therefore have the potential to detect QCD axion Dark Matter for masses between $10^{-7}\,\mathrm{eV}$ and $10^{-5}\,\mathrm{eV}$.

I Gravitational Waves and Spike Dependence

The phase difference between a vacuum inspiral and the one considered here is given by,

[TABLE]

where $\phi=-\frac{3}{4}\left(8\pi G_{N}\mathcal{M}_{c}f/c^{3}\right)^{-5/3}$ is the Newtonian vacuum phase evolution and $\tilde{\phi}$ is the phase evolution including dynamical friction, as given by Eq. (28) of Ref. Eda et al. (2015). The phase evolution of the gravitational wave signal provides a fundamental insight into the dynamics of the binary system. As shown in Fig. S1, the presence of a DM spike (with $\alpha>2.0$) produces a considerable phase shift when compared to the evolution of a vacuum inspiral. Again for $\alpha>2.0$, the specific phase evolution of any particular system can therefore be used to constrain $\alpha$ to high precision. For our baseline scenario of $\alpha=7/3$, constraints on $\alpha$ correspond to an error on the DM density of $\mathcal{O}(0.01\%)$ at $r\approx 1.3\times 10^{-8}\,\mathrm{pc}$. For $\alpha<2.0$, the phase difference becomes increasingly difficult to probe. Assuming the masses of the two objects can be independently measured to high precision, the constraint on $\alpha$ provides a direct constraint on the DM density local to the position of the NS Kavanagh et al. (2020). We do not account for errors on the overall normalisation of the DM density profile directly. The normalisation can also be measured, though it is degenerate with $M_{\mathrm{NS}}$ and $M_{\mathrm{BH}}$. To resolve the individual masses, higher order effects close to merger need to be accounted for. The errors associated with the individual mass determinations may dominate the error on the DM density normalisation at larger radii, but this is beyond the scope of the paper. We will address this in future work.

Figure S2 shows the constraint on the DM density from the GW signal, as described in the main text. The left and right panels show the constraint for $\alpha=2.0$ and $\alpha=2.5$, respectively. The error bars become larger for lower values of $\alpha$; this can easily be understood from Fig. S1. As the inspiral progresses, the GW frequency becomes larger (equivalently, the radius decreases). Similarly, the phase difference becomes ever smaller, gradually approaching the vacuum inspiral phase evolution and therefore providing no probe of the DM density. As $\alpha$ is increased, the large phase differences persist further into the inspiral, allowing one to probe the DM density closer to the IMBH.

Finally, in Fig. S3 we present the radio sensitivity for $\alpha=\{2.0,\,7/3,\,2.5\}$. As expected, the varying density, as shown in Fig. S2, amplifies or decreases the density close the IMBH. For $\alpha=2.0$, it is still possible to probe a small range of the QCD axion parameter space, although the constraint on the DM density becomes significantly worse (see left panel of Fig. S2). When $\alpha=2.5$, the density can be constrained extremely well down to small radii. The density is also increased by an order of magnitude compared to the $\alpha=7/3$ scenario, subsequently increasing the sensitivity by a similar amount.

II Dark Matter Velocity Distribution

We assume that the distribution of DM around the central black hole is spherically symmetric and that the velocity distribution of DM particles is isotropic. In this case, we can calculate the DM distribution function using Eddington’s Inversion Formula:

[TABLE]

Here, $\rho(r)$ is the density profile of the DM particles, while $\Psi(r)$ is the total gravitational potential, which in general includes a contribution from both the central mass and the mass enclosed in the DM halo. However, at small radii, the enclosed DM mass is small and we can typically neglect the contribution of the DM halo itself to the gravitational potential. We thus write $\Psi(r)=G_{N}\,M_{\mathrm{BH}}/r$ and re-express the density in terms of the potential:

[TABLE]

We therefore find:

[TABLE]

The DM speed distribution at a radius $r$ is then given by

[TABLE]

Here, we have defined $v_{\mathrm{max}}=v_{\mathrm{max}}(r)=\sqrt{2\Psi(r)}$, and we set the speed distribution to zero for $v>v_{\mathrm{max}}(r)$. With this definition, the speed distribution is normalised to one at any given radius:

[TABLE]

The peak velocity (far from the NS) is obtained by setting $\partial f(v|r)/\partial v=0$, giving:

[TABLE]

In Fig. S4, we show the DM speed distribution at several radii $r$, assuming $\alpha=7/3$ and $M_{\mathrm{IMBH}}=10^{4}\,M_{\odot}$. Solid lines show the speed distribution derived from a full numerical calculation of $f(\mathcal{E})$ (using Eq. (S4) and including self-consistently the potential due to the DM halo). Dashed lines show the approximate speed distribution given in Eq. (S5) (neglecting the potential of the DM halo itself). We see that in all cases of interest to us, $r\lesssim 10^{-8}\,\mathrm{pc}$, Eq. (S5) provides an excellent approximation to the full expression.

We note that as $r\rightarrow r_{\mathrm{ISCO}}$, the maximum DM speed tends towards the speed of light. As we discuss in the main text, the dominant effect from dynamical friction typically occurs at larger radii, where the DM speeds are lower and the non-relativistic formalism should apply. However, at a radius $r=3\times 10^{-9}\,\mathrm{pc}$, the maximum DM speed is $v_{\mathrm{max}}\approx 1.7\times 10^{5}\,\mathrm{km/s}\approx 0.56\,c$. Adding also the infall velocity toward the conversion radius, the DM particles are accelerated up to $\sim 0.8\,c$. This suggests that towards the end of the inspiral, our non-relativistic formalism would over-estimate the speeds which can be reached by the DM particles. However, even at $r=3\times 10^{-9}\,\mathrm{pc}$, the DM speeds are still only mildly relativistic ($\gamma\sim 1.7$), suggesting that this should be a small effect. We leave a more detailed analysis – including relativistic effects, boosting in the NS rest-frame and anisotropy of the infalling DM flux – to future work.

Finally, we have checked that the DM halo should survive the merger itself; the work done by dynamical friction during the five year inspiral is only a few percent of the total gravitational binding energy of the halo. A more detailed study of feedback on the DM halo in different systems is in preparation Kavanagh et al. (2020).

III Neutron Star Parameters

In the Goldreich-Julian model for the NS magnetosphere Goldreich and Julian (1969), the amplitude of the radio signal (produced by the resonant axion-photon conversion) is completely determined by the magnetic field strength at the NS poles $B_{0}$ and the spin period $P$. In particular, they determine the number density of charged particles around the NS (Eq. (7)), and consequently the plasma frequency (6) and the conversion radius (8). By plugging Eq. (10) into Eq. (9), we find the scaling relatation of the radiated power with respect to $B_{0}$ and $P$ to be

[TABLE]

where the quantities in the squared parentheses are almost independent of the NS parameters. Hence, for a given axion mass, the larger the NS magnetic field and spin period, the larger the radiated power. Furthermore, the radiated power can be significantly larger for $\cos\theta=1/\sqrt{3}$. We note that by neglecting the second term in the expression (12) for the velocity $v_{c}$, one can show that the minimum detectable axion-photon coupling from Eq. (14) scales with the axion mass as $g_{a\gamma\gamma}^{\rm min}\sim m_{a}^{-1/2}$.

In Fig. S5 we report the projected sensitivity curves of SKA2 for three different values for the NS magnetic field strength (left panel) and three different values for the NS spin period (right panel), while fixing $r=3\times 10^{-9}$ pc, $d=0.01$ Gpc, $\theta=90^{\circ}$ and $\alpha=7/3$. As expected, the larger the magnetic field and the spin period, the smaller the axion-photon coupling that can be probed by SKA2. Moreover, larger axion masses can be explored for larger magnetic fields or smaller spin periods. Increasing $B_{0}$, or decreasing $P$, causes the axion-photon conversion to occur at a larger radius (see Eq. (8)). The requirement that the the conversion radius is larger than the size of the NS ($r_{c}\geq r_{\rm NS}=10~{}{\rm km}$) is then satisfied for a wider range of axion masses.

We note that the values for $B_{0}$ and $P$ considered here cover almost all the possible properties of active NSs, according to the ATNF pulsar catalog Manchester et al. (2005). On the other hand, old dead NSs are expected to have low magnetic fields and large spin periods, providing weaker sensitivities. However, their properties are quite uncertain and model-dependent Safdi et al. (2018); Faucher-Giguere and Kaspi (2006); Popov et al. (2010). In all cases, we have also verified that the plasma remains bound to the NS. Even down to the innermost stable circular orbit, the forces from the NS magnetosphere dominate over the gravitational force from the BH by many orders of magnitude.

Finally, we conclude this appendix by discussing how the radio signal depends on the model for the NS magnetosphere. In deriving the expression for the conversion radius (8), we have assumed that the NS plasma is dominated by electron and positrons. The presence of ions would reduce the conversion radius once all the other quantities are fixed, since the plasma frequency (6) would be suppressed by the larger mass $m_{c}$. This would increase the radio signal according to Eq. (9) at the expense of probing smaller axion masses. Moreover, different analytic models and numerical simulations of the NS magnetosphere lead to different profiles for the plasma frequency Pétri (2016). This might increase or reduce the amplitude of the radio signal. For example, Ref. Safdi et al. (2018) pointed out that the electrosphere model instead of the Goldreich-Julian model for the NS plasma Krause-Polstorff and Michel (1985) provides larger or smaller radio signals depending on the polar angle and the misalignment between the magnetic dipole axis and the rotation axis.