Radio Line Properties of Axion Dark Matter Conversion in Neutron Stars

Axions are well-motivated candidates for dark matter. Recently, much interest has focused on the detection of photons produced by the resonant conversion of axion dark matter in neutron star magnetospheres. Various groups have begun to obtain radio data to search for the signal, however, more work is needed to obtain a robust theory prediction for the corresponding radio lines. In this work we derive detailed properties for the signal, obtaining both the line shape and time-dependence. The principal physical effects are from refraction in the plasma as well as from gravitation which together lead to substantial lensing which varies over the pulse period. The time-dependence from the co-rotation of the plasma with the pulsar distorts the frequencies leading to a Doppler broadened signal whose width varies in time. For our predictions, we trace curvilinear rays to the line of sight using the full set of equations from Hamiltonian optics for a dispersive medium in curved spacetime. Thus, for the first time, we describe the detailed shape of the line signal as well as its time dependence, which is more pronounced compared to earlier results. Our prediction of the features of the signal will be essential for this kind of dark matter search.


I. INTRODUCTION
Understanding the nature and origin of dark matter remains one of the greatest challenges in contemporary particle physics and cosmology. One of the most compelling explanations is that dark matter is of a particle nature consisting of some as yet undiscovered cold (non-relativistic) component in the present-day Universe. Given that to date, the simplest WIMP candidates for dark matter have not yet been observed, attention is beginning to focus on other dark matter scenarios. A particularly popular explanation is that dark matter may be composed of light (pseudo) scalar fields corresponding to axions-like particles (ALPs) [1][2][3]. These are attractive both for their simplicity and ubiquity in beyond Standard Model theories [4,5].
The axion was initially introduced to explain the absence of charge-parity (CP) violation 1 in Quantum Chromodynamics (QCD) [7][8][9][10][11] see also [12]. Regardless of its theoretical origin, one can still consider the possibility of dark matter consisting of generic ALPs with masses in the µeV range, interacting with photons via L aγγ = g aγγ aF µνF µν , where a is the axion field, and F µν andF µν are the photon field strength and its dual. * richard.battye@manchester.ac.uk † garbrecht@tum.de ‡ jamie.mcdonald@uclouvain.be § sankarshana.srinivasan@postgrad.manchester.ac.uk 1 Note that recently some authors have questioned whether the theory predicts CP-odd observables in the first place [6].
Historically, these ideas can be dated to a proposal by Pshirkov and Popov [34], while the more general question of mixing of axions with photons in neutron star magnetospheres was also considered in [35]. More recently, the subject has seen a renewed interest [16,[36][37][38][39][40][41].
One particularly useful aspect of radio detection of axion dark matter in NSs is that it has the potential to guide and complement existing haloscope experiments. Notable proposed and ongoing examples are MADMAX [42,43], HAYSTAC [44], ADMX [45], ORGAN [46] as well as many other proposed experiments involving novel condensed matter and metamaterial structures [47][48][49].
At present, a number of studies have begun to take data to constrain the pulsar signal for axion dark matter [50][51][52]. However, more work is required to properly characterise the shape and time-dependence of the radio line signal, including its (Doppler broadened) width, as investigated by us 2 in [16]. The authors of [41] developed a ray tracing procedure for deriving more accurate observational properties of the signal. In the present work, we extend this analysis to self-consistently account for an inhomogeneous and time-dependent magnetosphere as well as including the effects of gravity. This allows us to incorporate the bending of the rays due to varying refractive index and compute the Doppler broadening of the signal from the time-dependence of the magnetosphere.
The remainder of this paper is organised as follows. In sec. II and III we review axion conversion in neutron star magnetospheres and describe how to trace rays through the magnetosphere in a space and time dependent plasma. In sec. IV we present the signal properties resulting from our ray tracing analysis, including the shape of the Doppler broadened line signal and the effects of plasma in influencing the time-dependence of the signal. We finish in sec. V where we offer our conclusions and suggestions for future work.

II. AXION CONVERSION IN NEUTRON STAR MAGNETOSPHERES
We begin by introducing axion electrodynamics equations relevant for axion-photon mixing [16,36,53] a + m 2 a a = g aγγ E · B 0 , where B 0 is an external magnetic field, E is the electric field of the photon and a is the axion and ε is the permittivity. In general the solutions to these equations are complicated by the geometry of the magnetic field and plasma in relation to the propagation direction of the incoming axion, as discussed in the 2D simulations [16] and the de-phasing arguments of [53]. This suggests a more comprehensive treatment of mixing is warrented in future work. However, if one specialises to the "Planar case" case in which plasma gradients are aligned with the incident axion [36], and if one assumes the magnetic field is constant, the integration parameter z, then these of stationary phase to arrive at where B ⊥ = B 0 sinθ is the component of the magnetic field perpendicular to propagation and k γ = ω 2 − ω 2 p is the photon momentum and v em is the axion velocity at emission 4 . This leads to a conversion probability also quoted in ref [41], where ω p (x em ) =k em · ∇ω p is the projected plasma gradient onto the direction of propagation given by the unit vectork em . Note this incorporates the full angularθ dependence of the conversion probability, which is contained implicitly in B ⊥ and the directional derivative, which are computed implicitly in our code. At this point we make a comment about the form of the conversion probability, which results from a stationary phase approximation of the following integral where, Note thatk is the local average momentum associated to the average of the eigevalues of the mixing matrix (see appendix B in reference [16] for a more detailed explanation).
The method of stationary phase is essentially a Gaussian approximation of the exponent, where the width of the Gaussian is set by the second derivative of the phase ∆ p giving the ω p in the denominator of P a→γ . However, when rays are tangent to the critical surface, which is the case at the edge of the image, the projected derivativek em ·∇ω p vanishes, so that 1/ω p becomes singular. It might be tempting to interpret this as the limiting case of strongly adiabatic evolution, as discussed in [16]. However, such points actually represent a breakdown of the in strongly magnetised plasmas has been given in [55]. 4 Formally, in a strongly magnetised medium the dispersion rela- as pointed out in [53] (in contrast the dispersion used by some of the same author's in their earlier work [41] which used the isotropic dispersion relation above). This can introduce additional angular dependence as in [53]. These corrections can be incorporated into our pipeline in future work. However since the main goal of the present work is to understand the propagation of photons subsequent to conversion so as to allow a like-with-like comparison to the straight line rays of ref. [41] (re-produced in our fig. 3), we shall use their expression (5).
simplest treatment of stationary phase. Points where the first and second derivative of the phase vanish are known as degenerate stationary points and are the subject of study in their own right. The perturbative computation of the conversion probability results in the following integral Computing the integral in (6) for degenerate stationary points is a subtle topic intimately connected to the field of catastrophe theory [56]. We remark that if one expands the exponent to higher order, one obtains an exponent of the form ∼ z 2 ∆ γ + z 3 ∆ γ . Integrals with cubic exponents of this form are expressible in terms of Airy functions, which are again ubiquitous in optics applications of catastrophe theory. The stationary phase approximation for non-degenerate (∆ γ = 0) points, is valid so long as the first term is dominant over the second, this can be characterised by the dimensionless quantity (∆ γ ) 2 /(∆ γ ) 4/3 . This is in fact the only dimensionless quantity one can construct from ∆ γ and ∆ γ . In the present work, we take a precautionary approach, and excise any points from our analysis for which have (∆ γ ) 2 /(∆ γ ) 4/3 < 1. These typically occur near the edge of the image where rays become tangent to the critical surface. We leave a more careful treatment of degenerate stationary points for future work.
Returning to our discussion, we can now use eq. (5) to compute the radiant intensity I em at the point of emission from the critical surface: where 1/4π is the fraction of axions emitted in a particular direction for an isotropic axion distribution at the critical surface. This gives the energy flow per solid angle for a given point on the critical surface. Note the dark matter density is enhanced at the critical surface by a factor [36,41,57] where ρ ∞ DM is the dark matter mass-density at infinity, M is the mass of the neutron star, and v 0 is the velocity dispersion of dark matter.
Assume we observe the neutron star from a particular line of sight (θ, ϕ), defined relative to the center of the star, with the north pole aligned with the rotation axis. Consider the set of angles (θ obs , ϕ obs ) covering an observation region of angular size ∆Ω containing the image of the magnetosphere. The total collected flux F obs from this region is given by F obs (θ, ϕ) = ∆Ω dΩ(θ obs , ϕ obs ) I obs (θ obs , ϕ obs ) . (10) The area dA subtended by ∆Ω is where D is the distance between the observer and the star. We can divide this area into pixels characteristic size ∆b, labelled by i, with area dA i = (∆b) 2 and infinitesimal solid angle ∆Ω i so that dA i = (∆b) 2 = D 2 ∆Ω i . The integral (10) can then be approximated by This scenario is sketched in Fig. 1. Thus, to calculate the flux, we backtrace rays onto the critical surface. The specific intensity at the point of emission can then be related to that at detection (see appendix B or [58] and references therein) by using the property along rays where n is the refractive index. Hence we have the following relation where n em and n obs and ω em and ω obs etc. are measured in a coordinate system at rest with respect to the star. Combining this relation together with the definition (8) for the intensity at emission and using the fact that far from the star n obs 1 we arrive at wheref (r, r s ) = ω 3 obs /ω 3 em = (1 − r s /r em ) 3/2 1 is the gravitational red-shift factor. Hence in comparison to the formula in [41], we have an additional red-shift factor of order 1, and a factor of 1/n 2 em which compensates for the shrinking of the image in the image plane due to plasma lensing discussed in the next section. Finally we can relate the radiated power dP/dΩ to the flux F obs via

III. RAY TRACING
Ray tracing is a powerful technique for understanding the emission properties of astrophysical bodies and enables one to track the position, frequency and momentum of photons. The rays then contain all the information required to reconstruct the image of the object, as well as giving the angular power dependence, lensing effects and frequency distortions. Such techniques have been applied to both stars [59] and black holes [58,60,61], where in the latter case, the rays are geodesics of the spacetime metric rather than the plasma, but the principle is the same.
Ray tracing techniques were first applied to axion dark matter conversion in [41]. Here, the authors back-traced straight-line rays from the observer to the critical surface on which photons are produced, matching each ray onto its corresponding conversion amplitude. In this work, we examine how plasma affects ray tracing and the corresponding properties of the radio lines.
We shall see that the plasma has two very important effects which greatly modify signal properties in two important ways. The first is a new time-dependence resulting from the refraction of rays, which causes stronger pulsing of the signal. The plasma acts as a time-dependent lens, causing a variation in the number of rays reaching the observer, which fluctuates over the pulse period. The second effect involves the variation of frequency along rays, which allows us to determine the Doppler broadening of each photon, and by summing over all rays, derive the exact line shape of the signal in frequency space.
Our starting point for understanding the propagation is similar to the discussion found in [58], which is a generalisation of the flat space case discussed in [62] (and applied in refs. [63,64]) which describes the propagation of photons in an inhomogeneous and time-dependent plasma using Hamiltonian optics. We begin with a dispersion relation in a cold, isotropic plasma where g µν is the space time metric. Taking covariant derivatives of this equation, we arrive at where we used ∇ µ k ν = ∇ ν k µ since k µ = ∂ µ Θ is the derivative on the eikonal phase Θ in the relevant WKB approximation. We then define worldlines x µ (λ) associated to these rays satisfying where λ is an arbitrary worldline parameter. Putting this together we arrive at where Γ µ νρ are the Cristoffel symbols associated to the connection. We interpret the spatial components of the plasma derivatives on the right hand side of eq. (20) as an effective force, leading to the refraction of rays. Meanwhile, the temporal derivatives lead to frequency evolution along the worldline, as is apparent from eq. (18).
We choose a simple Schwarzschild metric to model the star's gravitational field where f (r) = 1−r s /r and r s = 2GM is the Schwarzschild radius for a neutron star of mass M . The refractive index n of the medium is defined by whereω is the frequency in a coordinate system at rest with respect to the neutron star and includes the gravitational red-shift.
Here ω = −k 0 is the co-moving frequency related to the temporal component of k µ . In the absence of timedependence in the plasma, ω is conserved along rays.

A. Doppler Broadening
When the plasma is time-dependent, ω evolves along rays according to eq. (18), which can be re-written in terms of coordinate time as The key point for this work, is to note that when the plasma background is time-dependent, as happens for the plasma around a neutron star, the photon frequency evolves according to eq. (24). We can obtain an estimate for the Doppler broadening which gives a frequency shift δω satisfying where ω is frequency before broadening, and x 0 gives the ray worldline. Formally we solve eq. (25) numerically for each ray, allowing us to build up the exact line shape for a given magnetosphere model.

B. Effect of Gravity
When taken in combination with refractive plasma effects, gravity plays an important role in influencing the characteristic size of the lensed image of the star. To understand how this happens, it is instructive to consider a simple spherically symmetric and stationary plasma with ω p = ω p (r). In this case, one has two conserved quantities: the frequency ω obs and angular momentum L which is related to the impact parameter b by L = ωb. One can derive a simple energy conservation equation for the radial coordinates [58] corresponding to motion in an effective potential Back-traced photons which are capable of reaching the critical surface, must have a distance of closest approach r min satisfying r min ≤ r c , where r c is the radius of the critical surface. Since r min is by definition a stationary point along the geodesic at r = r min we must have dr/dλ = 0. The maximum impact parameter b max , corresponds to those rays which just skim the critical surface. For these rays r min = r c . We therefore have the maximum impact parameter for rays which can reach the critical surface where we used the definition that at the critical surface, The key point to note is that b max sets the characteristic size of the image in the image plane, which in the toy example we describe here, is a circle radius b max .
We also know that since ω obs is the asymptotic frequency of photons, it satisfies ω 2 obs m 2 a (1 + v 2 0 ) where v 0 is the asymptotic velocity of the axion, set by the velocity dispersion of dark matter. Putting this together, we see that to leading order in v 0 and r s , we have that the characteristic size of the image is given by Let us now consider the scales at play. We can use a canonical model for the plasma density [65] (see eq. (41)) to estimate the size of r c by equating ω p (r c ) = m a . For the pulsar J0806.4-4123 used in this work and [41], we have B 0 = 2.5 × 10 13 G and P = 11.37 sec so that for a mass m a = 0.5 µeV we obtain a characters tic radius r c 5R. For a neutron star of mass M = M and radius R = 10km we have r s /R 0.3. Note this value is quite high, owing to neutron stars being very compact, and quite close to being black holes. Coming back to eq. (28), we see that for the NS values chosen above, r s /r c 0.06 v 2 0 ∼ 10 −6 so that gravity plays a vital role from a ray tracing perspective, in that it makes the area of the image A image ∼ r 2 image four orders of magnitude larger in comparison to a plasma analysis in flat space! Numerically, speaking, this greatly simplifies the task of locating and resolving the image of the critical surface in the image plane. Locating a larger region, is of course much easier than hunting for a highly lensed "pin prick". There is of course no such issue for the straight-line rays considered in [41] where the characteristic size of the image was just given by the geometric cross-section A vac image ∼ r 2 c . It is interesting to note that the authors of [41] claimed that gravity in the absence of plasma only produces a small percent-level correction to the total power. However, as explained above, when taken in combination with plasma, gravity in fact becomes an important component in making the problem numerically tractable by counter-balancing strong refraction from the magnetosphere. The relative image sizes with and without plasma can be seen by comparing the two panels in Fig. 2.
The interpretation of this tension between gravity and plasma is actually rather straightforward. Plasma is repulsive, and so tends to deflect rays away from the observer, which would otherwise reach the image plane in the vacuum case. Gravity meanwhile is attractive, counteracting the effects of plasma, leading to a larger image in the image plane. These distinctions are of course vitally important when it comes to accurate ray tracing, which is crucial both to derive the line shape and timedependence of the signal.
At this point, it is also interesting to relate the expression (27) to the refractive index appearing in eq. (15) which can be read off from eq. (22) and eq. (23) as where we have used the fact that ω 2 p = m 2 a at the critical surface. using ω 2 obs m 2 a (1 + v 2 0 ) we can write Note therefore that n em 1. It is interesting to relate this result to the effective size of the image as given by eqs. (27) and (28), from which we can now read off Hence although the image size is shrunk relative to the vacuum by a factor n 2 em 1, the rays are more intense, receiving a compensating enhancement 1/n 2 em . This is just the lensing principle: that focused light is brighter, so that one ends up with a smaller, but brighter image relative to the vacuum case. As a result, the total power in all directions is conserved when plasma effects are switched on and off but crucially the angular power distribution and time-dependence are different.
Underlying this is of course a conservation principle: so long as attenuation in the medium is neglected, the total power integrated over all emission directions from the magnetosphere must be the same with and without plasma effects. This is because refraction only bends rays, reassigning power to different outgoing directions relative to the vacuum case.

C. Ray Tracing Procedure
Our overall procedure adds many new features to the analysis performed in [41]. As emphasised elsewhere in the text, we include new refractive effects whereas [41] considered straight line rays in vacuum. We therefore must solve a differential equation for each ray which can be highly refracted in the magnetosphere affecting the time dependence of the signal.
Furthermore, we are able to quantitatively assess the Doppler broadening of rays for the first time. For each ray, we compute the Doppler broadening as follows. For a given frequency ω m a we first compute a ray backtraced from a given pixel. Rays from the same pixel with similar frequencies |ω − m a |/m a v 2 0 ∼ 10 −6 have the same world lines up to negligible corrections v 2 0 ∼ 10 −6 . Hence, for a given pixel, all frequencies of interest can be treated as being transported along the same ray to excellent approximation. All that then remains is to compute the frequency shift along that ray. Yet again, the different frequencies in this range experience the same frequency shift δω up to relative corrections O(v 2 0 ). This can be seen by examining the denominator of eq. (25). Consider a spread of frequencies at the point of emission ω em = m a + δω em where δω em v 2 0 m a . The frequency dependence appears in the denominator of the evolution equation, which can be Taylor expanded as 1/ω 1/m a (1 + δω/m a + · · · ) 1/m a + O(v 2 0 ) so that the frequency shift at the point of detection satisfies meaning that the Doppler broadening for a given ray is essentially achromatic within the relevant frequency range. Thus, the reason different rays experience different Doppler broadening is driven by their taking different paths x 0 (t) through the magnetosphere and ending on different points in the observing plane, rather than their different frequencies. Therefore, for a given ray, we can essentially apply an identical Doppler shift to all frequencies. As a result, there is no need to undertake an intractable back-scanning over a large number of observing frequency bins. Such an approach would only correct what is already seen in Fig. 5 to one part in a million. We also use an adaptive method to locate the image of the critical surface, which is much smaller than the unlensed version in [41] where such a procedure would not be necessary since straight line rays simply produce the geometric cross-section of the critical surface. Our explicit algorithm is as follows. Fig. 1. we begin by choosing a given observing direction (θ, ϕ). This gives the line of sight connecting the observer to the origin at the centre of the star. One then constructs an image plane perpendicular to this line of sight. This plane is then divided into pixels of side length ∆b. A light ray emanates from the centre of each. This light ray is then back-propagated from the image plane using eq. (20). As rays are back-propagated one either records a "hit" or a "miss" dependent on whether the ray reaches or is deflected from the critical surface.

As shown in
2. We then use the following adaptive method to resolve the image. First we perform a coarse-grained search with a larger pixel size to locate the disjoint regions in the image plane. Once found, these are then re-scanned with a smaller pixel size. This saves an inefficient high resolution scan over the entire image plane, most of which contains misses due to the image being lensed to a small, high intensity region, as can be seen in Fig. 2.
3. These rays can then be assigned the corresponding radiant intensity (8) at the point of emission on the critical surface. The total luminosity then follows from integrating over the image plane, given by summing up the intensity carried by each ray and multiplying by (∆b) 2 according to eq. (15).

4.
To compute the Doppler broadening, we then take these rays and evolve the frequency evolution along each ray according to eq. (24). This enables us to see how a line signal with frequency ω ∼ m a is broadened due to each ray i receiving a correction ω → ω + δω. The final frequencies at the point of detection can then be binned, so as to derive the exact shape of the signal in the frequency domain.
Our code is written in Mathematica using its in-built differential equation solvers and event locators to detect the critical surface. We also made use of stiff solver options which are required for especially oscillatory worldlines which are multiply reflected. Even though the solver uses an adaptive step size, one must also be sure to choose a maximum integration time-step along rays which is sufficiently small to achieve good convergence of results. In the future we hope to make our code publicly available.
We computed rays in parallel on 32-core cluster nodes. A full pulse profile corresponds to one of the curves in fig. 3. The time taken to compute one of these depends on the following factors (i) the resolution required, set by the total number of pixels: typically tens of thousands (ii) the number of time-steps sampled over the pulse period: more are required to resolve highly oscillatory pulse profiles (iii) the number of 32-core nodes available, which for our purposes was three (iv) the maximum integration time-step used along each ray. Another positive is that the performance of our procedure can be improved even further by implementing our algorithm in another language, e.g. Julia, C++ etc. where the in-built solvers are significantly faster.

IV. SIGNAL PROPERTIES
In order to make a phenomenological prediction for the signal, we take a simple Goldreich-Julian (GJ) model for the plasma density [65]. This begins with the magnetic field of an inclined rotating dipole [66] B r = B 0 R r 3 (cos χ cos θ + sin χ sin θ cos ψ) , (cos χ sin θ − sin χ cos θ cos ψ) , where χ is the angle between the rotation axis and magnetic dipole moment and ψ(t) = φ − Ω t. The neutron star has a surface magnetic field B 0 , radius R and rotational frequency Ω = 2π/P where P is the period. The GJ model gives the density of charge carriers where Ω = Ωẑ is the constant NS rotation vector. Neglecting the relativistic terms in the denominator, we ar-rive at [cos χ + 3 cos χ cos(2θ) + 3 sin χ cos ψ sin 2θ] .
The plasma frequency is then given by where α EM = e 2 /4π is the fine structure constant and m e is the electron mass 5 . We emphasise that our algorithm can be easily applied to any magnetosphere model in future work. We now go on to discuss two observationally relevant properties of the signal: Doppler broadening (which gives the frequency dependence of the signal) and the time-dependence of the signal amplitude as characterised by its pulse profile. Both these are relevant from an observational standpoint in that they play a role in determining sensitivity to the axion-photon coupling. 5 Note that here we assume two charge-separated regions consisting of positrons and electrons, as in [36]. However, as pointed out in [39], if the positively charged region consists of ions the plasma frequency experiences the replacement me → mp which can alter the location of the critical surface when equating ma = ωp.

A. Time dependence
We can characterise the time-dependence of the signal via the relative variance, defined by where · · · denotes pulse averaging.
The timedependence of the signal for various parameter choices is displayed in table I. Typically the time-dependence is stronger for lower axion masses, where it tends to also be larger than the variance reported in the vacuum case. This is also reflected in the comparison between the left and right columns in Fig. 3 which displays the pulse profile of the radiated power. We take as benchmark val-  I. Signal time-dependence. The relative timevariance of the signal σ = (dP/dΩ) 2 / dP/dΩ 2 − 1 where · · · denotes pulse averaging. The quantities σvac. and σ plas respectively give the variance with and without plasma refraction. We displace results for different observing angles θ and axion masses ma. The magnetic alignment angle was χ = 18 • with other model parameters as in previous plots. ues the three axion masses and observing angles used in [41]. Therefore, we see that plasma enhances the timedependence of the signal.

B. Doppler Broadening
One property of the signal which has remained less clear until now, is the precise width of the signal. We previously made an order of magnitude estimate for the Doppler broadening in [16]. Our treatment in this paper allows us to make the first systematic calculation of signal width and fully characterise the exact line shape of the signal. Understanding the signal width is important since the observation time required to achieve a given signal to noise ratio is given by the radiometer equation where T sys is the system temperature, S σ is the flux density noise level, A eff is the effective area of the telescope and ∆ω is the signal width. We therefore see that determining the signal width is crucial in order to obtain the observation time needed to reach a given sensitivity in g aγγ . In Fig. 4 we can see that photons emitted from different sections of the magnetosphere become blue or red-shifted according to eq. (24). We plot the line shape of the signal for the maximum and minimum masses of interest in Fig. 5. We also display the instantaneous shape of the line signal in fig. 7, from which we see that the Doppler broadening can be so sever that the unbroadened signal splits into several lines.These plots can be understood as follows. In both cases, the Doppler broadening is largest for those emission points which lie furthest from the centre of rotation, here the critical surface is moving with the greatest velocity and therefore imparts the most red/blue-shift (Fig. 4). In addition, at lower masses, the critical surface lies farther from the centre of rotation, and so Doppler broadening decreases with increasing mass.
To understand the different line shapes in Fig. 5 we note that for the lower mass m a = 0.5 µeV, the emission points are reasonably democratically spread over the critical surface (Figs. 4 and 6) resulting in a more top-hat like profile, with a large amount of power coming from weakly Doppler broadened points near the equator, resulting in a central spike. Meanwhile, for m a = 12.5µeV the critical surface is sunk mostly within the star, with only small polar cap regions protruding from the surface of the star. The photons here consist of a single blue-shifted and red-shifted patch of equal size, resulting in two distinct bumps in the line in the right panel of Fig. 5.
It is interesting to make an analytic estimate of the Doppler broadening. Using the GJ density in eq. (35), we can compute the frequency shift from eq. (25) to leading order in Ω which gives where ω 0 m a is the unperturbed frequency. The integral in general must be computed for each ray, and depends on the worldline (r 0 (t), θ 0 (t), ϕ 0 (t)). This is of course achieved in full by our numerical solutions. Nonetheless in order to make an order of magnitude estimate, we approximate a trajectory with θ 0 (t) θ c = constant and ϕ 0 (t) = ϕ c = constant. We also approximate dr/dt via the group velocity by taking ω 2 p (r)  Fig. 4. The left and right panels correspond to ma = 0.5µeV and ma = 12.5µeV, respectively. For the higher mass, the critical surface only extends in the small polar cap regions above the stellar surface, which we illustrated with a black sphere. Note the two images are not to the same scale with respect to one another. m 2 a (r c /r) 3 so that where r c is the characteristic radius of the critical surface which we define as This leads to .
(42) Using m a ω 0 and performing the integral, we obtain the characteristic size of the Doppler broadening This reproduces the characteristic size of Doppler broadening based on our preliminary estimate in [16]. From (41) we see that Doppler broadening grows with decreasing mass, where the radius of the critical surface is largest. From the upper axis in Fig. 5, we also see that eq. (43) gives a good order of magnitude estimate for the full numerical result. Note the frequency distribution of the un-broadened signal is inherited from the Maxwellian dark matter velocity distribution f DM (v) ∝ e −v 2 /v 2 0 , so that with frequencies ω = m 2 a (1 + v 2 ) the power spectrum of the unbroadened signal is Gaussian: where ∆ω DM = 1 2 v 2 0 m a and P tot. is the total power, integrated across all frequencies. Each infinitesimal frequency band centered on ω and contained within the original signal, there is a Doppler broadening where D(ω − ω ) is also normalised to 1. Note this is a one-to-many map, which takes an infinitesimal part of the original signal P 0 (ω) and maps it to a broader set of frequency bins. By linearity, we see that the original signal profile is mapped to a Doppler broadened signal In other words, the original signal becomes stretched by the Doppler broadening factor D(ω − ω ), with the final result given by convolving the original line shape with the Doppler shift function D(ω −ω ). The function D(ω −ω ) is of course computed by tracing the frequency evolution of rays, and binning the power at each frequency, with the final line shapes shown in Fig. 5.

V. DISCUSSION
In this work we have presented a framework in which to compute the detailed properties of radio lines resulting from the conversion of dark matter axions in the magnetospheres of pulsars. This reveals both an enhanced time-dependence relative to straight-line vacuum trajectories considered previously in [41] and allows us to rigorously compute both the shape and width of the signal in frequency space.
We also provided arguments based on elementary optics to underpin our results. We related the effective size of the image of the critical surface to the refractive index, which is a function both of the plasma density and gravity at the critical surface. This enables us to interpret the image size in terms of the competing lensing effects of gravity (attractive) and plasma (repulsive).
The Doppler broadening of the signal follows straightforwardly from the geodesic equations for rays propagating through a medium with a time-dependent refractive index. We also derived an order of magnitude estimate for the width of the signal, which supports the estimate obtained in our previous work [16].
Clearly what remains for future work is to provide a complete parameter scanning across axion masses, observing angles and magnetic alignment so as to produce the most conservative sensitivity curves for g aγγ . This could then be applied to astronomical data which would lead to the most robust constraints to date for the conversion of axion dark matter in neutron star magnetospheres.

Note Added
At the same time as our preprint for this paper was released, [53] also appeared which deals with similar questions addressed in this paper. Overall, both works reach the same broad conclusion that plasma-ray tracing is a necessary ingredient in accurately describing the signal, so that numerical ray-tracing should replace the naive modelling in [36] in all future work. While our approach chooses a given observing direction and back-traces rays onto the conversion surface (in the spirit of [41]), [53] uses a semi-stochastic approach, sampling over all-possible photon emissions directions at each point on the critical surface and forward-propagating rays to the sphere at infinity. For a large number of rays, this will then approach the true angular distribution of the signal. At this level, the difference is simply methodological and should not have and impact on observable results.
The first difference between the two works is that [53] includes anisostropic effects in the dispersion relation due the magnetic field. Another notable physical difference between the two papers is that our work incorporates gravity into the ray-tracing equations via minimal coupling. We have subsequently added a supplementary appendix A to this paper to illustrate the effects of modifying the strength of gravity. One can clearly see that the time-dependence of the signal is more pronounced as the effect of gravity is reduced, which seems to approach closer the results of [53] (no gravity) which reports greater time-variation of the signal than we do. This makes sense in light of the arguments of sec. III B, where we describe how, owing to a larger impact parameter induced by gravity a given observer backilluminates/samples a greater proportion of the emission surface. As a result, there is a more democratic spread of the power across the sky with increasing r s , resulting in smoother angular/time variation of the signal.
The authors of [53] also analyse many new effects at the level of the conversion probability itself which we do not include. Clearly in future work it would be interesting to combine all the novel effects considered in both references so as to incorporate the full array of corrections to previous work. One effect we note in passing, is the suggestion of de-phasing effects due to refraction causing the photon to go out of phase with the axion during the conversion process itself. Some preliminary estimates were made in [53] for the size of this effect. Clearly this is quite interesting, and the WKB estimates should be compared against full 3D simulations of the mixing equations (as hinted in [16]), perhaps along the lines of the comparison between [63,64] which compared WKB/raytracing against full electrodynamics simulations in a simple higher-dimensional geometry.
It would also be extremely interesting to carry out detailed benchmarking of both codes against exactly solvable analytic examples, though such an analysis represents a significant undertaking and clearly lies beyond the scope of the present paper.
Finally, we remark that we have carried out convergence tests on our results with the ray integration stepsize, number of rays and distance of the observing plane from the conversion surface. During development, we also tested our code against simple spherically symmetric plasma distributions and their analytic trajectories presented in [58]. We also tested our code in a toy-case for straight-line rays with a rotating ellipsoid with uniform surface luminosity, for which the pulse profile is simply given by the cross-sectional area of the ellipsoid.
In this appendix we briefly display the effect of modifying the strength of gravity in the ray-tracing equations by varying the dimensionless ratio r s /R in the ray tracing calculation. Note we keep the strength of gravity the same in determining the axion density via eq. (9), hence isolating the effect coming from the ray tracing.
In figs. 8 we present the time dependent profile of frequency integrated power for two values r s /R = 0.29, which is the standard value for the cases we have studied, and an artificially lower value of r s /R = 0.03, around a factor 10 lower. We see that there is a much larger time variation in the case where the effect of gravity is lowered. In particular, the peaks in the profile are much larger. Qualitatively, the effect that we see is similar to that seen in [53] who report a much larger variation as a function of t/P .
To investigate this further we present two plots in 9. The first is for an oblique rotator with χ = 18 • where we plot the fractional change in the power between the times when the power is at its maximum and minimum as in 8. The second is for an aligned rotator with χ = 0 • , but this time varying the observation angle and plotting the fractional difference between in the power at its maximum and minimum. Both show a characteristic increase as r s /R → 0 indicating that gravity is important. We note that the increase in the case of the aligned rotator can be estimated analytically.
In particular, we can understand this behaviour qualitatively in a simple analytic example where we consider the special case of a spherical plasma distribution discussed in sec. III B. This has the added advantage that it circumvents any complicated ray-tracing codes and is based on simple physical reasoning. In this case the critical surface is simply a sphere radius r c . The region illuminated in the image plane is a circle whose radius is described by some maximum impact parameter b max of rays which strike the critical surface. These rays emanate from within a circle drawn on the critical surface (shown by the dashed circle in fig 10). For simplicity, we take an isotropic surface emission. The luminosity is then given simply by integrating the surface emission contained within the dashed circle, which subtends an angle θ ap displayed in the figure. Hence, the larger the value of r c , the larger the impact parameter b max , and hence the greater the size of the circle on the critical surface, and the larger the value of θ ap . Therefore, varying θ ap is a proxy for varying r s . Of course the precise algebraic relation between θ ap an r s depends on the form of the spherical plasma distribution on r, and can be determined by solving the geodesic equations for θ evolution in spherical polar coordinates as in [58], but for now we remain general and choose θ ap as a proxy scaling parameter for r s . This is sufficient to build up physical intuition for the results presented in this appendix.
In the analytic example shown in the fig. 10, we consider a critical surface which has two luminous strips angular size ∼ ∆θ, with the surface luminosity being zero elsewhere. This is a proxy for the case of an aligned rotator in the GJ model where the strongest axion emission comes from similar strips above and below the equator near the line where n GJ 0 meets the stellar surface.
The key conclusion is that this toy example produces precisely the same qualitative behaviour as in fig. 9. It helps us to understand that this scaling behaviour with r s is a consequence of the fact that for a given observing direction, a larger proportion of the critical surface is probed, resulting in a more democratic spread of power across the sky. Overall, these examples suggest that gravity needs to be included in these calculations since it counteracts the effects of the plasma.

Appendix B: Covariant radiative transport in plasmas
We now derive the relation (14) used in the main text. By definition, the phase space distribution for photons satisfies where dN is the number of photons in the phase space volume d 3 xd 3 k. Formally, one must define phase space using a local coordinate system, since the background is inhomogeneous [67]. ωdN . Using the expression (B1) for dN , we arrive at We now re-write the momentum phase space element in terms of the refractive index n = k/ω, k = |k| .
Firstly we change to spherical coordinates where dΩ(k) = sin θ k dθ k dϕ k is the solid angle element in momentum space, defined with respect to some axis. Next we write dk in terms of the refractive index by using dk = d(nω) = dω n + ω dn dω .
The final term can be re-written in terms of the group velocity v g = n + ω dn dω Putting this together and combining eqs. (B2)-(B6) we can re-write the momentum volume element as This leads immediately to We again took a value ma = 7µeV. Here t max/min and θ max/min denote respectively the time and angle for which the power is maximum/minimum. Finally, we invoke Liouville's theorem, which states that by definition, the phase space density F is conserved along geodesics of Hamilton's equations. Hence we have where λ is the wordline parameter and x(λ), k(λ) satisfy eqs. (19) and (20).
To understand the power flow through a surface perpendicular to the wordlines, we introduce coordinates where x 1 ⊥ , x 2 ⊥ lie in the plane perpendicular to the direction of propagation, corresponding to x . The surface area element dA = dx 1 ⊥ dx 2 ⊥ . If the group velocity along the worldline is v g , then we have dx || = v g dt. In this coordinate system, the spatial volume element then becomes Hence in these coordinates, eq. (B8) reads This allows one to connect the distribution function F to the Poynting flux which is the energy per unit time per unit area: The specific intensity is then defined as where dΩ gives the solid angle in the direction of propagation. On comparing eq. (B2) with (B14), one obtains that and hence form eq. (B9) conservation of F along worldlines implies that along light rays. Hence we have the following relation between quantities at the point of emission on the critical surface, and points in the image plane: where ω obs and ω em etc. are measured in fixed (i.e. accelerated rather than freely falling) frames.