Adiabatic Axion-Photon Mixing Near Neutron Stars

One of the promising new proposals to search for axions in astrophysical environments is to look for narrow radio lines produced from the resonant conversion of axion dark matter falling through the magnetospheres of neutron stars. For sufficiently strong magnetic fields, axion masses in the $\mathcal{O}(10\mu{\rm eV)}$ range, and axion-photon couplings $g_{a\gamma} \gtrsim 10^{-12} \, {\rm GeV^{-1}}$, the conversion can become hyper-efficient, allowing axion-photon and photon-axion transitions to occur with $\mathcal{O}(1)$ probabilities. Despite the strong mixing between these particles, the observable radio flux emanating from the magnetosphere is expected to be heavily suppressed -- this is a consequence of the fact that photons sourced by infalling axions have a high probability of converting back into axions before escaping the magnetosphere. In this work, we study the evolution of the axion and photon phase space near the surface of highly magnetized neutron stars in the adiabatic regime, quantifying for the first time the properties of the radio flux that arise at high axion-photon couplings. We show that previous attempts to mimic the scaling in this regime have been overly conservative in their treatment, and that the suppression can be largely circumvented for radio observations targeting neutron star populations.

There are growing experimental efforts across the globe to search for dark matter axions using haloscopes [19], which typically attempt to measure the coupling of axions to photons, given by L aγ = − 1 4 g aγ aF µν F µν , where F µν is the photon field strength tensor, a is the axion field, and g aγ is a dimensionful coupling constant.The most successful approach to date involves constructing a small cavity whose electromagnetic modes can be tuned to match the frequency of the background axion field , however a variety of alternative ideas have also emerged which attempt to overcome the challenges of conventional cavity searches, allowing laboratory experiments to probe a broader range of axion masses and interactions (see e.g.[44][45][46][47][48][49][50]).Another approach is to look for signatures of axions in astrophysical environments (see e.g.[51][52][53] for recent reviews); these techniques are highly complementary to laboratory experiments since they are often capable of probing a wider range of axion masses, rely on different assumptions of the underlying distribution of axion dark matter, and can be used to break intrinsic degeneracies that arise in terrestrial searches.
Axions are generally thought of as feebly-interacting particles, implying that their interactions are, in most contexts, only expected to induce small perturbative effects on the systems of interest.For example, axion dark matter falling through a neutron star magnetosphere is typically expected to pass through the entire system unperturbed, only on rare occasions sourcing low-energy radio photons.Despite being a rare process, however, this signal can shine through astrophysical backgrounds thanks to: (i) the distinctive spectral shape of the radio signal, manifested as an extremely narrow spectral line (which sharply contrasts against smooth astrophysical backgrounds), and (ii) a large local axion number density-potentially exceeding ∼ 10 20 cm −3 -which can compensate for the inefficiency of axion-photon conversion.
Early observational campaigns [60,65,78,79] looking for radio lines produced from axion dark matter derived limits on the axion-photon coupling in this 'perturbative limit', i.e. they worked under the assumption that the axion-photon conversion probability was always small, P a→γ ≪ 1, implying that the radio luminosity scales as L radio ∝ g2 aγ .It was only recently pointed out in Ref. [62] that these assumptions can be strongly violated, particularly at large (but still viable) axion-photon couplings and for pulsars with strong magnetic fields.Instead of occasionally sourcing an on-shell photon, axions falling through the magnetosphere are expected to convert with O(1) probability.The story doesn't end there, however, as the newly produced photons will themselves encounter resonances 1 , converting back to axions with O(1) probability.In the large g aγ limit, the expectation is that photons typically convert back into axions before escaping the magnetosphere, resulting in highly suppressed radio luminosity.
As a first attempt to include these 'adiabatic conversions' into the calculation of the radio flux, Ref. [62] adopted the simplifying assumption that photon production at each resonance could be approximated by using a net effective conversion probability which is set by the product of the survival probability with the conversion probability, i.e.P eff a→γ = (1 − P a→γ ) × P a→γ .This approximation, which leads to an exponential suppression (see Eq. ( 4)) in the large coupling limit, is strictly speaking only valid when the axion-photon resonances take place on a spherical surface centered about the neutron star, and when the conversion probability is equivalent for all axions and all photons, i.e. it depends only on the radial distance from the neutron star.For realistic systems, neither of these assumptions hold, and it remains unclear how well this approximation reflects the true rate of photon production in the adiabatic regime.
In this manuscript, we develop an algorithm capable of carefully tracking the evolution of the axion and photon phase space around neutron stars, and characterize, for the first time, the scaling and properties of the radio flux produced from adiabatic resonant conversion of axion dark matter in the magnetospheres of neutron stars.For large axion-photon couplings and small axion masses, our algorithm recovers the approximate exponential suppression predicted in [62].This suppression, however, is only valid for a range of couplings -instead of being exponentially suppressed in the limit that g aγ → ∞, the radio flux instead asymptotes to a fixed finite value.
That being said, the suppression of the radio flux can be partially avoided at larger axion masses, where there is more room for axions to traverse the magnetosphere in such a way that they encounter only one level crossing (see top row of Fig. 1).These axions contribute to the radio flux, but are limited in number, implying the radio luminosity is phase-space suppressed relative to the naive perturbative approach in which re-conversions are neglected.We show that the observed suppression is crucially dependent on the geometry of the resonant surface around the neutron star, and provide approximate expressions which can be used to extrapolate the functional scaling of the radio luminosity from the non-adiabatic to the adiabatic regime, thereby evading the need for complex numerical analyses in the high coupling limit.This represents an important step in solidifying the limits derived in [62], and in establishing techniques that allow for future radio surveys to probe axions at large couplings.
The structure of this paper is as follows.In Sec.II we provide a general overview of mixing and propagation of axions and photons in magnetized plasmas.In Sec.III we discuss the algorithm that we develop to self-consistently track the worldlines and production probabilities of particles sourced by infalling axions.We implement these techniques in Sec.IV and examine the behavior and scaling of the radio flux for adiabatic axion-photon conversion.In Sec.V we give our conclusions.

II. AXION-PHOTON MIXING AROUND NEUTRON STARS
Axions falling through a neutron star magnetosphere can resonantly mix with low-energy electromagnetic modes when the four-momentum of the photon matches the four-momentum of the axion, i.e. k γ µ = k a µ .In the highly magnetized plasma found in the inner magnetosphere, the only super-luminous electromagnetic mode that can be excited is the Langmuir-O (LO) mode, whose dispersion relation is given by [63,64,80] 2 where ω p = 4παn e /m e is the plasma frequency of the medium3 , k = |k| is the modulus of the photon threemomentum and ϑ k is defined as the angle between the photon momentum and the magnetic field.Eq. ( 1) can be used to solve for the location of the resonances, which occur when [64] which reduces in the non-relativistic limit to ω p ≃ m a .The efficiency of resonant LO mode production has been computed analytically using various approximation schemes (always assuming that the conversion takes place in a very narrow region near the resonance itself, an assumption which is expected to hold to high degree in most contexts) [55,63,64,84], with the most recent calculation producing an axion-photon conversion probabil-ity given by [85] where v p = k/ω is the phase velocity of the photon at the resonance.This expression is only valid in the nonadiabatic limit (P a→γ ≪ 1), but is expected to generalize in the adiabatic limit (P a→γ ∼ 1) to the Landau-Zener formula [59] 4 where γ = P non-ad a→γ is the adiabaticity parameter.Once produced, photons are refracted away from the neutron star by the dense plasma.Owing to the highly non-linear trajectories of photons in this media, understanding the evolution and fate of the newly sourced photons requires dedicated ray tracing simulations (see e.g.[63,80,84]), which amounts to solving Hamilton's equations, given by where λ is the wordline of the photon, and the photon Hamiltonian in a magnetized plasma is given by Here, we have introduced where the '•' notation represents a contraction over the spatial indices, and the spatial dependence is understood to be implicitly embedded in ω p , the spacetime metric g µν and the magnetic 4-vector field B µ .Note that Eqns.( 5)-( 6) can also be used to solve for axion trajectories, but using the simpler Hamiltonian given by H a (x µ , k µ ) = g µν k µ k ν − m 2 a .When computing the evolution of axion and photon trajectories, we use the Schwarzschild metric, taking a characteristic neutron star mass M NS = 1 M ⊙ .For axions traversing the neutron star itself, we switch to the interior Schwarzschild metric [90] (which assumes a constant density on the interior of the star), adopting in this case a neutron star radius R NS = 10 km.These values are merely intended as 'ball-park' estimates, with measured systems suggesting typical neutron star masses closer to M NS ∼ 1.4 M ⊙ and R NS ∼ 10 − 14 km (see e.g.[91,92]); the impact of varying these parameters in the non-adiabatic limit has recently been discussed in [80], and alternative choices are not expected to qualitatively alter any of the conclusions drawn based on the rough estimates used here.
The main purpose of this paper is to study the nontrivial evolution of the axion and photon phase space in the adiabatic limit.This is accomplished by: (1) following infalling axion dark matter particles as they fall through the magnetosphere, (2) identifying all resonance points R i encountered during the infall, (3) assigning a phase space factor ξ i which accounts for the initial number density of infalling axions, the survival probability of the axion to reach resonance R i , and the probability of photon production (Eq.( 4)), and (4) iteratively repeating steps (2-4) with the newly produced photons until all possible resonances have been identified and all axion and photon trajectories have been traced to asymptotic distances.Owning to the large number of resonances that can be encountered, this procedure can become quite complex -as illustrated in Fig. 1, a single infalling axion can lead to anywhere between O(few) and O( 103 ) possible outgoing axions and photons.The radio signal can be computed by summing over the asymptotic position of photon trajectories, localized in some region in the sky, where each photon trajectory is appropriately weighted by the initial axions phase space and the probability that it was produced and survived [80].
In order to make concrete quantitative statements about the behavior of the radio flux in the adiabatic regime we adopt a fiducial model for the magnetosphere characterized by a dipolar magnetic field and a fully charge-separated Goldreich-Julian (GJ) charge density, which can be derived by searching for the minimal corotation charge density necessary to screen E • B in the magnetosphere.The GJ charge density yields a plasma frequency near the neutron star of [93] where Ω is the rotational frequency vector, m is the unit vector in the direction of the rotating magnetic dipole and ϑ m is angle between the two.The factor m • r = cos ϑ m cos ϑ + sin ϑ m sin ϑ cos(Ω t) encodes the angular factor between the magnetic axis and the vector r.The surface magnetic field strength is denoted B 0 .The fully charge-separated GJ model predicts small regions of vacuum, located at angles ϑ null , at the boundary of the charge separated regions 5 .While full-charge separation is expected to appear in dead neutron stars (see e.g.[58,94,95]), it is unclear the extent to which these features survive for more active pulsars.As we will show in the following sections, the presence of vacuum regions extending to the neutron star surface can play an important role in controlling the efficiency of radio production at high couplings.In order to make conservative statements about the adiabatic regime, we thus also consider the inclusion of a small 'boundary layer' of plasma around the neutron star, constructed in such a way that the large-scale features of the conversion surface (and plasma distribution) are unaltered, but the regions of vacuum near the neutron star are partially filled with a plasma density comparable to what is found at angles ϑ ∼ 0 and ϑ = π/2.Specifically, in the case of an active pulsar, we consider an additive boundary layer contribution ω p,BL to the plasma density of the form where ω p,0 is the GJ plasma frequency at the pole, r b = 0.3 × R max and δr = 0.1 × R max , where R max is the maximal radial extent of the conversion surface.The coefficients in r b and δ r have been chosen in such a manner that the filling of the null lines is significant, but the plasma at R max remains nearly unmodified.Other functional forms could also be adopted to perform the same function, however the conclusions will not be significantly altered so long as the large scale features of the conversion surface remain unaltered.
In the following we perform our analysis using both the GJ model, and the GJ + boundary layer -the former should be understood to be representative of a dead pulsar, and the latter as a conservative treatment of a more active pulsar 6 .

III. TRACING THE PHASE SPACE EVOLUTION
Here, we extend the forward ray tracing algorithm developed in [62,63,80] to self-consistently include the complete evolutionary tracks of all axions and photons sourced near the neutron star.The details of this algorithm are outlined below.
We begin by applying the Monte Carlo (MC) surface sampling algorithm developed in [63,80] to draw uniform samples from the resonant conversion surface, as defined in Eq. ( 2).An axion trajectory is initially backward propagated away from this initial condition R 0 = (x 0 , k 0 ) to an asymptotic distance, and a phase space factor ξ a i at each resonant point R i is recorded 7 , as indicated in Fig. 2.This factor accounts for the 6 Technically, the GJ charge distribution is only expected to be representative in the closed magnetic field lines, however the open field lines are volumetrically tiny in the region of interest and are thus not expected to play any important role in the evolution of these systems. 7In practice, the sampling scheme applies Liouville's theorem to relate the phase space at infinity to the phase space on the conversion surface.This is, however, not applicable in the strong coupling limit.Therefore, the final weight of the event has to be re-weighted by the probability that the infalling axion indeed survives its travel to the sampled conversion points.
asymptotic axion energy density which sourced the initial infalling trajectory, the effect of gravitational focusing, and the probability that the infalling trajectory leads to a photon at R 0 .In effect, this amounts to ξ a 0 ≡ 2n a,∞ / √ π (v 0 /v ∞ ) P 0 a→γ P i>0 a→a , where v 0 is the axion velocity at ξ 0 , n a,∞ and v ∞ are respectively the asymptotic axion number density and velocity, and P i>0 a→a represents the cumulative probability that the infalling axion is still an axion by the time it has reached the resonance R 0 8 .Starting from R 0 , we trace a photon trajectory out to infinity.Since the photon may hit one or more conversion surfaces [see e.g.Fig. 1] -potentially re-converting into an axion -all resonances encountered along its path must be tracked and assigned a weight ξ γ i .In turn, any newly sourced axions must also be propagated to infinity, and any conversion surfaces they encounter be assigned a weight ξ a i .This procedure is iterated until all possible resonances stemming for the original primary particle have been identified.Depending on the axion trajectory and the characteristic geometry of the conversion surface, each infalling axion trajectory can lead to anywhere from O(few) to O( 103 ) outgoing trajectories, making this a numerically challenging procedure.
The final radio flux at a given point on the sky (ϑ, ϕ) can be obtained by taking the collection of outgoing photon trajectories which end up within a small angular bin on the sky (at asymptotic distances), and summing over the weighted contributions of each of these photons, i.e. the power radiated in a region on the sky (ϑ 0 ± ε, ϕ 0 ± ε) is given by where N s are the number of samples drawn, and we have defined the photon weight function W i , which in the sampling procedure of [63,80] is given by9 Here, we have defined cos ϑ k∇E as the angle between k and ∇E γ , the pull-back metric h k which defines the momentum dependent conversion surface, the maximal number of resonant crossings per sample N max and the maximal radial distance used in the surface sampling algorithm R max , and the local axion number density n a,loc .The full tree, i.e. all possible outcomes, of the path of the outgoing photon should be considered, since even small axion-photon conversions may lead to detectable radio signals.The computational time for the full tree, however, is naively expected to scale as ∼ 2N res +1, where N res is the number of resonances encountered 10 , making this procedure significantly more computationally intensive than in the non-adiabatic limit 11 .In order to avoid 10 Note that the number of out-going trajectories Nout scales with the generation Ng as Nout = 2 Ng , and the number of resonances scales like Nres = 2 Ng −1 .These relations can be used to compute the total number of trajectories (including internal legs), which is given by Ntot = 1 + i≤Ng N out,i = 1 + 2 i≤Ng N res,i = 2Nres + 1 (where the factor of one comes from the initial trajectory). 11Note that the most computationally expensive part is the highly severe computational time for complicated trajectories (see e.g.bottom right panel of Fig. 1), we transition to a pure MC sampling after N = 5 conversion points have been encountered, implying that we consider (at most) 6 outgoing particles 12 -thus making sure that we include at least one outgoing photon in the event.In practice, this is achieved by always considering the branch with the highest weight, i.e. a particle is only propagated until it reaches a resonance, the outcomes stored in a pool, and the particle in the pool with the largest weight at any given time is propagated.
In most cases, the tree of possible outcomes is rather simplistic.For example, the average number of resonances encountered, N res considered in the next section are 3.1, 2.3, and 1.7, for the masses m a = 1.0×10 −5 , 1.0× 10 −5 and 2.6 × 10 −5 eV, respectively, in the GJ magnetosphere 13 .Despite many trajectories being simple, the MC selection process was triggered (i.e. more than 5 resonances encountered in the tree) in 13.3%, 6.2% and 0.1% of the trees, contributing to as much as 26.0%, 11.1% and 0.1% of the total flux at large couplings.

IV. RESULTS
Using the algorithm discussed above, we analyze the radio signal emanating from a neutron star with a dipolar magnetic field with surface strength B 0 = 10 14 G, a rotational period P = 2πs, and a misalignment angle ϑ m , which we set to zero for simplicity 14 .Owing to computational costs, we choose to keep these three parameters fixed throughout the analysis, varying only the axion mass, m a , and axion-photon coupling, g aγ .As we will show below, the scaling of the radio flux into the adiabatic regime is largely set by the geometry of the conversion surface -since shifting B 0 and P alter the geometry in a manner that is fully degenerate with a shift in the axion mass, and the role of ϑ m is at leading order to induce a small rotation in the conversion surface, we non-linear propagation near the resonances, implying that the computational time scales with the number of sub-branches in the tree, 2Nres + 1. 12 In addition, we include two stopping criteria to hinder potential rare semi-stable or complicated trajectories.First, we truncate the MC simulation if 50 resonances are encountered; such events are rare, occurring in only 13 of the ∼ 10 7 events included in the analysis.Second, the simulation is truncated when the simulated outcomes account for more than 1 − 10 −100 .This second threshold is overly conservative, and can be relaxed significantly to further reduce computation time. 13A number 1 means that only a single photon is forward propagated and a single axion backwards propagated as in Fig. 2. 14 We also take M NS = 1 M ⊙ and r NS = 10 km.The neutron star mass and radius are not expected to qualitatively change our conclusions (with the predominant effect being O(1) shifts in the overall flux [80]).Non-zero misalignment, on the other hand, has the dominant effect of smoothing out the differential power across a wider range of viewing angles (see e.g.[63]).
FIG. 3. Differential radiated power from dark matter axions converting into photons in the magnetosphere of a neutron star for axion masses ma = 1.0 × 10 −5 eV (first row), ma = 1.5 × 10 −5 eV (middle row) and ma = 2.6 × 10 −5 eV (bottom row), and for axion-photon couplings from gaγ = 3 × 10 −10 to gaγ = 10 −12 (see legend).The results are shown using the newly developed algorithm to compute the conversion trees from each infalling axion using either the GJ magnetosphere (left column), or the GJ magnetosphere with an additional boundary layer as described in Sec.II (middle column).In addition, the results with the conservative approximation Pa→γ = e −γ (1 − e −γ ) are shown in the right column for comparison.
believe the results identified here are quite general 15 .In Fig. 3, we plot the period-averaged differential power as viewed from an angle ϑ with respect to the axis of rotation, for three different axion masses and six choices of the axion-photon coupling (which smoothly extend the results from the non-adiabatic to adiabatic 15 The only notable subtlety is that the axion-photon coupling at which the adiabatic regime is encountered can shift to smaller or larger values, depending on the magnetic field strength and the plasma density in the magnetosphere.For this reason, our results should be interpreted qualitatively rather than quantitatively. regimes).Each of the differential power curves have been re-scaled by a factor of (g aγ /10 −10 GeV −1 ) −2 , so that the suppression of the power in the adiabatic regime can be more easily identified (in the non-adiabatic regime, this re-scaling causes all curves to lie on top of one another).We illustrate the evolution of the differential flux at large couplings using three distinct approaches.In the right column, we adopt the approximation scheme of [62], which amounts to assigning each photon an effective conversion probability P a→γ,eff = e −γ (1 − e −γ ).The left column, instead, shows a comparison with the full conversion tree as computed using the algorithm described in the preceding section.Finally, in the center column, In addition, we plot the results that would be obtained by assuming photons cannot re-covert to axions once produced -this is obtained by taking Pa→γ = 1 − e −γ , and effectively sets an upper limit on the flux (green dotted line).The small reduction in the green line at high couplings arises from photons which impact the NS surface (and are thus lost).we compute the full conversion tree including the 'boundary layer' contribution to the plasma frequency discussed in Sec.II.
A number of features can be readily appreciated from Fig. 3. First, the naive approximation scheme of [62] tends to consistently over suppress the flux in the adiabatic regime.Next, the suppression is largely, but not entirely, uniform across the sky -for small axion masses, the suppression is more apparent near the magnetic poles, but away from the poles the suppression is more uniform (note that in the case of misaligned rotators, this effect would be smeared across viewing angles).In addition, the suppression appears to be much more prominent for small axion masses, which corresponds to the scenario where the resonant conversion surface extends further from the neutron star surface (see Fig. 1).Finally, the existence of a small boundary layer of plasma around the neutron star tends to suppress the flux relative to the GJ magnetosphere, but not as much as the effective scheme adopted in [62].
These features can also be appreciated by looking at the sky-averaged flux as a function of the axion-photon coupling; Fig. 4 compares each of the three models for all three axion masses.Here, one can see that the radio flux is actually expected to plateau at sufficiently large axion-photon couplings, rather than become exponentially suppressed.The relative height of the plateau depends both on the axion mass and the on the existence of charge-separation in the magnetosphere.
Collectively, Figs. 3 and 4 lead to two significant conclusions: • Despite the fact particles, on average, encounter an even number of level crossings, the efficiency of these level crossings is not equivalent.As such, the approach of [62] naturally over-estimates the suppression of the radio flux in the adiabatic regime.
• Ref. [62] missed the importance of infalling axion trajectories which only encounter a single resonance, which despite often being uncommon can dominate the radio flux.Such trajectories occur when axions traverse the neutron star, entering or exiting near the charge separation boundary.For large axion masses, a larger fraction of the neutron star surface is 'exposed', allowing for a larger fraction of the infalling axion phase space to encounter single level crossings.

A. On the application to future searches
The computational cost of running the full conversion tree makes it difficult to fully embed within a more sophisticated analysis of radio data, such as the analysis performed in [62].As such, we attempt to develop below an approximate re-scaling technique that can be adopted in future work to approximate the suppression in the flux that arises in the adiabatic regime.
The simplest approximation that one can make to account for the transition from the non-adiabatic to the adiabatic regime, is that of [62], P a→γ = e −γ (1 − e −γ ).However, as we have shown, this approach is overly conservative in the adiabatic limit (cf.Fig. 4), and can lead to a significant underestimation of the total radio flux.On the other hand, the most optimistic approximation one can make is P a→γ = 1 − e −γ , i.e.only a single level crossing.
An alternative approach is to try and encode the suppression of the radio flux of each neutron star into an effective re-scaling parameter, which depends 16 on R max . 16The global suppression factor reflects the overall suppression for These suppression factors can then be included as 'effective' conversion probabilities which are turned on as the neutron star enter the adiabatic regime, P a→γ ≳ 0.2 [cf.Fig. 5].Using the data points in our sample we derive approximate re-scaling factors S = Φ(g aγ → ∞)/Φ(P a→γ = 1) for the sky-integrated flux -the suppression factors are shown as a function of R max in Fig. 6.We implement these suppression factors by adopting the adiabatic approximation at low couplings, P a→γ = e −γ (1 − e −γ ), and transitioning to our effective re-scaling P a→γ = S(R max )×(1−e −γ ) at larger values of g aγ -in the intermediate regime, we adopt the maximum of the two approaches.We illustrate the relative agreement of applying this approach to the sky-integrated flux in Fig. 7.

V. CONCLUSION
In this work, we have for the first time identified the behavior and scaling of radio emission produced from the resonant conversion of axion dark matter near neutron stars in the adiabatic (i.e.strong mixing) limit.We have done this by developing an MC sampling and ray tracing algorithm capable of carefully tracking the evolution of the axion and photon phase space.
Our results clearly indicate: (i) contrary to previous approximations, the radio flux is not exponentially suppressed at large axion-photon couplings, but rather plateaus to a fixed value, and (ii) the radio flux is not suppressed at all for small conversion surfaces (i.e.R max ∼ R NS ), which would be the case if there exist dead neutron stars in the field of view which support a maximal plasma density only slightly in excess of the axion mass.We further illustrate an approximate scaling relation which can be used to extrapolate radio observations into the adiabatic regime, circumventing the need to apply computationally expensive simulations -such as the one developed here -to large numbers of systems.
Our conclusions are based on a number of assumptions, most of which are thought to be well-justified; for the sake of clarity, we enumerate these assumptions below: 1. Axion-photon transitions are dominated by the resonant contribution, and it is valid to treat the resonance with the WKB approximation (i.e. that the background varies slowly relative to the axion wavelength).For axion masses capable of generating radio emission, this approximation is expected to be true over a majority of the magnetosphere, with the one exception perhaps being regions near the return currents and open field lines.
most viewing angles, with the exception near the poles, as can be deduced from Fig. 3.However, for slight misalignment angles, the in-homogeneity of the suppression is expected to be largely washed out.2. The adiabatic generalization of the non-adiabatic conversion probability is assumed to follow the Landau-Zener formula.This has been shown to be true in one-dimension and in an isotropic plasma (at least when the medium is smoothly and slowly varying) (see e.g.[59]), but has not be explicitly derived for a three dimensional anisotropic plasma.
3. Axions are assumed to be fully non-interacting away from the resonance, and the axion population is assumed to arise exclusively from either infalling axion dark matter, or from axions sourced from photons which themselves were sourced from axion dark matter (that is to say, local radiation from the magnetosphere is neglected).
4. No other exotic particle content is assumed to exist which could either alter the dispersion relations of these particles, or the assumption that their interactions with the ambient medium can be neglected.
5. The magnetosphere is assumed to be approximately characterized by a purely dipolar field and the Goldreich-Julian charge density.Higher-order magnetic multi-poles may exist near the star, but are not expected to have a large qualitatively impact on our results.For standard pulsars, deviations from the GJ model are expected along the open field lines and near the return currents, but the closed field lines (comprising nearly all of the near-field magnetosphere) are roughly expected to be well-characterized by the GJ values.The charge distribution may differ notably for millisecond pulsars, magnetars, and binary pulsar systems; however, given a model of any of these systems, the formalism developed here can be applied to make quantitative statements about each of these.
This work has important implications for the radio searches for axion dark matter (such as those performed in [60,62,65]), and will prove important as these observations are extended to other systems and to a broader range of frequencies.

FIG. 1 .
FIG.1.Illustration of possible final states arising from a single infalling axion (initial trajectory marked with red arrow).The infalling axion encounters resonances when ωp ≃ ma (green surface), potentially converting to a photon (orange); the photon in turn may encounter resonances during its propagation, potentially reconverting back into an axion (blue).This processes iterates until all possible outgoing trajectories have been identified.The collection of these final states, and their trajectories through the magnetosphere, are illustrated for two choices of axion mass (the upper panel illustrating the case of ma = 26 µeV, and the lower panels illustrating the case of ma = 10 µeV) and two differential sets of initial conditions.The top right panel illustrates the 'single level crossing' scenario discussed in Sec.IV, and the bottom right panel illustrates the potential complexity that can arise in these systems.

FIG. 2 .
FIG. 2. Similar to Fig. 1, but highlighting instead the implementation of the numerical sampling procedure, which only traces subsets of the full trees illustrated in Fig. 1.Here, an infalling axion is shown to cross several resonances (red points, labeled by Ri>0) before converting to a photon at the MC sampled conversion point, R0.The final rate is re-weighted by the probability that the axion survives the resonances and converts into a photon at R0, ξ γ 0 = 2na,∞/ √ π P 0 a→γ

35 P
FIG.4.Integrated radiated power from dark matter axions converting into photons in the magnetosphere of a neutron star for axion masses ma = 1.0 × 10 −5 eV (left plot), ma = 1.5 × 10 −5 eV (middle) and ma = 2.6 × 10 −5 eV (right).The results are shown for the conversion tree calculation with (blue solid line) and without (yellow dashed line) the boundary layer discussed in Sec.II.For comparison, the results with the conservative approximation Pa→γ = e −γ (1 − e −γ ) are plotted (red dashed dotted line).In addition, we plot the results that would be obtained by assuming photons cannot re-covert to axions once produced -this is obtained by taking Pa→γ = 1 − e −γ , and effectively sets an upper limit on the flux (green dotted line).The small reduction in the green line at high couplings arises from photons which impact the NS surface (and are thus lost).
FIG. 5. Distribution of the axion-photon conversion probability at the MC sample points x0.The distribution is shown for the three masses ma = 1.0 × 10 −5 eV (blue), ma = 1.5 × 10 −5 eV (yellow), and ma = 2.6 × 10 −5 eV (red).The solid lines indicate the GJ model, while the dashed includes the additional boundary layer.
FIG.7.Same as Fig.4, but including the re-scalings discussed in the main text.