Axion-Photon Conversion in Neutron Star Magnetospheres: The Role of the Plasma in the Goldreich-Julian Model

The most promising indirect search for the existence of axion dark matter uses radio telescopes to look for narrow spectral lines generated from the resonant conversion of axions in the magnetospheres of neutron stars. Unfortunately, a large list of theoretical uncertainties has prevented this search strategy from being fully accepted as robust. In this work we attempt to address major outstanding questions related to the role and impact of the plasma, including: $(i)$ does refraction and reflection of radio photons in the magnetosphere induce strong inhomogeneities in the flux, $(ii)$ can refraction induce premature axion-photon de-phasing, $(iii)$ to what extent do photon-plasma interactions induce a broadening of the spectral line, $(iv)$ does the flux have a strong time dependence, and $(v)$ can radio photons sourced by axions be absorbed by the plasma. We present an end-to-end analysis pipeline based on ray-tracing that exploits a state-of-the-art auto-differentiation algorithm to propagate photons from the conversion surface to asymptotically large distances. Adopting a charge symmetric Goldreich-Julian model for the magnetosphere, we show that for reasonable parameters one should expect a strong anisotropy of the signal, refraction induced axion-photon de-phasing, significant line-broadening, a variable time-dependence of the flux, and, for large enough magnetic fields, anisotropic absorption. Our simulation code is flexible enough to serve as the basis for follow-up studies with a large range of magnetosphere models.


Supplemental Material
Appendix A: Photon Dispersion Equation In this section we outline the derivation of the generalized dispersion relation relevant for the propagation of radio waves near the surface of neutron stars.First, we derive the dispersion relation for a cold and non-relativistic plasma (for a detailed overview, see e.g.[64] and references therein), and then generalize the derivation to the case where electrons and positrons have non-negligible boosts in the plasma rest frame.
We begin by writing down the linear equations of motion for non-relativistic charged particles in the presence of an external electromagnetic field, given by m s dṽ s dt = q s ( Ẽ + ṽs ⇥ B) , (A1) where we have introduced the subscript s to denote the species (e.g.electron, positron, proton, etc).The relevant Maxwell equations can be expressed as 17 Here the current is given by J ⌘ P s q s n s ṽs and n s is the number density.We can now assume an isotropic and homogeneous plasma, and take a plane wave ansatz for the fields as Ẽ = Ẽ1 e i( k•r !t) , B = B0 + B1 e i( k•r !t) , and ṽ = ṽ1 e i( k•r !t) , where B0 is the static background magnetic field.We furthermore take B0 to be in the ẑ direction and assume that the magnitude B 0 B 1 .Using Eqs.(A1-A3), one can solve for the various components of the current density Here, we have introduced the E ± notation to define E ± ⌘ E x ± iE y (and similarly for the velocity and current vectors).We have also defined the cyclotron frequency of species s as ⌦ s ⌘ q s B 0 /m s ; for the plasma of interest here, the dominant species are electrons and positrons, and thus ⌦ s = ⌥⌦ e .The dielectric tensor ✏ is generically defined as where the conductivity can be inferred from the current equation J = • Ẽ.With the current as defined in Eq. (A4), one can express the dielectric tensor as (A7) 17 In this section we will work with Gaussian units and take c = 1 (the latter being adopted throughout the text).
Maxwell's equations given in Eq. (A2) can be expressed in terms of the dielectric tensor as Combining these equations yields the so-called wave equation, given by where ñ = k/! is the index of refraction.Without loss of generality, we can take k to lie in the x-z plane; in this case, the dispersion relations of propagating electromagnetic modes can be directly determined by solving where we have defined ✓ as the angle between k and B0 .In the low frequency and large magnetic field limit, which equates to |⌦ e | !, !p , one finds the following three dispersion relations These correspond to the magnetosonic-t, Alfvén, and Langmuir-O mode, respectively (see e.g.[65]).
The dispersion relations above are derived assuming a non-relativistic plasma.In order to determine whether our calculations are robust against this assumption, we now derive generalized dispersion relations for charged species with arbitrary phase space distributions.This discussion largely follows the work of [65,[91][92][93] and the references therein.As before, we will work in the low frequency and strong magnetic field limit.In addition to taking |⌦ e | !, !p , the assumption of a strong magnetic field allows one to focus on plasma distribution functions that are one-dimensional (oriented along the magnetic field lines).
We begin by generalizing the equations of motion for charged particles in an external electromagnetic field.In full generality, the evolution of the distribution function of a species f s is given by the Vlasov equation [65,92,93] where ũ = ṽ is the proper velocity and = 1/ p 1 v 2 .Note that we still take B = B0 + B1 , with the background field B0 being oriented along the z-axis and much larger than B1 .The charge density and current density can be defined in terms of the distribution functions using In order to solve Maxwell's equations, it is useful to expand the phase space distribution in cylindrical velocity coordinates as By rewriting the Vlasov equation, this expansion allows one to write a single equation governing the evolution of each mode f s,n where the subscript s has been omitted for clarity.The operators L n and G n are given by In the above we have introduced the definitions ↵ = q/m, ⌦ = ⌦ s / , E = E x + i E y (the same relation holds for B and r), and we have defined One can also write the charge and current densities in terms of the moments f n , which yield In this notation, Maxwell's equations can be expressed as Up until now the analysis is still completely general and while the above equations are no easier to solve than the original versions, our notation has clarified that the current density is fully defined by knowing only f 0 and f .In order to identify these terms, we can exploit the low frequency approximation in which a small parameter ⇠ ⇠ !/|⌦ s | ⌧ 1 is introduced through the substitution ⌦ s !⌦ s /⇠.With this substitution, Eq. (A19) can be written as Due to the equilibrium distribution being gyrotropic [65], all f n with |n| > 1 should vanish if Ẽ !0, B ! 0 and r !0. Looking at the second of the above equations, this fact tells us that these f n cannot contain terms dependent on negative powers of ⇠.Furthermore, the overall form of Eq. (A31) shows that all terms in a given f n have to be at least one order of ⇠ higher than the lowest order term in f n 1 .As a result, it is possible to represent the phase space distribution in the following way One thus needs to start by determining f 0 , and all higher order f n can subsequently be found through Eq. (A31).We can further expand in the weak turbulence limit, ⌘ ⇠ E/B 0 ⌧ 1 [65], and keep only terms up to orders ⇠ 2 and ⌘, which is enough for our purposes.This allows us to write f 0 and f in the general forms Additionally, working up to order ⇠ 2 , and plugging Eq. (A32) into Eq.(A31) results in If we now switch to Fourier space and assume space-time dependencies given by exp(i k • r i!t), the above equations can be combined to derive explicit expressions for both f 0 and f .They are given by Here we used that B = k ⇥ Ẽ/! and introduced the parameters ⇣ = !k z v z and µ 0 = [⇣(@/@u ? ) + k z v ?(@/@u z )]/!.Notice that the artificial parameter ⇠ previously introduced to justify the expansion of Eq. (A31) has also been set to 1, justifying the validity of the above procedure.
At this point, one can directly compute the current density, and subsequently the dielectric tensor.Pulling out the density normalization from F s,0 , such that R d 3 uF s,0 = 1, the generalized dielectric tensor is given by FIG. 5. Comparison of the three electromagnetic modes present in the infinite magnetic field limit (note that the magnetosonict mode in this limit is synonymous with the free space mode, ! 2 = k 2 ).The Alfvén and Langmuir-O modes are shown for waves propagating at an angle ✓ = 0.2 (solid) and 0.5 (dashed) radians with respect to the magnetic field.Both modes are shown taking the limit of a non-relativistic plasma (red, purple) and for a Waterbag distribution with cut-o↵ velocity vc = 0.95 (see Eq. ( A42)).The vertical dotted line identifies the branch cut in the non-relativistic Langmuir-O mode arising at ! = !p.
Focusing on leading order terms in the low frequency limit, Eq. (A38) reduces to Some of these corrections contain thermal terms that can partially remove the suppression.Near the neutron star surface where the thermal corrections may be large, however, one expects the cyclotron frequency to be ⌦ e ⇠ 4 ⇥ 10 6 eV, which is about 10 orders of magnitude larger then the maximum frequency we are interested in.It therefore seems safe to neglect the higher order terms in Eq. (A39), even in the presence of large thermal corrections.This means that in the end we will only need to evaluate one term, namely To this end, we start by taking the thermal plasma distribution to be one-dimensional, i.e.F s,0 (u ?, u z ) = Fs,0 (u z ) (u ?)/u ?with normalization R du z Fs,0 = 1.Computing this remaining function yields where we have defined n k = k z /! = n cos ✓.To further simplify this expression one must define the shape of the one-dimensional distribution function.Since we are interested in studying the general behavior or the dispersion FIG. 6. Result of ray-tracing photon trajectories from source location at the conversion surface (black contour) of a neutron star, assuming a GJ model of the magnetosphere.Photons are sourced in the x z plane (with k = 10 3 ma r), allowed to propagate in three dimensions, and projected back into the x z plane.The three panels correspond to dispersion relations p , and Eq. ( 1).
relation, we adopt an intuitively straightforward distribution for both electrons and positrons called the Waterbag distribution [94].This assumes a linear flat distribution in velocities up to a cut-o↵ u c , i.e.
With this assumption, the integral in Eq. (A41) becomes with c and v c being the cut-o↵ gamma factor and velocity respectively.We can now directly solve for the dispersion relations using Eq.(A11).In this case, the magnetosonic-t modes remain unaltered by boosts, but the roots associated with the Alfvén and Langmuir-O are now obtained by solving the following One can nicely see that in the limit v c ! 0 and c ! 1, the non-relativistic dispersion relations for the Alfvén and Langmuir-O modes are recovered.An alternative limit of interest is that of a single relativistic species.In this case, the dispersion relation of the Langmuir-O mode is given by (A45) Notice that in the limit where cos ✓ ⌧ 1, photons follow the non-relativistic dispersion relation with the replacement !p ! !p / p c .In addition, in the limit ! 1, one recovers the free space distribution ! 2 = k 2 .This leads to the important observation that the inclusion of a boosted plasma is likely to reduce the importance of the various e↵ects studied in this work.
We demonstrate the behavior of the various modes that appear in the strong magnetized limit in Fig. 5. Here, we plot the spectral index as a function of the frequency (in units of !p ) for each of the modes at angles of ✓ = 0.2 and 0.5 radians with respect to the magnetic field.For both the Alfvén and Langmuir-O modes, we show the impact of taking a Waterbag distribution with cut-o↵ velocity v c = 0.95; the primary e↵ect is to lower (raise) the branch cut of the Langmuir-O (Alfvén) mode, such that a larger range of frequencies can propagate.The vertical dotted line highlights the branch cut that appears in the Langmuir-O mode at ! = !p .Photons sourced from axion conversion on the Langmuir-O mode trajectories will tend to travel to the right in the phase diagram as they escape the magnetosphere, and will quickly tend toward to the free space dispersion relation.Importantly, Fig. 5 illustrates that modes do not cross, and thus remain well-defined and isolated in the strong field limit.Since the rate is suppressed by e ⌧ , one can see that for these parameters the optical depth is beginning to suppress the flux coming from the lobes and torus, but leaves the flux coming from the throats e↵ectively untouched.
In order to illustrate the importance of adopting the correct dispersion relation -and in tracking the individual photon trajectories -we show in Fig. 6 the evolution of photons in the x-z plane propagating away from a neutron star conversion surface using either the free space dispersion relation , or the non-relativistic highly-magnetized plasma dispersion relation given in Eq. ( 1) (right).All photons are sourced with k = 10 3 m a r, where r is the unit vector directed from the origin toward the point of genesis.One can see that the angular inhomogeneities in the plasma frequency of the GJ model generate enormous anisotropic features in the photon trajectories.The impact of accounting for the magnetized nature of the plasma is more subtle, however; this e↵ect is most apparent near the magnetic poles, where photon trajectories can experience strong refraction.

Appendix B: Cyclotron Resonance
The dominant absorption process of low energy radio photons escaping the magnetosphere occurs via the cyclotron resonance, in which electrons in the ambient plasma are excited to higher Landau orbitals.For strong magnetic fields and relativistic plasmas, this occurs when [69-71] where k || = k cos ✓, v || is the velocity of the plasma in the direction of the magnetic field, and ⌦ e is the electron cyclotron frequency.Since the magnetic field falls o↵ as B / r 3 , one expects all photons to cross the cyclotron resonance -for non-relativistic plasmas this occurs at a fixed distance, while for relativistic distributions this may take place over an extended region.We show here that, in the simplified context of the GJ model, the cyclotron resonance is not necessarily negligible.We begin by noting that for non-relativistic plasmas the second term in Eq. (B1) can be neglected, and can be set to one.The optical depth of photons scattering of ambient electrons and positrons is given by where the cross section (assuming !⌧ m e ) for non-relativistic cyclotron absorption is given by [69] Here, r e ⌘ ↵/m e is the classical electron radius.Assuming the photon frequency is approximately constant, we find If one assumes outward radial trajectories this can be approximated as Here, `c is the location of the resonance along the path.Importantly, in the above we have assumed that the resonance condition is met within the light cone with radius R LC = 1/!NS , otherwise we set the optical depth to zero.We illustrate the optical depth ⌧ as a function of sky position in Fig. 7 for a strong magnetic field model with B s = 5 ⇥ 10 14 G.In this case, photons originating from either the bulge or the torus have an optical depth ⌧ ⇠ 1, and thus experience a suppression of the flux on the level of ⇠ 60%.This e↵ect can be far more extreme for larger magnetic fields and larger rotational speeds.

Appendix C: Revisiting the Conversion Probability
An important result of this paper is the correction of the conversion probability adopted in previous works (see main text).There are two problems which have gone overlooked: (i) the angular contribution to the conversion length, and (ii) the de-phasing of the photon and axion wave functions which results from non-linear photon propagation.In order to be clear on the origin of each e↵ect, we review the derivation of the conversion probability before highlighting where the novel e↵ects enter.The equations detailing the propagation of electromagnetic waves in the presence of axions, at leading order in the axion-photon coupling constant g a , are given by where m a is the axion mass, ✏ is the dielectric tensor, and B0 is the local magnetic field.These expressions are derived from Maxwell's equations coupled to axions, assuming a constant background magnetic field B0 .If we furthermore work in the high magnetization and low frequency limit, i.e. |⌦ e | !, !p , the dielectric tensor takes the form (see Appendix A for derivation) where ✏ zz = 1 ! 2 p /! 2 and the rotation matrix is given by Without loss of generality we have taken the external magnetic field to lie in the first quadrant of the x-z plane, at an angle ✓ from the z-axis.Moreover, we choose the photon and axion to travel along the z-axis.Notice that this is a di↵erent orientation relative to Appendix A, where previously we had chosen the magnetic field to be oriented along the z-axis and taken the photon momentum to be in the x-z plane.Nevertheless, the choice is irrelevant as our results only depend on ✓.After computing the dielectric tensor in Eq. (C2), it becomes clear that the photon polarization perpendicular to the external field (E y ) fully decouples.Furthermore, we find that the solution for E z can be re-expressed in terms of E x .Assuming an oscillatory time dependence with frequency !, one then obtains the following equation for the mixing between the axion and the E x component of the electric field Importantly, the general solution thus contains both a transverse and a longitudinal component, as expected for a Langmuir-O mode propagating at an oblique angle with respect to B0 .The longitudinal component will naturally be damped (evolving into a fully transverse O mode) as the wave propagates away from the neutron star.

↵
< l a t e x i t s h a 1 _ b a s e 6 4 = " a 7 Y q k B b A 1 x h F A P 7 R l P g y 5 w T l b < l a t e x i t s h a 1 _ b a s e 6 4 = " j t e / B B e W 4 f F r s R 6 t F + t 1 0 V q y l j O H 4 A e s t 0 / 4 n I 9 q < / l a t e x i t > a < l a t e x i t s h a 1 _ b a s e 6 4 = " I r e r C f r x X q 3 P h a t J a u Y q a I / s D 5 / A E z a l I 0 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D g < l a t e x i t s h a 1 _ b a s e 6 4 = " 6 m z Illustration of de-phasing introduced from non-linear axion trajectories.In the context of the one-dimensional formalism, the conversion length is computed from the de-phasing induced along straight-line trajectories.Should the axion trajectory deviate from this trajectory by an angle ↵, additional de-phasing will be introduced and the conversion length will be shortened.
After changing variables to the vector potential A x = E x /i!, it is possible to adopt a plane wave ansatz and express both fields as Notice here that the photon momentum can be inferred from the dispersion relation derived in Appendix A, and in the non-relativistic limit it is given by The seemingly arbitrary factor k 1/2 present in the photon amplitude emerges from the WKB approximation (see e.g.[64]).In the spirit of [95], we now expand the functions ⇠(z) and (z) in powers of the magnetic field as where we have accounted for the possibility of a strong positional dependence in the photon momentum.Employing both this expansion and the WKB approximation allows us to solve Eq. ( C4) for (1) , yielding where With this expression, one can write the photon field at a distance z (given a negligible initial field value) as Since the energy flux of a plane wave scalar field scales like f / k 2 a |a| 2 , with k a the momentum of the axion field, the conversion probability can be directly obtained from [42,96] 18 (C11) Notice that because conversion only takes place when k (z) ' k a , the pre-factors approximately cancel.Equation (C10) can be integrated using the stationary phase approximation, where the phase in the integral is expanded about @ /@z = 0. Keeping the leading order term in this Taylor expansion and performing both integrals in Eq. (C10), one arrives at where v c is the axion velocity at the conversion point.Importantly, we periodically find that conversion lengths of individual axions can be quite large -notice that this is actually quite problematic since our formalism is not valid in this regime.In order to avoid spurious features, we thus cut out photons whose conversion lengths exceed 1 km.This threshold is somewhat arbitrary, but any reasonable change in this threshold introduces e↵ectively no change in the physical observables.Now, as mentioned in the main text, the second derivative of the phase (and thus also the conversion probability) derived in [39,41,42] has been truncated at leading order in velocity.The derivative of ✓, however, only appears at next-to-leading order.If the axion speed at the conversion surface were always small this would not be a problem, but for conversions near the radius of the neutron star the speed can be as large as v c ⇠ 0.5.Consequently, the next-to-leading order correction can be rather significant, particularly for particles with non-radial orbits.Working up to second order in the velocity expansion, one finds @k @z ' This is related to the conversion length L c via Notice that if one keeps only the leading order in v c , considers perpendicular propagation (i.e.sin ✓ = 1), and adopts the radial trajectory approximation of @! p /@z, this result reproduces the conversion probabilities derived in Refs.[40][41][42].Here, we improve upon these approximations by directly calculating both directional derivatives using auto-di↵erentiation for each trajectory of interest, and find that deviations from the radial trajectory approximation can be significant.
The second concern that we must address is related to the fact that the three-dimensional mixing equations have not been solved (see e.g.[96] for recent progress in this direction).The derivation of the conversion probability had assumed a plane wave solution proportional to e ikz , which amounts to a one-dimensional simplification of the more general plane wave form e i k•r .If photons deviate strongly from this one-dimensional projection over distances r .L c , the axion and photon will no longer oscillate in phase, and the conversion probability will be markedly reduced.In order to address the importance of this e↵ect we begin by identifying the relationship between the conversion length and the phase di↵erence of the axion and photon in the limit where photons propagate along straight trajectories perpendicular to the external field.In this case, the phase overlap of the axion and photon is given by Eq. ( C12), but now with k = q ! 2 ! 2 p .If one directly computes this phase from photon trajectories generated with the

FIG. 11.
Comparison of the time-integrated di↵erential power as approximated in [39] (black), with a run in which (i) absorption and premature de-phasing is neglected, (ii) the radial derivative approximation @ `k !3ma/(2rcvc) is applied universally, (iii) photons are assumed to travel radially outward, (iv) the size of the surface element at each point is re-scaled to be consistent with the radial approximation, and (v) the factor of |v • n| has been removed.The above assumptions should be roughly consistent with the approximation of [39].Results are shown for the fiducial model used throughout this text (left) and the magnetar PSR J1745-2900 (right), but taking ✓m = 0.01 and ma = 10 6 eV (left) and 10 5 eV (right), and ga = 10 12 GeV 1 .
non-magnetized plasma dispersion relation (i.e.! 2 = k 2 + ! 2 p ), one finds that the axion-photon phase di↵erence for all photons (which have not undergone refraction) at the conversion surface is (L c ) ⇠ ⇡/2.As can be seen from the derivation above, the generalization to non-perpendicular directions of the phase involves the replacement of k by that given in Eq. (C6).Nevertheless, this equation still assumes that photons travel on straight line trajectories (notice that this is apparent in the ansatz in Eq. ( C5)) -additional de-phasing may enter as photons refract.In order to account for this e↵ect we approximate the de-phasing by defining a new phase overlap where ↵ is the angle between the initial axion momentum and the photon momentum after some time (and thus depends on time implicitly), and we have replaced the one-dimensional distance z by the path length `.In our analysis we attempt to account for this additional de-phasing by identifying the distance each photon has traveled when ⇤ = ⇡/2; we associate this path length with the 'corrected' conversion length L 0 c .Since the conversion probability is / L 2 c , we re-weight the MC weight by a factor of (L 0 c /L c ) 2 .This procedure is illustrated schematically in Fig. 8.While a rigorous treatment of the three-dimensional mixing is required in order to robustly assess the impact of this e↵ect, we believe that our prescription o↵ers a reasonable estimate of its magnitude.

Appendix D: Parameter Dependence of Sky Maps
The results presented in the main text are shown for a fiducial set of parameters of the GJ model which are expected to be representative of neutron stars in the Milky Way.The purpose of this section is to address the sensitivity of the results presented to the choices of these parameters.In particular, we focus on the parameters likely to induce the largest e↵ects; these include: the misalignment angle ✓ m , the rotational frequency of the neutron star !NS , and the magnetic field strength at the surface B s .We also discuss implications of changing the axion mass.
We begin in Fig. 9 by considering the impact of the misalignment angle on the isotropy (top) and line dispersion (bottom).Specifically we consider three choices of ✓ m corresponding to 0.05 (left), 0.2 (center), and 0.6 (right) radians.The trend in the both cases is immediately apparent -larger values of ✓ m induce stronger anisotropies, in particular with respect to photons produced in the throat, and larger line widths.Furthermore, large values of ✓ m imply less sensitivity to the viewing angle, since the throats more broadly sweep through large areas of the sky.
Since we have discussed the impact of varying the above parameters on line width in the main text, we focus next instead on their impact on the radiated power.We illustrate in Fig. 10 the time-averaged di↵erential power per viewing angle dP/d✓, as a function of viewing angle for various parameter choices.We plot for comparison the predicted power using the formalism of Ref. [39] assuming ✓ m = 0 (all other parameters are either set to the fiducial values used in the main text or to those of the magnetar).We highlight two important results of these figures.Firstly, the flux is strongly reduced relative to previous predictions.Secondly, the sensitivity to the viewing angle is strongly increased (spanning, in some cases, nearly 5 orders of magnitude).
In order to ensure the code developed for this work is reliable, we also present a consistency check in which we attempt to reproduce the results of Ref.
[39] by systematically removing the e↵ects included here.We do this by (i) removing all absorption, (ii) removing the de-phasing cut on L c , (iii) taking the radial derivative approximation @ `k !3m a /(2r c v c ), (iv) assuming photons sourced at some location ri travel radially outward thereafter, (v) rescaling the size of the surface element at each point by a factor of |n • r| (this comes from the fact that [39] has a self contradictory assumption of radial trajectories that are perpendicular to the normal), and (vi) removing the factor of |v • n|.As shown in Fig. 11, this procedure perfectly reproduces the results of Ref. [39].The panels display results for the fiducial model with ✓ m = 0.01 and m a = 10 6 eV (left) and m a = 10 5 eV (right) (and taking g a = 10 12 GeV 1 ).

Appendix E: Time Dependence of Sky Maps
An important feature of the signal discussed in this work is the expected time-dependence of the flux -strong time variations allow for the signal to be more easily disentangled from spurious line signals.As mentioned in the main text, the expected time-dependence of the signal can be read o↵ the flux maps by drawing horizontal lines and tracing the flux over a period of rotation.A single example of this procedure is illustrated in Fig. 2. In this section we present two additional illustrations to highlight the time dependence of the flux for a variety of viewing angles and two misalignment angles.
In each case we bin the flux at each value of over a narrow range of viewing angles (defined by angles of constant ✓) in order to extract the flux as a function of time.Specifically, we take viewing angles of ✓ = 0.3, 0.6, 1, 1.3, and 1.6 radians (these regions are highlighted in the left panel of Fig. 12), and a width in ✓ of 0.02 radians.We show the projected flux (in arbitrary units) for each slice in the right panels of Fig. 12. Depending on the viewing angle, the

FIG. 7 .
FIG.7.Optical depth ⌧ for neutron star with Bs = 5 ⇥ 10 14 G, !NS = 1 s 1 , ✓m = 0.2, and for an axion mass ma = 10 6 eV.Since the rate is suppressed by e ⌧ , one can see that for these parameters the optical depth is beginning to suppress the flux coming from the lobes and torus, but leaves the flux coming from the throats e↵ectively untouched.

FIG. 9 .
FIG. 9. Fractional rate (top) and line width (bottom) in each pixel at distances r ⇠ RLC .Figures from left to right illustrate ✓m = 0.05, 0.2, and 0.6 radians (note that scales are not equivalent).