Generalized Ray Tracing for Axions in Astrophysical Plasmas

Ray tracing plays a vital role in black hole imaging, modeling the emission mechanisms of pulsars, and deriving signatures from physics beyond the Standard Model. In this work we focus on one specific application of ray tracing, namely, predicting radio signals generated from the resonant conversion of axion dark matter in the strongly magnetized plasma surrounding neutron stars. The production and propagation of low-energy photons in these environments are sensitive to both the anisotropic response of the background plasma and curved spacetime; here, we employ a fully covariant framework capable of treating both effects. We implement this both via forward and backward ray tracing. In forward ray tracing, photons are sampled at the point of emission and propagated to infinity, whilst in the backward-tracing approach, photons are traced backwards from an image plane to the point of production. We explore various approximations adopted in prior work, quantifying the importance of gravity, plasma anisotropy, the neutron star mass and radius, and imposing the proper kinematic matching of the resonance. Finally, using a more realistic model for the charge distribution of magnetar magnetospheres, we revisit the sensitivity of current and future radio and sub-mm telescopes to spectral lines emanating from the Galactic Center Magnetar, showing such observations may extend sensitivity to axion masses $m_a \sim \mathcal{O}({\rm few}) \times 10^{-3}$ eV, potentially even probing parameter space of the QCD axion.


I. INTRODUCTION
Astronomy is the art of inferring details about astrophysical environments through indirect measurements on the messengers they emit, be they photons, gravitational waves or neutrinos.A basic problem in astronomy is therefore to model the production and propagation of these messengers from their point of emission to the moment of detection at the observatory.Taking the case of electromagnetic signals, this entails calculating the photon production mechanism and the subsequent evolution of photons as they pass through astrophysical media.In general this requires tracking, amongst other things, their intensity, frequency, polarisation and refraction.In turn these features affect the power, directional dependence, time-variation and spectral morphology of the signal.
Axions were originally introduced many decades ago to address one of the major outstanding problems in highenergy theory, the so-called Strong CP Problem [26][27][28][29]; this is effectively the question of why QCD seems to conserve charge-parity symmetry, or equivalently, why the electric dipole moment of the neutron appears to be so unnaturally small.Today, the term axion is typically used, however, to refer to the broader class of light pseudoscalars, regardless of whether they solve the strong CP problem; such particles are nevertheless well-motivated candidates for new fundamental physics, as they generically arise in well-motivated high energy theories such as String Theory [30][31][32][33][34].
One of the recent and more compelling proposals to indirectly search for axions in astrophysical environments involves looking for radio photons generated from axionphoton mixing in the magnetospheres of neutron stars, where the large magnetic fields and ambient plasma serve to resonantly amplify the mixing process.The presence of axions in these environments can give rise to a num-ber of distinctive signatures at radio energies, including narrow spectral lines from the conversion of axion dark matter [21,22,[35][36][37][38][39][40][41][42][43][44][45][46][47][48][49][50], and large broadband emission generated from axions locally sourced in the magnetosphere itself [25,51,52].In recent years, ray tracing has played an increasingly prominent role in understanding of the inhomogeneity, spectrum, and temporal evolution of these radio signals, and has proven to be a fundamental tool in accurately interpreting radio searches for axions [22,25,44,49,50].
The goal of this manuscript is to develop a generalized ray tracing framework capable of treating the production of radio photons from axions in astrophysical plasma, with a particular focus on the treatment of ray tracing through a highly magnetized plasma in curved spacetime, as is required for axion searches near neutron stars.To date, the ray tracing algorithms used in this field have only incorporated a subset of the relevant effects, including either the presence of an anisotropic plasma [22] or curved spacetime effects in an isotropic plasma [23].Moreover, the methodology of these two algorithms differ markedly, with one approach back-propagating photons from an asymptotic observer, and the other forward propagating photons from the point of production (see Fig. 1 for an illustrative example of the two setups).This work serves to unite these frameworks, highlighting a variety of important subtitles in the field which have thus far gone overlooked.While the focus of this paper will be directed toward applications of ray tracing for axion searches, the formalism is sufficiently general to be applicable, and of interest for, the broader astronomy and astrophysics communities.
The structure of this paper is organized as follows.In Sec.II we present an overview of geometric ray tracing.In particular, we derive a general formalism in which asymptotic observables (such as the radiated power in a given emission direction from a source) can be reconstructed by either (i) backward propagating rays from an infinitesimal surface element located far away from the emitting region, or (ii) forward propagating rays from the emitting region to asymptotic distances.This formalism can be straightforwardly applied to arbitrary spacetime metrics and any background medium -the generalization of these techniques to photons and axions in magnetized plasmas and to curved spacetime are the focus of Secs.III and IV.In Sec.V, we apply the formalism developed in the preceding sections to the specific problem of searching for radio spectral lines produced from the resonant conversion of axion dark matter in the magnetospheres of neutron stars.We investigate the interplay of a number of important effects, including, e.g. , the impact of plasma anisotropy, strong gravity, multiply reflected photons, and a proper kinematic matching of the resonance.We also revisit the extent to which radio observations can be used to search for axions near the galactic center magnetar SGR J1745-2900; here, we apply an improved modeling of the charge distribution in the magnetosphere, showing that expected e ± densities found in the closed magnetic field lines of magnetars shift the expected signal to higher frequencies, potentially generating signals in the O(100) GHz -THz regime.Using rough rough estimates of the magnetic field and charge normalization, we show that an improved understanding of the properties of SGR J1745-2900 allow current and future radio and sub-mm telescopes to probe unexplored regions of the axion parameter space.We conclude in Sec.VII.

II. RAY TRACING
In this section, we outline the basics of ray tracing required to describe radiative transport in arbitrary media.We introduce Hamilton's equations, which allow one to identify the wordlines of photons, and discuss distribution functions of photons and their phase-space integrals, which allow for the computation of observables such as energy flux or radiant intensity.We use two distinct methods to numerically calculate observational signatures relevant for astronomy.The first approach uses image-based, backward ray tracing from the observer [23] and the second implements a forward ray tracing routine which directly samples the phase-space of photons produced at source, and reconstructs asymptotic observables by propagating these photons to large distances [22].
When discussing radiative transport, we are typically interested in radiation whose wavelength is much smaller than characteristic variational scales of the medium through which the radiation propagates.In that case, geometric optics applies, and one can describe photons in terms of their local phase-space coordinates given by position x µ momentum k µ .
The dispersion relation of these photons is then described by setting a Hamiltonian H = H(k, x) to zero H = 0, which gives a relation k 0 = k 0 (k, x, t) describing the energy of the mode.This defines a family of curves (k µ (λ), x µ (λ)) in phase-space, where λ is the wordlineparameter on such curves, along which H is everywhere vanishing: i.e. dH(k(λ), x(λ))/dλ = 0.By applying the chain-rule, one can see that which is satisfied provided x µ and k µ obey Eqs. ( 2) are the well-known Hamilton's equations, which allow one to determine the trajectory of photons from the source to the observer.In order to connect the properties of the source (e.g.emissivity) with the asymptotic observables (e.g.observed flux), we can introduce Boltzmann's equation [53,54]  where f = f (k, x) is the phase-space distribution of photons and S[k, x] is a source term whose precise structure depends on the emission process.Here, one can see that Hamilton's equations (Eq.( 2)) give the characteristics (or orbits) of the Liouville-Vlasov operator appearing on the left-hand side of Eq. (3).Thus, along worldlines, we have in the absence of collisions (S = 0), this implies that f γ is conserved along rays, i.e.
which is equivalent to Liouville's theorem [53].Equation ( 5) is the key equation in numerical ray tracing routines; it allows the asymptotic distribution of photons to be reconstructed by saturating the space with rays, computing the value of f γ at source, and propagating this conserved quantity along rays.In other words, ray tracing effectively amounts to solving the 8-dimensional equation (3), which describes the spacetime dependence of the photon distribution function f γ .By integrating Eq. ( 3), one can also derive a continuity equation for energy [54].Assuming a (quasi) stationary background (such that ∂ t H = 0), Eq. ( 3) can be premultiplied by k 0 , placed on-shell, and integrated over 3momentum d 3 k and a spatial volume d 3 x = dV.This procedure yields where dA is the surface element of V (which we assume lies outside the source), ω is the (on-shell) photon frequency, v g is the photon group-velocity, and Q is defined by where C = S/(∂ k0 H) is a renormalised collision kernel whose denominator arises from dividing both sides of Eq. ( 3) by ∂ k0 H before integrating over phase space to obtain Eq. ( 6).For a stationary solution, the first term on the left-hand side of Eq. ( 6) vanishes, and the total power emanating from the source is Eqs. ( 7) and ( 8) are essential for ray tracing, as they allow for the construction of observables from the individual rays.As we will show in the following subsections, Eq. ( 7) lies at the heart of the forward ray tracing procedure, while Eq. ( 8) is key to understanding backward tracing; Eq. ( 6) serves to unite these frameworks, since for stationary solutions these two quantities are equal.

A. Backward ray tracing and Imaging
Backward ray tracing is the traditional approach used to infer the emission properties in different directions from a source.This approach proceeds by constructing an infinitesimal surface far from the object (oriented with a surface normal parallel to the line of sight), and tracing parallel rays from the surface 'backwards' in time until they encounter the source itself.The differential power flowing the infinitesimal surface is given by where θ g is the angle between the photon group velocity and the surface.Note that if the infinitesimal surface is sufficiently far from the source then one can safely assume the photon to be in vacuum, implying v g = 1, cos θ g = 1, and d 3 k = dΩ ω 2 dω.Following [21,23], the power per solid angle per unit frequency is then given by a summation over each of the rays i, appropriately weighted by their contribution to the surface area dA i , their energy ω i , and their phase space factor f i γ , i.e.
The only unknown quantity in Eq. ( 9) is f i γ , which is determined by defining f γ at the source, and using the fact that f γ is conserved along rays.
In order to implement this procedure numerically, we define the asymptotic surface to be a plane consisting of square pixels with side length ∆b.Rays are sourced through the center of each pixel, labelled by (i, j), in a direction perpendicular to the surface itself.For simplicity, we assume photons are monochromatic, which allows one to evade sampling over energy.This procedure generates a differential power given by dP(θ, φ) dΩdω = i,j Note, however, that non-uniform sampling of the asymptotic surface can, in general, dramatically expedite the numerical calculation.

B. Forward ray tracing
An alternative approach to backward ray tracing is to directly simulate photon production from the collision integral Eq. ( 7), and trace the rays forward to all parts of the sky surrounding the source.For a (quasi) stationary background, the power flowing at infinity is simply equal to the power produced at the source, i.e.
Rather than sampling rays from an asymptotic surface (as in backward ray tracing), forward ray tracing works by stochastically sampling the photon phase space at the source, i.e. it uses Monte Carlo integration to directly compute the right hand side of Eq. (12).Photons are then propagated away from the source and binned on a sphere at infinity; the angular power distribution in a direction (θ, ϕ) is then reconstructed by summing over the weighted rays which end up in a small bin centered about (θ, ϕ).
At this point, one may be concerned by the fact that Eq. ( 12) relates the integrated power locally to the integrated power at infinity, while the forward ray tracing procedure described above relies on the fact that this connection between local and asymptotic power also holds at the differential level.Notice that from Liouville's theorem, it follows that the number of particles in a phase-space element is conserved along rays, so that dN = d 3 kd 3 xf is constant along rays.More explicitly, we can write this as dN = d 3 k cos θ g v g dtdA, so that the power flowing through an infinitesimal surface at the point of emission, is equal to the power flowing out of another surface at infinity (so long as these points are connected by a ray).It is ultimately for this reason that the forward propagation approach is valid, since all rays in a given part of the sky conserve the phase-space element.An equivalent way of viewing this is via the conservation of etendue.This means that not only is the total integrated power conserved (and invariant under a re-partitioning of the Monte-Carlo integration), but any two sub-surfaces connected by rays, also have conserved power flowing between them, and are therefore also invariant under a re-partitioning of the Monte Carlo integration.
In order to provide an illustrative example of the forward ray tracing procedure, consider the case of uniform isotropic emission from a finite volume V C .By drawing N s uniform samples over d 3 x and d 3 k, one can write the differential power flowing through a patch of the sky (θ 0 ± ϵ θ , ϕ 0 ± ϵ ϕ ) as where θ f,i and ϕ f,i are the final angular coordinates of the photon after propagation, and we have defined the function As in the case of backward propagation, the choice of uniform sampling may not be optimal, and one may instead prefer to implement importance sampling.As illustrated by this example, the forward-tracing method is fully general and can be used to solve the radiative transport problem in any setup.Ultimately, however, our purpose is to study the production of photons from axions; this problem is more subtle in that axionphoton mixing is a one-to-one process, meaning that it only occurs on particular surfaces in phase-space (namely FIG. 2. Photon Emission Points .Subset of photon emission points from rays which propagate to an observing direction (θ, ϕ) = (35 • , 0 • ), projected onto the point of photon production, for the case of backward ray tracing (left) and forward ray tracing (right).Gray contours illustrate the 2D surface around the neutron star defined by ma = ωp.Results are shown for the case of an anisotropic dispersion relation in curved spacetime, with MNS = 1 and ma = 10 −5 eV.
at locations where the dispersion relations of axions and photons become degenerate).This condition collapses Eq. ( 12) to a momentum-weighted sum over kinematic surfaces rather than volumes.In the following sections we provide a more detailed discussion on how to generalize the Monte Carlo integration procedure to the case of axions.

C. Comparing Forward and Backward Ray Tracing
Having summarized the general approach to forward and backward ray tracing, we outline here the potential advantages and disadvantages of each.
• Forward propagation inevitably generates samples across the entire sky (i.e. the sphere at infinity surrounding the source).In the event that the observing angle is known, i.e. one is only interested in the power radiated in a particular direction on the sky, the forward propagation approach clearly suffers from oversampling, since only a fraction of the samples are actually used in the computation of the relevant observable.In the event that this sampling is uniform, the over sampling may be severe, however there exist many forms of adaptive sampling algorithms which can be used to improve the sampling efficiency.This is not a problem for backward propagation, since by construction all samples originate from the angle of interest.The issue of oversampling can be seen in the comparison provided in Fig. 2; in this example, low-energy photons are assumed to be sourced from axions near a neutron star (see Sec. III), and their point of origin for a particular viewing angle (taken here to be (θ, ϕ) = (35 • , 0 • )) is reconstructed using both backward and forward ray tracing techniques.While the rays originate from the same location near the neutron star, the density of samples in the forward ray tracing approach is significantly reduced with respect to the backward ray tracing example.
• Conversely, generating full sky distributions of the flux is more complicated in the context of backward propagation, as one must scan over (and interpolate between) many viewing angles (whereas this is a natural output of the forward sampling procedure).Full sky distributions can prove useful for understanding the fundamental physics and observables in the problem at hand, extracting the time profile of the signal, and in some cases, full sky distributions are required in order to marginalize over uncertainties associated with an unknown viewing angle.This could be circumvented by sampling the observing directions and positions on the observing plane (see Ref. [23]) stochastically (i.e. , via Monte Carlo integration) rather than deterministically.
• Backward propagation also allows one to reconstruct physical images -see Fig 3.In this work, we are concerned principally with neutron stars, which are too small to resolve and so the precise structure of the image plays no role from a data analysis perspective in that context, since at low resolution, one is sensitive only to the integrated power over the image plane.While image reconstruction could become important in other contexts, it is unclear whether any additional information is added that cannot be directly obtained from the photon source locations (which is information that is available in both approaches).
• As will be shown shortly, divergences arise in the conversion probability for axion-photon mixing.These divergences are naturally regulated by the phase-space measure (see discussion around Eq. ( 20)).The removal of these divergences in the forward tracing is therefore somewhat tautologous, as the cancellation occurs as soon as the integral Eq. ( 20) is written down.This integral (which is then Monte Carlo sampled in the forward-tracing numerical routine) is therefore trivially convergent.However, in backward ray tracing, the conversion probability leads to quantities which diverge at an individual ray-level; in principle, these divergences should be regulated by the phase-space measure in the image plane.Backward ray tracing therefore provides a powerful consistency check for the treatment of phase space, kinetic theory and the conversion probability in this work and [54], and indeed the equivalence and finiteness of forward and backward results is ensured by the expression Eq. ( 6).
The existence of two independent numerical routines has therefore proved extremely fruitful in developing both numerical and analytic understanding of the problem at hand.
• Gravity plays a stabilizing role in backward ray tracing, since rays must be back-propagating from regions of high refractive index, to regions where rays are almost evanescent (ω p ≃ ω).This has a tendency to deflect rays away from the production surface.In the absence of gravity, ultra-high numerical precision is therefore needed for rays to converge to the point of resonance where ω p ∼ m a ≃ ω.
In the case of forward propagation this problem is evaded, allowing one to independently study the impact of gravity on ray propagation.

D. Revisions, Theory and Code Testing
In this section, we briefly outline some of the changes which have been made both to our presentation and understanding of the theory concerning ray tracing methods since previous work [22,23].We also describe revisions and improvements to the underlying codes.The backward ray tracing code is implemented in Mathematica, whilst the forward ray tracing code uses Julia, with both requiring parallelization on computing clusters.
We have carried out extensive tests of the methods used in this paper.We have verified at a ray-by-ray level that the two codes produce the same photon worldlines according to Eqs. (2) by backward propagating rays with initial conditions inferred from forward propagation from the source.We carried out this comparison for two equivalent forms of the dispersion relation presented in [56] and [57].Note that two important realizations emerged in attempting to generate agreement at the ray-by-ray level.First, fundamental constants (such as the speed of flight) need to be defined equivalently to a high level of precision, and second, the image plane in the backward ray tracing approach needs to be placed at sufficiently large distances such that the vacuum approximation is valid (which is an order of magnitude larger than that used in the original work of [23]).In addition, highprecision ODE solvers must be applied in both cases.
We checked that rays which escape to a particular direction (θ, φ) in the sky emanate from the same sourcing regions as those given by back-propagating from that same part of the sky (see Fig. 2).We confirmed that the total power computed from forward propagation (which is trivially consistent with Eq. ( 12)) is equivalent to the power inferred by backward tracing to within a few percent.This was also done across a full range of observing angles, and for different photon frequencies.
We have also verified mathematically, using analytic expressions found in [23,58], that, for an isotropic plasma and a spherical emission surface, the analytic expressions Eq. ( 11) for the power from backward tracing, and Eq. ( 12) for forward tracing, are equal.Note we even verified this in curved spacetime for a Schwarzschild metric.
We also now implement correct kinematic matching of axions and photons at k a = k γ , and incorporate a variety of curved-spacetime effects in each code.In addition, we make use of the new conversion probability derived in Ref. [54].In particular, the methods implemented here provide a much deeper understanding of phase space, and verify numerically that the divergence in the probability is regulated by phase space measure.
In the backward ray tracing code, we have also implemented an improved event-location routine to identify where photons are produced.Additionally, we now choose to describe ray tracing in terms of the photon distribution function f γ , rather than the radiant intensity I γ of the photon (as was done in previous work [23]).These are equivalent and related by f γ = I γ /(ω 3 n r ), where n r is the ray-refractive index described in [59].In particular, our discussion of phase-space in this work and Ref. [54] now allows us to properly understand the connection between the forward and backward ray tracing methods of [23] and [22].As discussed in the sections below, the original work of [23] did not include the possibility of photon reflection -ray tracing was terminated when the first stopping condition, m a = ω p , was encountered.This can prove problematic for photons sourced near the neutron star surface, since a sizable fraction of these photons undergo rapid oscillations between large plasma gradients.This has now been corrected by tracing photons for a much longer period of time, ensuring they have entered, and then escaped, the magnetosphere.Finally, in the backward ray tracing code we now implement a coupled set of first order ODEs in Eq. ( 2) rather than a single second order ODE given by Eq. ( 20) of Ref. [23], which is only possible in an isotropic medium.
The forward ray tracing code of [22] also included one notable update of the surface sampling procedure, which now allows one to appropriately treat the resonance condition (defined by the location where k γ µ = k a µ , rather than ω p = m a , with the plasma frequency defined as ω 2 p ≡ 4πn e /m e ).Note that the modified resonance condition implies that photons are not sourced by a single surface, but rather a foliation of surfaces, each being locally defined by the relative orientation of the axion momentum and the magnetic field; as such, the new Monte Carlo sampling procedure picks out a single surface within the foliation by first defining the orientation of ka , and then proceeding to uniformly sample the associated two-dimensional surface (as outlined in [22]).

III. AXION-PHOTON CONVERSION AND MAGNETIZED PLASMAS
We turn now to the application of ray tracing to the production of photons from axions in magnetized plasmas.In Ref. [54], it was shown that the axion collision integral is equivalent to where |M a→γ | 2 is an effective squared-matrix element for axion-photon conversion, f a is the phase-space distribution of axions, and the axion energy E a is defined by E 2 a = k 2 + m 2 a .Note the 3-momentum integral can be performed trivially so that upon substituting Eq. (15) into Eq.( 12) we find the power is defined as Using the integral identity for any smooth function G where dΣ is the area element on the surface defined by G(z) = 0, the integration of the delta function in Eq. ( 16) can then be performed with respect to the spatial integral, giving where Σ k is a surface defined by the set of points so that over a continuum of k values, the Σ k define a foliation of distinct surfaces.As a special case, note that for an isotropic (unmagnetized) cold plasma, we have and p so that the surface on which the energies become degenerate is given simply by ω p = m a , in which case the foliation collapses to a single surface Σ, independent of k.In general, however, the dispersion relation (and hence E γ ) can also depend on the direction of k, generating a distinct surface for each value of the momentum; as described below, this is the case in a strongly magnetized plasma, as found near the surfaces of neutron stars.
Returning to the expression Eq. ( 18) for the axion collision integral, we can then insert a factor of v a and cos θ n in the denominator and the numerator, where cos θ n is the cosine of the angle between the phase velocity v p = k/E γ and the normal to the surface surface Σ k .Here v p is the phase velocity of the axion, equal to that of the photon at the resonant surface.We then also use the relation This allows us to write where is a conversion probability, which, as shown in Ref. [54], defines the ratio of axion to photon phase-space densities at the resonance: The general form of P aγ of course depends on the medium in question.For strongly magnetized plasmas (as is relevant for this work), the photon mode of interest is the Langmuir-Ordinary (LO) mode (see [22]); in non-relativistic plasmas, the energy of the LO mode satisfies where θ B is the angle between B and k.Using this expression for the photon energy, one finds a conversion probability given by [54] Eq. ( 20) thus relates the power P flowing through a bounding surface A to a term involving an integration over axion phase-space at the source.The surface A can then be chosen to be the sphere at infinity, thus relating the power at source to the power measured by observers far away.
As a reminder, the general idea behind the forward propagation approach is to directly compute the right hand side of Eq. ( 20) locally at the source using Monte Carlo integration, and then use ray tracing to relate the local properties of emission to the asymptotic properties by exploiting the fact that f γ is conserved along trajectories once outside the source (Eq.( 5)).Alternatively, the backward ray tracing approach tracks the conserved phase-space density f γ along photon world-lines until k γ = k a , at which point one assigns a value f γ = P aγ ×f a .Note that in the case of backward ray tracing, trajectories must be traced until they have passed through and fully escaped the gravitational potential of the neutron stars -this is a consequence of the fact that each photon worldline can encounter multiple level crossings, and the total photon phase space is the weighted sum of each of these.

IV. GENERALIZATION TO CURVED SPACETIME
In this section, we discuss how to generalize the results presented in previous sections to curved spacetime.We focus in particular on the case of photon production via axions in strongly magnetized plasmas in the presence of gravitational fields.This is particularly relevant in the vicinity of neutron stars, which are the most compact stars in existence, with the ratio of the Schwarzschild radius to the neutron star radius being r s /R ∼ 0.3.
The importance of incorporating gravitational effects is multi-faceted.Firstly, gravity alters photon trajectories relative to flat spacetime, acting as an attractive force which counteracts repulsive refraction due to plasma gradients [23,53,57,58].Next, axions fall towards the star along geodesics of the stellar spacetime metric.Energymomentum conservation demands that axions and photons should have matching 4-momenta at the point of creation, k a µ = k γ µ .For kinematic self-consistency, one should therefore demand that gravitational effects are incorporated self-consistently in the axion and photon dispersion relations.By doing so, one changes the geometry of the foliation Σ k in Eq. ( 20) relative to flat space, with the innermost and outermost surfaces in the foliation separated by a greater distance than in flat space.In addition, metric contributions will enter the infinitesimal area dΣ k appearing in the collision integral.
A related consequence is the fact that gravity increases the density of axion near the star [36] -this arises because geodesic bundles of axions become focused in regions of strong gravity (causing a greater number of axions to cross the resonant conversion region per unit time), and can lead to a sizable enhancement of the amplitude of the radio flux.This effect has been included in previous work via a renormalization of f a , but in generalization must be done self-consistently with the aforementioned effects.
Turning to the production process itself, the analytic expression for the conversion probability in flat space is intimately connected to phase-space, kinematics, dispersion relations, and worldlines of photons and axions [54].Crucially, divergences emerging in the conversion probability are shown in flat space to be canceled by compensating terms in the phase-space integration in Eq. ( 20).As a result, a fully self-consistent generalization of kinetic theory and ray tracing to curved spacetime is also relevant for properly describing the efficiency with which photons are produced.
Finally, as pointed out in [44], strong magnetic fields and larger values of g aγγ can lead to O(1) axion-photon transitions.Computing the radio signals in this scenario requires self-consistently tracking axion and photon trajectories through all resonant crossings (as well as the trajectories of axions and photons sourced from those resonances); although we do not go beyond perturbative production in this work, the techniques developed here lay the groundwork for such follow-up analyses.
In the remainder of this section, we outline how to incorporate curved spacetime effects into each step of the problem.

A. Magnetized Plasmas in Curved Spacetime
Let us begin by defining the generalized Hamiltonian that describes the covariant treatment of waves in a magnetized plasma.The covariant treatment of waves in a magnetized plasma on curved spacetime can be found in [56,57,[60][61][62] (the less general case of an isotropic plasma is studied in [53,58]).For a non-relativistic plasma, with fluid velocity u µ , phase space coordinates (x µ , k µ ) and background electromagnetic field strength tensor F µν and spacetime metric g µν , we find the Hamil-tonian is given by where and with where B αβ = h αµ h βν F µν , ε µναβ is the Levi-Civita tensor, h µν is the projection operator onto the rest-frame of the plasma, ω L is the Larmour frequency, and ω and p µ are the frequency and 4-momentum in the rest-frame of the plasma.We can also define the magnetic field via [57] Using these relations, we have where the first two identities follow from the antisymmetric properties of the Levi-Civita tensor when contracted with multiple u µ , which also implies u•ω L = 0.It therefore follows that the Hamiltonian can be re-written as We can now go to the strong magnetization limit relevant for neutron star magnetospheres.This is characterized by the limit ω L ≫ ω, ω p .In that limit, the dispersion relation H = 0 can be equivalently realized by the Hamiltonian where B, which is identical to the dispersion relation in [57] in the limit of a non-relativistic plasma 1 .
We can also define an angle between B µ and k µ where For stationary plasmas and spacetimes, we can choose the fluid to be instantaneously at rest with respect to the space-time metric.For example, for a Schwarzschild metric where A ≡ (1 − r s /r) and r s ≡ 2GM , we can define generalized phase-space coordinates (t, r, θ, φ) and (k t , k r , k θ , k φ ) and a fluid velocity In this case, 2 , with |k| 2 = g ij k i k j and where i, j ∈ (r, θ, φ) label the spatial components of tensors.We also have B 0 = 0. Hence in this frame, we can write the dispersion relation as which gives with the lower and upper signs corresponding respectively to the Alfén and LO mode of Eq. (23).Hence, to construct rays in curved spacetime, we solve Hamilton's equation Eq. ( 2) with Eq. ( 34).In the next subsection, we describe how to generalize forward ray tracing to curved spacetime, focusing in particular on how to generalize the area measure Σ k appearing in Eq. (20).

B. Radiative Transport in Curved Spacetime
In this section, we discuss how to generalize the kinetic theory for the production of photons to curved spacetime.In particular, generalizing forward ray tracing to curved spacetime necessitates generalizing the expression Eq. ( 12), or, more specifically to axions, the surface integrals in Eq. (20).In flat space the infinitesimal power flowing per unit time dt can be written as where we have explicitly written To generalize this to curved spacetime, we recognize that dtdΣ k is to be interpreted as the spacetime volume element on a 2+1 dimensional submanifold.For a metric of the form Eq. ( 36), the temporal piece of the volume element is generalized by taking The generalization of dΣ k to curved spacetime is given by taking where |h k | is the determinant of the pull-back metric h k corresponding to embedding the surfaces in Eq. ( 19) in the background spacetime of g µν .The metric h corresponds to a 2D spatial surface, which we can label by their coordinates (θ, φ) with the radial coordinate of the surface corresponding given by where the k subscript reminds us that these surfaces form a foliation parametrized by k.The h k = h k (θ, φ) can be read off from the spatial line-element ds 2 of the Schwarzschild metric, defined by Using the chain-rule we can project this line element onto the surfaces corresponding to Eq. ( 43), giving Inserting this result into Eq.( 44), one can directly read of the elements of h αβ ; the determinant is then given by where, for compactness we have dropped the k subscript on r.This provides the generalization of dΣ k to curved spacetime.Notice that the expression in Eq. ( 47) gives precisely the gravitational redshift factor that accounts for energy loss of photons propagating out of the gravitational potential.
Putting everything together, we find the power integral in curved spacetime can be expressed as Notice that here d 3 k corresponds to the volume element of any orthonormal tetrad on the manifold used to define momentum space locally.It can be rotated into any other orthonormal tetrad basis using a local SO(3) rotation, which does not change the momentum volume measure, because the determinant of the SO(3) matrix is 1.The angle θ n should be defined covariantly via where n µ is the unit normal to the surface Σ k .

C. Gravitational Focusing of Axions
In general, the asymptotic axion distribution is expected to follow a Maxwell-Boltzmann distribution with a non-relativistic momentum dispersion k 0 = m a v 0 with v 0 ∼ 220 km/s, so that where we have assumed for simplicity that the pulsar is at rest with respect to the rest frame of the asymptotic axion distribution.
To compute the phase-space density near the star, we first make use of conservation of energy, where |k ∞ | is the asymptotic momentum of the axion at infinity and s is an axion worldline parameter.Note that knowledge of the asymptotic solutions, Liouville's theorem, and kinematic constraints like Eq. ( 51), are in general not sufficient to solve for the functional form of the distribution function globally -that is precisely why ray tracing is employed.However, in the present context, owing to the uniform and isotropic boundary conditions, we can construct solutions using simple conservation arguments, which effectively allow us to circumvent ray tracing for axions, allowing us to determine their distribution analytically.We know, by Liouville's theorem that and using the energy conservation equation (i.e.Eq. ( 51)) on the right-hand side, and Liovuille's theorem on the left hand-side, one has Since this holds for any ray worldlines k(s) and x(s), this must therefore give the global solution The reason we are able to determine the distribution of axions given boundary conditions on a surface at infinity, without needing to ray-trace, is a special case arising from the high degree of symmetry in the problem.Indeed, the solution we have constructed is nothing more than f a ∝ exp(−H a /k2 0 ) where H a (x, k) = k 2 /(2m)−2GM m a / |x|, with the normalisation fixed by boundary conditions.This is a trivial solution to the collisionless, stationary, axion Boltzmann equation (55) where the left-hand side gives the Poisson bracket.
Clearly the solution we have constructed satisfies this equation since f a = f a (H a ) gives a trivially vanishing poisson bracket.For more complicated distributions, such as, e.g. , the scenario in which axions are predominantly confined to gravitationally bound miniclusters, these simple solutions do not apply, and one must also apply ray tracing to the in-fall of axions, see, e.g. , [49].
Returning to our main discussion, we see the distribution in Eq. ( 50) gives axions with characteristic asymptotic momentum O(k 0 ), for simplicity, we therefore choose to evaluate the ray tracing algorithms by isolating a single momentum k = k 0 , for axions, i.e. , lim that is, a monochromatic, zero-width distribution with characteristic momentum k = k 0 .This procedure can be shown to yield nearly an identical density enhancement as the Maxwellian distribution 2 , but significantly reduces the sampling required in the backward ray tracing approach since only a single photon frequency must be sampled in the image plane 3 .We return in the later sections to illustrate the impact of this simplifying assumption.Following similar arguments to above, we determine the global distribution to be In particular, integrating this equation over momentum space, gives where k 2 c = k 2 0 + 2GM m 2 a /|x| so that the axion number density is enhanced by the ratio of the asymptotic and escape velocities, due to focusing of the gravitational field.Putting this together, we can express the axion distribution function in terms of the local number density of axions: where ω c = m 2 a + k 2 c and n c and v a are the axion number density and phase velocity, respectively, both evaluated at source.
One might wonder why it is that this gravitational focusing of axions is not reversed when photons climb out of the gravitational potential.The reason for this, is that photons do not experience the same refractive index as axions as they exit the potential -they are in plasma, not vacuum, and asymptote to a refractive index n γ → 1 as they move away from the star, whilst axions approach n a → v 0 .This is especially clear in the isotropic plasma, where the distribution function f = I/(ω 3 n 2 r ) is conserved along rays [23,53] -here, we have defined n r , the refractive index, and I is the radiant intensity.Using conservation of f thus implies where in the second to third equations, we used f γ = P aγ f a .Putting this together, we see firstly that a so that gravitational focusing increases the intensity of axions.Secondly, we have where we used, n ∞ r,γ = 1 and n ∞ r,a = v 0 .We see that (in this illustrative isotropic case) the extent to which the gravitational focusing of axions is undone as photons exit the plasma potential, is precisely captured by the ratio of the asymptotic refractive indices of the two species.We notice that by equating each of the equations Eq. ( 60) and Eq. ( 61), the volume element change [23] from gravitational infall of axions, ∝ (ω c /ω ∞ ) 3 ≃ (1 − r s /r) −3/2 , is undone when photons exit the plasma, since this term is purely gravitational.However, the ratio of the refractive indices in Eq. ( 61) quantifies the difference of the two potentials through which axions enter and photons exit the plasma.Similar logic holds for the anisotropic medium.

D. Conversion Probability in Curved Spacetime
The conversion probability for axion-photon conversion in 3D magnetized plasmas has recently been computed in flat space [54], yielding the result given in Eq. ( 24) (an expression which is valid in a (quasi) stationary background).Formally, given we are including curved spacetime effects in our analysis, self-consistency implies we should also generalize the production rate in Eq. ( 24) to curved spacetime.
One important reason for doing so is that divergences in the conversion probability should be regulated by the phase-space measure.More specifically, in flat space the FIG. 4. Equivalence of Forward and Backward Ray Tracing.Period-averaged differential power as would be seen by an observer viewing at an angle θ with respect to the rotational axis, computed using forward (solid) and backward (dashed, points) ray tracing.Results are shown for the isotropic ('iso') and anisotropic ('aniso') dispersion relations and an axion mass of ma = 10 −5 eV.
surface normal of Σ k is parallel to ∇ x E γ , so that divergences occurring when the phase-velocity v a is perpendicular to ∇ x E γ are regulated by the flux-like projection v p • Σ k , where Σ k is the directed surface element.
The generalization of the conversion probability to the case of an isotropic plasma in curved spacetime is straightforward.The dispersion relation is given by which gives From Eq. ( 19), this implies that family of surfaces Σ k collapses to a single emission surface Σ, on which which is independent of k.In that case, the unit normal to surface Σ is given by and the angle in Eq. ( 49) between this normal and the phase velocity, is given by We also generalize the conversion probability in the isotropic case by modifying the flat space results according to where in the isotropic case we have Collectively, this yields a conversion probability given by If this is inserted into Eq.( 48), the divergence occurring in the conversion probability Eq. ( 69) (arising when k • ∂ω 2 p → 0), is regulated by the cosine angle in the measure of Eq. ( 48), such that where p ).This is manifestly convergent.Generalising Eq. ( 69) to the anisotropic case in curved spacetime is non-trivial, as this involves formulating kinetic theory for photons in curved spacetime [63] -note that this is similar to what is done for the case of scalars in Ref. [64].Such a generalization is clearly important, but lies beyond the scope of the present work.Instead, for the purpose of making progress, we employ a phenomenological generalization of the conversion probability similar to what is presented above, taking as before the substitution to arrive at an ansatz for the anisotropic conversion probability in curved spacetime of the form Similarly, we generalize the angle in curved space by taking which guarantees that the differential power computed using the forward propagation approach is convergent.However, since this is only a partial generalization to curved spacetime, divergences are not canceled in the backward ray tracing approach (see Appendix A for a more detailed discussion of the origin of these divergences).In order to regulate these divergences, we introduce an IR cutoff in the backward ray tracing approach on the effective conversion length, defined by L c ≡ k • ∂E γ , which serves to exclude any contribution with L c > 1km4 .This should not be confused with any type of formal regulatory procedure; rather it is a means to excise points which lead to divergences that would otherwise be regulated with a formal generalization of the conversion probability to curved space time.

V. AXION DARK MATTER DETECTION WITH NEUTRON STARS
Having outlined the generalized ray tracing procedure, we can now return to the problem at hand -namely, the application of these approaches to radio searches for axion dark matter near neutron stars.We choose to focus in particular on spectral lines arising from a smooth background distribution of axion dark matter, which are among the most well-studied indirect probes of axions in these environments [21, 22, 35, 36, 38, 39, 41, 43-45, 48, 50, 65].One should bear in mind, however, that these tools are more broadly applicable to similar searches, such as those looking for transient lines from the encounters of miniclusters and axion stars with neutron stars [40,42,46,47,49], and broadband radio searches generated from a locally sourced population of axions [25,51,52].
We begin by establishing a set of fiducial parameters which define our baseline model, given by an axion mass of m a = 10 −5 eV, an axion-photon coupling g aγγ = 10 −12 GeV −1 , a neutron star mass of M NS = 1M ⊙ , a neutron star radius of r NS = 10 km, a surface dipolar magnetic field of B 0 = 10 14 G, a neutron star rotational frequency of Ω NS = 1Hz, and a misalignment angle θ m = 0 radians.In what follows we normalize the axion distribution to asymptotic number density n a,∞ = 0.45 × m −1 a GeV cm −3 .Deviations from these parameters below are always explicitly stated.
Each of the examples in this section is illustrated assuming a perfectly dipolar magnetic field and a GJ charge density, defined by a e ± charge density where θ m is the misalignment angle between the rotational and magnetic axis, and the time dependence of the plasma due to rotation has been embedded in the term m • r = cos θ m cos θ + sin θ m sin θ cos(Ωt), where θ is a polar angle given by taking rotational axis as the north pole.Throughout this paper, we choose coordinates such that θ also gives the angle between line of sight to an observer and the rotational axis.The GJ charge distribution is expected to be a reasonable approximation across most of the closed field lines of standard pulsars, and near the surfaces of dead neutron stars [38,66,67].We revisit this assumption when applying our results to the galactic center magnetar, as the charge densities near such objects are expected to be highly enhanced with respect to the GJ model [68].Before illustrating how each assumption and free parameter impacts the projected radio signal, we begin by illustrating the agreement between forward and backward ray tracing, using the formalism developed in the preceding sections.In Fig. 4 we show the differential power (averaged over the rotational period of the neutron star) generated from resonant axion-photon transitions, and computed using either an anisotropic, or an isotropic, dispersion relation in curved spacetime.We show results for an axion mass of 10 −5 eV, but have also confirmed agreement at other masses.These numerical results confirm the equivalence of forward and backward ray tracing, in accordance with the theory laid out in Sec.II.This demonstrates explicitly the equivalence of sampling the collision integral and explicitly evaluating the asymptotic flux via back-tracing, with each scenario manifested in the right-and left-hand side of Eq. ( 6), respectively.
Before continuing, it is worth highlighting that achieving such a high level of agreement requires not only a unified formalism (see Sec. II), but also high-precision numerics, with the results being strongly sensitive to a variety of different factors, such as: the initial conditions of the photon in the backward ray tracing procedure, the value of fundamental constants (we find maintaining consistency across many decimal places is typically required), the precision of the ODE solver, etc.We now go on to illustrate how varying the size of different physical phenomena discussed in this and proceeding sections affect observational signatures from axion dark matter conversion in neutron stars.

A. Physics of ray tracing
The goal of this subsection is to illustrate the importance of the various physical effects within our models, which until now, have not yet been self consistently incorporated.This includes the combined impact of plasma anisotropy and curved space on photon propagation, the effect of varying the neutron star mass and radius, the importance of including multiply reflected photons, and the impact of simplifications to the asymptotic velocity distribution of axions prior to in-fall.In Appendix B we also examine the extent to which gravity can be encoded in the initial conditions of the photons, and the importance of imposing the proper kinematic matching condition 5 .In each case, we analyze the impact of including FIG. 5. Effects of Plasma Anisotropy.Period-averaged differential power as would be seen by an observer viewing at an angle θ with respect to the rotational axis.Results are shown for two choices of the axion mass (ma = 10 −5 and 10 −6 eV).The isotropic ('iso') case combines a dispersion relation g µν kµkν + ω 2 p = 0 and conversion probability P iso aγ of Eq. ( 69).The anisotropic ('aniso') uses the dispersion relation of Eq. ( 39) and the conversion probability P aniso aγ in Eq. ( 72).All results include gravity.Plots are generated using the fiducial parameters of the main text.FIG. 6. Neutron Star Mass (Renormalized).Same as Fig. 7, but setting the axion-photon conversion probability to 1 for all rays, and re-normalizing the weights of each ray by the factor Ri given in Eq. 75 (this re-normalization ensures the sky-integrated power is fixed across all neutron star masses).Result is shown for an axion mass of ma = 10 −5 eV (left) and 10 −6 eV (right).and/or neglecting an effect by computing the rotationperiod averaged differential power ⟨dP/dθ⟩ as a function of the viewing angle, θ.We emphasize that it is not always straightforward to isolate individual effects, since a self-consistent treatment usually involves not just one, but many modifications.As an example, one can consider that gravity plays a role not only in modifying the trajectories of individual rays, but also enters the photon initial conditions (which in turn enter the conversion probability, producing an effect both on the weighting of photons and their propagation) and axion number density.We clarify below when isolating an effect breaks self-consistency.
In the next two subsections, we begin by addressing the question of the relative importance of treating the interpretation is something of a subtlety that could cause confusion for less familiar readers.
anisotropy of the plasma and curve spacetime, the former having been treated in [22] and latter in [23].

Plasma Anisotropies
We now examine differential power for an axion mass of 10 −5 eV and 10 −6 eV, assuming either an isotropic or anisotropic plasma, which correspond to the dispersion relations of Eqs. ( 62) and (39), respectively.In each case, we also adopt a conversion probability that is self-consistent with the choice of dispersion relation; namely, we use P iso aγ (Eq.( 69)) for the isotropic plasma, and P aniso aγ (Eq.( 72)) for the anisotropic plasma.Modifying the photon dispersion relation only redistributes outgoing photons in phase space, keeping modifying how power is distributed on the sky while maintaining a fixed sky-integrate power (assuming photon absorption can be neglected).Meanwhile, changing the conversion probability between the isotropic and anisotropic cases affects the overall normalization of the power, in accordance with Eq. ( 20).In particular, this leads to a larger overall power output for an anisotropic plasma.For comparison, we also plot a comparison between the anisotropic and isotropic scenarios, but setting the conversion probability in both cases to 1 for all rays -this allows us to isolate the impact of plasma anisotropy on the evolution of photon trajectories.The effect of gravity is included in all cases.
The results for isotropic and anisotropic plasmas are shown in Fig. 5 (the left panel includes the conversion probability, while in the right panel it is set to unity).We see that for the parameters chosen, the axion mass plays the dominant role in changing both the morphology and the normalization of the power profile, with the plasma anisotropy driving a more subtle, although nonnegligible, change of shifting power away from the poles and equatorial plane -this latter effect arises from the θ B dependence in the photon dispersion relation.At larger masses (m a = 10 −5 eV), we see that differential power varies by an order of magnitude across θ ∈ (0, π), while the effect of assuming an isotropic plasma can induce variations (for fixed values of θ) at the order of magnitude level, although it is worth noting that the average difference is not so pronounced.For an axion mass of 10 −6 eV (where the conversion surface extends to larger characteristic radii), the differential power across the sky can vary by more than three orders of magnitude, with the effect of plasma anisotropy introducing as much as a one order of magnitude shift in the power (although only a certain viewing angles).
A range of approximations [21,22,36,41,44,50,69] have been used throughout the literature to model plasma physics in the context of axion searches, both at the level of dispersion relations and production itself.In this work, we are now in a position to collate these various approximations and make some assessment of their relative importance from a observational point of view.More explicitly, given the results reported in Fig. 5, the question arises as to the significance of the differences between simplified isotropic scenarios and more complete anisotropic plasma descriptions, the latter combining anisotropic effects both in dispersion relations and conversion probabilities.
One should bear in mind that the results of Fig. 5 concern aligned rotators (θ m = 0), whilst in general neutron stars are non-aligned with θ m ̸ = 0.In the case of a frequency domain analysis (in which one uses timeintegrated observations), the primary effect of inducing a misaligned rotation axis is to isotropize the periodaveraged flux (see e.g. the differential power curves in [22]).

Curved Spacetime Plasma Effects
Having discussed the effect of plasma anisotropy, we turn now to the role of curved spacetime effects.We begin by focusing on the role of gravity in altering the photon dispersion relation (and hence ray-propagation) [70], and then return to the impact of gravity on the local intensity emitted at the conversion surface.
In order to understand the effect of gravity we begin by running the ray tracing analysis using a variety of different neutron star masses ranging from M NS = 0.5 M ⊙ to M NS = 2.2 M ⊙ .We begin by isolating the effect of gravity on the propagation of photons by re-scaling the photon weights in such a way that the sky-integrated asymptotic power remains fixed.This is done by setting the conversion probability of each ray to one, and re-scaling each photon by a factor where r i is the initial radius of the photon at the resonance.Note that the choice of normalization of the axion number density and photon momentum are arbitrary, and thus one should not attempt to interpret the normalization of the differential power as being physical.The result of this procedure is shown for two choices of axion mass in Fig. 6.For m a = 10 −5 eV (i.e. when the conversion surface is close to the neutron star), the effect is most pronounced near the magnetic poles, leading to a variation in the differential power up to a factor of ∼ 5 -at other viewing angles, however, the effect is typically no more than factor of two.As expected, the effect of increasing the neutron star mass leads to an isotropization of the radio flux [23].For smaller axion masses, the conversion surface shifts away from the neutron star and the effect of gravity on the propagation of rays is suppressed; specifically, for the case of m a = 10 −6 eV, we find gravity only modifies ray propagation at the O(10%) level.In the Appendix, we attempt to understand the extent to which previous approximations adopted in [22] are valid; in particular, we show that embedding the effect of gravity into the initial conditions, but neglecting gravity during the propagation, tends to induce negligibly small errors in the differential power (we caution, however, that this is an empirical result, and is not guaranteed to hold in all contexts).We now return to analyzing the full effect of gravity, which enters not only the propagation of rays but also the local power emitted from the resonant conversion surface.We plot in Fig. 7 the impact of increasing the neutron star mass from 1 M ⊙ to the maximally allowed value of ∼ 2.2 M ⊙ for two choices of axion masses (unlike Fig. 6, here we maintain the conversion probability and the appropriate weights of the rays).Fig. 7 illustrates that the neutron star mass plays a non-negligible role in the prediction of the radio flux, and should likely be included in future modeling.
In Fig. 8 we illustrate the impact of varying the neutron star radius, assuming a neutron star mass of either 1 M ⊙ or 2.2 M ⊙ .Here, we see that the characteristic neutron star radius adopted in previous works tends to lead to an underestimation of the radio flux by a factor of ∼ 2 -this effect is nearly uniform over the sky, and arises primarily from the fact that the increase of radius (at fixed surface magnetic field strength) tends to increase the net surface area over which resonant axion-photon transitions take place.
The initial ray tracing performed in Refs.[22] and [23] had observed notable deviations in the inferred time profiles induced by resonant axion-photon mixing (note that time variation is directly related to the variation of the differential power with viewing angle, θ).Ref. [23] had attributed this difference to the effect of gravity -it was argued that gravity tends to isotropize the signal, washing out the strong variation observed in the time profiles of [22].Here, we show that this conclusion is not correct; instead, the variation observed in Ref. [22] was larger due to three contributing factors: (i) [22] focused on lower mass axions, which have larger time variation due to the larger characteristic conversion surface, (ii) Ref. [23] had not included multiply reflected photons, which tends to increase the time variation (see Fig. 9, and the following subsection), and (iii) the anisotropy of the plasma at low masses can further enhance the variation in the time-domain.

Multiply Reflected Photons
The backward ray tracing method [23] is essentially a standard approach to the problem of radiative transfer [59].Photon production is calculated via an integral along the photon worldline, as depicted in Eq. ( 4).As the photon worldline is back-propagated, it can encounter multiple sites of photon production from axions.This happens wherever axion and photon dispersion relations become degenerate, corresponding to a family of discrete values λ i of the worldline parameter at which photons are produced.In the original work [23] we only considered the first level crossing 6 and λ = λ 1 .In the present work, we include contributions to the asymptotic value of the photon distribution f γ (λ → ∞) from all production sites along the photon worldline, so that This effect is especially pronounced within the throats of the magnetosphere, where the worldines of photons experience multiple reflections, causing them to repeatedly encounter level-crossings.
In Fig. 9, we illustrate the impact of including these multiply reflected photons for two different axion masses.Here, one can see that accounting for multiply reflected photons tends to induce a factor of ∼ 2 enhancement in the differential power for viewing angles that are roughly aligned with the throats of the magnetosphere.This effect is enhanced for smaller axion masses, where the throats become more prominent.

Asymptotic Velocity Distribution
The examples provided thus far have assumed for simplicity that the asymptotic axion speed distribution can be treated as a delta function fixed to |⃗ v| ∼ 220 km/s, corresponding to a monochromatic photon signal.This is the assumption adopted in [23,39,50], while [22,44] treated the full Maxwellian distribution.In Fig. 10, we show that for both axion masses of interest, the simplification of neglecting the width of the asymptotic energy distribution tends to induce negligible variations in the inferred power.Importantly, however, this statement is only valid when the neutron star is assumed to be at rest with respect to the galaxy, and may not hold for strongly boosted neutron stars.

VI. GALACTIC CENTER MAGNETAR
The Galactic Center magnetar SGR J1745-2900 has an inferred dipolar field strength of 1.6 × 10 14 G, a rotational period P ∼ 3.76s, and has a two-dimensional projected distance from the Galactic Center of ∼ 0.17pc [71], making it a promising target for axion searches.SGR J1745-2900 is in fact frequently adopted as a benchmark for developing projected sensitivities for future observations [21,22,36,39,44], and various groups have attempted to place constraints using existing radio observations (see e.g.[48]).All analyses to date have assumed that the charge distribution can be described by the GJ model, which is derived by determining the minimal corotating charge density needed to screen the product of Magnetars, however, are expected to have magnetospheres which differ markedly from the standard pulsar population, with charge densities greatly exceeding the minimal co-rotating density.Bearing in mind that the current understanding of magnetar magnetospheres is far from complete (see e.g.[72,73] for recent reviews outlining the current understanding of magnetars), we describe below the expected properties of these systems, and use the techniques outlined in the previous section to revisit sensitivity estimates that could be achieved with existing and future data.Needless to say, the validity of these projections hinges upon a variety of rough approximations, and thus they should be treated with some skepticism; nevertheless, these projections serve to answer the important question of whether magnetars could be useful in the future to probe regions of axion parameter space which are conventionally inaccessible to other indirect axion searches.
Magnetars are very young and highly magnetized objects whose bright X-ray emission is powered magnetically, rather than rotationally (implying the energy losses exceed the spin-down power of the neutron star), see e.g. , [73].These objects are strongly variable and exhibit a broad range of phenomenon ranging from X-ray bursts, glitches and antiglitches, non-uniform spin-down, giant flares, and even fast radio bursts (see, e.g. , [74] for a recent association between fast radio bursts and a nearby magnetar).
It is believed that the variable magnetar activity is driven by the evolution of the ultra-strong magnetic field -a shifting in the structure of the magnetic field internal to the neutron star deforms the crust, inducing a strong shearing of the crust which subsequently drives electric currents into the magnetosphere.The net result is a twisted, nearly force-free, magnetosphere threaded by strong electric currents ⃗ j ∼ ∇ × ⃗ B. Phenomena like flaring can then be explained, e.g. , via the presence of instabilities that appear as the magnetic twist exceeds a critical threshold (see e.g. , [68,72,73,75]).
In the context of axions, there are two important dis- tinctions between the magnetospheres of magnetars and those of standard pulsars.First, the assumption of a purely dipolar field is broken -the twisted magnetosphere induces a complex topology to the magnetic structure, and the magnetic field strength may easily exceed the inferred dipolar value [76][77][78][79][80].The second important distinction comes from the fact that the minimum charge density flowing from the twisted field configuration greatly exceeds the minimum co-rotational charge density identified in the GJ Model.This can be seen by noting that the force-free condition, ⃗ E • ⃗ B = 0, implies a minimal charge density near the neutron star given by [75] where the first term is the GJ charge density, and the second ρ tws is the minimal charge density required to support a twisted field configuration; note that for nontwisted field configurations (as, e.g., is the case with a standard dipolar field) ∇ × ⃗ B = 0, and thus ρ ∼ ρ GJ (at least in the closed zone near the neutron star, where the force free condition is expected to be satisfied).The twisted field configuration expected to arise in the magnetospheres of magnetars is supported and stabilized by the presence of a strong electromagnetic current |j| ∼ ∇ × ⃗ B ≫ ρ GJ , which is sourced from e ± pair production processes near the neutron star.The electron/positron number density n e ± can be derived by writing the local charge and current density as where n ± and v ± are the number density and velocity of the e ± , and then directly solving for n + + n − ; for a semi-relativistic plasma [68] with |j| ≫ ρ this quantity is roughly bounded to be n ± ≳ 2| ⃗ j|/e [75,81].
Ref. [75] investigated the charge distribution in the context of a globally twisted dipole configuration, showing this leads to a characteristic current on the order of (see also [81]) Introducing a charge multiplicity factor λ (defined with respect to the minimal charge density ∼ 2| ⃗ j|/e) and assuming the magnetic field strength is purely a function of radius, decaying as a dipolar magnetic field |B| ∼ B 0 (r NS /r) by a factor n e /n GJ e ∼ 10 3 × λ, which is in agreement with inferences of the charge densities obtained from observations of the resonant cyclotron absorption of X-rays [82].Notice that Eq. ( 80) implies the magnetospheres of magnetars may support resonances for axions with masses near, and above, 10 −2 eV.
One of the fundamental questions sitting at the forefront of the field for many years is how such strong currents are sustained in magnetar magnetospheres.Recent work on the electrodynamics in super-QED field strengths has shown that the hard X-ray spectrum extending to energies E ≳ 10 keV observed in many magnetars can arise from a highly collisional semi-relativistic plasma with a characteristic density 10 − 20 times larger than the minimal current density given in Eq. ( 79).These enhanced densities are sustained via a combination of ohmic heating and pair creation, and may be necessary in order to explain a number of observed phenomena including rapid X-ray brightening, concentrated thermal hotspots, and thermal X-ray emission [68].
The goal of this section is not to provide an accurate description of axion conversion in magnetar magnetospheres, but merely to point out the extent to which the axion searches performed in, e.g. , Ref. [48], are modified when more realistic assumptions are adopted.In this vein, we take four fiducial models, which, in spite of their simplicity, are expected to give some rough estimation of the sensitivity that radio experiments could have to axion-photon conversion in these systems.These models consider two distinct values of magnetic fields, one with the minimal value B 0 = 1.6 × 10 14 G (corresponding to the pure dipolar magnetic field) and the other with B 0 = 4 × 10 14 G (while this number is somewhat ad hoc, intended to show the impact of a moderate magnetic enhancement although an order one enhancement, we note that it could be seen as roughly consistent with the twist factor inferred in [83]), and two distinct values of charge multiplicity parameter λ (specifically, we take a uniform value of λ = 1, and λ = 20); note that the λ = 1 is intended to represent a minimal lower bound, with λ ∼ O(10) being closer to the value predicted by [68].We therefore consider four magnetar models, M1: λ = 1, B = 1.6 × 10 14 G, M2: λ = 20, B = 1.6 × 10 14 G, M3: λ = 1, B = 4 × 10 14 G, M4: λ = 20, B = 4 × 10 14 G.
In order to simplify the analysis, we treat the magnetic field as being purely dipolar (despite this assumption being inconsistent with the adopted charge density magnetic field structures); in general, the magnetic field structure should be obtained by solving the Grad-Shafranov equation (see e.g.[84]), however including such an effect is beyond the scope of this work.We note, however, that there are three effects which are expected to arise as one includes the geometric effects appearing in twisted configurations.First, O(1) angular factors shift the photon production efficiency, and therefore induce comparable shifts in the photon anisotropy (as compared with a dipolar field configuration).Next, the radial dependence of twisted magnetic fields is modified with respect to the dipolar case (falling as r −2+p , with p < 1 for a twisted field and p = 1 for the dipolar field) [75,85].Finally, and most importantly, the optical depth for highly twisted fields can be increased.For the moment, we estimate the optical depth using the dipolar configuration but caution that a more careful assessment of this effect may be needed in the future.
In order to be conservative, we choose to remove axionphoton conversion in open magnetic field lines, since active pair production and current flows in these regions are expected to require a more sophisticated level of modeling.These excised regions correspond to the white hatching in Fig. 11.Using Eq. ( 80) we see that field lines are characterized by the curves r/ sin 2 θ = L, where L gives the maximal radial distance of the field line from the neutron star.Open field lines are those which extend beyond the light-cylinder at r LC = Ω −1 NS , i.e., they satisfy L ≥ Ω −1 NS .When considering axion-photon conversion, we do not include photon production occurring on open magnetic field lines.
As a word of caution, dense return currents running along the closed field line could produce secondary effects not included by this 'cutting' procedure, such as the redirection and funneling photons near the conversion surface.These effects are not included here, but will be investigated in future work.
Resonant cyclotron absorption, occurring when the frequency of the radiation matches the cyclotron frequency ω = ω c , can be increasingly important for lowfrequency radiation emitted near highly magnetized objects (Ref.[22] had shown using the GJ model that the optical depth can be O(1) for the Galactic Center magnetar).Assuming the cyclotron resonance occurs at large distances from the magnetar where the trajectories can be approximated as radial, the optical depth is roughly given by [22] τ ∼ π 3 where it is understood that all quantities are evaluated at the point of resonance.
The characteristic charge densities spanned by Eq. ( 80) suggest that magnetars will be efficiently producing radiation across frequencies from O(1) GHz -O(5) THz, and thus we use a combination of sub-mm telescopes to develop our sensitivity projections.In particular, we adopt projections for current telescopes based on observations by the Green Bank Telescope (GBT), the Atacama Large Millimeter Array (ALMA), and the Stratospheric Observatory for Infrared Astronomy (SOFIA), which have broad bandwidth coverage over the O(10)GHz-THz regime.In addition, we include a separate set of projections for the Square Kilometer Array (SKA), which will cover frequencies from 50 MHz to 24 GHz.We compute the radio spectrum at seven fixed axion masses, corresponding to observing frequencies of ∼ 2.4, 20, 35, 100 GHz, 500 GHz, 1 THz, and 4 THz.SKA observations will cover the lowest two frequency bins, and assuming a system equivalent flux density SEFD = 0.098 Jy, a bandwidth of 10 −5 × m a , and an observing time of 10 hours, SKA will have sensitivity (at the 95% confidence level) to radio lines at each observing frequency of 5 and 2 µJy [107].We use the GBT telescope to establish sensitivity in the lowest three frequencies -here, we adopt a sensitivity (across all frequencies) consistent with the quoted 95% confidence upper limit used in the analysis of [48], S lim ∼ 0.3 mJy (computed using a bin width Left panel corresponds to sensitivity that could be achieved using current telescopes (namely a combination of GBT, ALMA, and SOFIA), while the right panel includes projected sensitivity from SKA.The vertical bars on each point reflect the 1σ variation in the inferred limits which are obtained by randomly sampling the orientation of Earth as defined with respect to the rotational axis.Shown for reference are current the QCD axion band (purple) [86], and constraints from globular clusters (gray) [87], CAST (dark blue) [88], axion haloscopes (gold) [89][90][91][92][93][94][95][96][97][98][99][100][101][102][103][104][105][106], pulsar polar cap cascades (light blue) [25], and GBT observations of the Galactic Center (green) [44].
of δf ∼ 28 MHz).At 100 7 and 500 GHz, we adopt a sensitivity for ALMA consistent with the quoted capability for 60 seconds of observations and a line width of 8 δf /f ∼ 10 −6 , which corresponds to S lim ∼ 5 mJy and 25 mJy for 100 and 500 GHz respectively [108].Note that GC magnetar has been observed up to frequencies of a few hundred GHz, with the observed flux density sitting below the 10 mJy level [109].The sensitivity of the two highest frequency bins is set to 100 mJy, which is roughly the 4σ line sensitivity for 900 seconds of observation estimated in [110].In general, one would either marginalize over the unknown parameters, or take the value which give the most conservative constraints [48,50] -since our goal, however, is only to provide an indicative idea of rough sensitivity, we simply fix θ m = 0, d = 8.3 kpc, M NS = 1.4 M ⊙ , and r NS = 12 km.Note that since we do not vary θ m , which plays an important role in determining the line width, we make the simplifying assumption 7 GBT can also observe at 100 GHz -using the online sensitivity calculator for GBT, we find a couple hours of observation tends to give comparable sensitivity estimates. 8The line studied here is expected to be slightly wider than this level, however the sensitivity scales weakly with bandwidth, and can easily be compensated for using additional observing time.
that the entirety of the signal is contained within a single frequency bin -this bin is assumed to have a value of 10 −5 × m a for all telescopes except GBT, where we take the bin width used in the observations of Ref. [48].
The projected sensitivity to the axion-photon coupling in each of our four magnetar models M1-M4 motivated above, is shown in Fig. 12.In this figure we show the flux predicted at a typical viewing angle, defined by generating 10 3 samples and selecting the median value, and the ±1σ variations about the median value.Fig. 12 illustrates that SGR 1745-2900 may produce observable radio emission up to ∼ 4THz, however only if there exists a sizable non-dipolar contribution to the magnetic field and the charge density exceeds the minimal expected value by an order 10 value -nevertheless, emission up to ∼ 500 GHz is still expected across all models, potentially allowing parts of the QCD axion band to be explored using existing instruments.These sensitivities should still be interpreted with caution, as systematic uncertainties have not been folded in, and the impact of non-dipolar field modeling has not yet been explored.

VII. CONCLUSIONS
In this work, we have constructed a generalized ray tracing framework capable of analyzing radio signals sourced from resonant axion-photon mixing in astrophysical plasmas, focusing in particular on the treatment of these interactions in highly magnetized plasmas and curved spacetime.We have explicitly shown how these calculations can be self-consistently embedded using either a 'forward ray tracing' approach (in which one samples from the photon phase-space at production, propagates the photons to far distances, and reconstructs observables from the final photon distribution) or a 'backward ray tracing' approach (the more conventional ray tracing approach, in which one propagates rays from an observing plane far away from the source to the point of production); while these approaches use different methodology, we have shown using detailed phase-space arguments (Sec.II-III) that these must yield identical results.We then demonstrated this spectacular agreement explicitly through extensive numerical comparison of the two codes.Note that this is a highly non-trivial result, as even small deviations in the definitions of fundamental constants or accumulated errors in the ODE solvers can generate sizable effects.
Previous work on ray tracing in astrophysical axion searches have included only a subset of the effects studied here, focusing either on the propagation of photons through a magnetized plasma in flat space [22] (using forward propagation) or through an isotropic plasma in curved space [23] (using backward propagation).This work unites these frameworks and allows for a thorough investigation of each of the assumptions adopted in the literature thus far.Our primary conclusions are as follows: • In curved space, the anisotropy of the plasma tends to squeeze the radiated power to small angles (but away from the magnetic pole).This can cause the observed power to deviate from the isotropic scenario by potentially an order of magnitude or so, depending on the axion mass and the viewing angle.
• For large neutron star and axion masses, gravity can induce sizable modifications to photon trajectories, and tends to isotropize the radiated flux.
For small neutron star and axion masses this effect becomes negligible.Interestingly, we find that previous approaches which had included the effect of gravity in the initial conditions but not in the propagation of photons are extremely accurate, despite not being self-consistent.
• Varying the neutron radius within a range of values permitted by the equation of state can lead to a factor of ∼ 2 shift in the sky-integrated power.This effect is predominantly driven by the change in the resonant surface area.
• The improved backward ray tracing algorithm now accounts for multiple photon production sites along photon worldlines (see Sec. V A 3). Including this effect can enhance the total power by a factor of a few (the effect being more pronounced for smaller axion masses).Crucially though, the inclusion of these effects is needed to have agreement between the two ray tracing approaches used in this paper.
In this work (and Ref. [54]), we have gone to great lengths to fully develop kinetic theory in anisotropic media, which has additional complications relative to isotropic backgrounds.In particular, generalizing the dispersion relation for photons in anisotropic plasmas to curved spacetime [57,[60][61][62], is somewhat more involved that an isotropic plasma [53,58].Similarly, the geometry of the conversion surface is also more complicated (both in flat and curved spacetime), and consists of a foliation of multiple production surfaces.In turn, these correspond to more complicated kinematic matching of axions and photons at the conversion surface.
Related to this discussion of anisotropic media, we have paid particular attention to including curved spacetime effects in our routines, which are required to selfconsistently incorporate gravity across the full range of physical effects.For self-consistency, gravity should also be incorporated into the conversion probability itself, such that when integrating over phase-space (see Sec. III) the resulting power is convergent.In this work, and Ref. [54], we have shown that in flat space, both the isotropic and anisotropic conversion probability gives finite results.In Sec.IV, we also offered arguments for generalising the conversion probability in isotropic plasmas to curved spacetime, showing that this generalized form of the conversion probability leads to convergent results.One of the remaining open problems, however, remains a full generalization of the anisotropic conversion probability to curved spacetime.Presumably the answer lies somewhere in generalising the phase-space and kinetic theory arguments of the present work using techniques laid out in Refs.[63,64], though we leave such derivations for future work.
With a view to observations, in Sec.V, we have used our newly developed ray tracing framework to revisit sensitivity estimates of radio and microwave telescopes to spectral lines emanating from Galactic Center magnetar SGR J1745-2900.Here, we provided an extensive discussion on the state-of-the-art knowledge of charge distributions in magnetars, showing that previous approximations relaying on the Goldreich-Julian charge distribution have likely underestimated the characteristic plasma frequency near the surface of the star.Using four distinct models, we show that the high plasma densities near the surface of the magnetar are capable of generating electromagnetic signatures up to the O(THz) regime, with current and future infrastructure potentially covering significant unexplored regions of axion parameter space.This work provides greater motivation for understanding magnetar charge densities, as such searches may provide one of the unique avenues for indirectly probing this wellmotivated region of axion parameter space.
Ray tracing has proven to be an invaluable tool in astronomy and astrophysics, and in recent years has emerged an increasingly important approach in the indirect search for axions.This paper has developed the fundamentals necessary to incorporate ray tracing into astrophysical axion searches, and has for the first time investigated and quantified the validity of a variety of different assumptions adopted in previous applications and searches.While there still exist open questions which need to be addressed, such as the impact of uncertainties in the charge distribution and near-field magnetic field configuration, and how axions and photons mix in high magnetized environments, the framework developed here lays the much needed groundwork for the future development of indirect axion searches in neutron star mag-netospheres.taking M NS → 0 (since the initial conditions for the photon momentum here are fixed in both the curved and flat space analysis), and thus we do not expect the results to mimic those of Fig. 6.
The second approximation, which was adopted in [22,41,49], assumes the resonance occurs along a single two-dimensional surface appearing at m a ≃ ω p .As mentioned in the main text (see Sec. IV B), resonances occur across a foilation of surfaces, defined by k a µ = k γ µ , which can extend to slightly larger and smaller radii than what one would infer by applying the former approximation.In Fig. 14, we show the relative importance of including, or neglecting, the proper kinematic matching condition.Note that in the latter case, in order to ensure the photon in on-shell, we set initial conditions of photons by taking ω γ = ω a and the unit 3-momentum vectors to satisfy ka = kγ .The normalization |k| of the photon 3-momentum is then inferred from Eq. (38).The radial width of the foliation of surfaces defined by the appropriate resonance condition k a µ = k γ µ tends to be below the 10% level, and thus the effect is not expected to be large; this expectation is confirmed in Fig. 14, which shows that the differential power is only slightly modified for a narrow region of viewing angles.

1 FIG. 1 .
FIG.1.Demonstration of Forward and Backward Ray Tracing.Backward ray tracing (left) and forward ray tracing (right) procedures.We show the Langmuir-Ordinary (LO) modes for a Goldreich-Julian plasma density in a strong magnetic field plasma with dispersion relation of Eq. (23).Rays emanate from a plasma isosurface ωp = 10 −5 eV.In the case of backward ray tracing, we show an image plane at θ = 1.2 whilst in forward-tracing gray trajectories show all rays sampled from the isosurface, whilst red rays denote those binned into a viewing angle of θ = 0.7.Other values were B = 10 14 G, ma = 10 −5 eV, α = 0, P = 2π, MNS = 1M⊙.

FIG. 7 .
FIG. 7. Neutron Star Mass.Same as Fig. 5, but showing the impact of varying the neutron star mass from 1 M⊙ to 2.2 M⊙.Result is shown for an axion mass of ma = 10 −5 eV (left) and 10 −6 eV (right).Bottom panel in each case shows the relative difference (in percent) with respect to the fiducial neutron star mass of 1 M⊙.

FIG. 8 .
FIG. 8. Neutron Star Radius.Same as left panel of Fig. 5 but varying the neutron star radius and neutron star mass between rNS ∈ [10, 12] km and MNS ∈ [1, 2.2] M⊙.As before we show in the bottom panel the relative difference (in percentage) with respect to the fiducial models (taken to be those with rNS = 10 km).

FIG. 9 .
FIG.9.Multiple Photon Production.The impact of including multiple resonant photon productions along the photon worldline in backward ray tracing.For ma = 10µeV and the same fiducial NS parameters used in previous plots.

FIG. 11 .
FIG. 11.Magnetar Plasma Distributions.Left: Log-10 of the plasma frequency in the GJ model for an aligned rotator with B0 = 1.6 × 10 14 G and P = 3.76 s.The vacuum regions defining charge separation in the GJ model can clearly be seen extending across the diagonals.Center: Log-10 of the plasma frequency as computed using Eq.(80) for the same parameters.Right: Log-10 of the ratio of the plasma frequency in the preceding models.In all cases the open field lines have been excised (shown with white hatched region) as the charge densities are expected to differ notably from the minimal force-free values in this region.

FIG. 13 .
FIG. 13.Artificial Gravity.Result of adopting initial conditions consistent with photon production from axions near a neutron star with mass MNS = 1 M⊙ (left) and 2.2 M⊙ (right), but setting MNS = 0 in the ray tracing procedure.