Axion-photon conversion in 3D media and astrophysical plasmas

With axions now a primary candidate for dark matter, understanding their indirect astrophysical signatures is of paramount importance. Key to this is the production of photons from axions in magnetised astrophysical plasmas. While simple formulae for axion-photon mixing in 1D have been sketched several decades ago, there has recently been renewed interest in robust calculations for this process in arbitrary 3D plasmas. These calculations are vital for understanding, amongst other things, the radio production from axion dark matter conversion in neutron stars, which may lead to indirect axion dark matter detection with current telescopes or future searches, e.g., by the SKA. In this paper, we derive the relevant transport equations in magnetised plasmas. These equations describe both the production and propagation of photons in an arbitrary 3D medium due to the resonant conversion of axions into photons. They also fully incorporate the refraction of photons, and we find no evidence for a conjectured phenomenon of dephasing. Our result is free of divergences which plagued previous calculations, and our kinetic theory description provides a direct link between ray tracing and the production mechanism. These results mark an important step toward solving one of the major open questions concerning indirect searches of axions in recent years, namely how to compute the photon production rate from axions in arbitrary 3D plasmas.


Introduction
Astrophysical environments offer one of the most exciting arenas in which to probe physics beyond the Standard Model, providing extreme conditions in which to search for new laws of nature [1].Amongst these, observations of compact objects, namely stars and black holes, are opening new windows onto the Universe.Neutron stars in particular have become a source of increasing focus for the dark matter community, since they could provide indirect evidence for the existence of dark matter axions.The question of what constitutes dark matter remains one of the most challenging open questions in basic science.Axions [2][3][4][5][6][7][8][9][10][11] are now becoming one of the most popular explanations.Indeed, axion dark matter has been the target of past, present and future laboratory searches .While we await the development of the next generation of axion dark matter detectors, neutron stars continue to offer a powerful indirect probe in the axion mass range m ϕ ≃ µeV .Dark matter axions can convert into radio photons in the strong magnetic fields of neutron stars [43][44][45][46].The conversion occurs in those regions where the axion and photon dispersion relations become degenerate.This occurs when ω p ≃ m ϕ , where ω p is the plasma mass.This mechanism has been the subject of extensive study [47][48][49][50][51][52], spanning a range of radio observations.Efforts have also been made on the theoretical side to model the signal [53][54][55][56][57].Even so, one of the outstanding theoretical questions has been how to calculate the production process itself.The particular challenge is to understand the conversion from axions to photons in the limit in which the wavelengths of the two fields are much smaller than characteristic variational scales of the background plasma, so that they can be treated locally as particle states in the spirit of the WKB approximation.Whilst the classic reference [58] has become the go-to work for describing axion-photon mixing in this limit, it deals only with mixing between axions and photons propagating in a single direction, making it suitable only for describing 1D problems.Reference [44] employed a heuristic description of resonant axion-photon mixing in plasmas, but this too was essentially a 1D description.It was soon noted [47] that 3D effects would be important in describing the production.
There are a number of challenges when undertaking the calculation in 3D.First, the large hierarchy of scales in the WKB limit makes solving the problem numerically difficult.Second, even if the problem could be solved numerically, in order to implement ray-tracing strategies [49,54], the conversion probability must be computed millions of times for each parameter choice of neutron star environment and axion mass.It is therefore crucial to derive an analytic expression for the production of photons from axions in arbitrary 3D media.
To derive an analytic expression for the photon production rate, one must address a number of issues.First, due to the vector nature of photons, multiple polarisations of the electric field couple to one another in Maxwell's equations in anisotropic media.Furthermore, one has to recover a particle and phase-space picture, starting from a set of classical wave equations, namely those of axion electrodynamics.This requires one to somehow identify worldlines of photons and integrate along these to obtain the intensity of photons due to axion conversion.Next, one must have a systematic calculation scheme for exploiting the hierarchy of scales between axion and photon wavelengths, and much larger background scales.This entails performing an expansion in gradients along the lines of a WKB approximation.An additional challenge is that axions and photons in general have different dispersion relations away from the resonance, which complicates finding classical solutions in terms of plane waves.This issue was raised already in Ref. [47].
One also has to deal with the fact that plasma modes are not vacuum photons, but rather collective plasma excitations of charge carriers and the electromagnetic field.This must be incorporated into any photon production rate.In addition, owing to the 3D nature of astrophysical plasmas, one has to incorporate the fact that photons will be refracted relative to axions, and therefore the two particles move on different paths as they enter and exit the resonant point.This fact was noted in Ref. [54], which queried the extent to which photon refraction may affect the conversion probability.Finally, one should ensure that any divergences appearing in the conversion probability are regulated by appropriate physical quantities appearing in the phase-space integration over axion and photon momenta, without the need to impose arbitrary regulatory cutoffs.
Reference [56] approached the calculation of axion-photon conversion in 3D magnetised plasmas by directly solving the equations of axion electrodynamics, seeking plane-wave solutions in the WKB limit to derive a first-order transport-like equation for the electric field amplitude sourced by axions.This equation was then solved using a stationary phase approximation at the resonance to derive the outgoing electric field of the photon, which was subsequently used to define a rate for axion-photon conversion.
In this work, rather than attempting to solve classical field equations directly, we show instead that kinetic theory provides a simple, elegant and powerful framework, which offers a consistent language to address all the issues mentioned above.In this formalism, the calculation can be carried out in a controlled way by deriving a Boltzmann-like transport equation for the phase-space distribution of photons and axions, denoted f γ and f ϕ , respectively.The pertinent Boltzmann-like equation is where H is the Hamiltonian describing the propagation of photons in the medium in question, and E ϕ and E γ are the axion and photon energies, respectively.Here, the delta function in the collision term on the right-hand side imposes the kinematic condition for the conversion to take place.The remaining terms in the collisional piece are the magnetic field B ext , the unit photon polarisation 3-vector ε and the axion-photon coupling g aγγ .Equation (1.1) describes the production of photons due to axions and their subsequent propagation along photon worldlines.By solving the associated characteristics of the Liouville-Vlasov operator that appears on the left-hand side, which are nothing more than Hamilton's equations one can compute the solution for the phase-space distribution of photons f γ along these integral curves, thereby obtaining the ratio at the point of conversion.In general, this expression depends on the kinematics of the photon and axion at that point of conversion.One medium of particular relevance is the strongly magnetised plasma, where the cyclotron frequency greatly exceeds the photon and plasma frequencies.For this case, we obtain an explicit expression . (1.4) The formula (1.4) is the central result of our paper, describing the production of photons from axion dark matter in the plasma surrounding neutron stars.We stress, however, that the description above applies equally to mixing between arbitrary particles in any medium.
A large part of this paper is therefore dedicated to justifying, on the basis of a first-principles calculation that the form of Eq. (1.1) is correct.Our scheme also allows one to make direct contact with ray-tracing routines [54,55,59,60].Not only does this make clear the physical interpretation at every step of the calculation, but the final result is well-behaved when integrated over the photon emission surface, and free from divergences in a formal sense, so long as perturbation theory remains valid.To derive a consistent set of transport equations from which Eq. (1.1) emerges, we use a first-principles approach formulated within the Schwinger-Keldysh closed time path formalism [61,62].This framework provides a path-integral representation for the generating functional of statistical expectation values and, in this way, provides a means for us to extract the correlation functions in the plasma.The Schwinger-Dyson equations for these Green's functions can then be derived from the corresponding effective action, in particular, the socalled two-particle irreducible effective action [63] and its generalisation in the closed time path formalism [64].From there, and following approaches similar to those introduced by Kadanoff and Baym [65], and Danielewicz [66], the Schwinger-Dyson equation can be used to derive evolution equations for phase-space distribution functions, and it is these that will allow us to describe the axion-photon mixing in the magnetised plasma.The above approaches to non-equilibrium quantum field theory and transport phenomena have found wide application, particularly in the context of leptogenesis (for reviews, see Refs.[67,68] and references therein) and baryogenesis [69,70] (both describing the generation of the baryon asymmetry of the universe), and in analyses of neutrino quantum kinetics [71,72] (e.g., in core collapse supernovae).
By projecting the resulting transport equations onto photon worldlines, we are able to derive a simple expression for the conversion probability of axions and photons in terms of a differential equation in a single integration parameter along photon trajectories.In this way, we obtain a description of photon production from axions in a magnetised plasma, taking us closer to a solution of one of the major open theoretical questions in indirect radio searches for axions.
The structure of this paper is as follows.In section 2, we use first-principles techniques from non-equilibrium quantum field theory to derive transport equations for photons that describe their production and propagation in arbitrary media, and then use these equations to derive the explicit form of the photon source term arising from axions.In section 3, we then employ a method of characteristics to solve the transport equation along the worldines of photons, which allows us to derive a conversion probability for the production of photons from axions.Finally, we apply these results to a strongly magnetised plasma, thereby obtaining the central result of our paper.This is, in fact, extendable to mixing with photons and other particles in general 3D media, as demonstrated in section 5 for the example of a parity-even scalar.Finally, in section 6, we offer our conclusions.Further technical details are provided in the appendices.

Kinetic Equations
To arrive at kinetic equations that describe the evolution of the photon phase-space densities f c (with c labelling different polarisations), we are primarily interested in the Schwinger-Dyson equations that determine the photon Wightman propagators D > µν (x, y) = ⟨A µ (x)A ν (y)⟩ ρ and D < µν (x, y) = ⟨A ν (y)A µ (x)⟩ ρ , as expressed in the Heisenberg picture.The subscript ρ on the expectation values indicates crucially that these are statistical expectation values of the vector fields, weighted by the density operator ρ, which, in the Heisenberg picture, is timeindependent and encodes the initial state of the ensemble.The relevant Schwinger-Dyson equations can be derived from first principles within the so-called Schwinger-Kelydsh closed time path formalism by means of the two-particle irreducible effective action.The technical details of these frameworks can be found in Refs.[73,74].However, it is sufficient for our purposes to begin with the Schwinger-Dyson equation for the Wightman function [74,75] where D A νσ (z, y) is the advanced photon propagator and whose partial derivatives are understood to act on the first spacetime coordinate x of the Wightman propagators D <,> νσ (x, y).The projection operator P µα is defined by P µα = u µ u α − η µα , where u µ is the 4-velocity of the plasma and η µα is the Minkowski metric.The term proportional to 1/ξ is a covariant gauge-dependent term.Π µν R (x, z) is the retarded photon selfenergy, and Π <,>µν (x, z) are the Wightman self-energies (which are related to the absorptive parts of the usual Feynman self-energy).These self-energies contain contributions from both axions and the background plasma.In the present setup, the contribution arising from the plasma will give the dominant contribution to the photon dispersion relation, while the piece arising from axions will give a production term for photons.
For a given two-point function or self-energy B = D, Π, it is helpful to define a "Hermitian" function where the subscripts R and A refer to the corresponding retarded and advanced functions, respectively.Noting that we can then write and recast Eq. (2.1) in the form of a Kadanoff-Baym equation wherein we can already identify the collision terms on the right-hand side.
Ultimately, we would like to derive a governing equation for photon distribution functions f γ , which are densities in phase space.It therefore proves convenient to move from configuration space (x, y) to Wigner space (k, X), wherein we Fourier transform with respect to a relative coordinate s µ = x µ − y µ to introduce the four-momentum k and leave behind a dependence on the central coordinate X µ = (x µ + y µ )/2.In doing so, the convolution integral over the spacetime coordinate z in Eq. (2.6) is traded for an infinite series of derivatives in the Wigner-space coordinates (k, X).Further details are provided in Appendix A.
When the wavelength of quasiparticle excitations is short compared to background gradients, such that ∂ X • ∂ k ≪ 1, we can truncate this series of gradients to first order1 to obtain wherein with (2.9) The terms at first-order in gradients involve Poisson brackets of the form The kinetic equation is obtained from the anti-Hermitian part of Eq. (2.7), leading to where and we have defined (2.12) The self-energy Π H is the Hermitian part of the retarded self-energy, which can be obtained according to Eq. (2. 3), and it can be read off from the (loss-free) classical photon polarisation tensor.We note that we have neglected gradient corrections to the collision terms on the righthand side of the kinetic equation (2.11), which do not contribute to the outgoing photon flux at second order in the axion-photon coupling, as described further in Appendix A.
To leading order in gradients, Eq. (2.12) defines the inverse of the retarded Wigner function D R , such that Of particular importance are the eigenvectors ε b and eigenvalues H b of D, where Latin characters a, b, c, . . .index the eigenbasis.We choose the label ϕ for axions to avoid confusion with the label a.The eigenvectors satisfy The ε c give the physical polarisations of the eigenmodes, and setting H c (k, X) = 0 gives the dispersion relation of the mode.The notation H c for the eigenvalues is suggestive, since, as we will see in subsequent sections, they will become the Hamiltonians that describe the propagation of each eigenmode.
By projecting Eq. (2.11) onto the polarisation basis vectors, we can derive a kinetic equation that describes the transport of the plasma mode associated to ε a .The full details are given in Appendix A. Explicitly, we find where H.c. indicates the addition of the Hermitian-conjugate term.Equation (2.15) generalises the results of Sec.2.3.2 of Ref. [74] for scalars to Abelian gauge fields (see also Ref. [78]).
As explained in greater detail in Appendix A, we only consider dispersive effects from the plasma, as embodied by the Hermitian part of the photon response function Π H in Eqs.(2.12) and (2.15).This gives the mass-shell condition for photons and captures refractive effects represented by the Poisson bracket {D, D < }.We do not consider absorptive plasma contributions to photons arising from Π <,> , which correspond to scattering, production or absorption of photons due to real processes that involve scatterings with charge carriers in the plasma.Hence, we only consider axion contributions to the gain and loss terms Π <,> for photons appearing on the right-hand side of Eq. (2.15).
Since we are interested here in the production of photons due to axions, this approximation is justified provided the mean free path of photons is long compared to the length scales over which axions resonantly convert to photons.This is equivalent to neglecting the attenuation coefficients in the standard equations for radiative transfer [79,80].Of course, in the context of neutron stars, attenuation via the cyclotron resonance can be significant [54] for more strongly magnetised stars.However, this occurs in regions far from where axionphoton conversion takes place.We therefore do not consider absorptive plasma effects in this work, since they are not relevant for perturbative axion-photon production in the local environments that we have in mind.
To proceed, we require a solution for D < that involves the photon distribution functions.We begin by noting that, since D is self-adjoint, we can express D R in terms of the eigenvectors and their inverse eigenvalues via Eqs.(2.13) and (2.14) as The pole prescription iησ(k 0 ) involves the sign function σ(k 0 ), which ensures the pole structure corresponding to physical states is consistent with the causal properties of the retarded function.Any solutions D <,> for the Wightman functions must satisfy where ρ is the spectral function.We can then use Eqs.(2.16) and (2.17) to express the spectral function as We are now in a position to seek Ansätze for the Wightman functions D ≷ , which must respect Eqs.(2.17) and (2.18).We consider solutions of the form Evidently, Eq. (2.19) satisfies Eq. (2.17).The f c = f c (k, X) are the distribution functions of quasiparticles (e.g., plasma modes) within the medium.Note that, in order for the f c to be interpreted as distribution functions, each H c ≡ H c (k 0 , k) must correspond to a single dispersion relation, so that, from H c (k 0 , k) = 0, one obtains an equation k 2 0 = E(k) 2 , where E(k) gives the energy of a physical eigenmode in the system.This can be obtained from the Hermitian part of Eq. (2.7), as shown in Appendix A. For our purposes, we will need only the positive roots k 0 > 0 to which the analysis that follows will be restricted.Were any of the H c to correspond to multiple physical states, this would spoil the interpretation of the f c as corresponding to distinct eigenmodes within the system.We do not offer a general proof that this is always true for the decomposition in Eq. (2.16) (though on physical grounds it seems such arguments must surely exist).Instead, we have verified, for instance, for a strongly magnetised plasma, that there are precisely three physical modes, and that each one corresponds to the zeros of a different H c .This is explained in Appendix B, where we analyse the spectral structure of the operator (2.12) in greater detail.Similarly for an isotropic, cold, unmagnetised plasma, we find that there are modes k 0 = |k| and k 0 = |k| 2 + ω 2 p , and each of these corresponds to zeros of distinct H c .
We remark in passing that the summands in Eq. (2.19) are diagonal in polarisation labels.Formally, one could imagine a more general Ansatz that incorporates quantum correlations between different polarisation eigenstates, which would schematically involve a sum of the form Here, the f cc ′ with c ̸ = c ′ would quantify the size of off-diagonal quantum correlations between different polarisations.Similarly, H cc ′ gives an effective mass-shell condition of off-diagonal correlations.This typically gives a mass-shell condition that is an averaging of the mass-shell conditions for the diagonal states c and c ′ , as was shown in Ref. [47,69,81].This form is analogous to the generalised density matrix describing quantum correlations between different neutrino flavours (see, e.g., Refs.[82,83]).
To generate off-diagonal correlations, symmetries in the theory must be broken in such a way that polarisation eigenstates (or flavour, in the case of neutrinos) are not conserved.For photons in a magnetised plasma, since different polarisations have different dispersion relations, mixing is in general kinematically forbidden.The exception would be if a background field can inject momentum transfer between two different modes, or if there is a level-crossing in which the dispersion relations of different polarisations become degenerate, allowing for resonant transitions to occur.The latter is discussed briefly in Appendix B after Eqs.(B.13) and (B.14).Assuming that none of these special cases occur at the point of axion-photon conversion, off-diagonal correlations will not be generated and our diagonal Ansatz (2. 19) is sufficient to characterise the problem at hand.This should fully capture classical physics encoded in axion electrodynamics.We therefore leave a deeper treatment of off-diagonal quantum-correlations for future work.
Returning to the derivation of the transport equations, we treat the quasiparticles as being stable, with real mass-shell relations given by the zeros of H c .Inserting then the Ansatz (2.19) into Eq.(2.15), we arrive at After some relatively straightforward manipulations of the basis vectors and their derivatives (see Appendix C), the left-hand side of this equation can be simplified to give The left-hand term is now recognisable as the standard Vlasov operator [74] as appears in the collisionless Boltzmann equation.Furthermore, the notation for H c now becomes apparent: it is nothing more than the Hamiltonian for the modes.Finally, the terms on the right-hand side of Eq. (2.21) are the collision terms, with the first and second terms corresponding to the gain and loss processes for photons, respectively, as is apparent from the quantum statistics factors f c and (1 + f c ), with the latter corresponding to Bose enhancement of the production process.This is precisely a quantum Boltzmann or transport equation, which governs the evolution of the phase-space distribution f c of each quasiparticle photon mode.
The interpretation of the ε a ≡ ε a (k, X) as polarisation vectors readily follows by considering the equation of motion for the one-point function or more explicitly, from Eqs. (2.1)-(2.2), In the linear regime, wherein we can neglect any terms in Π µν H that depend on the onepoint function, we see that A µ is the zero mode of the same operator D that appears in the Schwinger-Dyson equation for the Green's functions.In the homogeneous limit Π H (x, z) = Π hom H (x − z), the classical solution A µ can therefore be decomposed in momentum space in terms of the homogeneous polarisation vectors ε a (k) that satisfy the Fourier transformed eigenvalue problem Thus, in the weak gradient expansion, i.e., the WKB limit, the response function Π H (k, X) should locally resemble that of an infinite homogeneous medium Π hom H (k) with constant properties, e.g., constant plasma frequency ω p , and magnetic field B. This means that we can identify Π R (k, X) = Π hom R (k) B→B(X),ωp→ωp(X) and promote the homogeneous polarisation vectors ε a (k) to inhomogeneous polarisation vectors ε a (k, X) that satisfy in accordance with the Ansatz for the retarded propagator in Eq. (2.16).
To make the physical properties of our transport equation more manifest, it is informative to integrate Eq. (2.21) over positive frequencies by acting with ∞ 0 dk 0 .This leads to an explicit form where we have dropped the subscript c labelling the photon mode and have written the photon distribution function for a given mode as f γ to distinguish it from the axion distribution function f ϕ to be introduced later.Note that we use the shorthand notation and the chain-rule identities (see, e.g., Ref. [84]) In addition, we used the fact that the denominator |∂ k 0 H| can be written as sign(∂ k 0 H)∂ k 0 H, with the sign function cancelling on the left-and right-hand sides of Eq. (2.27).We identify v g as the group velocity of the mode and ∇ X E as the standard force term.The term ∂ t E gives a generalisation of the force to time-dependent backgrounds.Note that for an anisotropic medium, the group velocity is, in general, not parallel to k (see, for instance, Ref. [80]).

Gauge Invariance and Physical Electromagnetic Fields
The left-hand side of Eq. (2.27) is expressed in terms of physical quantities.We would also like to express the collision term on the right-hand side using gauge-invariant objects involving the electromagnetic fields.We begin by noting that the terms ε * µ ε ν Π ≶ µν and ∂ k 0 H are not separately gauge invariant.The reason for this is that the polarisation vectors ε µ are normalised to unity with η µν ε * µ ε ν = 1, and the covariant normalisation is not conserved under Lorentz transformations.However, the ratio of these two quantities is Lorentz invariant, and can therefore be recast in terms of physical fields.To show this, we start by writing where now E is an arbitrarily normalised vector with E µ ∝ ε µ , whose norm is given by Under gauge transformations, the gauge-field A µ transforms as for some function Λ.In momentum (Wigner) space, the polarisation vectors E therefore transform as It then follows that, under gauge transformations, where the terms involving contractions of k µ with Π ≶ µν all vanish on-shell, since is slightly more involved.We give the full details in Appendix D.
We now show further how to express E 2 ∂ k 0 H in terms of the energy density of the plasma mode.Specifically, it is straightforward to show, starting from the identity E * • D • E = HE 2 , dividing by the photon frequency ω, differentiating with respect to ∂ k 0 and putting everything on-shell (i.e., taking As shown in Appendix D, since the right-hand side of Eq. (2.34) is gauge invariant, it follows that the left-hand side can be evaluated in any gauge to derive an expression for E 2 ∂ k 0 H.We choose temporal gauge in which E 0 = 0 (and drop the covariant gauge term by taking ξ → ∞ in Eq. (2.12)).In this case, the electric field satisfies , where ϵ ij and Π ij H are the spatial parts of the dielectric and polarisation tensors, respectively.Using these relations and the explicit form of Eq. (2.12) allows us to write Eq. (2.34) as We can then use Faraday's law in momentum space.This can be inserted into Eq.(2.35), giving which, using the temporal gauge identity E 2 = − |E| 2 /ω 2 , can be written as (2.37) Herein, we have defined as the ratio between the electric energy density U E given by the numerator, and total energy density U in the denominator [56,85].Relations similar to (2.36), which relate derivatives of the eigenvalues of the wave operator (2.12) to the ratio of stored energies, have been discussed at length in Ref. [86] and Sec.2.3.6 of Ref. [87].Finally, since E * • Π ≶ • E is also gauge invariant, it too can be evaluated in temporal gauge, to give Here, we have again used as the unit 3-polarisation of the electric field for the given eigenmode.Putting all these steps together, we can insert Eqs.(2.37) and (2.39) into Eq.(2.29) to get This can then be inserted into the right-hand side of Eq. (2.27) to yield a form of the Boltzmann equation in terms of electromagnetic fields Note that in the equation above, and henceforth in the main text, we shall simply write lower-case x instead of X for the central coordinate -this is partly for ease of notation and to make its interpretation as a position phase-space coordinate more readily apparent.The appearance of the stored energy in the plasma, as captured by the factor R = U E /U defined in Eq. (2.38), has also been identified in scattering processes in plasmas [86] and was discussed in Ref. [56].

Axion Source Term
Having derived a consistent set of transport equations that determine the evolution of the photon distribution functions, we now turn our attention to isolating the photon production process via axions.To this end, we assume that the axions and photons are not in equilibrium with each other.We also neglect the emission of photons directly from the plasma, and focus on the regime where the mean free path of photons is much larger than the spatial extent of the plasma, such that absorption processes can be neglected.Formally, this corresponds to the decomposition of the photon Wightman self-energy into where Π ≶ pl and Π ≶ ax correspond to the production/absorption of photons due to plasma processes and the axion background, respectively.The approximation amounts to isolating the contribution to photon production from Π ≶ ax and neglecting any effects arising form Π ≶ pl .Working then in the limit that f γ ≪ 1, we can take the Born approximation in the collision terms and only the gain term remains on the right-hand side of the Boltzmann equation.In this regime, putting Eq. (2.21) on-shell gives the following equation for the photon distribution functions: We reiterate that the polarisation eigenlabel c on f c γ , ε c and H c has been suppressed; it is understood that we deal with production of a particular photon eigenmode, whose distribution function we label simply f γ .
We recall that the object H = H(k, x) is the Hamiltonian of the mode, with k µ and x µ giving the 4-momentum and coordinates of the photon, respectively.Setting H(k, x) = 0 gives the dispersion relation for that mode.For example, in an unmagnetised cold plasma, one has H = k µ k µ − ω 2 p .The object Π <µν ax (k, x) is the Wigner transform of the photon Wightman self-energy, and corresponds to production due to axions.It is defined by Here, j µ ϕ is the effective current arising from the axion-photon interaction where F µν ext is the dual external field strength tensor, originating, for instance, from a background magnetic field.The form of the production term on the right-hand side of Eq. (2.44) has the form of the standard photon-production term for gauge theories [88][89][90][91], except that here it applies to more general polarisations relevant for arbitrary inhomogeneous and anisotropic media.It is straightforward to show that, to leading order in g aγγ (see Appendix E), the Wightman self-energy contribution from axions is given by where f ϕ (k, x) is the phase-space density of axions.The result (2.47) is generated from the diagram in Fig. 1.Furthermore, we should remember that the k 0 appearing in Eq. (2.47) is on-shell with respect to the photon dispersion relation, i.e., k 0 = E γ (k, x).Since the axion energy satisfies ϕ , it follows that the delta function can be written as where the dots on the right-hand side denote contraction of Lorentz indices.Note that production is on a mode-by-mode basis, i.e., this equation describes conversion from axions into a given photon eigenmode.We emphasise again that the polarisation 4-vector ε is computed in-medium and is not simply the vacuum polarisation vector.In summary, we see formally that, by solving Eq. (2.49) for f γ , we directly obtain the asymptotic distribution of photons as they escape from magnetised plasma regions and become transversely polarised photons in vacuum.This recasts the problem as one of radiative transfer.

Photon Production and Propagation
The result (2.49) gives the governing equation for the phase-space density of photons.By solving this 8-dimensional phase-space equation, one could in principle completely reconstruct the phase-space density of photons.However, it is also natural to use a ray-tracing prescription in which we track trajectories of individual particles.This is equivalent to solving this equation through a method of characteristics.Indeed, ray-tracing approaches have been used extensively [49,[51][52][53][54][55] to derive astrophysical signals of axions.For this reason, it is more useful to solve these equations for f γ computed along the worldline of photons, whose trajectories can then be used to construct observables, such as the radio flux from an astrophysical environment.
With this in mind, we construct solutions using a method of characteristics, which consists of identifying those integral curves (x µ (λ), k µ (λ)) in phase space that satisfy where λ is a worldline parameter.From the chain rule, it follows that along these curves, Eq. (2.49) reads By solving this equation and propagating the solutions out to infinity, one immediately obtains the asymptotic distribution of photons measured by an external observer.This situation is illustrated in Fig. 2. Note also that E γ,ϕ (λ) is a shorthand for E γ,ϕ (k(λ), x(λ)) and that the axion energy function is evaluated along the photon worldline.Working, as before, in temporal gauge with ε 0 = 0 and integrating Eq. (3.2) gives for λ > λ c , where λ c is the value of the worldline parameter for which the energies satisfy All the quantities on the right-hand side of Eq. (3.3) are understood to be evaluated at the resonance λ = λ c .Primes in the denominator of Eq. (3.3) denote differentiation with respect to λ.Now, for a stationary background, E γ is conserved, so that E ′ γ (λ) = 0.By using where, in the final equality, we used the chain-rule identity (2.28).We emphasise again that k(λ) is the 3-momentum along the photon worldline.Defining a "conversion probability" P aγ via we can read off where v p = k/E γ is the phase velocity of the photon, equal to the phase velocity of the axion at the resonance by energy-momentum conservation.We can use results from the previous section (especially Eq. (2.36)) to re-express this as where ε is the unit polarisation 3-vector of the electric field and B ext is the external magnetic field.It is important to remark that this equation is valid for any quasi-stationary medium, so long as the energy (i.e., dispersion relation) and polarisation vectors are known.Notice in particular that the gradient term in the denominator is projected along the axion worldline, which is in the direction of k, rather than the photon worldline, which is in the direction of v g .In an anisotropic medium, v g and k are in general not parallel [80].
The arguments laid out in this and previous sections have provided a rigorous means both for deriving and solving transport equations for photons to compute the production due to axions.Nonetheless, it is interesting to note that the same equation (3.8) emerges from more heuristic arguments involving naive second quantisation of plasma modes.To see this, one can take the classical modes in the medium and apply second quantisation to these quasiparticles so that the gauge field can be written as [92] A where c c † k and c c k are creation and annihilation operators for the plasma eigenmode c.We emphasise again that these create collective plasma excitations, not vacuum photons, such that ε c are the (local) in-medium polarisation vectors of the inhomogeneous plasma.In this case, one can simply "write down" a Boltzmann equation where the right-hand side is a collision integral containing the matrix element for axion-photon conversion into a particular photon mode.The latter can be read off from Eqs. (2.46) and (3.9), giving Upon performing the 3-momentum integration in the collision integral, we get We can again project this equation onto worldlines, as was done previously, to arrive at Integrating this equation over λ and using Eq.(3.5), we arrive again at Eq. (3.7), justifying our conversion probability in terms of more familiar quantum field theory arguments involving Boltzmann equations and collision integrals.

Total Power
It may look like the denominator in the conversion probability (Eqs.(3.7) and (3.8)) gives uncontrolled divergences when v p is perpendicular to ∇ x E γ .However, when integrated over phase space, this would-be divergence is cancelled by the phase-space measure, as we now explain.
To show this, we first return to the Boltzmann equation (2.49), which can be multiplied by k 0 to arrive at Note that the left-hand side can be written as Next, we can integrate this equation over a finite spatial 3-volume V, whose bounding surface has area element dA, and a region in 4-momentum space d 4 k running over all momenta.
Upon performing these integrations, we arrive at where we have used To derive the second term, we applied the divergence theorem with respect to the spatial integral.We also discarded the surface terms arising from the The source term Q is defined as The integral identity in Eq. (3.16) is, of course, nothing more than a continuity equation for the photon flux.The first term gives the rate of change of the total energy in the volume V, which changes due to flux through the boundary (left-hand side, second term), energy shifts in photons due to a non-stationary background (left-hand side, third term), or due to production (right-hand side).
Next, we can make use of the integral identity where dΣ is the area element on the surface defined by G(x) = 0. Applying this to the spatial integration on the right-hand side of Eq. (3.16), by treating δ E γ (k, x) 2 − E ϕ (k) 2 as a function of x for fixed k, we arrive at where Σ k is the spatial surface on which E γ (k, x) = E ϕ (k).We can leave the integrand in Eq. (3.19) unchanged by multiplying the numerator and denominator by the same factor v p cos θ n , where v p is the magnitude of the phase velocity and θ n is the angle between the phase velocity and the surface normal of Σ k denoted by n.By doing this, we have where dΣ k = ndΣ k is the directed surface element, so that n • vp = cos θ n .Hence, we have Inserting this into Eq.(3.16) gives (for a stationary Hamiltonian) From the right-hand side of Eq. (3.22), we immediately see that the would-be divergence in the probability, arising when v p is perpendicular to ∇ x E γ , is cancelled by the flux-like projection dΣ k • v p in the phase-space measure.Note, however, that the limit ∇ x E γ → 0 does represent a divergence, corresponding to strongly adiabatic conversion.This case is more complicated and will be discussed in the conclusions.

Magnetised Plasmas
We now come to one of the main physical situations in which this conversion can be active: magnetised astrophysical plasmas, as are encountered in the magnetospheres of neutron stars.

Weakly Magnetised Plasmas
In a weakly magnetised, cold plasma, ω c ≪ ω, ω p , where ω = E γ is the photon-frequency (in natural units), ω c = e |B| /m e is the cyclotron frequency, and ω p = e 2 n e /m e is the plasmafrequency with n e being the electron number density.In that limit, the photon dispersion relation becomes In addition, v p • ∇E γ = (ω p /E γ )v p • ∇ω p , and the polarisation vectors are transverse to the direction of propagation, so that , where θ is the angle between k and B ext .The ratio of the total and electric energy in a weakly magnetised plasma turns out to be U/U E = 2.All together, using Eq.(3.8), this yields This reproduces the 1D formula for propagation perpendicular to B (θ = π/2) in the nonrelativistic limit E γ /ω p ≃ 1, as derived in Refs.[47,53].

Strongly Magnetised Plasmas
For a strongly magnetised plasma, we can use the relation between the permittivity and the polarisation tensor We choose coordinates in which k = (0, 0, k) points in the z-direction, and B ext lies in the y-z plane, with B ext = |B ext | (− sin θŷ + cos θẑ), so that θ is the angle between k and B. For a strongly magnetised plasma, neglecting relativistic effects2 , we find the permittivity 3-tensor is given by [93] From this, we can read off, using the equation for the electric field inferred from Eqs. (2.14), the following dispersion relations: ) which follow by setting the H c = 0.These are the magnetosonic-t, Langmuir-O (LO) and Alfvén modes, respectively.To leading order in g aγγ , the LO mode is the only one which can both be produced on resonance and propagate out of the magnetosphere as ω p → 0. The LO mode has electric polarisation unit 3-vector It is also straightforward to show that for the LO mode, U/U E = 2.We verified this by explicit calculation using Eqs.(4.8) and (2.38).This can also be seen by noting that, since the permittivity takes the form ϵ = I−A/ω 2 for some matrix A, it follows that , but the first term vanishes due to the equations of motion, leaving U = 2 |E| 2 = 2U E , as required.
From Eqs. (3.8) and (4.8), we therefore read off This result reproduces the isotropic result (4.2) in the limit θ → π/2 in which the LO mode becomes ordinary and behaves as though in an isotropic plasma.Equation (4.9) is the central result of our paper, providing the conversion probability for axions to photons in a strongly magnetised plasma.Note that we have also verified, via an alternative and more explicit calculation, the correctness of the form of Eq. (4.9).This we did in covariant gauge using the covariant form (3.7) and deriving an explicit form of the eigenvalue derivative ∂ k 0 H.This calculation is given in the Appendix F. We also state here the explicit form of the energy gradient, which corresponds to taking spatial derivatives while keeping k constant.This gives where

Neutron Stars
It is interesting to plot the conversion probabilities in a realistic astrophysical setup.Neutron stars have been the subject of considerable attention as targets for indirect detection of axions [43-45, 47, 48, 50-55, 94, 95].Here, the plasma profile can be approximated by the Goldreich-Julian model [96].This model consists of an oblique rotating magnetic dipole [97], whose components, in polar coordinates, read where α is the angle between the rotation axis and magnetic dipole moment and ψ(t) = φ NS − Ω t, where (θ NS , φ NS ) are polar coordinates defined with respect to the rotational axis.The neutron star has a surface magnetic field B 0 , radius R and rotational frequency Ω = 2π/P , where P is the period.The Goldreich-Julian model gives the density of charge carriers as where Ω = Ωẑ is the constant neutron star rotation vector.Neglecting the relativistic terms in the denominator, we arrive at Next, we assume an electrosphere model, in which the magnetosphere is separated into positive and negative charges with mass m e , so that the plasma frequency is given by ω p = e 2 |n GJ | /m e .For a fixed frequency ω, the conversion surface on which k µ a = k µ γ is then defined by where θ is the angle between B and k, not to be confused with the polar angle θ NS .The final approximation in Eq. (4.15) holds for non-relativistic axions with k ≪ m a .In this limit, all the kinematic surfaces in the foliation Σ k in the right-hand side of Eq. (3.22) become independent of k and collapse to a single surface Σ, given by ω p = m a .It is useful to define a coordinate system for the vector k given by where x, ŷ and ẑ are defined by the Cartesian axis related to the polar coordinates (r, θ NS , φ NS ) and (θ k , ϕ k ) give the direction of propagation of the axion.We can then define an averaged probability over all k directions by where dΩ k is the solid angle appearing in the momentum phase-space measure d 3 k = dkk 2 dΩ k , and θ n is the angle between the phase velocity v p = k/ω and the unit normal n to the conversion surface Σ k .The normal n is parallel to ∇ x E γ .The factor cos θ n in the definition (4.17) has been pulled out from the phase-space measure in Eq. (3.22) to regulate P aγ .The resulting equation (4.17) is therefore a function of the remaining phase-space coordinate x (as well as k = |k|, but this is fixed by kinematics once we specify a frequency ω).The radial coordinate r can also be eliminated by using the resonant condition ω p = m a , which defines a critical surface r = r c (θ NS , φ NS ).Hence, on this surface, ⟨P aγ ⟩ is a function of the polar coordinates θ NS and φ NS .This function is displayed as a function of polar angle θ NS in Fig. 3. Notice that, for certain angles (θ NS , φ NS ), the conversion surface penetrates the star, and no production occurs, as shown by the excised gray regions.In Fig. 4, we also show how divergences in the conversion probability P aγ of Eq. (4.9) are regulated by the flux angle cos θ n , which appears naturally in the phase-space measure, as discussed in Sec.3.1.Note that due to geometry, varying θ k between (0, π) does not necessarily mean cos θ n fully ranges over (0, 1).
The results presented here should be compared with those of Ref. [56], wherein the Authors analyse the equations of axion electrodynamics directly.The resulting expression for the conversion probability differs from our result (4.9), as obtained from the transport equation for the photon distribution functions.Understanding the origins of the discrepancies between these approaches, and any qualitative and quantitative differences in the conversion probabilities, is clearly of great interest, not least from the perspective of analysing how kinetic theory yields results consistent with the classical wave equations.We leave this for future work.We remark, however, that the expression for the conversion probability in Ref. [56] requires the introduction of an arbitrary cutoff (see Fig. 3 of Ref. [56]) to regulate divergences.By contrast, our expression (4.9) does not lead to divergent results, and the conversion probability is rendered finite (see Sec. , v a = 220km/sec and g aγγ = 10 −14 GeV.We display in blue the 1D isotropic conversion probability (4.2) and in green, the new result in Eq. (4.9) of this work.Gray regions show where the production surface would be inside the star, where no photons are produced from axions.

Parity-Even Couplings
Before concluding, we remark on the ease by which this approach can be generalised to couplings of parity-even scalars to the Maxwell term.Such couplings often arise via triangle diagrams in theories with gauge-singlet scalars that couple directly to the Standard Model fermions or mix with the Standard Model Higgs, or in certain scalar-tensor theories of gravity that involve additional gauge-singlet scalar fields that are non-minimally coupled to the Ricci scalar.In the latter case, the (inverse) Primakoff effects have been used directly to constrain such models [98,99].
In these cases, the relevant interaction takes the form In the presence of a background magnetic field B ext , this leads to wherein we have neglected terms that involve gradients of the external magnetic field.This is justified since we have in mind fields ϕ whose wavelength is much shorter than that of the background field, thereby allowing us to neglect ∂ i B ≪ ∂ i ϕ.This leads to a contribution to the photon self-energy given by  4.9) (green), as a function of polar momentum angle θ k of the axion.We fixed ϕ k = 0.2 and neutron star angles θ NS = 0.5 and φ NS = 0.6π.Other values are as in Fig. 3.The blue curve shows how the divergences in the probability are regulated by the phase-space measure, as discussed in Sec.3.1.The gray dashed line indicates the pole in the conversion probability at cos θ n = 0, when the incoming axion momentum is tangent to the conversion surface.
in place of the expression in Eq. (2.47) for the parity-odd scalar coupling.This leads to a conversion probability (5.4) Applications of the present approach for constraining the couplings of parity-even scalars to the Maxwell term are left for future work.

Conclusions
In this work, we have provided a robust derivation for the conversion of axions to photons in a strongly magnetised plasma with arbitrary 3D geometry.In fact, our calculation is valid in any medium for which the dispersion relations and polarisations of photons are known, and where the photon mean free path is large compared to the spatial extent of the production region.This calculation has been wanting for many years, since the mechanism of resonant photon production in neutron stars was first re-popularised some half a decade ago [44,45].Indeed, since then, there has been an active programme of both theoretical and observational work carried out to understand and search for these signals [47, 48, 50-55, 94, 95].It is striking that despite this activity, a reliable calculation of the most fundamental part of this story -the production mechanism of photons themselves -has remained elusive.What has been lacking until now is the identification of an appropriate theoretical language to describe axion-photon conversion in such a way that the complexities of 3D geometries, multiple polarisations, background inhomogeneities and differing dispersion relations of axions and photons (amongst other problems) can be dealt with.In the present work, we have shown how kinetic theory -that is, a theory describing the phase-space of a stochastic ensemble of particles and their collisions -is a powerful formalism to describe axion-photon mixing.
To derive a set of transport equations (also called, Boltzmann, kinetic or radiative transfer equations), we deployed powerful techniques from non-equilibrium quantum field theory.This allowed us to derive the transport equations (see Eqs. (2.44) and (2.49)), which describe the propagation and production of photons directly.Without the tools of non-equilibrium quantum field theory at our disposal, the form of these equations would simply have to be intuited on the basis of generalising other formulae for bare photons in vacuum.Furthermore, deriving the correct form of the source term arising from axion production would have been challenging without these techniques.That being said, rather pleasingly, the final answer, having been derived on the basis of a rigorous calculation, then has an immediate physical interpretation in terms of a classical Boltzmann equation for particles with a classical Hamiltonian H.The source term (collisional piece) resembles the standard source term for gauge-boson production [88][89][90][91].
These equations give a full 3D description of the system and completely encode the kinematics and dispersion relations for axions and photons, which are in general distinct away from the resonance.At first sight, solving the transport equations in 8-dimensional phase space looks unwieldy.However, by recasting these equations in terms of characteristics (i.e., the photon rays arising from a system of Hamilton's equations), we were immediately able to probe the 3D structure of the system, capture the refraction of photons and, at the same time, the conversion from axions to photons within a 3D environment.By evolving the system along photon worldines, we were then able to solve for the outgoing phase-space density of Langmuir-O modes and derive a "conversion probability", which gives the ratio of photon to axion distribution functions.Rather satisfyingly, the conversion probability reproduces the previously derived limit for axion-photon conversion in an isotropic plasma, and the generalisation to 3D has a similar structure when written in terms of polarisation vectors.Our final answer is free from divergences within the limits of perturbation theory, without the need to impose arbitrary cutoffs.
As discussed in the main text, the expression for the conversion probability (4.9) obtained in this work from kinetic theory differs from that obtained in an existing analysis of the classical wave equations of axion electrodynamics, see Ref. [56].Moreover, we have seen that the divergences appearing in our conversion probability density are rendered finite by the phase-space measure in the integrated conversion probability; specifically, by the flux angle between the incoming axion momentum and the conversion surface Σ k .This occurs without the imposition of a separate cutoff, as was necessary in the analysis of Ref. [56].Notwithstanding this difference, it is clear that a detailed and critical comparison of these two approaches and the resulting expressions for the conversion probability is warranted in order to reach a consensus on the correct description of resonant photon production.This will require further analytical and numerical studies beyond the scope of the present work, and it may be presented elsewhere.
We now discuss some other conclusions of this work, which relate to questions raised elsewhere in the literature.
• "Dephasing" does not occur Previous works [54] have considered whether the conversion probability will experience suppression due to photon refraction relative to the axion which sources it, owing to the two particles having different dispersion relations away from the resonance.This speculative effect was termed dephasing [54] and was parametrised by an ad hoc phenomenological factor to suppress the probability according to the amount of refraction.This was done to make conservative estimates for the conversion probability lest the 1D answer overestimate the true conversion probability.
In this work, we were able, for the first time, to account fully for the refraction of photons when computing the conversion probability.Indeed, since our kinetic equations apply to any 3D environment, and the method of characteristics we used to solve them precisely captures the refracted photon worldlines, the final conversion probability must contain (subject to the approximations used) all the relevant information concerning 3D effects.Equation (3.8) should therefore fully account for such effects.How then is photon refraction contained in this formula?
The answer is that it is already contained in the factor 1/|∇ x E γ | appearing in the probability (see Eq. (4.9)).On the one hand, this term can be thought of as capturing the physical width of the resonance, in that it quantifies the size of spatial gradients in the background.In this sense, it defines an effective conversion length over which the resonance persists.However, this term also completely captures the refraction of the photon, and occurs directly from the force term ∂ x H in the Boltzmann equation.Thus, 1/|∇ x E γ | 1/2 can also be thought of as capturing a generalised notion of a radius of curvature R ∝ 1/|∇ x E γ | 1/2 of the photon due to refraction, which quantifies the degree to which the photon is accelerated as it passes through the resonance.This is true whether the photon is accelerated parallel or transverse to the sourcing axion, so that it is more helpful to think of the conversion probability as scaling as P aγ ∝ R 2 with 1/R proportional to the acceleration experienced by the photon.
Hence, the size of background gradients describes the physical width of the resonance but, at the same time, quantifies the amount of refraction (acceleration) experienced by the photon.The two effects are one and the same, simultaneously encoded in the factor 1/|∇ x E γ |.Wider resonances correspond to more weakly refracted photons and vice versa.
Finally, one might question whether some of the approximations that we have made could be the reason why such dephasing is not present in our result.Let us therefore catalogue the main approximations made, of which there are three.They are: (i) the quasiparticle approximation in which decay or absorption of photons is neglected in the left-hand side of the Boltzmann equation; (ii) the truncation of the gradient expansion, i.e., the WKB approximation; and (iii) the Born approximation, in which we work to leading order in g aγγ .Clearly (i) has nothing whatsoever to do with refraction.Similarly (iii) is not relevant, since dephasing was claimed to happen to leading order in g aγγ .Whether or not dephasing occurs at higher order in perturbation remains to be seen.
Thus, the only tentative avenue for question would be our truncation of the gradient expansion (ii), employed in Eq. (A.7).However, truncation at this order reproduces all the required classical equations of kinetic theory.Indeed, this truncation is precisely the WKB approximation, in which particles emerge as a consequence of assuming fluctuations in fields occur at wavelengths much smaller than the characteristic scales of background inhomogeneities.The higher-order terms in the gradient expansion are therefore highly suppressed by 1/(kL) ≪ 1, where L is the variational scale of the background and k is a characteristic momentum of the particles.Hence, any such dephasing cannot meaningfully be hidden in higher-order gradients, since, if such a truncation is forbidden, the WKB approximation is not valid, and it makes no sense to talk about localised particles in phase space in the first place.The whole concept of relative refraction between individual particles then breaks down.We therefore conclude that, at leading order in g aγγ , our calculation is robust, and all relevant effects of refraction are captured in our formula and encoded in the term 1/|∇ x E γ |.
• Resonance width and wave equations versus S-matrix elements and kinetic equations It is interesting to restate the explicit form of the Boltzmann equation quoted in Eq. (2.49) of the main text From the delta function on the right-hand side, it naively appears that the resonance happens instantaneously when E γ = E ϕ .However, if one were to solve classical wave equations for the mixing between the photon field A µ and the axion field ϕ, as was done in, e.g., Refs.[44,47], one would observe that the photon field takes a finite amount of integration time to grow across the resonance.The typical scale over which the photon is resonantly driven is set by some length, L c ∝ 1/(k • ∇E γ ).How then, should we reconcile these two pictures?The answer lies in the fact that ultimately we are only interested in the asymptotic axion and photon states.That is, we require the asymptotic value of f γ along photon worldlines in the kinetic equations, or the asymptotic values of the electromagnetic fields in the classical wave equations.The conversion probability inferred from the wave equations is then captured by the ratio of the amplitudes of the asymptotic ingoing axion field and outgoing photon field.However, the asymptotics of classical field equations are precisely captured by the S-matrix of tree-level quantum field theory (see Fig. 1).Furthermore, the Schwinger-Dyson equations introduced in Sec. 2 directly generate the collision integrals, and therefore provide a means to read off these S-matrix elements.The physical width of the resonance is then determined by finding the roots of the delta function in Eq. (6.1), which gives factors proportional to 1/|∇ x E γ |.Hence, although the kinetic equation does not contain information about the transient evolution of the wave equations through the resonance, it does fully capture the asymptotic behaviour of solutions, and therefore provides a powerful and elegant route to computing the conversion probability that would be derived if a correct solution to the classical equations could be sought.
• Divergences Our final result is free from divergences, which were thought to arise as a consequence either of a breakdown of stationary phase or as a consequence of a breakdown of WKB [56].Indeed, in a previous work [55], some of us argued that the naive divergence observed when k is perpendicular to ∇ x E γ arose from a breakdown of stationary phase due to so-called degenerate stationary points [100].Such divergences occur because the conversion probability P aγ is proportional to 1/|k • ∇E γ |.However, we have now shown that such divergences are naturally regulated by the phase-space measure, as explained in the paragraph immediately preceding Sec. 4. Furthermore, our calculation does not actually involve any stationary phase approximation at all.Hence, these "perpendicular" divergences do not require arbitrary regulation, but are naturally regulated by the phase-space measure.The exception is when |∇ x E γ | itself becomes small.However, this represents a breakdown of perturbation theory, as we discuss below.
• Breakdown of perturbation theory and the adiabatic limit The one divergence that persists in our result occurs when |∇ x E γ | becomes small compared to terms in the numerator.Here, there is no regulatory behaviour from the phase-space measure.Instead, this divergence represents background gradients becoming large, i.e., the width of the resonance becoming large.As such, the resonant interaction can no-longer be treated as a short-lived perturbation, and the conversion becomes adiabatic, as discussed in Ref. [47].In this case, perturbation theory breaks schematically when In this adiabatic limit, one must derive a resummed conversion probability involving an expression to all orders in g aγγ .In the case of 1D conversion, some of us showed in Ref. [47] that the conversion probability for axion-photon mixing in the adiabatic limit was given by resumming the perturbative result P aγ → 1 − exp(−P aγ ).This expression is the well-known Landau-Zener result [101].Physically, this effect corresponds to a process in which the resonance persists for so long that, inside the resonant region, axions convert into photons and then convert back into axions again, with this process repeating ad infinitum to higher and higher orders in g aγγ .As pointed out in Ref. [94], such effects can be important in strongly magnetised media, such as occur around magnetars.In that instance, one would have to derive a generalised Landau-Zener result for 3D mixing valid to all orders in g aγγ .This would require finding a non-perturbative solution to our kinetic equations, which generalises the Landau-Zener result to 3D.Undertaking such a calculation lies beyond the scope of the present work, but it is of importance if one is to place limits for larger values of g aγγ using strongly magnetised stars.
In summary, we have provided a robust calculation for resonant mixing of photons with axions in 3D media.This can then be directly inserted into numerical routines or analytic calculations.This may help to settle one of the major open theoretical questions in indirect searches for axions in astrophysical environments, leaving mainly astrophysical uncertainties (e.g., magnetosphere modelling) as the remaining area for improvement in the robustness of future indirect axion searches.
The resulting transforms are functions of the phase-space variables (k, X); X can be thought of as the position and k the momentum.The Wigner transform of convolutions of two-point functions A and B satisfies where X = (u + v)/2 is the average coordinate and the diamond operator is defined by: The superscripts (1) and ( 2) indicate action on the first and second argument in the round brackets.Thus, if we were to perform a Wigner transform of Eq. (2.1) or (2.6), we would arrive at an infinite series of derivatives arising from background gradients.Note that, when expanding to first order in gradients, we can write where we made use of the Poisson bracket defined in Eq. (2.10), which we will use frequently in the subsequent approximations.We further note that the Wigner transform of the kinetic operator D 0 in Eq. (2.2) can be written as with where the subscript X in Eq. (A.8) indicates the presence of X-derivatives that act on D <,> (k, X).With these specifications, the Wigner transform of the Kadanoff-Baym equation (2.6) reads where the operator ♢ is defined in Eq. (A.6).Following Refs.[69,70,102,103], to obtain kinetic and constraint equations, we take respectively the anti-Hermitian and Hermitian parts of Eq. (A.10), which yields .11)where square brackets on the indices k and X indicate antisymmetrisation defined by refraction of photons in combination with the right-hand side of Eq. (A.11), which corresponds to their production.Retention of gradients of Π H therefore allows us to fully capture refraction of photons within the production mechanism, addressing issues of so-called "dephasing" raised in Ref. [54].
Taking the constraint equations (lower sign in Eq. (A.11)) and neglecting the gradients and the collision terms on the right-hand side, we obtain For general D <,> , i.e., with possible correlations between the physical polarisation states, the commutators do not vanish.They only do so when D is polarisation-diagonal, as we shall assume eventually.The form (A.14) of the constraint equation can be viewed as an inhomogeous equation for D <,> sourced by [Π <,> , D H ] + .Formally then, the most general solution consists of the sum of the general solution to the homogeneous and a particular solution satisfying the inhomogeneous equation.The source term involves the Hermitian function D H that can be derived, in turn, from the causal, i.e., retarded and advanced propagators D R,A via Eq.(2.3).The equations for the causal propagators on the closed time path involve the spectral selfenergy Π A = Π > − Π < .Now, the Π <,> are dominated by absorptive plasma effects Π <,> pl that we do not consider in the present work.
We can nonetheless solve the inhomogeneous equation in the zero-width approximation without fully specifying Π A .For that purpose, we note that Π µν A must be consistent with the boundary conditions for the causal propagators, i.e., one may approximate its leading effect through an appropriate pole prescription.Further details are not needed because we know that any special solution to the inhomogeneous equation must be consistent with the spectral function (A.15) The latter can here suitably be obtained in the zero-width approximation.Namely, we can solve the equation for the retarded and advanced propagators in Wigner space and in the approximation of zero width (i.e., Π µν A → 0) where the eigenstates and eigenvalues of D are defined in Eq. (2.14), σ is the signum function, and η is positive and infinitesimal.Substituting into Eq.(A.15), this leads to Now, we turn to the solutions of the homogeneous part of the constraint equation, i.e., the possible contributions to D <,> that satisfy This equation relates the asymptotic flux of the contributions fc to a (surface) source term on the right-hand side.However, the surface term on the right-hand side vanishes if the surface is chosen to be in an asymptotic region outside compactly supported sources from axions.In other words, the volume can be chosen such that Π <,> = 0 on its bounding surface, A, so that the right-hand side vanishes.This means that the flux contribution from fc through an asymptotic enclosing surface vanishes, i.e., c Further, since each integral and integrand in Eq. (A.36) is positive, this sum vanishes if and only if fc vanishes everywhere on the surface.We conclude, therefore, that fc vanishes outside the region of axion sources and therefore does not contribute to asymptotic solutions.Hence, from an observational perspective, although the Poisson bracket terms on the right-hand side of Eq. (A.28) locally give a contribution to the phase-space density of photons, these terms do not contribute to the asymptotic solutions.In particular, these solutions are off-shell (as can be seen from Eq. (A.27), which contains no mass-shell condition for the photon) and are active only within the axion source.This presumably corresponds to offshell transient solutions in the classical axion electrodynamics equations, but these do not propagate to infinity away from axion sources.We emphasise that these arguments are a consequence of the Born approximation used throughout this paper, i.e., they only remain true to leading order in g aγγ .On-shell and off-shell solutions could of course mix at higher order in g aγγ .Furthermore, the discussion above means that great care should be taken when identifying asymptotic states when comparing classical wave equations to the kinetic equations used here.
Since we are interested only in those contributions to the asymptotic photon solutions away from the site of axion conversion, we can effectively isolate those parts of Eq. (A.28) which determine f 0 c , defined in Eq. (A.29), as the piece that survives into asymptotic regions away from production.This is tantamount to using an effective kinetic equation This equation therefore captures those on-shell solutions that can propagate away from the axion source, which is the main result appearing in Eq. (2.15) of Sec. 2.

B Eigenmodes in Strongly Magnetised Plasmas and the Spectral Function
In this appendix, we analyse the spectral structure of the operator (2.12), defined by which inverts the retarded Green's function D R according to (see Eqs. (2.13) and (2.12) of the main text) In Sec. 2, we showed that these relations give rise to the spectral decomposition presented in Eq. (2.16) of the main text.Our purpose here is to demonstrate that, for a strongly magnetised plasma, three of the four eigenvalues H c have zeros, which correspond to distinct physical states in the plasma.Importantly, we show that two different eigenvalues do not share a common physical zero.Throughout, we use the same coordinates appearing in Eq. (F.3).
The eigenvalues H c of D for an arbitrary R ξ gauge cannot be written down analytically.However, their product, which gives the determinant of the matrix, has a simple expression where i labels the eigen-energies ω i given by Notice that the prefactor 1/ξ in Eq. (B.3) prevents the matrix from being singular.This is because removing the gauge-fixing term by taking ξ → ∞ allows for a trivial zero eigenvalue solution D µν k ν = 0.This follows from Π µν H k ν = 0 and the projection piece in the first two terms of Eq. (B.1) annihilating k µ .We see then that setting the determinant to zero gives the physical eigenmodes.Importantly, these are independent of ξ, as required by the fact that the physical shell structure should be gauge independent.
Recall that, in Eq. (2.19) of the main text, we used the spectral decomposition of the Wightman functions where we interpreted the f c as corresponding to the phase-space density of the plasma eigenmode labelled by c.To justify this, we need to ensure that the zeros of each H a correspond to only one physical state, otherwise different terms in this sum in Eq. (B.7) will mix different states, spoiling the interpretation of the f c as corresponding to the occupancy of different physical modes.To do this, we must ensure that the eigenmodes are in one-to-one correspondence with the (zeros of) H a , so that no two H a share a common root.Formally, we should not have H a (ω i , k) = 0 and H b (ω i , k) = 0 for a ̸ = b for any ω i .
To prove this is the case for the system above, suppose that we have identified a root k 0 = ω i such that, for some H a , we have The question then arises as to whether or not there is another H b̸ =a for which this is also a root.Consider the sum of all possible products of eigenvalue triples We know that any term in the sum containing H a vanishes when placed on the ω i −shell.
Since there are four eigenvalues, there is only one term left on shell, which is the product of the other three eigenvalues, i.e., If this term is non-vanishing, we know that no other eigenvalues have a common root with H a .However, we also know that products of eigenvalues can be obtained from the moments of the characteristic polynomial of D defined by It follows that by differentiating with respect to s and setting s = 0 Hence, if we can verify that Eq. (B.12) is non-vanishing when k 0 = ω i for each ω i , we have shown that no two eigenvalues H c share a common root.We find, using the same coordinates as Eq.(F.3) and the calculations given in Appendix F, that the eigenmodes satisfy These expressions are clearly non-zero in general, meaning that the H c do not share common roots.They only vanish for particular values of ω p and θ.Clearly, they vanish if ω p = 0, since all modes become degenerate in vacuum.Then there is the special case θ = 0, which corresponds to a level crossing between the mt and the other two modes, where k LO,Alfén = ω.Barring level crossings, the eigenvalues H c have distinct roots, with each corresponding to one of the mass-shell conditions above.This can be seen in Fig. 5, where we see that the zeros of the H c do not overlap and that one of the eigenvalues (which we label H ξ ) is never zero.The latter arises due to the presence of the gauge term involving ξ, which regulates the determinant, preventing it from vanishing off-shell.Hence, there are only three zeros, corresponding to the three physical propagating states in the medium.We have therefore confirmed that the zeros of the H c are disjoint, justifying the physical interpretation that each f c in Eq. (B.7) corresponds to the distribution functions of different physical modes.

C Poisson Brackets
In this appendix, we derive the Boltzmann equation (2.21) quoted in the main text.To do so, we show that inserting the spectral Ansatz Eq. (2.19) for the photon Wightman function gives rise to Eq. (2.21).Substituting Eq. (2.19) into Eq.(C.1) gives Eq. (2.20) quoted in the main text where, as in earlier appendices, we allow ourselves to use a as a polarisation label, since here it cannot be confused with axions.Expanding the Poisson bracket by the product rule and using ε a * • ε c = δ ac , we find To simplify this expression, we take which can be differentiated to yield From this, the following three properties are obtained: We prove identities I-III at the end of this appendix but, for now, we take them as read and proceed with the calculation.Using these identities, we see that any terms involving derivatives of the polarisation vectors vanish in Eq. (C.3), so that where in the last equality we have used the product rule for the Poisson bracket and the trivial fact that {H a , δ(H a )} = 0. Inserting the result (C.9) into Eq.(C.2) gives This is precisely the expression (2.21) of the main text.As promised, in the remainder of this appendix, we prove the identities I-III.

Identity I
To prove this identity, we begin by multiplying Eq. (C.5) from the right by ε a and use (C.4) to obtain The second terms on each side cancel, leaving one with as required.

Identity II
We first multiply Eq. (C.5) by a different eigenvector ε c from the right, which gives Here, we used the orthonormality condition ε a * • ε c = δ ac in deriving the first term on the right-hand side, and also the identity D • ε c = H c ε c to obtain the second term on the left-hand side.We can then multiply this by δ(H c ), which eliminates the second term on the left-hand side, yielding Finally, we can multiply this by ε where, in the first term on the right-hand side, we used δ ac to set a = c in ε a • ∂ x ε c * .The first term on the right-hand side of Eq. (C.15) vanishes, because ε a • ∂ x ε a * = ∂ x (ε a • ε a * )/2 = 0, since ε a • ε a * = 1.This leaves However, the second term on the right-hand side is clearly symmetric under interchange of k and x.As such, it vanishes by virtue of antisymmetrisation, so that This completes the proof of identity II.

Identity III
We begin by contracting Eq. (C.5) with ∂ x ε a , which gives Clearly, the first term on the right-hand side vanishes due again to ε a * •∂ x ε a = ∂ x (ε a * •ε a )/2 = 0, from which it follows that Similarly to the proof of the previous identity, the two terms on the right-hand side are symmetric under interchange of k and x, which follows in the second term from the Hermitian nature of D. Hence, upon antisymmetrisation, they both vanish, and we have completing the last of the three derivations for the identities above.

D Gauge Invariance of the Stored Energy Term
In this appendix, we show that E 2 ∂ k 0 H is gauge invariant.This result was required to derive the result (2.42), which recasts the Boltzmann equation in terms of physical electromagnetic fields.To establish the gauge invariance of E 2 ∂ k 0 H, we begin with the identity I proved in Appendix C, namely where we recall, from the discussion preceding Eq. (2.29), that E has arbitrary normalisation.Now, consider a gauge transformation (2.31) to a new polarisation vector The left-hand side of Eq. (D.1) then reads We can rewrite the second term on the right-hand using where η is the Minkowski metric, carrying the spare index associated to the derivative.The second term on the right-hand side of the equation above vanishes on-shell, since E * • D = 0.
Next, since we wish to allow for arbitrary gauges, we should remove the covariant gauge term by taking ξ → ∞ in Eq. (2.12).We then have that, even for off-shell

E Axion Collision Term
The axion contribution to the photon Wightman self-energy Π < , appearing in Fig. 1, can be read off from the interaction vertex (2.46).It is From Eq. (A.3), we have so that the Wigner transform sends The derivatives ∂ X are by assumption sub-dominant in the gradient expansion and can be neglected.Similarly, by Taylor expanding about s µ = 0, we have where the second term is again gradient suppressed, since, in Wigner space, they are O(∂ k • ∂ X ) = O(1/(kL)) ≪ 1, where 1/L sets the size of background gradients and k is a typical momentum scale.Similarly, we have F µν ext (y) ≃ F ext µν (X).Next, using the Wightman function for a spin-zero field [74] G < (k, X) = 2πδ(k 2 − m 2 ϕ ) {f ϕ (k, X)θ(k 0 ) + θ(−k 0 ) [1 + f ϕ (−k, X)]} , (E.5) where f ϕ is the axion distribution function, we arrive at the following expression for the photon Wightman self-energy in Wigner space, quoted in Eq. (2.47): (E.6)

F Conversion Probability in Covariant Gauge
In this appendix, we reproduce the same expression for the conversion probability (4.9), using covariant gauge.Our starting point is the covariant form of the conversion probability in Eq. (3.7) To compute the polarisation 4-vector ε, we require the general polarisation 4-tensor appearing in Eq. (2.12).This can be inferred from the components Π ij H , which can be read off from the classical theory via Eqs.(4.3) and (4.4).The remaining components Π 0i H and Π 00 H can then be read off [87] by using the transverse condition We can now use this to solve the following eigenvalue equation in a covariant gauge: Solving this equation on-shell, we find the polarisation 4-vector of the LO mode to be (up to an overall normalisation) ε LO µ = sin θ cos θ, 0, It is straightforward to verify, using E i = F 0i = ∂ 0 A i − ∂ i A 0 , with A µ ∝ ε LO µ , that this reproduces the expression for the physical electric field polarisation (4.8).Using Eq. (F.6), we compute the matrix-element-like term in the numerator of Eq. (F.1), viz., where εLO ν is the unit 4-vector defined by εLO ν = ε LO ν /(ε LO * ρ η ρσ ε LO σ ) 1/2 .Next, we require the expression for ∂ k 0 H, which can be written on-shell as Upon computing the right-hand side of (F.8), we find Note that Eqs.(F.7) and (F.9) are independent of the gauge parameter ξ.We now divide Eq. (F.7) by Eq. (F.9) and insert the result into Eq.(F.1) to obtain This reproduces exactly Eq. (4.9).Thus, we have shown, via a different calculational route and in a Lorentz covariant gauge, the same form of the probability in Eq. (4.9), which was derived using temporal gauge and physical electromagnetic fields.Note that the approach in this appendix also circumvents any physical discussion about the stored energy in the modes, which came from substituting ∂ k 0 H for the expression (2.36) involving the stored energy in the plasma.Instead, we have worked directly with the eigenvalues of the matrix in Eq. (F.5).

Figure 1 .
Figure 1.Contribution to the photon Wightman self-energy Π < ax (x, y) due to axion interactions.The optical cut through the axion line gives the standard inverse Primakoff process a → γ.

Figure 2 .
Figure 2. Worldline method for computing axion-photon conversion.The phase-space density f γ is integrated along the photon worldline (red) with affine parameter λ.The axion worldline is shown in blue.Production occurs at the point E γ k(λ), x(λ) = E a k(λ) on the kinematic surface Σ k with unit normal n ∝ ∇E γ .The axion and photon 4-momenta are equal at the resonance k µ a = k µ γ , but their group velocities, v γ g and v a ∝ k, which are tangent to worldlines, are different, since in an anisotropic medium, the photon 3-momentum and group velocity are not parallel [80].

1 Figure 3 .
Figure3.The integrated conversion probability ⟨P aγ ⟩ = dΩ k v p cos θ n P aγ evaluated on the critical surface r = r c (θ NS , φ NS ) given by m a ≃ ω p .We choseφ NS = π, θ m = 0.4, t = 0, m a = 10µeV and B = 10 14 G, Ω = 1Hz, R = 10km with ω = m a (1 + v 2 a ) 1/2, v a = 220km/sec and g aγγ = 10 −14 GeV.We display in blue the 1D isotropic conversion probability (4.2) and in green, the new result in Eq. (4.9) of this work.Gray regions show where the production surface would be inside the star, where no photons are produced from axions.

1 Figure 4 .
Figure 4. Angular dependence of the anisotropic conversion probability in Eq. (4.9) (green), as a function of polar momentum angle θ k of the axion.We fixed ϕ k = 0.2 and neutron star angles θ NS = 0.5 and φ NS = 0.6π.Other values are as in Fig.3.The blue curve shows how the divergences in the probability are regulated by the phase-space measure, as discussed in Sec.3.1.The gray dashed line indicates the pole in the conversion probability at cos θ n = 0, when the incoming axion momentum is tangent to the conversion surface.
k µ , D(ξ → ∞) • k = 0. (D.5)From this, it follows that the first term on the right-hand side of Eq. (D.4) also vanishes, leaving us withE ′ * • (∂ k D) • k = 0, (D.6)establishing that the second term on the right-hand side of Eq. (D.3) vanishes.Returning now to the last term on the right-hand side of Eq. (D.3), we see that it can be written ask • (∂ k D) • k = k • ∂ k (D • k) − k • D • η = 0, (D.7)where both terms vanish by making use again of the relation (D.5).Thus, having established Eqs.(D.6) and (D.7), we see that both the second and third terms on the right-hand side of Eq. (D.3) vanish, leading toE * • (∂ k D) • E = E ′ * • (∂ k D) • E ′ .(D.8)This establishes the gauge invariance of E * •(∂ k D)•E.Furthermore, since, from Eq. (D.1), this is equal to ∂ k HE 2 , it follows that ∂ k 0 HE 2 is also gauge invariant, proving the result quoted in Sec.2.1.