Radio Signals from Axion Dark Matter Conversion in Neutron Star Magnetospheres

We show that axion dark matter (DM) may be detectable through narrow radio lines emitted from neutron stars. The neutron star magnetosphere hosts a strong magnetic field and a plasma frequency that increases towards the neutron star surface. As the axions pass through the magnetosphere, they can resonantly convert into radio photons in a narrow region around the radius at which the plasma frequency equals the axion mass. The bandwidth of the signal is set by the small DM velocity dispersion far away from the neutron star. We solve the axion-photon mixing equations, including a full treatment of the magnetized plasma and associated anisotropic dielectric tensor, to obtain the conversion probability. We discuss possible neutron-star targets and how they may probe the QCD axion parameter space in the mass range of ~0.2-40 $\mu$eV.

The QCD axion is one of the best-motivated dark matter (DM) candidates: in addition to explaining the observed abundance of DM [1][2][3], the axion may also resolve the strong CP problem [4][5][6][7]. However, testing the axion DM hypothesis is notoriously difficult, due to the fact that the axion is predicted to interact only very weakly with ordinary matter. In particular, axions interact with photons through L aγ = − 1 4 g aγγ aE · B, where g aγγ ∼ 10 −15 GeV −1 for an axion of mass m a = 10 −6 eV. Direct detection experiments such as ADMX [8][9][10] are just beginning to probe small regions of the QCD axion parameter space by exploiting this small coupling between the axion and electromagnetism. In this Letter, we discuss a new avenue for testing the axion DM hypothesis: axion DM indirect detection. We show that in the presence of axion DM, monochromatic radio signals are emitted from neutron stars (NSs) due to axion-photon conversion within the high magnetic field regions in the NS magnetosphere.
As we point out in this Letter, and was noted in [11], the finite electron density in the plasma within the NS magnetosphere gives the photon a mass m γ . 1 This is a key ingredient for obtaining resonant conversion, which can only take place when energy and momentum conservation are both satisfied. This photon mass is expected to fall off monotonically with distance from the NS surface. Thus for any axion mass smaller than the plasma mass at the surface, there exists a radius r c , the conversion radius, at which the axion mass equals the plasma mass and the axion-photon conversion process takes place resonantly. However, the plasma around a NS is highly magnetized, with the cyclotron frequency Ω c much larger than the plasma frequency ω p , and the effective photon mass is anisotropic with strong angular dependences. In this Letter, we solve the axion-photon equations of motion within the NS magnetosphere, including a full treatment of the magnetized plasma, to calculate the axionphoton conversion probability. We find that if the axion mass equals the photon mass at r c , then conversion takes place resonantly over a distance L ∼ r c v c /m a , where v c is the axion velocity at r c .
The process of axion-photon conversion in neutron-star magnetospheres is related to the detection mechanism utilized by axion helioscopes, such as the CAST experiment, to search for relativistic axions produced inside of the Sun which propagate to Earth [12][13][14]. Solar axions have a thermal spectrum, E ∼ keV, and can convert to photons within CAST's static magnetic field B. When the momentum transfer q = m 2 a /(2E) is much less than the length of the B-field region L (qL 1), then conversion takes place with probability P aγ ∼ g 2 aγγ B 2 L 2 . However, for large m a the conversion rate drops rapidly because it is no longer possible to satisfy both energy and momentum conservation simultaneously; the axion and photon have different dispersion relations, which becomes more apparent for larger values of m a and longer distances L. To circumvent this issue and maintain sensitivity to high-m a axions, CAST fills the B-field region with 4 He and 3 He to give the photon a mass so that the axion and photon have the same dispersion relation. By varying the pressure of the gas, the plasma frequency can be adjusted to scan over a range of m a values. With this technique, CAST has set some of the strongest limits on axion-like particles in the mass range m a ∼ 10 −4 -10 0 eV [13,14].
Most previous efforts to detect axion DM have focused on direct detection; see [15] for a detailed review.
The majority of such experiments utilize the coupling of the axion to electromagnetic fields, though some experiments, such as CASPEr [16], directly use the axion coupling to nucleons. For example, ADMX utilizes a microwave cavity to induce DM axion-photon conversion in the presence of an external B-field. ADMX has already constrained axion DM scenarios in a narrow mass range around 10 −6 eV, and future runs of ADMX should expand the sensitivity to axion masses in the range m a ∼ 10 −6 -10 −5 eV [10,17]. The HAYSTAC collaboration will try to push the reach to masses as high as 10 −4 eV using similar technology [18][19][20], while MAD-MAX [21] will probe a similar mass range with dielectric haloscopes [22]. A separate set of experiments, such as ABRACADABRA [23,24] and DM-Radio [25,26], are working to test the axion DM hypothesis at lower masses, potentially down to 10 −9 eV, by exploiting the coupling a E · B in the limit where the axion wavelength is much larger than the size of the experiment. Our work complements these approaches by providing an avenue for indirect detection of QCD axion DM in the mass range ∼ 0.2 − 40 µeV utilizing existing and planned radio telescopes.
NSs have long been recognized as promising targets for axion searches due to their strong magnetic fields. Previous efforts have focused on either photon-axion conversion leading to spectral distortions in the outgoing electromagnetic emission [27] or the conversion of thermal axions, produced inside of the NS, into photons in the NS magnetosphere [28]. However, neither of these processes require the axion to be DM, and also neither are sensitive enough to probe the QCD axion. In particular, thermal axions are ultrarelativistic and hence cannot undergo resonant conversion in the magnetosphere [29], but DM axions are only mildly relativistic and resonant conversion is obtained over a broad range of parameters. Our work builds upon this previous work by calculating the radio flux from axion DM conversion into photons within the magnetosphere. Assuming that the axion makes up all of the DM, we show that radio searches may be sensitive to QCD-axion-strength couplings g aγγ . 2 Neutron Star Magnetosphere. -The magnetic field in the vicinity of the NS surface is thought to be well described by a dipole configuration, with an axism that is misaligned from the rotation axis (which we take to be the z-axis) by an angle θ m . Charged particles are stripped from the surface of the NS at the magnetic poles and accelerated along open field lines, producing the non-thermal, pulsed radio and gamma-ray emission seen from pulsars. These regions near the magnetic poles are characterized by a high-density, boosted plasma. On the other hand, the NS "lobes" consist of closed magnetic field lines and likely much more modest plasma densities. We will take a simplistic approximation to the neutronstar magnetosphere, inspired by Goldreich and Julian (GJ) [30]. The GJ model gives the minimum plasma density necessary in the presence of the NS magnetic field, by finding a self-consistent solution to Maxwell's equations when particles on magnetic field lines corotate with the star.
Though originally proposed for aligned NSs with θ m = 0, the GJ derivation applies equally well to misaligned NSs, and gives a charge density where Ω = 2π/P with P the NS spin period, and θ is the polar angle with respect to the rotation axis. We will take the charge density as a rough estimate of the electron number density: n e = |n c |. 3 The plasma frequency is ω p ≈ 4παn e /m e , so that within the GJ model where is the component of the magnetic field along theẑ direction. Note that m ·r = cos θ m cos θ + sin θ m sin θ cos(Ωt) (4) depends on time due to the rotation of the NS. In (2) we have neglected the relativistic correction in the denominator of (1), which can be important for millisecond pulsars but is typically a percent-level correction for the pulsars with large P that we will be concerned with in this analysis. In practice, the true plasma density is likely more complicated than the simple GJ model. In particular, there could be non-trivial time dependence and boosts within the plasma. However, the GJ model provides a straightforward starting point for this analysis, which we hope can be improved in future work with more realistic models for the NS magnetosphere. In this analysis, we focus only on the region of closed field lines where the plasma is expected to be nonrelativistic, leaving the complications of boosted plasma near the magnetic poles to future work. As we show below, the axion-photon conversion occurs resonantly within the vicinity of the conversion radius r c , defined to be the radius at which the plasma frequency equals the axion mass. Using the expressions above, we find that in the GJ model the time-dependent r c is given by Conversion Probability from Mixing Equations. -Since the axion DM starts out non-relativistic far away from the NS and is accelerated to semi-relativistic velocities at radius r c , we can approximate the axion trajectories as radial. In the Supplementary Material (SM) we give a set of physical arguments that may be used to understand the parametric dependence of the axion-photon conversion probability. Here, we calculate the conversion probability by solving the coupled wave equations for the axion-photon system in the presence of the interaction term − 1 4 g aγγ aE · B in the Lagrangian, which leads to mixing between the axion a and the component of the photon vector potential A that is transverse to the axion's motion but coplanar with the magnetic field.
Following [29], we assume radial plane wave solutions of the form a(r, t) = ie iωt−ikrã (r) and A (r, t) = e iωt−ikrÃ (r), where k 2 = ω 2 − m 2 a . As we will show, the resonant conversion takes place in a narrow enough region around r c that we may neglect the r dependence of ω. Similarly, while the dispersion relation for k holds for both the axion and the photon at r c , the photon dispersion changes away from the conversion radius due to the continuously varying plasma mass. We account for both of these effects in turn. The analytic arguments presented below are supported by a full numerical analysis in the SM, where we also derive the equations of motion for the coupled axion-photon system in the plasma.
Near r c , we may use the WKB approximation |Ã (r)| k|Ã (r)| and |ã (r)| k|a (r)|. The mixing equations reduce to the first-order ordinary differential equation where ξ = sin 2θ ω = m a 1 + v 2 c , and k = m a v c . Above, we have defined θ to be the angle between the propagation directionr and the magnetic field B.
For r r c , the axion-photon system no longer strongly mixes, but the amplitude of A modulates due to the varying plasma frequency of the medium. This effect is familiar from wave mechanics as it is exactly analogous to the increasing amplitude of ocean waves as they approach the shore (though in our analysis, we are considering the opposite case of waves leaving the shore). The net effect is a suppression of the outgoing electromagnetic wave by a factor of √ v c asymptotically far away from the NS, To calculate the energy flux in electromagnetic radiation asymptotically far from the NS, we may use the formalism of transition amplitudes by analogy to timedependent perturbation theory in the Schrödinger equation, working to first order in ∆ B [29]. Taking initial conditionsÃ (r 0 ) = 0 andã(r 0 ) = a 0 , and neglecting the modulation of the outgoing electromagnetic wave for now, (6) gives This expression represents |A (r)| 2 /a 2 0 and may be interpreted classically as the ratio of energy density in the electromagnetic field to the energy density in the axion field at a radius r. Taking r → ∞ and including the amplitude modulation of the outgoing electromagnetic field, we evaluate (8) by the method of stationary phase to obtain the energy transfer function with L = 2πr c v c /(3m a ) and p ∞ aγ ≡ v c lim r→∞ p aγ (r). Note that L may be interpreted as the distance over which the resonant conversion takes place atθ = π/2. While derived forθ = π/2, the expression in (9) holds for genericθ to leading order in v c .
Radiated Power. -Next, we calculate the electromagnetic power emitted from the NS by axions converting into photons. We first note that since the NS plasma is optically thin, Thomson scattering of photons is negligible for the long-period NSs under consideration [29], and thus outgoing photons do not scatter. Because m −1 a r c for m a in the MHz-GHz range, L is parametrically smaller than r c , and thus conversion takes place in a small region around r c . We thus estimate the radiated power P by multiplying the flux of DM through a surface subtending a solid angle dΩ at r c by the energy transfer function: where ρ rc DM is the DM mass density at r = r c and v c is the DM velocity at r = r c . We note that all quantities on the right-hand side of (10) depend on θ, θ m , and t through their dependence on r c (see (5)). The factor of two comes from the fact that the DM may convert into photons either on its way in to the NS or out of the NS; if it is converted on the way in, then the photon is reflected back out, since the higher-density plasma acts as a mirror to photons of frequency less than the plasma frequency.
Let us assume that asymptotically far away from the NS, the DM has density ρ ∞ DM and is described by a Maxwell-Boltzmann velocity distribution with where v 0 ∼ 10 −3 is the DM virial velocity. We can calculate v c by conservation of energy: v 2 where v ∞ is the DM speed asymptotically far away from the NS, which is typically much smaller than the escape velocity 2GM NS /r c . Liouville's theorem maps the phase-space distribution from asymptotic infinity to r c : where v ∞ [v] denotes the velocity at asymptotic infinity that gives velocity v at radius r c . Integrating (12) and expanding in the limit v 2 This then leads to an expression for the radiated power: with Both (14) and (15) are formally only valid so long as r c > r 0 ; no resonant conversion takes place inside the NS. This regulates the otherwise-divergent denominator in (15), and also gives a strong angular and time dependence to the signal. Radio Telescope Sensitivity. -The radio flux at Earth is given by F (θ, θ m , t) = dP(θ, θ m , t)/dΩ/d 2 , where d is the distance from us to the NS. However, the quantity that is more relevant for radio observations is the flux density, defined as S = F/B, where B is the bandwidth. Energy conservation implies that the expected bandwidth of the signal is B ∼ (v 0 /c) 2 m a /(2π). This is because the frequency of the emitted photon asymptotically far away from the NS must be equal to the energy of the infalling DM particle that created it, regardless of the radius at which the DM was converted. In other words, the kinetic energy gained by the infalling DM is exactly canceled by the gravitational redshift of the outgoing photon, and the bandwidth of the signal is determined by the velocity dispersion in the asymptotic DM distribution. This leads to the estimate for the flux density This should be compared to the minimum detectable flux at a radio telescope, which is given by where SNR min is the minimum signal-to-noise ratio, SEFD is the system-equivalent flux density, n pol is the number of polarizations (we will take n pol = 2), B is the bandwidth, and ∆t obs is the observation time. As an example, the Arecibo Telescope has SEFD ∼ 2 Jy. Eq. (16) holds for sources whose bandwidth is wider than the intrinsic frequency resolution of the telescope, which we will assume is always the case. In the SM we discuss the optimization of the bandwidth for axion DM radio signals and how to account for the non-trivial time-dependence of the light curve in our sensitivity estimates. Neutron star targets. -As pointed out in [11], nearby isolated neutron stars (INS) make excellent targets for radio signals from axion DM conversion. This group of ∼7 NSs is characterized by their proximity to Earth ( < ∼ 500 pc), strong magnetic fields (∼10 13 G), long spin periods (∼5 s), and lack of observed pulsed, non-thermal emission (see, e.g., [31][32][33]). Importantly, since these NSs do not exhibit radio emission, we can estimate the sensitivity to axion DM assuming that we are limited by thermal noise in the radio telescope rather than background radiation from the NS. Additionally, the lack of non-thermal emission suggests that pair production at the NS surface is inefficient and that the GJ model for the plasma density may hold throughout the magnetosphere [11].
In Fig. 1 we show the sensitivity in g aγγ from 100 hrs of observation of one of the isolated NSs, J0806.4-4123. This NS has a period P ≈ 11.37 s, magnetic field B 0 ≈ 2.5 × 10 13 G, and is at a distance d ≈ 250 pc from Earth [33]. We also assume M NS = 1M and r 0 = 10 km. The QCD axion is predicted to lie within the band, while limits from CAST and ADMX (current and projected) are indicated. We have taken θm = 10 • and the solid (dashed) curves assume θ = 90 • (θ = 120.5 • ). The lower mass cutoff is set by the lowest available frequency of current radio telescopes, while the high-mass cutoff comes from requiring the conversion radius to be outside the NS radius.
We take SEFD ≈ 2 Jy for our estimates, though this may be improved with future instruments such as the Square Kilometer Array. Our sensitivity curves are defined by 1σ significance, as discussed in the SM. We show sensitivities calculated for two pulsar geometries. The solid curve takes a generic value θ = 90 • for the polar angle of Earth in the NS frame, while the dashed curve is tuned to θ = 120.5 • , which gives the near maximal signal at low masses and which is also highly pulsed from emission at orientations where r c → r 0 . In both cases we take a generic misalignment angle θ m = 10 • . Note that the low-mass cutoff is set to the m a /(2π) = 50 MHz threshold of typical radio telescopes, while the high-mass cutoff is determined by the maximum mass for which the conversion radius is outside the NS radius. Another class of targets are NSs that occupy regions of high DM density and/or low velocity dispersion. For example, consider the magnetar SGR J1745-2900, which is located R ≈ 0.1 pc away from the Galactic Center [34][35][36][37]. This magnetar has a magnetic field B 0 ≈ 1.6 × 10 14 G and a period P ∼ 3.76 s [34,35]. While the magnetar was first discovered in X-rays, a highly pulsed and variable ∼mJy radio signal has been observed from the magnetar (see, e.g., [36,37]). The DM density in this region is highly uncertain. Using the NFW and Burkert fits from [38], we find that the DM density at R = 0.1 pc is enhanced by a factor of 2 × 10 5 for the best-fit NFW profile, relative to the local density, but only a factor ∼4 for the best-fit cored Burkert profile. The cored profiles, however, may be in tension with new data from the Galactic bulge (see, for example, [39,40]). On the other hand, if the DM distribution is described by a generalized NFW profile with an index γ = 1.5, which is allowed by the kinematic data available, then the enhancement would be ∼10 7 . There is also the possibility of a DM density spike near Sgr A*, the supermassive black hole at the center of the Galaxy. With the density spike, the DM density at R = 0.1 pc could be enhanced by a factor ∼10 9 relative to the local density [41].
In Fig. 1, we show the projected sensitivity from 100 hrs observation of SGR J1745-2900, assuming both the NFW DM profile (blue) and spike profile (purple), with solid and dashed lines representing the two geometries θ = 90 • and θ = 120.5 • . We take v 0 = 200 km/s and d = 8.5 kpc for the distance to the Galactic Center and assume M NS = 1M and r 0 = 10 km as before. Despite the fact that pulsed radio emission has been observed from this magnetar, we have made these sensitivity estimates under the assumption that the dominant noise source is the thermal noise in the telescope. Since the non-thermal radio emission is pulsed, non-pulsed (or pulsed but out-of-phase) DM-induced flux would likely still be dominated by telescope noise. Interestingly, as seen in Fig. 1, observations of SGR J1745-2900 could be sensitive to the QCD axion over multiple orders of magnitude in m a , depending on the DM density profile. However, we stress that this sensitivity estimate relies on the GJ model, which may not apply to the magnetar.
Alternatively, one could consider isolated NSs within dwarf galaxies. In the Sagittarius dwarf galaxy, the central DM density is enhanced by a factor ∼5 × 10 5 compared to the local density in the solar neighborhood, and the velocity dispersion is low, v 0 ∼ 10 km/s [42]. For this estimate we have taken the cored DM density profile from [42]. The globular cluster M54 appears to be coincident with the center of the Sagittarius dwarf galaxy, with the cluster having a core radius ∼1 pc, a mass ∼2 × 10 6 M , and a distance of around ∼20 kpc from Earth [43,44]. Given the mass of M54, there are likely many hundreds of NSs within the central core [45]. Assuming that just one of these NSs has the properties of J0806.4-4123, we would obtain the sensitivity to g aγγ shown in Fig. 1 (labeled INS in M54). If there are N such INS's in the field of view, then we may expect the sensitivity to improve as 1/ √ N . The fact that all NSs radiate at the same frequency from axion DM could make even more distant galaxies promising targets.
A narrow radio line from a NS target could provide a striking signature of axion DM. On the other hand, in the absence of a signal, it will be difficult to set a robust limit on g aγγ because of challenges in understanding confidently the plasma density and time-dependent dynamics in the inner regions of the magnetosphere. Towards that end, it would be useful to incorporate the physics of axion-photon conversion into NS simulations [46]. Such work should lead to more precise predictions for the radio-line signal.
Note added -While this work was in the final stages of preparation, Ref. [47] appeared, which addresses similar questions. Our work differs in several respects, but importantly where we do overlap we disagree in detail with their results for the conversion probability, radio flux, and projected sensitivity.

Radio Signals from Axion Dark Matter Conversion in Neutron Star Magnetospheres Supplementary Material
Anson Hook, Yonatan Kahn, Benjamin R. Safdi, Zhiquan Sun This Supplementary Material contains additional calculations and examples that support the results presented in the main Letter and further illustrate the phenomenology of a hypothetical signal. We begin by giving a heuristic derivation of the conversion probability that does not rely on having to solve the coupled differential equations but rather is based on more general physics arguments. We also give an extended derivation of the mixing equations and their solution, including the effects of strong magnetic fields in the plasma, along with numerical examples that support the analytic solution for the conversion probability used in the main Letter. Next, we discuss the non-trivial light curves and polarization profiles that might be expected from an axion signal. Finally, we describe how to estimate the radio sensitivity, accounting for the optimal bandwidth for an analysis of the radio data given knowledge of the DM velocity distribution and also knowledge of the non-trivial light curves, and highlight the strong angular dependence of the signal.

CONVERSION PROBABILITY PARAMETRICS
We first demonstrate how to understand the parametric dependence of the axion-photon conversion rate. Parametrically, the axion photon conversion rate is of the form P (a → γ) ∼ sin 2 Θ sin 2 (∆kL) (S1) where tan Θ ∼ Bg aγγ ω/(m 2 a − m 2 γ ) is the mixing between the axion and the photon, ∆k is the difference in momentum between an axion and photon of the same energy, and L is the distance over which the conversion occurs. As is well known, relativistic axions converting in a magnetic field in vacuum have m γ = 0 and ∆k ∼ m 2 a /ω so that Θ 1 and the conversion probability scales as This gives the conversion probability familiar from experiments such as CAST.
In this work, we consider the very different regime where axion-photon mixing is maximal due to m γ ∼ m a , but the axion is non-relativistic. As a direct consequence of the assumption that the photon effective mass is the same as the axion mass, we have Θ ∼ O(1). When mixing is maximal, the difference in momentum between the two propagating eigenstates is ∆k ∼ Bωgaγγ k , where we are using (6) (or (S9) below) and ignoring anyθ dependence for simplicity. Thus the conversion rate is given by Note that while (S3) is similar to (S2), the derivation and region of validity is completely different. For example, if the axions were assumed to be relativistic instead of non-relativistic, the conversion probability would instead be highly suppressed by the large magnetic fields of the NS [29]. Because NSs are macroscopic objects with large magnetic fields that extend over macroscopic distances, the length scale L is not determined by the change in the magnetic field, but instead by the change in the plasma frequency of the photon. Eq. (S3) assumes coherent conversion, i.e. that photons which were generated from conversion at r and r + L add coherently. If the plasma mass changes over this distance, then photons generated at r 1 can instead interfere with photons generated at a different radius r 2 . L is roughly the distance over which photons generated at the beginning start to interfere with photons generated at the end, i.e. L ∼ 1/δk, where δk is the difference in momentum at different locations. kδk ∼ m γ δm γ and due to the power-law dependence of the photon effective mass on r, we also have δm γ /m γ ∼ L/r c . Thus we have the estimate that We note that our result is parametrically different from the estimate in [11], which assumes L ∼ r c . The final axion photon conversion rate is In the main Letter, and in more detail below, we carefully calculate the axion-photon conversion probability from the mixing equations, but the parametric scaling can be understood by the arguments just put forth.

AXION-PHOTON CONVERSION
In this section, we give additional details of the equations of motion of the axion-photon system and the approximate solution to these equations for axion DM in the NS magnetosphere.

Equations of motion
The plasma surrounding a NS is a cold, highly-magnetized plasma and due to the large magnetic field it is strongly birefringent. The effect of the magnetic field and plasma is taken into account by introducing a dielectric tensor. Because the photon has both transverse and longitudinal modes, it is convenient to work directly with the electric field. The propagation of electromagnetic waves and axions in a plasma is determined by where the magnetic field B is taken to be the external magnetic field due to the NS. We work over a small enough region of space where the change in gravitational potential may be neglected. The electric displacement field D is given by [48] where the magnetic field is taken to be at an angleθ from the z-axis in the positive y − z quadrant, and R yz θ is the rotation matrix byθ in the yz-plane. The coefficients in the dielectric tensor are given by Taking the Fourier transform in time, and going to the high-magnetization limit (Ω c ω, ω p ), which sends → 1 and g → 0, E x decouples from the equations, and E z does not propagate and can be solved for algebraically. The mixing matrix simplifies to where B t = |B| sinθ is the component of the B-field transverse to the direction of motion. Equation (6) is derived from a WKB approximation of this second-order differential equation, assuming outgoing plane waves in the radial direction with a local coordinate system defined byr =ẑ and defining A = E y /iω. Note that in (S9), we have neglected the vacuum birefringence term from strong-field QED. Ref. [29] points out that this term can drastically affect the axion-photon conversion probability; however, this term is negligible for non-relativistic DM axions. Vacuum birefringence changes the condition for resonance from ω p = m a to Ref. [29] points out that for sufficiently large κ and ultrarelativistic axions with ω m a , the maximal mixing condition (S10) can never be satisfied because the left-hand side becomes negative. Fortunately, for nonrelativistic or mildly relativistic axions with ω = O(1) × m a , this is never a problem. Indeed, in this case we have (S13) Even for the largest NS magnetic fields considered in the Letter, the effects of the vacuum birefringence term are a percent-level perturbation and can be neglected for semi-relativistic axions, and ω p = m a can always be satisfied. More specifically, as long as B < ∼ 10 15 G (the size of the largest known magnetar fields), the maximal mixing conditions can always be satisfied, but for these large fields one would need to consider the effect of κ when solving for r c .

Damping of the photon wave
As mentioned in the main Letter, it is important to take into account the damping of the outgoing electromagnetic wave in calculating the energy transfer from the axion field to the outgoing electromagnetic radiation. To isolate this effect, we neglect the mixing terms and consider the photon wave equation with a spatially-dependent plasma mass: We now define k(r) = ω 2 − ω 2 p (r) and use a second WKB approximation, performing an expansion in the small dimensionless parameter ε = 1 k 2 dk dr . We have verified that this parameter is small for r > r c . This gives the solution (for r > r c ) Solving (6) for the photon field leads to the expression in (9).
The axion-photon evolution at finite r In the main Letter, we presented an approximate expression for p aγ valid at r → ∞. In this subsection, we give the result at finite r. First, we considerθ = π/2, whereθ is the angle between the B-field andr, and then we discuss the generalization to arbitraryθ. Takingθ = π/2 and performing the integral in (9) at finite r, we find with L = 2πr c v c /(3m a ) and the function G(x) defined by in terms of the Fresnel C and S integrals. Note that lim x→∞ G(x) = 1 and that G r−rc L rises quickly from 0 at r < r c to values near unity at r ∼ r c , over a distance of order L, while the amplitude modulation of the outgoing wave is encapsulated by the second line in (S16).
At arbitraryθ, the length L, which is found from the stationary phase approximation to the integral in (8), is modified to and we should also make the substitution in (S16) . (S19) Moreover, the conversion radius r c is modified to (S20) Note that B(r c ) 2 L 2 , the quantity which appears in p aγ , has no explicit dependence onθ up to O(v 2 c ) corrections. In Fig. S1, we compare our analytic estimate (S16) to numerical solutions to the full 2 nd -order mixing equations in (S9), with boundary conditions imposed at r = 223.05 km (r = 223.5 km) for the blue (black and red) curves. We take the frequency ω ≈ m a (1 where v c is the DM velocity at the conversion radius. In this example, we take r 0 = 10 km, NS mass M NS = 1 M , B 0 = 10 14 G, θ m = 0, P = 1 sec, m a = 1 GHz, and g aγγ = 10 −12 GeV −1 . This implies that the transition radius is at r c ≈ 224 km and that at this radius the axion velocity is v c ≈ 0.11, in natural units, where we neglect the O(10 −3 ) initial axion velocity v 0 asymptotically far away from the NS. We illustrate three different angles, θ = 90 • , 72 • , and 45 • ; for each angle, the angleθ betweenr and the magnetic field is indicated in the figure. As we change θ, we change the conversion radius r c given in (5) because of the angular dependence of the plasma frequency. However, in order to highlight the differences between the differentθ we chose to fix the normalization of the plasma frequency to that found at θ = 90 • . That is, at the level of the equations of motion in (S9), we keep the implicit dependence of ω p on θ fixed and only vary the explicit dependence ofθ. In the left panel, the solid curves show the numerical solutions, while the dashed show the analytic approximations, which are seen to agree well. Indeed, we see that both the peak values and the asymptotic values for both the analytical and numerical solutions are independent ofθ (and hence θ). This is further illustrated in the right panel, where we compare the numerical solutions to the asymptotic value for p ∞ aγ (black, dashed), which is the same for allθ to leading order in v c .
We note that there are small discrepancies between the analytic approximation (9) and the numerical solutions in these examples. These may be due in part to the necessity of setting the boundary conditions very close to r c rather than at r 0 to avoid contamination of the solution by the spurious exponentially-growing mode when ω < ω p . Still, the difference between the analytic and numerical results is less than ∼ 10% at large r.

LIGHT CURVES AND POLARIZATIONS
The misalignment between the pulsar's magnetic axis and rotation axis leads to non-trivial light curves. Fig. S2 shows the change in the conversion radius and the radiated power over the pulsar's rotation period as a function of the rotation phase for an example NS.
The misaligned NS also leads to non-trivial polarization structure. This is because the electric field of the radio emission is aligned with the tangential component B t of the NS's magnetic field, but the direction of B t changes as a function of Ωt. A more careful analysis is needed to verify if the polarization structure survives Faraday rotation induced from the magnetosphere, but we neglect this effect for now and study the polarization of the photons as they are emitted. We define the basis vectorsˆ 1 =θ andˆ 2 =φ for the polarization, whereθ andφ are the polar and azimuthal unit vectors in the NS's frame. Without loss of generality, we may consider the Earth to be at φ = 0 in the frame of the NS, such thatm ·r is as given in (4). The time-dependent polarization vectorn(t) of the radio emission is given byn cos θ m sin θ − sin θ m cos θ cos(Ωt) ˆ 1 − sin θ m sin(Ωt)ˆ 2 cos θ m sin θ − sin θ m cos θ cos(Ωt) 2 + sin 2 θ m sin 2 (Ωt) .
In Fig. S3, we illustrate the linear polarization profiles, neglecting the possible effect of Faraday rotation, in theˆ 1 andˆ 2 directions over a pulsar period. Note that for θ = 90 • and a small misalignment angle θ m = 15 • , the outgoing wave is almost completely polarized over the whole rotation period.

RADIO SENSITIVITY ESTIMATES
In this section we give additional details for how we make the radio sensitivity estimates in the main Letter.  Figure S3. Theˆ 1 (black) andˆ 2 (blue) polarization components of the electric field as a function of the phase Ωt for various θ, with θm = 15 • . Note that these polarization curves are valid for any NS with these alignments, though depending on ma and the NS properties, the NS may not radiate over the whole period since rc may drop below the conversion radius (see Fig. S2).

Bandwidth optimization
Considering that the asymptotic DM velocities, which determine the spread of photon frequencies, have a dispersion ∼v 0 , it follows that we expect the bandwidth B of the radio signal to be B ∼ m a v 2 0 /(2π), where the 2π simply arises from the translation from angular frequency. However, it is worth calculating the optimal bandwidth for an analysis of the radio-telescope data given that we are trying to optimize SNR Having a smaller bandwidth increase our signal-to-noise ratio because of the √ B in the denominator, but on the other hand this also means that we have less flux F , which acts to decrease the signal. In this section, we compute the optimal bandwidth by maximizing F/ √ B. Let us define where I(ω 1 , ω 2 ) is the function to maximize and where f ω (ω) is the distribution of radio frequencies in the lab-frame for an axion signal with mass m a . In the DM frame, the DM velocity distribution is given by (11). Let us assume that the NS is boosted with a velocity v b with respect to the DM frame. Then, the distribution of frequencies in the NS frame, but asymptotically far away from the NS, is given by where we have definedω ≡ ω/m a . This implies that with x ≡ v b /v 0 and ω 1,2 = m a + m a δω 1,2 v 2 0 . We have also defined B ≡ (ω 2 − ω 1 )/2π andB ≡ 2πB/(m a v 2 0 ). For x = 0, we find the maximum atδω 1 = 0.013 andB * ≈ 1.12, with value I(w 1 , w 2 ) m a v 2 0 /(2π) ≈ 0.74. For x = 1, we find the maximum atδω 1 = 0.032 andB * ≈ 1.97, with value I(w 1 , w 2 ) m a v 2 0 /(2π) ≈ 0.58. Some explanation is warranted here. First, we have not accounted for the gravitational potential of the NS for the reason mention in the main Letter: a DM wave with energy ω asymptotically far away from the NS will produce an electromagnetic wave of frequency ω when measured asymptotically far from the NS, since both the axion-photon mixing equations and the NS's gravitational potential conserve energy. Second, we have not accounted for the boost of the lab frame with respect to the NS because this boost only affects the radio emission; such small boosts are not important for relativistic particles. On the other hand, the NS boost with respect to the DM distribution is important because DM is non-relativistic. withm ·r and r c (t) functions of time. Then, we see that

Time-dependent light curves
where ∆t obs is the total time the radio telescope observed, and φ is the azimuthal coordinate on the NS. That is, m ·r = cos θ m cos θ + sin θ m sin θ cos(φ) , r c (φ) = 224 km × r 0 10 km where we can also identify φ = Ωt. Evaluating (S32) and solving for g aγγ such that χ 2 Asimov = 1 leads to the projected sensitivities in Fig. 1.

Angular dependence of sensitivity
In Fig. 1 we showed the projected sensitivity to g aγγ as a function of mass for two different angles θ, with θ m = 10 • . However, it is important to understand how the sensitivity changes a function of θ, given that in practice this will be a random angle that depends on the relative orientation of the Earth and the NS being observed. In Fig. S4 we show how the sensitivity changes as a function of θ for the two cases θ m = 10 • , 30 • . We illustrate two different masses, m a = 2 × 10 −7 eV and m a = 2 × 10 −6 eV, and normalize the sensitivity to the maximum (worst) value, g max aγγ , found over all θ. Note that the large enhancement due to the strong beaming is found over a relatively wide range of θ, so we may expect such an enhancement for generic NS targets.