Radio Signal of Axion Photon Conversion in Neutron Stars: A Ray Tracing Analysis

Axion dark matter can resonantly convert into photons in the magnetospheres of neutron stars (NSs). It has recently been shown that radio observations of nearby NSs can therefore provide a highly sensitive probe of the axion parameter space. Here we extend existing calculations by performing the first three-dimensional computation of the photon flux, taking into account the isotropic phase-space distribution of axions and the structure of the NS magnetosphere. In particular, we study the overall magnitude of the flux and its possible time variation. We find that overall signal strength is robust to our more realistic analysis. In addition, we find that the variance of the signal with respect to the NS rotation is washed out by the additional trajectories in our treatment. Nevertheless, we show that SKA observations towards J0806.4-4123 are sensitive to $g_{a\gamma\gamma}\sim 3\times10^{-13}\mathrm{\,GeV^{-1}}$ at $m_a\sim 7\times10^{-6}\mathrm{\,eV}$, even when accounting for Doppler broadening. Finally, we provide the necessary code to calculate the photon flux for any given NS system.


I. INTRODUCTION
The QCD (quantum chromodynamics) axion was first introduced in 1977 by Peccei and Quinn as a solution to the strong CP problem of the QCD sector [1][2][3][4]. Depending on the production mechanisms at play in the early Universe, the QCD axion can behave like cold collisionless matter, therefore allowing it to account for a fraction, or all, of the dark matter (DM) in the Universe [5][6][7]. The QCD axion is therefore considered one of the most well-motivated DM candidates to date. The assumption that the Pecci-Quinn (PQ) symmetry is broken after inflation leads to the axion mass of m a ¼ 26.2 AE 3.4 μeV [8], though larger uncertainties may arise from the contribution of topological defects [9][10][11]. If, on the other hand, the PQ symmetry is broken before inflation this mass constraint can be relaxed to give the classical window of 10 −7 eV ≲ m a ≲ 10 −4 eV [12][13][14][15][16][17].
A large range of observational strategies to directly detect axion particles now exist [18][19][20][21][22][23][24][25][26][27][28][29][30][31][32][33][34] (see Ref. [35] for a recent review). Many of these experiments exploit the axion's coupling to electromagnetism L ¼ −ð1=4Þg aγγ F μνF μν a ¼ g aγγ E · Ba to induce axion-photon conversion in the presence of magnetic fields. For the QCD axion, the axion mass and axion-photon coupling strength are proportional, m a ∝ g aγγ [36][37][38][39]. On the other hand, axionlike particles (ALPs) do not have the same coupling/mass relation and can therefore take on a wider variety of parameter combinations. Although not connected to the strong CP problem, ALPs are a generic prediction from the spontaneous breaking of approximate global symmetries in beyond the standard model physics as well as compactifications of higher dimensions in string theory [40,41]. ALPs therefore represent a prime target for searches of new physics.
As well as terrestrial direct searches, a variety of astrophysical observations can be used to constrain the ALP parameter space. For example, Refs. [42][43][44] examined the possibility of stimulated axion decay for various astrophysical targets, showing that dedicated observations of dwarf spheroidal galaxies can potentially improve upon current constraints by an Oð1-10Þ factor. Observations of stellar lifetimes also serve as a sensitive probe of ALPs. The hot plasma in the interior of a star readily produces low mass particles which allow for energy transport out of the stellar environment and have the potential to change a stars normal evolution [45][46][47].
Reference [48] recently suggested that the magnetospheres of neutron stars (NSs) can induce enough axionphoton conversion to be subsequently observed by radio telescopes, potentially probing QCD axions. 1 Importantly, the finite electron density of the plasma induces an effective photon mass [50] which, when equal to the axion mass, allows for the conversion to become resonant. Since, in the simplest scenarios, the photon mass monotonically decreases with the distance to the NS surface, there exists a continuum of resonant conversion surfaces corresponding to different axion masses. The conversion process also conserves energy (up to Doppler broadening [51]), which allows for an inference of the axion mass directly from the radio signal. This is particularly valuable since most terrestrial experiments must fine tune their experimental setups to gain sensitivity to particular axion masses. The two approaches of direct and indirect detection are therefore highly complementary.
Radio signals from a collection of NSs, such as the bulge population in the center of our galaxy [52], were investigated in Ref. [53]. Multimessenger signals (radio and gravitational waves) from black hole-neutron star inspirals were studied in Ref. [54]. A detailed study of the mixing equations used to describe axion-photon conversion is provided in Ref. [51,55].
In this work we present a more complete treatment of the radio signal calculation for individual NSs. To this end, we extend the recent analytic treatment presented in Ref. [48] by performing a numerical ray-tracing computation of the conversion process. This allows us to fully account for the isotropic phase-space distribution of axions in the vicinity of the NS. We find significant qualitative and quantitative differences with respect Ref. [48], and provide the necessary code to calculate the flux for various parameter combinations. Our results impact the overall reach of future radio searches for QCD axions and ALPs, as well as the optimization of realistic search strategies.
The paper is organized as follows. In Sec. II we describe the signal calculation, paying particular attention to the raytracing algorithm. In Sec. III, we report our results for the NS J0806.4-412 and calculate the sensitivity of nextgeneration radio telescopes to the radio signal and its potential time variability. Finally, we conclude in Sec. IV.

II. SIGNAL CALCULATION
Here we present the formalism and assumptions behind the ray-tracing method. In particular, we discuss the axionphoton conversion probability, the dark matter distribution, the neutron star's magnetosphere, and the photon flux seen on Earth. Where appropriate, we follow the calculations from Ref. [48].

A. Conversion probability
Photons in a plasma acquire an effective mass ("plasma mass") through interactions with free charges, which is given by [49] where n c is the charge carrier number density, m c the charge carrier particle mass, and α is the fine-structure constant. Resonant conversion can occur when the photon plasma mass approximately matches the axion mass m a ≃ ω p . This condition singles out a resonant conversion shell which is dependent on the axion mass. Using the WKB and stationary phase approximations, one can show that the probability of an axion converting into a photon while traversing a resonant conversion region is given by Here, g aγγ denotes the axion-photon-photon coupling mentioned above and B ⊥ is the strength of the NS magnetic field perpendicular to the axion trajectory. The plasma mass derivative ω 0 p denotes the derivative along the axion trajectory at the point of resonant conversion. For the special case of radial trajectories, it is given by jω 0 p j ¼ jdω p =drj ¼ 3 2 m a =r (as in [48]). Further details about the calculation of the conversion probability can be found in Appendix, including a critical comparison with earlier literature.

B. Phase space distribution of dark matter at the neutron star's surface
The phase-space distribution (PSD), fðr; vÞ, describes the statistical properties (spatial positions, r, and velocities, v) of a group of particles. Assuming that the PSD is stationary gives Liouville's theorem [56] from which we can see that the PSD is conserved along the trajectories of the system. We can therefore equate the PSD at infinity to the PSD at the NS surface, where the subscript infinity refers to quantities far from the NS. The distribution of DM is assumed to be isotropic in the rest frame of the galaxy. The standard halo model predicts an isotropic Maxwellian distribution far from the NS [57,58], given by where ϱ DM ∞ is the DM density and v 0 is the spread of the distribution (discussed below). Equations (3) and (4) show that the PSD close to the NS surface will be isotropic as long as v 2 ∞ does not depend on the direction of v. We can see this independence directly from energy conservation which allows us to relate the velocity at infinity to the physical velocity seen by a local observer at radius r (the radial coordinate of the Schwarzschild metric) from the NS center (in the limit v 2 where G is Newton's constant and M NS is the mass of the NS. The isotropy of the PSD at the NS surface holds as long as the NS velocity is small with respect to the galactic rest frame. 3 Equations (10) and (11) of Ref. [59] directly yield the local DM PSD: Note that due to energy conservation, the minimum velocity at a given radius is v min ðrÞ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2GM NS =r p . Integrating Eq. (6) with respect to v such that jvj ≥ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2GM NS =r p , we obtain: The DM velocity dispersion v 0 is small enough that, in practice, we can take the limit x ≫ 1 which allows us to neglect the right-hand side (rhs) of Eq. (7) and work with the simpler expression: Furthermore, conservation of energy for infalling axions further yields the DM velocity In addition to the overall normalisation of the signal, we must worry about the width of the line. We consider two contributions to this width: firstly from the intrinsic velocity distribution of DM far from the neutron star and secondly from the overall Doppler broadening due to the collective motion of the magnetosphere associated with the NS spin. The former is given by the Maxwell-Boltzmann distribution which leads to B ∼ ðv 0 =cÞ 2 m a =ð2πÞ where B is the bandwidth (discussed below). The latter is described in Ref. [51], where we take a conservative approach by setting the conversion surface to be the maximum conversion radius for any given axion mass. 4 It is not clear that our approach is a good description of the true line width since each pixel contributes some fraction of the total observed flux. In addition, a full calculation of the width would need to consider reflection and transmission components of the signal separately and further investigate the turbulence in the plasma at the scale of the photon wavelength. We leave this to future work and instead show optimistic and conservative sensitivity curves from the two broadening effects.

C. Neutron star magnetosphere
As a concrete example, we use the Goldreich and Julian (GJ) model [60] as a description of the magnetosphere around the NS (for a review on pulsar magnetosphere see Ref. [61]). The GJ model assumes that the magnetosphere is corotating with the NS and provides the following analytic expression for the charge number density: where ω ¼ ð2π=PÞẑ is the constant NS rotation vector with P being the NS spin period, θ is the polar angle, and B is the local magnetic field. The latter is assumed to be in a dipole configuration with axis along the directionm. In spherical coordinates ðr; θ; ϕÞ, we have: ðcos θ m cos θ − sin θ m cos θ cos ψÞ 2 This follows from gravitational redshift where E ∞ and E r are the axion energy at infinity and radius r, respectively. 3 Note that this is not necessarily a good assumption since the NS and the DM are moving nonrelativistically. We will extend this formalism to account for boosts into the NS's reference frame in future work. 4 where θ m is the misalignment angle betweenẑ andm, ψðtÞ ¼ ϕ − ωt, B 0 is the magnetic field strength at the NS poles, and r NS denotes the NS's radial size. See Fig. 1 for a visualization of the system with the relevant angles labeled.
The plasma mass (1) is computed by taking n c ¼ jn GJ j and assuming the presence of electrons and positrons only (m c ¼ m e ). By neglecting relativistic corrections [second factor in Eq. (10)], we get where the time dependence resides in the term m ·r ¼ cos θ m cos θ þ sin θ m sin θ cos ψðtÞ: ð13Þ The resonant axion-photon conversion region is identified by the matching relation ω p ðr; θ; ϕ; tÞ ¼ m a . For a given angle θ, we would expect the radio signal to display some time variation as the NS rotates, due to varying resonant conversion surface observed from Earth. This variation is discussed in Ref. [48] where they assumes radial trajectories for the infalling axions. This assumption implies that the radio signal arises from the specific direction θ connecting the NS to the Earth. For an isotropic PSD, we expect that the absence of any preferred axion trajectory will suppress the time-dependence of the radio signal as it is instead given by the cumulative flux from all possible axion trajectories (as discussed below). Before concluding, we note that the NS magnetosphere model plays an important role in defining the local properties of the NS plasma. In the present paper, we consider the GJ model with electrons and positrons to make a direct comparison with the analytic results discussed in Ref. [48]. Importantly, the GJ model does not describe inhomogeneities in the plasma structure which, in real NSs, may exist. These inhomogeneities may significantly affect the signal and must be accounted for in future work [62]. In addition, the GJ charge density may significantly differ from the true charge density-this difference is typically described through the multiplicity which can vary significantly [63]. However, our ray-tracing algorithm is implemented in such a way that the magnetosphere model can be straightforwardly modified to account for more complicated and realistic scenarios capturing local perturbations [61]. The impact of different models (such as the electrosphere model [64] or numerical simulations of the NS magnetosphere [65][66][67]) on the signal prediction is left to future work.

D. Calculating the photon flux
We assume that axions can convert into photons only at points, denoted r c , where the plasma mass equals the axion mass. The conversion probability P a→γ ðr; uÞ depends on the local magnetic field strength (and therefore the position) as well as the direction of the axion at that point, given by u with respect to the radial trajectory at r c . Photons entering regions with a larger plasma mass will be reflected, providing a factor of two greater flux when accounting for trajectories parallel and antiparallel to the line of sight (as mentioned in the next subsection). 5 Reflected waves will be will be also Doppler broadened [51], which we will further discuss below. Note that we neglect multiply scattered photons.
The total flux expected from a neutron star can be written in terms of the intensity as where the angular integral is over a region ΔΩ which covers the neutron star. We first consider the emission through a planar conversion region (for a → γ) that is taken to be perpendicular to the line-of-sight toward the NS, and at a distance D from the observer. Its physical size is taken to be infinitesimally small, dA ¼ D 2 dΩ. We assume that the conversion region is immersed in an isotropic distribution of axions with number density n a and velocity v a . The current of axions that traverses the conversion region in direction dΩ 0 is given by In general, resonant conversion (associated with k 2 a ¼ ω 2 − ω 2 p ðrÞ) and reflection (associated with m a ¼ ω p ðrÞ) do not happen at the same place. For radial trajectories, one can show that the difference between both points is given by Δr ¼ r c 3 k 2 c =m 2 a (which follows from ω p ∝ r −3=2 ). The size of the resonant conversion region is given by L 2 ¼ 2π=3 · r c k c =m 2 a , from which follows that L ≃ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi 2πΔr=k c p . We have separation when λ c ≪ Δr, where λ c is the de Broglie wavelength at resonance.
where α is the angle between the line-of-sight from the neutron star to the observer and the direction of the axions. The first two factors give the size of the conversion region projected onto the axion direction. The third and fourth factors together denote the number density of axions moving in direction dΩ 0 , and the last factor is their velocity. The signal intensity is given by As usual dΩ is the angular size of the observed region, and dA 0 is the detector area. The second equality holds since we can exchange the role of solid angle and observed area (using dΩ ¼ dA=D 2 and dΩ 0 ¼ dA 0 =D 2 ). In the last step, we set cosðαÞ ≃ 1. The photon flux from a NS (still assuming a perpendicular conversion region) can therefore be written as where we took P a→γ to be the probability of axion-photon conversion upon traversal of the conversion region, and we split up the integral into sums over different line-of-sights i, which each contribute ΔΩ i ¼ Δb 2 =D 2 to the integral (and P i ΔΩ i ¼ ΔΩ). Note that these line-of-sights are very close to parallel at the scale of the neutron star.
In order to generalize to nonperpendicular emission planes, the integral over dΩ would have to be, in principle, replaced by an integral over the 2-dim sub-manifold C. However, we can simply parametrize this submanifold by the angular direction seen from the observer, Ω. This is possible since we only observe (from Earth) the first crossing of the conversion region (everything else is absorbed). In that case, the area of this submanifold can be calculated as where α is the angle between the normal of the submanifold C at each point and the line of sight. The last factor accounts for the deprojection of the submanifold when integrating over dΩ. Interestingly, this factor cancels the cos α that we obtained in Eq. (15). As a result, the equation for the flux F in terms of sums over line-of-sights, Eq. (17), remains the same for nonplanar emission regions, provided that the quantities in the parentheses are evaluated at the point of the conversion.

E. Computational approach
The total radio flux is computed through a ray-tracing algorithm defined by the following steps: (1) We define the region of interest (ROI) as a planar surface ΔA perpendicular to the line-of-sight toward the NS at a distance D from the Earth and a distance d from the NS center. The latter is chosen to be the maximum distance at which the resonant axion-photon conversion occurs. Typically, we have d ¼ Oð100 kmÞ for the minimum axion mass considered. The ROI is divided into square pixels of size Δb, whose centers identify a specific photon trajectory i. (2) For each pixel, we back-propagate the photon by numerically computing its geodesics starting from the center of the corresponding pixel and taking the initial velocity to be perpendicular to the ROI. where we require that the resonant conversion occurs outside the NS surface. The factor of two takes into account the reflection of photons on the way toward the neutron star. We do not consider the contribution of nonresonant conversion since it is generally subdominant. (5) The total radiated power is obtained by summing the contributions from all the pixels and the total radio flux is simply given by This expression matches the one reported in Eq. (17) multiplied by the axion mass. Throughout our analysis we consider a flat metric (classical approximation)-all photons propagate in parallel straight lines. This is a good approximation for the case at hand. The Schwarzschild metric differs from the Minkowski one by terms proportional to εðrÞ ¼ r s =r with r s being the Schwarzschild radius. Such a ratio reaches its maximum value at the NS surface where ε ¼ Oð0.1Þ. We checked that the total radiated power changes by no more than a few percent when considering the Schwarzschild metric. 6 There are also additional effects associated with the NS's spin and described by the Kerr metric. This Kerr metric reduces to the Schwarzschild one in the limit of ε 0 ðrÞ ¼ 2πI=Pr ≪ 1 where I and P are the moment of inertia and the spin period of the neutron star, respectively. For the isolated NS J0806.4-4123 analyzed in the next section, the maximum value of ε 0 ðrÞ at the NS surface is of the order of 10 −5 for reasonable NS moments of inertia [68].
Throughout, we also neglect any general relativistic corrections to our results. Since the Lorentz factors that we encounter for axion velocities can be as large as γ ∼ 1.2, this means that our results are in general only accurate to within AE20% (although the details depend on where most of the observable conversion happens, and we expect much higher accuracy in most cases that are of interest here).
We further checked that the numerical calculation of the radiated power is converged with respect to: (i) the spatial resolution along each trajectory that is used to identify the resonant conversion region; (ii) the pixel resolution of the ROI that sets the number of trajectories. We find that the contribution of each trajectory typically varies within 0.1% for a spacial resolution of 1 m. The pixel resolution is instead fixed by the requirement that the total radio flux does not vary by more than 1% when increasing the number of trajectories.

III. RESULTS
Up to now we have not fixed our calculations to a specific system apart from assuming that the magnetosphere is described by the GJ model. The signal is highly dependent on the system being observed. We therefore fix the reference model to be the isolated neutron star J0806. 4-4123, which has a period P ¼ 11.37 s, a magnetic field B 0 ¼ 2.5 × 10 13 G (1 G ¼ 10 −4 T), and is located at a distance D ≈ 250 pc from the Earth [69]. Moreover, we take M NS ¼ 1 M ⊙ and r NS ¼ 10 km for the mass and the size of the neutron star as well as ϱ ∞ DM ¼ 0.3 GeV=cm 3 and v 0 ¼ 200 km=s for the local density and velocity dispersion of DM particles, respectively. In the following, we report the predictions for the radio flux for different angular configurations of the system, and discuss the sensitivity of next-generation radio telescopes to the signal and its time variability. Note that we expect the results of J0806.4-4123 to be qualitatively similar to other NSs. The code we provide is flexible and can be used for any isolated NS system.

A. Flux predictions
In Fig. 2 we show the radiated power of the reference NS from each pixel of the ROI, which corresponds to a planar surface of 100 km × 100 km with 6277 2 pixels (Δb ¼ 30 m). We checked that this choice of resolution indeed meets our convergence criterion. Moreover, we consider the misalignment angle θ m ¼ 15°, the direction θ ¼ 58.31°and ϕ ¼ 0. We take the axion-photon coupling g aγγ ¼ 1.0 × 10 −12 GeV −1 and the axion mass m a ¼ 0.5 μeV. As can be seen from the plot, there exist very bright pixels corresponding to trajectories for which the conversion region is very close to the NS surface or the crossing with the resonant conversion surface is almost tangential. The former implies that the resonant conversion takes place in regions with very high magnetic fields. The latter, instead, implies that the derivative of the plasma mass along the trajectory is almost zero, therefore strongly enhancing the axion-photon conversion probability. The position of these bright pixels changes during the NS rotation as shown in the video (link is provided in the caption of the figure). In Fig. 3 we report the total radiated power (summing the contribution of all the pixels) as a function of the instantaneous phase during a whole NS period for a few benchmark cases. In particular, the plots from left to right correspond to an increasing axion mass. The solid, dashed, and dotted lines represent the results for a polar angle θ of 36°, 54°, and 72°, while the azimuthal angle is fixed to ϕ ¼ 0. Most importantly, our numerical flux predictions (shown in light blue) are compared with the analytic calculations (shown in black) in which the axion PSD is completely radial [48]. It is clear that taking into account the contribution of all the trajectories (isotropic PSD) has significant implications. First, we find that an isotropic PSD provides an averaged normalization of the radiated power which is typically up to an order of magnitude smaller than the estimated power from a radial PSD. Second, the ray-tracing calculation does not provide a sharp cut-off for the radio signal at large axion masses. The analytic prediction requires the axion-photon conversion to occur outside the NS surface and therefore implies an upper value for the axion mass that can produce a radio signal, as obtained by setting r c ¼ r NS in Eq. (12). This can be seen by the fact that the black lines are nonzero only in a specific time window during the NS period and, remarkably, are absent in the last plot with m a ¼ 12.5 μeV. Our numerical calculation therefore allows one to extend the radio sensitivity curves to larger axion masses. Thirdly, the isotropic PSD erases the majority of the time variability of the radio signal. The light blue lines are indeed almost flat, while the black curves have large peaks for specific instantaneous phases, especially for low axion masses (see first plot). function of the polar angle θ. The three plots correspond to three different misalignment angles θ m , while the values of the axion mass are represented by different colors. We note that the averaged radiated power increases as one considers larger axion masses up to a certain value where it starts to decrease, as can been seen for m a ¼ 12.5 μeV (yellow line). For axion masses larger than a certain threshold, fewer trajectories have a crossing point outside the NS surface and, consequently, cannot contribute to the total radio signal. However, we highlight once again that this is not a sharp change as occurs in the analytic derivation where the radio signal suddenly becomes zero [48]. Figure 5 shows that the relative variance of the radio signal significantly depends on the misalignment angle, reaching 30% for the extreme value θ m ¼ 90°(last plot). For reasonably small values of θ m , the relative variance is instead practically equal to one, implying an almost negligible time variability of the signal. Figure 5 shows that the variability of the signal as a function of time is small for realistic values of the misalignment and viewing angles. As discussed above, the time variation of the signal has been reduced (in the majority of cases) to Oð0-5Þ%. We therefore forecast the sensitivity of radio telescopes to a line detection, although we discuss how well these searches can find the remaining time variability.

B. Radio sensitivity
Assuming the thermal noise of the radio telescope is Gaussian, the signal-to-noise ratio (SNR) is given by where S a→γ is the flux density of the source, SEFD is the system equivalent flux density, B bandwidth of the signal, τ obs the observation time, and the factor of two simply accounts for the number of polarizations. The flux density is given by S a→γ ¼ F=B. The bandwidth is described in Sec. II where we consider two contributions to the broadening of the line. The sensitivity to both scenarios is shown in Fig. 6 where the solid red line corresponds to broadening from the DM velocity distribution only and the red dashed line accounts for Doppler broadening from the rotation of the magnetosphere as well.
Realistically, the radio sensitivity will be limited by a telescopes ability to remove false positives. For example, interference from radio sources on Earth or satellites in the field of view can produce radio line emission. FIG. 6. Projected sensitivity to the axion-photon coupling from radio observations. We consider the isolated NS J0806.4-412 and assume τ obs ¼ 100 hrs. The two red lines correspond to the sensitivity limit for two line broadening scenarios as described in the text. The red solid line only accounts for the DM velocity distribution far from the NS where as the red dashed line also accounts for Doppler broadening from the rotation of the NS magnetosphere. The red band shows the minimum coupling required to detect the time variation of the signal (here we neglect Doppler broadening).
A confirmation of the signal therefore requires an interferometer capable of taking ON/OFF spectra of the source in radio quiet regions. We therefore consider the future telescope SKA which is projected to have SEFD ∼ 0.098 Jy [70] and assume an observation time of τ obs ¼ 100 hrs. Figure 6 shows the 2σ sensitivity 7 to the radio line for various masses. In addition, we show the region of parameter space that could correspond to the QCD axion (blue), the current and future ADMX sensitivity in dark and light grey respectively [18,26], and the CAST sensitivity in orange [24].
Although reduced, the time variability of the radio line would provide a striking confirmation of its astrophysical origin. We therefore estimate the scaling from a 2σ measurement of the line to a 5σ detection of the variability. We therefore make the substitution where σ T denotes the variance of the signal over a NS rotation. From Fig. 5 we see that, for θ m ¼ 18°, we have 0.0102 ≲ σ T ≲ 0.0172 depending on the axion mass. We can then scale line detection sensitivity to the smallest coupling with a detectable time variability as where SNR T is the signal-to-noise ratio for a detection of the time variation which we set to 5. The red band in Fig. 6 shows this scaling for the range of values for σ T . In particular, we can see that a significantly larger coupling (and subsequently larger flux) is required to detect the small variation in the signal. For larger values of θ m we see greater variation in the signal, making it easier to detect, although this effect is still greatly suppressed with respect to Ref. [48].

IV. DISCUSSION AND CONCLUSIONS
In this work we have calculated the expected flux from axion-photon conversion in a NS magnetosphere. To correctly account for the isotropic phase space distribution of the DM at the NS surface we performed a ray-tracing procedure, back propagating photons to the conversion surface to find the total flux at Earth. Our work builds upon Ref. [48] with two primary results: (i) The predicted flux is large enough to potentially detects ALPs down to g aγγ ∼ 10 −12 GeV −1 ; (ii) The time variation of the signal is small for realistic values of the misalignment angle, 0°≲ θ m ≲ 30°, where we see, at maximum, an Oð5%Þ modulation of the signal. Both of these effects can be seen in Fig. 3 in which the blue lines indicate this work and the black lines show the radial approximation of Ref. [48] where the time variation is greatly enhanced. We note however that our analysis currently neglects two important physical effects: reflection and refraction of photons escaping after the conversion process. Both of these effects could change the magnitude and time variability of the final signal. Nevertheless, we have performed the most detailed calculation of the signal to date and show that the flux is still observable for a variety of angular configurations, as shown in Fig. 4. Figure 5 shows the variance of the signal for different viewing and misalignment angles from which it is clear that only for unrealistically large misalignment angles, θ m ≳ 30°, do we see an appreciable modulation of the signal.
In Fig. 6 we show that observations of NS targets, such as J0806.4-4123, with future telescopes such as SKA will probe unexplored regions of the ALP parameter space. Although the time variability of the signal is small compared to previous estimates, a measurement of this variability would provide a striking signature of the signal's astrophysical origin. We therefore estimate the coupling required to make a 5σ detection of the time variation for θ m ¼ 18°and θ ¼ 54°, showing that observations J0806. 4-4123 could still detect this variation for couplings below the CAST limit.
By accounting for the isotropic phase space distribution of DM we are able to extend the sensitivity of SKA to higher axion masses than Ref. [48]. The high mass cut off of the sensitivity is set by the requirement that the conversion process occurs outside the NS interior which, in our setup, occurs at a different mass for each pixel. At high axion masses we therefore retain a fraction of the overall flux induced by nonradial trajectories. This effect is also reflected by the reduction of sensitivity at high masses m a ≳ 10 −5 eV, as seen in Fig. 6.
Although we have made a crucial step toward calculating the true signal, there are a number of caveats that should be addressed in future work. First, we neglect the boost of the NS with respect to the galactic rest frame. In practice this boost would mean that the NS sees a prevailing DM wind, similar to the wind studied in direct detection experiments [71]. Accounting for this effect is relatively simple if the boost is known but we leave this to future work. Second, we assume that the GJ model is a good approximation of the NS magnetosphere which, for realistic NSs, may not be the case. Future work should systematically understand how realistic NS magnetosphere models can effect both the magnitude and width of the radio line. Although this may only be possible with full 3D simulations, a systematic study of different analytic models would provide valuable information for the space of possible signatures. As mentioned in Sec. II C, inhomogeneities in the plasma can significantly affect the signal calculation and should also be accounted for in future work. In addition, the multiplicity of the charge density can vary significantly from the GJ value, potentially changing the conversion region's size and its distance from the NS. Reference [51] recently studied the Doppler broadening of the line due to the motion of the magnetosphere, showing that there is potentially a significant contribution to the overall width. Moving forward, it is important that the precise width and its evolution in time is computed more accurately, accounting for turbulence in the plasma as well as the transmission and reflection components of the signal. In particular, multiply scattered photons have to be taken into account consistently. This is especially important for conversion in the throat of the NS magnetosphere as the photon production can be greatly enhanced in this region but may not escape to infinity without scattering. We leave this to future work. Overall, we have taken a step toward understanding the true signal of axion-photon conversion from a NS-future work will build upon our framework by incorporating many of the physical effects mentioned above and reducing the number of assumptions we made in this work.
Finally, we emphasise the complementarity between indirect and direct searches for axion DM. Given a detection of a radio line from a NS it would be easy to confirm its DM nature through measurements with cavity searches such as ADMX. Importantly, the frequency of the radio signature would allow for an inference of the axion mass and subsequently reduce the frequency range through which direct experiments would need to search (a primary issue in direct searches for axions in resonant cavities). Both direct and indirect approaches therefore represent fundamental tools in the search for axion dark matter. Code used for the calculations throughout this work can be found at https:// github.com/mikaelLEROY/AxionNS_RayTracing.
If we assume that the magnetic field varies only slowly around the conversion region, the integrals can be evaluated analytically and one obtains (in the limit η → ∞) where all quantities are evaluated at the conversion point. Equation (A7) is equivalent to the expressions obtained in [51], and differs from Ref. [48] by an additional factor 1=v a (which enhances to overall emission). Note that throughout we neglect corrections of the order of the Lorentz factor and always approximate γ ≃ 1 (which is accurate to within ≲20% in our scenarios).