Axion-photon conversion in strongly magnetised plasmas

Axion dark matter can resonantly convert to photons in the magnetosphere of neutron stars, possibly giving rise to radio signals observable on Earth. This method for the indirect detection of axion dark matter has recently received significant attention in the literature. The calculation of the radio signal is complicated by a number of effects; most importantly, the gravitational infall of the axions onto the neutron star accelerates them to semi-relativistic speed, and the neutron star magnetosphere is highly anisotropic. Both of these factors complicate the calculation of the conversion of axions to photons. In this work, we present the first fully three-dimensional calculation of the axion-photon conversion in highly magnetised anisotropic media. Depending on the axion trajectory, this calculation leads to orders-of-magnitude differences in the conversion compared to the simplified one-dimensional calculation used so far in the literature, altering the directionality of the produced photons. Our results will have important implications for the radio signal one would observe in a telescope.


Introduction
Axions comprise a class of hypothetical particles that were first conceived as a consequence of the 'Peccei-Quinn' solution to the strong CP-problem in particle physics [1][2][3][4]. Axions are naturally produced in the early universe, and provide an increasingly popular candidate for explaining dark matter [5][6][7][8][9][10][11][12]. 1 Testing the axion dark matter hypothesis is a central goal of contemporary fundamental physics.
Any search for axion dark matter relies on the transfer of energy from the astrophysical axions into Standard Model particles. One promising route to achieve this is through the axion-photon coupling, which receives 'universal' contributions from axion-pion mixing, thus providing a rather robust and model-independent target. Most experimental and observational ideas to search for axion dark matter are based on the axion-photon coupling.
Astronomical searches for axion dark matter via the axion-photon coupling seek to leverage strong astrophysical magnetic fields to effectively transfer energy from the axion field into photons (or vice versa). Particularly interesting are environments where plasma effects generate an effective mass for the photon that allows the axion and photon dispersion relations to intersect. In such a degenerate situation, dark matter axions can resonantly convert into photons, enhancing the photon signal and possibly leading to signals detectable with telescopes on Earth.
Neutron stars are ideal astrophysical environments for converting dark matter axions into photons [13,14]: Neutron stars can host extremely strong magnetic fields (up to ∼ 10 15 G), and are surrounded by a plasma that admits resonant axion-photon conversion near the surface of the star. Due to the plasma densities in the neutron star magnetospheres, resonant conversion is most promising for axions in the neV-meV mass range, including the perhaps best motivated µeV-meV mass range for QCD axion dark matter. The corresponding photon signal would be a narrow spectral line in the MHz-THz (radio) frequency range, with its frequency set by the axion mass. Radio telescopes combine excellent spatial and frequency resolution with sensitivity to very weak radio signals. For these reasons, axionphoton conversion in neutron star magnetospheres has recently become a very active topic of research, and significant effort has been devoted to characterising the possible signal [15][16][17][18][19][20][21][22][23][24], and searching for it in radio data [25][26][27][28].
In order to calculate the radio signal from the conversion of dark matter axions into photons in the magnetosphere of neutron stars, one must model a number of effects. First, the phase space of the axion dark matter evolves under the gravitational influence of the neutron star. Second, in the region where the effective photon mass approximately matches the axion mass, one needs to solve the system of equations describing the axion-electrodynamics to compute axion-photon conversion. Finally, one must propagate the electromagnetic modes excited by the axion through the inhomogeneous magnetosphere to the observer far away from the neutron star. Regarding the first effect, previous work has either assumed that axions are on purely radial orbits (in the frame on the neutron star) as they cross the zone of resonant conversion, or that the axion velocity distribution is isotropic at the conversion zone; neither of these assumptions is correct. We discuss how to compute the appropriate axion phasespace density under the influence of the gravitational field of the neutron star for a neutron star with speed relative to the rest frame of the dark matter distribution comparable to the dark matter velocity dispersion, as one expects for a realistic neutron star. The main work in this paper is to present a calculation of the resonant axion-photon conversion in the neutron star magnetosphere (i.e., the second effect above) which fully accounts for the anisotropic plasma density in the magnetosphere. This fully three-dimensional calculation improves over the simplified one-dimensional calculation previously used in the literature, and we will discuss it in more detail in the following paragraph. We do not perform any calculation of the third effect, the propagation of the photons from the conversion zone to the observer at infinity, in this work. This effect can be dealt with by 'ray-tracing' calculations, following the electromagnetic modes through the anisotropic plasma of the neutron star magnetosphere, see Refs. [20,23,24].
Let us now discuss the core of this work, the calculation of resonant axion-photon conversion in highly magnetised anisotropic plasmas, e.g., in the magnetosphere of a neutron star. A central assumption of essentially all previous studies of this problem is the use of an effectively one-dimensional formalism for calculating the conversion rate of axions into photons, based on the well known calculation for transverse modes [29]. 2 The one-dimensional, linearised system can be expressed as a 'Schrödinger-like' equation, with time replaced by the spatial coordinate along the direction of propagation. The 'conversion probability' of axions into photons calculated from this equation can then be interpreted as the flux transfer: i.e. the ratio of energy densities of the transverse mode of the electric field and the axion. However, in this paper we show that the intuition built around transverse modes breaks down in highly aniosotropic media, such as the magnetosphere of a neutron star. The anisotropy causes a strong mixing between different photon polarisations, so that propagating states are no longer purely transverse or longitudinal. 3 It follows from this that the naive 'conversion probability' no longer describes the physically interesting flux transfer.
We revisit the classical problem of axion-photon mixing in anisotropic media, including relativistic plasma effects, and derive a new 'Schrödinger-like' equation governing the evolution of the Langmuir-O (LO) mode, which is the propagating mode in the plasma that is excited by axion-photon conversion. Our calculation has several important implications. First, the anisotropy of the neutron star magnetosphere results in the photon mode evolving along a different direction than that given by its momentum. This invalidates the onedimensional analysis used in previous work, and can lead to orders-of-magnitude change in the axion-photon mixing for certain axion trajectories. Second, the 'Schrödinger-like equation' now includes damping/growth terms, that can invalidate the unitary evolution of the system. Third, the energy stored in the LO modes is not given by the transverse portion of the field only, and the full Hamiltonian must be used to calculate the flux transfer. Consequently, the flux transfer is no longer identical to the 'conversion probability' calculated from the Schrödinger-like equation, so that a naive over-reliance on the quantum mechanical analogue would lead to incorrect conclusions about axion-photon mixing.
The remainder of this paper is organised as follows. We begin by discussing axionphoton conversion in an anistropic plasma and a magnetic field in section 2, and derive the new form of the 'Schrödinger-like' equation. In section 3, we find the relation between the 'axion-photon conversion probability' and the flux transfer. In section 4 we discuss the calculation of the axion phase space in the neutron star magnetosphere for an neutron star moving with arbitrary speed relative to the rest frame of the dark matter distribution. In section 5 we present numerical results for our conversion calculation and compare them to the results one would obtain in the simpler one-dimensional calculation. As we will see, the conversion probability can differ by orders of magnitude between our calculation that fully accounts for the anisotropy of the plasma and the simplified one-dimensional calculation. We conclude in section 6. Further technical details of our calculation of axion-photon conversion can be found in the appendices.

Axion-photon conversion in aniostropic plasmas
We consider dark matter axions that are gravitationally accelerated towards a neutron star and may convert to photons in its magnetosphere. Due to the large gravitational field of the neutron star, the axions will generally be mildly relativistic, a particularly complicated regime from the point of view of calculating the axion-photon conversion. In such a regime, medium effects cannot be ignored, however the axion also has significant momentum and so the conversion has a strong directional dependence. Our calculation will focus on the regime where the WKB approximation can be applied, i.e. where the axion has momentum larger than the first derivatives of the E-fields. As long as this condition holds, our calculation will apply regardless of the relativistic (or not) nature of the axion. Axion-photon mixing is maximised in regions satisfying the resonance condition, which is approximately given when the plasma frequency is equal to the axion mass ω p m a . Since this condition is typically fulfilled only in small regions, we solve the mixing problem in the flat-space approximation, neglecting gravity. This approximation is expected to be accurate for most (but not all) possible trajectories of infalling axions.
The classical theory of axions coupled to electromagnetism is governed by the Klein-Gordon equation and Maxwell's equations, modified to include the axion-photon mixing. Following common practice (see, for example, Refs. [14,30,31]) we can linearise the equations of motion in a background with the magnetic field B NS , and restrict (without loss of generality) to electromagnetic modes that oscillate with frequency ω: where we have set the relative permeability to µ = 1. The electric field is denoted by E, and the displacement field by D. We would now like to solve these equations for dark matter axions traversing a neutron star magnetosphere. We must first express the equations (2.1) and (2.2) in a convenient coordinate system for any given axion. Two basis vectors can be defined respectively by the direction of motion, and the external magnetic field, which generically are linearly independent. We denote the axion's direction of motion byẑ, so that k = kẑ. This will also be the initial direction of any produced photons, though the changing medium will result in the photon momentum evolving over longer distance scales. We take B NS to lie in the (y, z)-plane at an angle θ(z) with the positive z-axis. In the following, we assume θ = {0, π} unless otherwise specified. This definesŷ(z) and consequently alsox(z) byx(z) =ŷ(z) ×ẑ. We will see that the conversion generally occurs over a small region, and so we will neglect gravity and assume that the axion travels along a straight line and d dzẑ = 0. It then follows that d dzx ∝ŷ and d dzŷ ∝x. The expressions for the differential operators ∇ 2 and ∇(∇ · E) in the (x,ŷ,ẑ) basis are then identical to those in a fixed Cartesian coordinate system. In general, θ varies due to gravitational bending, the resulting photon refracting in the medium and the fact that the magnetic field itself changes. Again as the region where axion-photon conversion occurs should be small, we take θ to be almost constant. The electric displacement field for a highly magnetised plasma is given by the constitutive relation [32] where ij is the dielectric tensor with the plasma frequency ω p . One can see that the anisotropy of the dielectric tensor mixes E y and E z , so one is not in general left with purely transverse of longitudinal propagating states. Here, we have assumed that the particles forming the plasma are moving along the field lines with speed v || in the neutron star frame, and We have also defined where F i || is the distribution function of the species i in the plasma. Note that for an ultrarelativistic plasma, γ −3 γ −1 [33]. In Refs. [14,19,20], the electrons in the plasma were treated as non-relativistic, i.e. γ = 1. While we allow for relativistic electrons, we note that infalling axions will be at most mildly relativistic. For mildly relativistic axions, one would expectω ω, however, in principle the distinction betweenω and ω can lead to differences for high momenta axions traveling through relativistic plasmas. For dark matter falling into a neutron star, this will simply change the effective ω p by a factor ofω/ω which is trivial compared to the significant modelling uncertainties of neutron stars (discussed in section 5.1).
When the axions are decoupled from the system, one can solve for the propagating electromagnetic states of the medium. Note that such states only apply locally: as one moves through the medium, the eigenstates deform adiabatically. As noted in [23], the relevant electromagnetic mode (i.e. the mode excited by axion-photon conversions) is the so-called Langmuir-O (LO) mode, whose polarization has both longitudinal and transverse components. As we discuss in Appendix A, in principle, the axion can also excite purely longitudinal modes. However, such purely longitudinal modes are non-propagating and are hence irrelevant for the electromagnetic signal observed far from the neutron star. The LO mode, on the other hand, evolves adiabatically into a purely transverse vacuum mode as it exits the neutron star atmosphere [33]. The specific form of the LO mode is found in Appendix A to beÊ where we have defined an effective relativistic plasma frequencyω 2 p = ω 2 ω 2 p γ −3 /ω 2 . In the non-relativistic limit, as has been mostly considered in the literature, there is no distinction betweenω p and ω p . Our task is to calculate the conversion of axions into the LO mode. While the modification to the dispersion relation was noted in Ref. [23], the effects of the mixing of transverse and longitudinal modes on the conversion calculation and the implications for how one defines a photon were not explored. As we will see, such effects can change the conversion probability for a given axion trajectory by orders of magnitude.
We can now write out the x, y, z components of equation (2.2). We will be considering the case of axion photon mixing in a slowly varying plasma, so we can neglect second order derivatives that do not involve the z-direction. In addition, E x does not directly couple to the axion as, by constructionx⊥B NS . While E x can be generated by derivative terms, it is suppressed when compared to E y,z so we can also neglect it. At lowest order this leaves the following coupled equations: ∂ 2 E y ∂z∂y (ω 2 −ω 2 p cos 2 θ)E z +ω 2 p cos θ sin θE y + ω 2 g aγ aB NS cos θ . Previous studies of axion-photon conversion near neutron stars have all neglected the lefthand side of these equations, featuring mixed y-and z-derivatives. 4 Neglecting mixed deriva-tives ignores the mixing between E y and E z and is equivalent to assuming that E y only evolves along z. This does not hold for general trajectories, and is only a good approximation if θ π/2, i.e., when the axion happens to be traveling perpendicularly to B NS . To simplify matters, we can eliminate E z from the equations. This will allow us to solve directly for the propagating mode, and neglect non-propagating contributions to E z , which are discussed in Appendix A. By substituting equation (2.8) into (2.7) and neglecting any higher order derivatives that follow, we end up with a single partial differential equation for E y and a: We neglect also the derivatives of the axion driving term, as they will only have a small impact on the evolution of E y .
To attempt an analytic solution, we assume that the change inω p , and so k, is small and thus the direction of propagation stays the same. Of course, this is not true in general: photon trajectories will bend in a changing medium. Assuming small changes inω p , k necessarily restricts us to solving over some finite region where the photon path is approximately constant. The region over which the majority of the generation of an axion-induced E-field occurs is referred to as the conversion zone. While one can estimate the radius of curvature for the changing direction of the photon, such as was done in Ref. [23], our focus in this work is not on the propagation of the converted photons. As most of the generation of the propagating LO mode via the axion should be localised to a smaller region, separating propagation and conversion allows for an analytic approximation of the fields generated by the axion, which could then be propagated numerically via ray tracing code such as in Refs. [23,24].
As long as we are treating the medium as slowly varying, we can write the fields as an envelope term, denoted by a tilde, multiplying a plane wave, To move further, we will take the WKB approximation and assume that kẼ y (y, z) ∂Ẽ y (y, z)/∂y, ∂Ẽ y (y, z)/∂z ∂ 2Ẽ y (y, z)/∂z 2 . As we will see, the smaller the derivatives the longer the length scales over which axions can resonantly convert, leading to higher conversion probabilities. Thus a breakdown of the WKB approximation, which occurs when the first derivatives of theω p vanish, corresponds to a region of enhanced axion-photon conversion. The WKB approximation allows us to then simplify Eq. (2.9) to where we have used the on-shell relation for the axion, m 2 a = ω 2 − k 2 , and defined (2.12) The damping/growth term D appears as an artifact of assuming that the photon does not experience any curvature, i.e., that it is well described as only having momentum in the z-direction, and is given by The back-reaction on to the axion field should be small, so we can treatã as constant.
To progress, we can write a new differential operator (2.14) This operator allows us to finally write the axion-photon mixing equations as a Schrödingerlike equation As ∂/∂s is a function ofω p and θ, and in generalω p (s), any integration over s will in principle be curvilinear. However, if we are only concerned with small regions, we can assume thatω p is only slowly varying and so s is constant with respect to y, z. Assuming a constantω p we can writeŝ 1 1 +ω Comparing with Eq. (2.6) we see thatŝ is actually perpendicular toÊ LO . In other words, the propagating mode E LO is transverse with respect to the direction it evolves in, rather than the direction of propagation. While the mode is not transverse, it evolves in a similar manner to a transverse wave. Note that ∂/∂s is not defined with unit norm for simplicity, however that means one has to be slightly careful when talking about length scales associated with s. For non-relativistic axions (ω p = m a ), we can write down simple expressionŝ s cos θŷ + sin θẑ . (2.17c) In this limit we can simply show the relationship between B NS , z, E LO and s, depicted in Fig. 1. We can see in this limit that B NS ||E LO : the LO mode couples perfectly to the axion via E · B. For ultrarelativsitic axions with ω ω p , we see insteadÊ −ŷ,ŝ ẑ and ξ sin 2 θ giving back the usual mixing of an axion with a transverse photon [29].
Ignoring the overall phase, as shown in Appendix B, Eq. (2.15) can be solved in the stationary phase approximation to givẽ whereω p = ∂ω p /∂s, we have allowed for a small change in θ given by θ = ∂θ/∂s and s 0 is given by solving for the resonance condition (2.19) Figure 1. We show the various axes associated with axion-photon conversion in an anisotropic medium. The axion propagates in theẑ direction, which is at an angle θ from the magnetic field B NS . The magnetic field is defined to lie in theẑ,ŷ plane. In the non-relativistic limit, the E-field associated with the propagating LO mode,Ê LO , is given byB NS , i.e.,Ê LO is parallel to the B-field. The magnitude ofÊ LO evolves along the perpendicular direction, denoted byŝ.
One can see how the LO mode extrapolates between longitudinal and transverse behaviour as one varies θ. If θ = π/2, the dispersion relation is that of a regular transverse mode, with a resonance given byω p = m a . If instead θ = 0, then the resonance condition is that for an axion converting to a longitudinal plasmon with negligible spatial dispersion,ω p = ω, as for example studied in Ref. [34]. In the non-relativistic limit, ω → m a , and the resonance condition simply becomesω p = m a for any θ. The conversion length can be read off from the exponent in the stationary phase approximation (similar to Ref. [19]), giving where the latter statement holds in the non-relativistic limit. To evaluate the damping term in Eq. (2.18), we note that by assumption, the equation is only valid in a small region around the conversion zone where higher order derivatives and the bending of photons can be neglected. Assuming that ∂ω p /∂y = O(∂ω p /∂s), one can see that over the conversion zone (s 0 ± L/2) thus, for simplicity, we will neglect this damping in the probability of conversion. We can also include the longitudinal component of the LO mode, giving Actually, as discussed in Appendix A, there are two contributions to E z . As we are concerned with propagating modes that may exit the neutron star, here we only consider the longitudinal component of the LO mode. There is also a term that directly comes from the axion-E z mixing, however, it only persists in the region where the axion mixing is strong, and follows the phase of the axion, rather than the LO mode. Thus, after a length L, this axion-mixing term will dephase with the LO mode and not contribute to the E-fields exiting the neutron star. If we then takeω p = ω = m a , we get the simple expressions We kept the factor of k 2 in the denominator of E z to avoid an unphysical divergence at θ → {0, π}. Note that the E z component can actually dominate the E-field: which has a maximum at |E z /E y | m a /2k for k m a . Even for large k, E y and E z are generically of similar magnitude.

Probability of conversion
When expressing the classical problem of axion-photon mixing as a Schrödinger-like equation (with time replaced by the spatial coordinate z or s), it is natural to make use of the language and techniques originally developed within quantum mechanics to simplify calculations. The most natural object to compute from the Schrödinger equation is the 'conversion probability' (|A y /a| 2 ), i.e. the norm square of the transition amplitude of going from a pure axion state (a = 1) into a pure photon state (A y ). However, it is important to remember that this is not actually a quantum mechanical calculation, for which one would have to use the correct Schrödinger equation (with time derivatives and = 0), and carefully defined wavefunctions. In the vaccuum case, the 'conversion probability' calculated from the classical Schrödingerlike equation is of direct physical interest as it gives the ratio of the energy densities stored in axion and electromagnetic fields [29]; this definition of conversion probability has been used in Refs. [14,19,20,23] in the context of axion-photon conversion near neutron stars. However, as we discuss in this section, when taking into account the 3D complexity of the neutron star magnetosphere and the correct definition of propagating modes in the plasma, the naive 'conversion probability' is no longer the physically relevant object to compute. In essence, such a definition does not include the fact that the Hamiltonian is modified by the presence of matter, so would miscount the number of photons generated if interpreted as an actual conversion probability.
The key point of this section is that the (time-averaged) energy density stored in the propagating LO mode is no longer simply 1/4(1 +ω 2 p /ω 2 )|E y | 2 , as is the case for an isotropic non-relativistic plasma. To calculate the energy stored in the LO mode, we must consider that the medium has both temporal dispersion, and is anisotropic, giving [35,36] While the H-field terms are suppressed at order k 2 /ω 2 relative to the E-field terms, near the neutron star surface, k can be relatively large compared to ω, meaning that significant energy can be stored in the H-fields. For a propagating state given by E LO , we find Equation (3.2) reduces to 1/4(1+ω 2 p /ω 2 )|E y | 2 in the non-relativistic limit and when θ = π/2, and reduces further to |E y | 2 /2 for ω =ω p , and it is only in this case that the naive conversion probability (|A y /a| 2 ) gives the ratio of energy densities stored in the electromagnetic and axion fields. Notably, the energy stored in the propagating mode is enhanced for B-fields almost aligned with the direction of propagation. In the absence of medium losses, the energy stored in the photon field should be conserved: when propagated through a slowly changing medium the propagation state will adiabatically deform, but maintain stored energy. Thus, taking only the transverse component of the stored energy into account, as done in previous studies, neglects a significant part of the stored energy, or in other words, miscounts the number of photons that are converted.
For a simple expression, we can finally write the ratio of the energy densities of the electromagnetic propagating mode and the axion field in the non-relativistic limit (m a k) as As the system is considered here to be time independent, a single axion should convert to a single photon with energy ω and so the ratio of energy densities should correspond to the ratio between the axion and photon fluxes. Note that this flux transfer is distinct from the 'conversion probability' (|A y /a| 2 ) calculated from the Schrödinger-like Eq. (2.15), essentially reweighting it via the Hamilitonian. However, if one did a full quantum mechanical calculation, the expectation value of the conversion probability of axions to LO photons should agree with this flux transfer [37]. Unlike the result of Ref. [14], no propagating photon is generated from a non-relativistic axion in a longitudinal B-field (i.e., for B NS ẑ). However, we should stress that θ = {0, π} is a special condition and the order of limits matters. For purely longitudinal conversion one should perform a dedicated calculation (i.e., solve Eq. (2.8) directly or perform an analysis in the vein of Ref. [34]). We should note that Eq. (3.3) breaks down under more extreme conditions. If |ω p | = 0, in other words, ifω p is constant along s, the conversion length becomes infinite. To regulate such divergences, one must cut-off L when it is competitive compared to other limiting scales, such as the radius of curvature of the photon trajectory or that over which second derivatives are relevant [23], cf. Appendix B.
Our ratio can be contrasted with the 'conversion probability' calculated when the nontrivial Hamiltonian and the mixing of transverse and longitudinal modes are neglected, If we take θ = π/2 so thatŝ =ẑ, reducing back to a purely one dimensional problem, then we see that both expressions agree, as one would expect We have now derived in Eq. (3.3) an analytic expression for axion-photon conversion in highly aniosotropic systems. In cases where the axion is traveling at a non-trivial angle to the external B-field there can be significant alterations to the flux transfer. To see what impact the 3D calculation has in a neutron star, we must consider a specific neutron star model.

Phase space distribution
In this section, we discuss the phase space distribution of the axions at the conversion surface. An excellent discussion of the change of the phase space distribution of collisionless particles under the influence of the gravitational field of an astrophysical object such as a neutron star can be found in Ref. [38]. We will assume that far away from the neutron star, the velocity distribution of the axions in the galactic rest frame is given by a Maxwell-Boltzmann distribution with velocity dispersion σ v truncated at the Galactic escape velocity v esc as in the Standard Halo Model [39][40][41], where H(x) is the Heaviside step function, and the normalization factor is given by We can shift this distribution into the frame of the neutron star, where v NS is the velocity of the neutron star relative to the galactic rest frame. The sixdimensional phase space density far away from the neutron star is then where n ∞ a (ρ ∞ a ) are the axion number (mass) density far away from the neutron star. Liouville's theorem states that the phase space density is conserved along trajectories. Hence, at some point r NS in spherical coordinates centered on the neutron star, the phase space density is given by where v ∞ = v ∞ (r NS , v) is the velocity at infinity for an orbit with velocity v at r NS . To distinguish from our axion-centric frame as used in the previous sections, we use NS as a subscript; in other words, the position vector is given byr NS = (r NS , θ NS , φ NS ) in spherical coordinates. Let us stress here that Liouville's theorem states only that the phase space density is conserved along trajectories, it does not make any immediate statements about the relation of the functional form of the phase space density at different r. In order to obtain that functional form, one must compute v ∞ (r NS , v). For a Newtonian potential, energy, angular momentum, and the Laplace-Runge-Lenz vector of a particle moving in the potential are conserved, such that [38] v with v ∞ = v 2 − 2GM NS /r NS , the gravitational constant G, the mass of the neutron star M NS , andr NS = r NS /r NS . It is interesting to note that due to the symmetry of the system, if the velocity and mass distribution of the axions is homogeneous and isotropic (in the frame of the neutron star) far away from the neutron star, i.e. v NS = 0 in Eq. (4.3), the velocity distribution will be isotropic at any r NS . In this case, one can immediately compute [20] With Eqs. (4.4)-(4.6), the dark matter number density at r NS can be obtained by integrating over all kinematically allowed velocities, v For an homogeneous and isotropic (in the frame of the neutron star) axion phase space distribution far away from the neutron star and v esc σ v , one can find an approximate closed solution for n a (r NS ), see Ref. [20]. We stress that their solution [and Eq. (4.7)] is only applicable if, far away from the neutron star, both the axion density and velocity distributions are both homogeneous and isotropic (in the frame of the neutron star) and if v esc σ v . A typical neutron star moves with speed v NS of order σ v relative to the Galactic rest frame, breaking the assumption of an isotropic velocity distribution in the frame of the neutron star.
While the full six dimensional phase space contains all relevant information, for any realistic conversion location, 2GM NS /r NS {σ v , v esc }. Hence, almost the entirety of the axion's speed comes from the infall. Thus, the distribution over v can approximately be treated as a delta-function δ(v − 2GM NS /r NS ). In other words, we care about the distribution of incoming angles, rather than the distribution over incoming speeds. This allows us to define Thus armed, we can calculate the flux transfer for general axion distributions, including non-vanishing speeds of the neutron star relative to the galaxy.

Comparison with 1D calculations in the Goldreich-Julian Model
It is clear that for specific values θ there can be significant variation between the 3D and 1D versions of the axion-photon conversion. As the plasma is generated and shaped by the neutron star magnetic field, the changes in ω p are driven by the changing B-field, and so may also be correlated with θ.
Here we explore what influence a more complete axionphoton conversion calculation has under realistic conditions. After a brief discussion of the Goldreich-Julian neutron star model in section 5.1, we will compare results for the axionphoton conversion in our full 3D calculation with results in the conventional 1D calculation in sections 5.2-5.3.

The Goldreich-Julian Model
To consider a realistic scenario, we turn to an analytic description of a neutron star. One commonly used model is the Goldreich-Julian (GJ) model [42]. The GJ model gives a simple analytic solution by requiring a charge density to satisfy the condition E·B = 0 at the surface of the neutron star in the presence of a strong, rotating dipole magnetic field given by [cos θ m cos θ NS + sin θ m sin θ NS cos(Ωt)] , where B 0 is the magnetic field strength at the surface of the neutron star (which is at located at r 0 ), and θ m is the misalignment angle between B NS and Ω. The rotation vector of the neutron star is given by Ω = Ωẑ NS with Ω = 2π/T and T the neutron star spin period. We denote the polar angle of a given location with respect to the rotation axis z NS to be θ NS . The charge density is given by While we consider all charges to be electrons for simplicity, such that n GJ simply gives the number of electrons, our conversion calculation still applies to more complicated plasmas. The plasma frequency of an electron plasma is ω p 4παn e /m e , giving Similarly to Ref. [14] we have also defined m ·r NS = cos θ m cos θ NS + sin θ m sin θ NS cos(Ωt) .

(5.5)
It should be noted that it is unlikely that the GJ model holds over the entirety (or any) of the neutron star. The nature of neutron star atmospheres is highly uncertain, due to both the difficulty in modeling such large electromagnetic systems and the lack of observational data [43][44][45][46]. Asides from possible higher multipole components of the magnetic field [47], pair creation cascades at the poles may lead to significantly higher electron densities, as well as highly boosted electrons [45,48,49]. Our formalism allows for both γ 1 and n e = n GJ , however, for specific examples we will remain within the GJ model.

Directionality of emitted photons
To see the impact of a 3D formalism, we can take the simplest example and ask: how does the flux transfer change for different incoming axion trajectories? As we are considering the GJ model, for all cases γ = 1, meaning thatω p = ω p . In section 3 we derived the flux transfer R in a frame centered around the axion's velocity (ẑ) and the direction of the B-field. However, as discussed above, the properties of the plasma are usually described in the frame of the neutron star, typically in spherical coordinates. To write down ourx andŷ directions in terms of the incoming axion directionẑ and magnetic fieldB we can usê For an example, we consider a neutron star with B 0 = 10 14 G, θ m = 0.4 and T = 1 s. Taking some point on the sphere, θ NS = 0.5, φ NS = 0.6π, we can compare the flux transfer at the axion conversion surface as a function of the incoming axion angles, denoted by θ v , φ v in the neutron star coordinate frame. Due to the dependence of the LO mode dispersion relation on θ and k, Eq. (2.19) indicates that the radius of conversion actually is slightly different for different incoming axion angles. However, the change in flux transfer due to small changes in conversion radius is subdominant compared to the impact of including aniosotropic effects of the medium. For a benchmark, we take the axion to have mass m a = 25 µeV and coupling g aγ = 10 −14 GeV −1 . Raising or lowering m a changes where the conversion zone is located: higher mass axions convert closer to the neutron star where magnetic fields are larger, usually leading to greater flux transfers. If the axion mass is too large, however, it is greater than the plasma frequency everywhere outside the surface of the neutron star and no resonant conversion can occur. As the plasma frequency is set by the magnetic field strength of the neutron star, for each star, there will be a maximum m a which will allow for resonant axion-photon conversion.
We plot the results of the 3D calculation, [the relativistic version of Eq. (3. 3)], and of the 1D calculation, Eq. (3.4), in the right and left panels of Fig. 2, respectively. To avoid divergences when the WKB approximation breaks down, we use a simple cut-off estimated by the typical lengths scales over which the second order derivatives should come into play, derived in Appendix C We see that the number of photons generated on a given trajectory can differ between the 3D and the 1D formalism by several orders of magnitude. Furthermore, regions that were highly disfavoured in the 1D calculation can turn out to have maximal photon production when the aniosotropic nature of the medium is taken into account.
To understand better the underlying physics behind the regions of enhanced and suppressed conversion, we can take a single slice of θ v , φ v space. In Fig. 3, we show a slice of constant θ v , in particular, θ v = π/3. The top panel of Fig. 3 shows the 'conversion probability' for the 1D calculation and the flux transfer calculated with 3D effects, with the regulating cut-off flux transfer [Eq. (5.7)] shown by the gray dashed line. In the middle panel, we show the ratio R/P 1D a→γ : up to three orders of magnitude enhancement can occur for specific axion trajectories. The bottom panel shows the relevant derivatives for each calculation, ∂ω p /∂z in the 1D case and ∂ω p /∂s for the 3D calculation. The order-of-magnitude discrepancies between the 1D and 3D calculations of the flux transfer occur when ∂ω p /∂z = ∂ω p /∂s. While there are still differences caused by correctly calculating the energy stored in the LO mode and the correct angular dependence, the primary difference between calculations comes from the difference betweenŝ andẑ: when ∂ω p /∂s ∂ω p /∂z, the answers are similar up to the different angular dependence. The divergences in R, P a→γ (which are regulated by the cutoff) occur when ∂ω p /∂(s, z) = 0, as indicated by the gray dashed line in the bottom panel of Fig. 3. In these examples, θ is negligible and so the cut-off is a stronger regulator. Note that the relationship betweenŝ andẑ does not depend on m a : even for very non-relativistic axions which could convert far out from the NS there would still be significant differences in the flux transfer in a given direction.
We can see that treating the LO mode as essentially a transverse mode (possibly with different dispersion, as in [23]) leads to a significant miscalculation of the axion-photon flux transfer. To have a reliable calculation of the axion-photon flux transfer for a given axion trajectory, one must take into account both the full energy stored in the LO mode, and, more importantly, the distinction between the direction of the axion's momentum,ẑ, and the direction over which the energy density evolves,ŝ.

Integrated power
Asides from the directional dependence of the photons emitted from a given point, one may also be concerned with the total power emitted over some portion of the neutron star. Studies of ray-tracing suggest that the photons bend relatively quickly as plasma frequency drops [23,24]. The momentum acquired by the photon as it transitions to vacuum must come from the gradient of the medium, so photons that begin mildly relativistic would become approximately parallel to the change inω p . Because of this, one can make a rough estimate that the total power exiting the neutron star will end up on the same trajectory, determined by ∇ω p . We stress however that the axion can be reasonably relativistic, and so one should do a full ray-tracing calculation for reliable results.  Figure 4. Axion photon conversion integrated over axion phase space, R W , as a function of incoming neutron star polar angle θ NS for 3D (green) and 1D (blue) calculations. The 3D calculation uses the full flux transfer of axions to photons, whereas the 1D case only considers the 'conversion probability' interpreted from a 'Schrödinger-like' equation. We choose a neutron star with B 0 = 10 14 G, θ m = 0.4 and T = 1 s. The functions are calculated at the conversion zone located at an azimuthal angle φ NS = π in the neutron star frame. We assume an axion mass of 25 µeV with a coupling strength g aγ = 10 −14 GeV −1 as well as an isotropic velocity distribution of the axions. We mark the region where the conversion radius would be less than the neutron star radius r 0 with the red dashed lines labeled 'Inside NS'.
To calculate the power radiated from a given point on the neutron star r NS , we must integrate over the incoming axion angles θ v , φ v , weighted with an appropriate phase space distribution. To focus on the difference caused by properly treating the aniostropy of the medium, we will neglect the small differences in conversion radii for different axion trajectories and use a single radius of conversion given by the transverse photon dispersion, m a = ω p . The variation due to choosing different radii seems to be at most a factor of 2 when one considers the longitudinal photon dispersion (ω = ω p ) instead, and does not change the general relationship between the two calculations. Similarly, we neglect contributions from θ , which again gives a subdominant contribution for the cases we consider. The distributionweighted axion-photon flux transfer R W is given by where f 5 is the five dimensional phase space defined in Eq. (4.9). For the 1D case we use the 1D 'conversion probability', Eq. (3.4) instead of R.
We can now show R W for a slice along the neutron star conversion surface assuming an isotropic axion distribution. We keep φ NS fixed at π and vary θ NS to derive Fig. 4. The structure generated in a slice through θ NS is in general behaviour independent of the specifically chosen φ NS . Integrating over all photon directions generally smooths out differences between the 3D and 1D calculations. However, there are still noticeable differences. Overall, the 3D calculation leads to more energy stored in the EM fields, leading to more photons radiated from the neutron star. As the 1D calculation neglected energy stored in the non-transverse components of the Hamiltonian, we would naturally expect to see additional 'conversion' just from correctly counting the number of photons generated. This enhancement is particularly noticeable around the neutron star 'throat', i.e. the region where the conversion radius vanishes as Ω · B NS → 0, just outside the red boarders in Fig. 4. Inside the neutron star itself, the plasma density no longer follows the GJ model; given the extreme densities inside a neutron star the plasma frequency is much larger than the axion masses of interest everywhere inside a neutron star. Thus, no propagating photons can be produced from axion-photon conversion inside a neutron star. One would expect that issues of photon curvature, as explored in Ref. [23], as well any concerns about gravity bending the axion path would be most pronounced in the throats, making this region generally unreliable. One would likely need to do numerical calculations for a robust solution including a gravitationally bending axion trajectory.
Having a similar total conversion rate for isotropic dark matter infall is not surprising, as ultimately the length scales and magnetic fields that give the general order of conversion are the same. However, where the axions are converted, and in what directions the photons emerge, is significantly impacted. The plasma effectively lenses the outgoing photons, which when combined with the limited line of sight and Doppler shifting of light can give a complicated, time dependent signal. Thus any aspect which changes the location and direction of conversion will impact axion searches significantly. Because of this we focus on the changes to the power and directionality of a given point on the neutron star, as the total power radiated from a neutron star is not usually relevant to an individual observation of an individual neutron star. In addition, for compact clumps of axions, such as axion stars [50][51][52][53][54][55] or miniclusters [56][57][58][59][60][61], one would expect a much stronger dependence in the resulting signal depending at what angle the compact object hits the neutron star. Thus, changes to the production location and direction can have a disproportionate effect on observational constraints.

Conclusions
The axion is one of the premier candidates for dark matter. It displays unique wavelike phenomenology requiring both novel techniques in the lab and astrophysical searches. Resonant conversion to photons in the atmospheres of neutron stars is one such possible method for discovering the axion. Neutron stars are complicated systems and therefore calculating the projected signal is extremely challenging. Here we have focused on one piece of the larger puzzle, and together with recent advances in ray-tracing help lay the road to a robust prediction for what might be observed on Earth.
To this end, we have calculated the production of photons via axions in strongly magnetised, anisotropic environments. The prototypical example of such an environment is a neutron star, but our calculation is generic. Previous calculations have not explored the anisotropic nature of the medium and have focused only on the transverse part of photon modes. We show, however, that the mixing of transverse and longitudinal modes into Langmuir-O modes causes the direction over which the amplitude of the electromagnetic fields evolve to be different from the direction of the axion. We find that rather than evolving in the direction of propagation, the waves evolve also in the transverse direction. This directional difference can lead to orders of magnitude change in the flux transfer R in a given direction.
Furthermore, taking only the transverse component of the stored energy into account neglects a significant part of the stored energy, or in other words, miscounts the number of photons that are converted. Taking only the transverse component essentially confuses a classical, Schrödinger-like solution with the quantum mechanical conversion probability. We show that taking both features into account significantly modifies the probability of conversion over a given axion-trajectory. While a truly complete solution would require numerically solving the equations of motion in a varying medium, we provide simple analytical formulae which can be used to estimate the radio signals generated in neutron stars by axions. This will allow the robustness of axion limits to be increased in the case of non-observation, such as Refs. [25][26][27][28], or allow one to estimate the neutron star and axion properties in the case of a discovery.
The first term is just the longitudinal component of the LO mode, with the second corresponding to direct axion mixing. In the limit of a perfectly aligned magnetic field with θ = {0, π} Eq. (A.4) reduces to the expression for a pure longitudinal excitation in a lossless plasma [31,62] (A.5) We can see that the two different terms in E z have very different propagation behaviour. After some distance (given by the conversion length L)Ẽ y will be out of phase withã. Notice that the terms proportional toã are not generally propagating: as soon as the axion mixing goes away, so to does this contribution to E z . Thus, for the purpose of solving for an outgoing radio signal from the neutron star, the mode that the axion couples to is the LO mode, which will adiabatically transform into a transverse mode as it exits the system [33].

B Stationary phase approximation
To provide an analytic expression for axion-photon conversion, in Eq. (2.15) in the main text we derived a Schrödinger-like equation Rather than a numeric solution, we can use the stationary phase approximation to find an analytic solution. The premise of the stationary phase approximation is that highly oscillatory integrands tend to cancel when integrated over many oscillation periods, and that the dominant contribution therefore often comes from stationary points, say s 0 , where the phase, f (s), is locally constant, i.e. f (s 0 ) = 0. In the standard form of the stationary phase approximation, the non-oscillatory factor of the integrand, say g(s), is approximated as constant around s 0 , and the phase is expanded to quadratic order: f (s) ≈ f (s 0 )+ 1 2 f (s 0 )(s− s 0 ) 2 , and the limits of integration are taken to ±∞. For s > s 0 , this gives the approximation, 3) The quality of the stationary phase approximation of course depends on the accuracy of the assumptions that it relies on. It's convenient to define a (squared) 'conversion length' as In the limit where |f (s 0 )| is larger than all other scales in the problem, the conversion length is small, and the stationary phase approximation is expected to be very accurate. However, in many physical situations there are no large hierarchies of scales, and one needs to consider corrections to the standard stationary phase formula. Moreover, extending the integration range to ±∞ is not always justified, and the full integral that one wants to evaluate may only be well-defined in a small region around the stationary point, such as for our Eq. (B.2).
In the context of relativistic axion-photon conversion, these issues were recently discussed in Ref. [63], where corrections to Eq. (B.3) were derived from considering a finite integration range, and accounting for contributions neglected in the standard stationary phase approximation. For the purpose of our analysis, the main point is that when one considers the integration range [s 0 − ∆s, s 0 + ∆s], the stationary phase approximation gives where Err denotes the error function. The norm of the error function grows linearly with ∆s L for small arguments, and performs slowly damped oscillations around unity for ∆s L > 1. This in particular implies that the standard formula only applies when it is justified to take ∆s L. Equation (B.6) applies for any ratio of ∆s/L, but simplifies considerably when ∆z < L, and the integral becomes up to corrections to cubic order in ∆z/L. This implies that the formula (B.3) can be regulated and modified in a very simple way, by replacing L → √ 2∆z, when it is justified to consider ∆z < L, for example when the WKB approximation breaks down, or the assumption of constant θ fails. For a more detailed discussion of how to evaluate also the contributions from |z| > ∆z, see Ref. [63].
The most important factor comes from f (s 0 ). While we have so far considered θ to be constant, if it is slowly changing the largest impact will be felt in the conversion length [23] as the impact of dephasing the axion and photon is much greater than the small changes to the definition ofŷ. We can then see that where we have neglected an overall phase and whereω p = ∂ω p /∂s. This is the analytical form of axion-induced E y that we will use in our subsequent analysis. The stationary point condition (cf. f (s 0 ) = 0) translates into the 'resonance condition': For most axion trajectories, this conversion length is smaller than the length scales associated with the gravitational curvature of the axion trajectory, the curvature of the photon path due to the anisotropic medium, and the scale over which the magnetic field changes. When this is the case, the standard stationary phase approximation applies. However, this is not always the case, and in some cases the conversion length can get significantly enhanced, as we now discuss.

C Regulating divergences
The flux transfer R, given by Eq. (3.3), was calculated assuming the WKB approximation holds: i.e., that the axion momentum is much larger than the first derivatives ofẼ y , which are larger than the second derivatives. However, R diverges in the limit that the first derivative vanishes (i.e., the conversion length becomes infinite). Such a divergence must be regulated for two reasons. The first is the breakdown of the WKB approximation, and the second is the breakdown of the assumption that the conversion happens over length scales smaller than other relevant ones in the neutron star, for example the radius of photon bending due to the changing media.
In the first case, it is evidently untrue that the vanishing first derivatives are larger than the second derivatives. Rather than attempting to solve the full 3D equations by hand with higher derivatives, as E x , E y and E z all appear with their second derivatives, we can use a simple estimate to institute a cut off on the diverging flux transfers. In the 1D case, it was shown in Ref. [14] that an order-of-magnitude estimate could be obtained by using the simple formula where L is the conversion length. The assumption of coherent conversion requires that the change in phase between two locations is small. Thus, we can estimate L by checking the distance over which k changes. We can approximate k changing over some length L by writing k(L) = k 0 + k 1 L + k 2 L 2 , (C.2) where destructive interference occurs when δk = k − k 0 satisfies δkL π (as the phase enters via e i(kL−ωt) . As we are concerned with the case where k 1 L k 2 L 2 , we can neglect the linear term and see that L 3 π/k 2 .

(C.3)
As the only change in k comes from the changing plasma frequencyω p , in the limit of vanishing first derivatives kk 2 =ω pωp,2 , where similarlȳ ω p (L) ω p,0 +ω p,2 L 2 . (C.4) To get a conservative estimate for L, we can use that at second order, all directions may play a role. Thus while the second derivatives with respect to some directions might be small or vanishing, we know that the radial direction will have a power law dependence. Thus, assuming that ther direction is the relevant one will give a fairly conservative estimate for the total conversion probability. We stress that this is not a full solution of the 3D Maxwell equations, but rather a conservative way of regulating our WKB formalism. Taking the power law dependence ofω p , we then find that L = πv c r 2 c ω p 1/3 .

(C.5)
Putting in some typical numbers, we see that the typical conversion length should increase by about an order of magnitude, when compared to purely radial estimates of Ref. [14], leading to a hundredfold increase in the conversion probability. The second case relates to the overall trajectory of photons in the medium, for example, that the photon bends and so dephases with the axion, or other curvature effects of the medium [23]. While such a constraint may provide a more stringent limit, particularly in the neutron star throat, to calculate such effects in full requires ray tracing calculations [20,23,24], which are beyond the scope of this work. When these effects are considered, one should use the stricter of the two cut-offs. As we are only interested in demonstrating the influence of a 3D, aniosotropic calculation of the flux transfer rather than a start to finish signal calculation for observational purposes, we can simplify matters and just use a cut-off based on the second derivatives of the E-field,