The Fourier formalism for relativistic axion-photon conversion, with astrophysical applications

We study the weak mixing of photons and relativistic axion-like particles (axions) in plasmas with background magnetic fields, ${\bf B}$. We show that, to leading order in the axion-photon coupling, the conversion probability, $P_{\gamma \to a}$, is given by the one-dimensional power spectrum of the magnetic field components perpendicular to the particle trajectory. Equivalently, we express $P_{\gamma \to a}$ as the Fourier transform of the magnetic field autocorrelation function, and establish a dictionary between properties of the real-space magnetic field and the energy-dependent conversion probability. For axions more massive than the plasma frequency, ($m_a>\omega_{\rm pl}$), we use this formalism to analytically solve the problem of perturbative axion-photon mixing in a general magnetic field. In the general case where $\omega_{\rm pl}/m_a$ varies arbitrarily along the trajectory, we show that a naive application of the standard formalism for 'resonant' conversion can give highly inaccurate results, and that a careful calculation generically gives non-resonant contributions at least as large as the resonant contribution. Furthermore, we demonstrate how techniques based on the Fast Fourier Transform provide a new, highly efficient numerical method for calculating axion-photon mixing. We briefly discuss magnetic field modelling in galaxy clusters in the light of our results and argue, in particular, that a recently proposed 'regular' model used for studying axion-photon mixing (specifically applied to the Perseus cluster) is inconsistent with observations. Our formalism suggest new methods to search for imprints of axions, and will be important for spectrographs with percent level sensitivity, which includes existing X-ray observations by Chandra as well as the upcoming Athena mission.

Determining the elementary particle content beyond the established Standard Model is a central goal of contemporary high-energy physics. The QCD axion (cf. [1]) and axion-like particles (axions) comprise a wellmotivated class of hypothetical particles that frequently appear in extensions of the Standard Model, including effective theories derived from string theory [2]. Both the QCD axion and axions can be understood as pseudo-Nambu-Goldstone bosons of broken, approximate symmetries, and the QCD axion provides the leading candidate solution to the strong CP-problem. In this paper we will simply use 'axions' to refer to the QCD axion and axion-like particles, as our results apply equally to these particles.
Axions can naturally be very light, with feeble couplings to matter and radiation, and provide an increasingly popular candidate for explaining the nature of dark matter [3][4][5]. Characteristic of axions is their coupling to electromagnetism through the Lagrangian term L int = 1 4 g aγ aF µνF µν where g aγ denotes the axion-photon coupling and the axion, a, is assumed to have a mass m a . This term induces a mixing between the photon and the axion in backgrounds with non-vanishing electromagnetic fields, opening the possibility to interconvert axions and photons. Such interconversion underpins the majority of the experimental and observational efforts to search for axions (cf. [6][7][8][9][10]). Particularly powerful are searches for axion-induced distortions in the spectra of luminous X-ray and gamma-ray sources located in galaxy clusters . State-of-theart analyses using high-quality data from the Chandra X-ray observatory have bounded axion-induced spectral distortions from the AGN at the centre of the Perseus cluster to the few percent level, leading to some of the strongest limits on the axion-photon coupling to date (g aγ ≤ 7.9 × 10 −13 GeV −1 for m a < 10 −12 eV [19]). The observational prospects are very good for even more sensitive searches for axions using existing and planned X-ray and gamma-ray telescopes, such as Athena [37] and the Cerenkov Telescope Array [38]. However, such studies will be limited by modelling uncertainties affecting the axion-photon mixing. The conversion probability is in general a complicated function of the mode energy, the axion parameters m a and g aγ , as well as the plasma density and magnetic field along the particle trajectory. The robustness of the predictions to astrophysical modelling uncertainties has been investigated by several groups over the past decade [22,31,[39][40][41], typically finding that limits on g aγ can change by a factor of a few depending on the magnetic field model employed. A recent, unusual contribution to these studies is reference [41], which investigated the robustness of certain gamma-ray constraints assuming -as a limiting but ostensibly observationally consistent case -that the cluster magnetic field in Perseus is highly regular, finding much weaker limits than those of [24]. 1 Common to all these studies is the reliance on very simple magnetic field models to generate physical intuition (typically, a constant magnetic field), and numerical simulations to examine more involved models. This has left room for some confusion as to what properties of the magnetic field really drive the appearance of features in the conversion probability, and how robust the predictions really are.
In this paper, we revisit the theory of axion-photon conversion, and develop a powerful new method calculating and interpreting axion-photon mixing. Our results are equally valid for applications in the laboratory and space, but our main focus is on astrophysical applications involving conversion of high-energy photons into relativistic axions. Our approach is inspired by -and will be particularly useful for -the emerging sub-field of precision X-ray searches for axions [14][15][16][17][18][19]36], which presently constrain the mixing probability to be no larger than a few percent. For such weak mixing, the system is in the perturbative regime, and the goal of this paper is to demonstrate the significant conceptual, calculational and methodological advances that can be gained using perturbation theory.
It is well-known that the classical, linearised axionphoton system can be written as a Schrödinger-like equation, with time replaced by the spatial coordinate, say z, along the direction of propagation [42]. In quantum mechanics, the asymptotic, perturbative transition amplitude can be expressed as a Fourier transform of the interaction Hamiltonian. In this paper, we show that a similar, but more subtle, statement also holds for classical axion-photon mixing.
In the simplest case, which is directly relevant for gamma-ray searches, the mass of the axion is always larger than the plasma frequency. We find that, to leading order in g aγ , the axion-photon transition amplitude involving such massive axions is given by a sum of Fourier cosine and sine transforms of the (relevant component of the) magnetic field, B. The conversion probability is then given by the power spectrum of B. We derive a version of the Wiener-Khintchine theorem that shows that the conversion probability is equal to the Fourier cosine transform of the magnetic autocorrelation function. Strikingly, this means that questions about the spectrum of oscillations induced by axion-photon con-version map directly onto questions about the real-space properties of the magnetic autocorrelation. This constructively answers the questions of what properties of the magnetic field are reflected in the conversion probability. We demonstrate how to apply this formalism in a series of examples, ranging from the simple to the general. In particular, we analytically calculate the conversion probability and the magnetic autocorrelation function for a magnetic field expressed as a Fourier series. Since any physically relevant magnetic field can be expressed in such a way, this explicitly solves the problem of weak axion-photon mixing, for sufficiently massive axions.
Somewhat more subtle is the case of massless axions. In this case, the transition amplitude is no longer given by a Fourier transform obtained by integrating over the spatial coordinate z, as in the massive case. However, the amplitude can be written as a Fourier transform obtained by integrating over the variable Again we find that the transition amplitude is given by cosine and sine transforms, but this time of the function G = B/ω 2 pl . The conversion probability is given by the power spectrum of G, or equivalently by our version of the Wiener-Khintchine theorem, as the cosine transform of the autocorrelation function of G (in ϕ-space).
In general, there may be some regions where m a < ω pl , and others where m a > ω pl . In this case, it is not possible to express the transition amplitude as a single Fourier transform. However, we find that the transition amplitude reduces to a sum of 'non-resonant' contributions that are given by Fourier transforms of G (one for each region where m a − ω pl has a definite sign), and a sum of 'resonant' contributions (from points where m a = ω pl ).
An often used method to analytically evaluate resonant axion-photon mixing is the stationary phase approximation. We show that, in the relativistic limit, the resonant amplitude calculated from the stationary phase approximation is enhanced by a factor of ω/m a , which is very large for light axions emerging from X-rays or gamma rays. However, we also show that a naive application of the stationary phase approximation to the resonant conversion gives a result that is highly inaccurate. We address this issue by deriving a modified form of the stationary phase approximation that is relevant for relativistic axion-photon conversion: with this new formula, we see that in careful calculations of the amplitude, the nonresonant contributions are generically at least as large as the resonant contributions.
Fast Fourier Transforms (FFTs), and related methods, have revolutionised digital signal processing by providing highly effective ways to numerically evaluate the discrete Fourier transform. We show that with our formalism, axion-photon conversion can be evaluated using methods based on FFTs. This drastically reduces the computational effort in determining the effects of axion-photon mixing, and allows for effective marginalisation or Monte Carlo over astrophysical magnetic fields and plasma densities. See also [43] for a discussion of FFT techniques applied to axion searches.
Since magnetic auto-correlations determine the predictions of perturbative axion-photon conversion, it is important to re-examine what is known about astrophysical magnetic fields in the environments most promising for axion searches. To this end, we discuss the expected properties of magnetic fields in galaxy clusters. We review the arguments for turbulence in the intracluster medium (ICM), and the classes of models used in astrophysical axion searches. In particular, we critically examine the recently proposed 'regular model' of [41], finding it at odds with observations.
We suggest that future studies of the magnetic field autocorrelation function in state-of-the-art magnetohydrodynamic simulations of the ICM will be very useful in improving the sensitivity of axion searches. Moreover, our formalism could lead to new methodologies for axion searches, e.g. by generating the relevant mixing probabilities directly from an observationally inferred class of autocorrelation functions, without ever explicitly solving the Schrödinger equation.
This paper is organised as follows: in section II we review the classical theory of axion-photon conversion, and how it can be understood through non-relativistic quantum mechanics. We analytically develop the new formalism in sections III-V, for the cases of a comparatively heavy axion (section III), a very light axion (section IV), and the intermediate, general case (section V). We then test relevant aspects numerically in section VI, and discuss numerical implementations using discrete, Fourier-like transforms. In section VII, we discuss what is known about magnetic autocorrelations in the relevant astrophysical environments, and we indicate new possible directions for the methodology of axion searches. We conclude in section VIII.
The most important new results of this paper are (in no particular order) i) the relation between the non-resonant conversion probability and the magnetic autocorrelation function given by equations (30), (79) and (115); ii) accounting carefully for both resonant and non-resonant contributions when ω pl = m a at one or more points along the trajectory, cf. section V B; iii) the demonstration that one can use highly efficient FFT methods to calculate the conversion probabilities, cf. section VI B; iv) the identification of new, promising methods that are made possible by the new formalism, and which can open up new directions for axion searches, cf. section VII C.
to a Schrödinger-like equation from which the transition amplitudes can be calculated order-by-order in perturbation theory [42], in direct analogy with time-dependent perturbation theory in quantum mechanics. The axionphoton interaction is described by the following Lagrangian where a is the axion field, A µ is the photon field, F µν is the electromagnetic tensor,F µν = 1 2 µνλρ F λρ . We have neglected the effects of Faraday rotations and QED birefringence at low energy. 2 This equation of motion are the Klein-Gordon equation for the axion field, and an extension of Maxwell's equations. The linearised equations around a static background magnetic field B 0 are given by where we included the plasma frequency ω pl . Most works on axion-photon mixing (from [42] to the more recent X-ray and gamma-ray studies, see e.g. [46]) further simplify this system by assuming that all background quantities only vary along the z-direction, which is taken to be the axion or photon propagation direction. 3 We consider a right-moving plane wave in the z-direction, and write the components of the background magnetic field respectively as B x and B y . Equations (2) now become In the relativistic limit, the equations of motion can be reduced to first order by appealing to the rotating wave approximation: ω 2 +∂ 2 z = (ω +i∂ z )(ω −i∂ z ) 2ω(ω −i∂ z ). 4 After a shift of the axion field a → −ia, the equations of motion are 2 The Faraday effect is inversely proportional to the energy, and completely negligible at X-ray energies or higher. Moreover, at X-ray energies, the QED birefringence contribution would in a typical galaxy cluster environments with µG magnetic fields and plasma frequencies of the order of ω pl ∼ 10 −12 eV lead to a fractional correction to ∆γ , that we define in (8), of order 10 −16 . The QED contribution scales like E 2 B 2 , and is less suppressed, and sometimes even non-negligible, in environments with strong magnetic fields probed at very high energies (cf. [45] for a recent example). 3 For a recent discussions of extensions to anisotropic plasmas, see [44,47]. 4 An alternative method to arrive at the first-order equations is to use the WKB approximation, which is applicable also in the non-relativistic case [44].
These classical mixing equations are now of the form of the Schrödinger equation, with time replaced by the spatial coordinate z [42]: Here the axion field and the components of the vector potential are components of the 'Schrödinger-picture state vector' where we have suppressed the dependence on the mode energy ω. The basis vectors are assumed to be zindependent, but the coefficients of the state vector evolve with z. The Hamiltonian is decomposed into free and interaction parts and where j = x, y. 5 The formal solution to the Schrödingerlike equation is and where P z denotes the path-ordering operator. The z-evolution operator U (z, 0) is implicitly energy dependent. Numerical solutions can be found by discretising the z-direction into a sufficiently large number of 'cells' so that the Hamiltonian is approximately constant in each cell, and the total evolution operator U is the product of the evolution operators for all cells, appropriately ordered. Such a numerical approach is very common in astrophysical searches for axions, but becomes computationally costly when the mixing occurs over large distances with non-trivially varying magnetic fields, when the energy resolution needs to be finely sampled, and when a large number of possible magnetic field configurations must be considered. However, the structure of equation (5) is suggestive of different approach: perturbation theory in direct analogy with time-dependent perturbation theory in quantum mechanics. This is most conveniently studied in the interaction picture, where the equations of motion are given by where The general solution can be expressed as an expansion in the interaction Hamiltonian, which equivalently corresponds to an expansion in the coupling constant, g aγ . To linear order in perturbation theory, the state in the interaction picture is given by In the original Schrödinger-picture basis, this corresponds to The perturbative expansion can be expressed using Feynman diagrams, where each order in H I corresponds to an additional axion-photon mixing vertex, and where the zeroth-order evolution operator U 0 propagates the state between the vertices.
We are interested in the transition amplitude of an initial photon state emerging as an axion. This corresponds to the initial condition Ψ(0) = (1, 0, 0) T if the photon is linearly x-polarised, and Ψ(0) = (0, 1, 0) T if its y-polarised. With j = x, y the leading order (LO) transition amplitude is then given by [42] The conversion probability is obtained by squaring this transition amplitude: P γj →a = |A γj →a | 2 . This conversion probability calculated from the (classical) Schrödinger-like equation gives the ratio of the squared axion field to the squared electric field, which also corresponds to the observationally relevant flux ratio of axion and photons.
Next-to-leading-order (NLO) corrections appear at order O(g 3 aγ ) in the amplitude from the 3-vertex Feynman diagram: 'photon → axion → photon → axion . The LO conversion probability appear at O(g 2 aγ ), with higher-order corrections appearing at O(g 4 aγ ). Thus, the perturbative expansion is generically a good approximation when the conversion probability is small, with fractional corrections at the same order as the conversion probability (i.e. a 1% perturbative conversion probability generically receives corrections at the order of 0.01%, cf. section VI D for a detailed discussion).
The above equations lead to unitary time evolution, and so neglects possible photon absorption. This is a good approximation when the medium is optically thin over the region at which axion-photon mixing occurs, which is the case for many astrophysical applications of practical relevance, cf. e.g. [16]. However, including photon absorption is straightforward and leads to a modified Schrödinger-like equation that features a non-Hermitian Hamiltonian: encodes the damping. The the real function β depends on the absorption cross-section and astrophysical parameters, such as the free electron density. Since [D(z), H 0 (z )] = 0, it is straightforward to include damping as a zeroth-order modification to the perturbative expansion: equation (10) in the interaction picture still holds, but with the unitary operator U 0 now replaced by the non-unitary transfer matrix so that Ψ int (z) = T −1 0 (z, 0)Ψ(0) and H int = T −1 0 H I T 0 . This way, photon absorption modifies the zeroth-order propagator, and to linear order the Schrödinger-picture state is given by: In the following, we will neglect photon absorption.
In this paper, we will almost exclusively focus on polarised transition probabilities, involving a single component of the magnetic field. However, many bright astrophysical sources that are suitable for searches for axions are unpolarised. The conversion probability for an upolarised source of photons is P γ→a = 1 2 P γx→a + P γy→a , and the survival fraction of the unpolarised flux is In the following, we will primarily consider the conversion probability for a linearly polarised photon (with respect to a fixed, Cartesian coordinate system), with the understanding that the unpolarised probability can easily be obtained using formula (16). For a discussion on polarised signals from axion-photon conversion, see [48].
Finally, and for context, we note that solving the Schrödinger-like equation for axion photon mixing is equivalent to solving the von Neumann equation for the corresponding density matrix, ρ, as is done in much of the literature. The equation of motion of ρ is, and the diagonal elements of ρ are interpreted as the square of the wavefunctions of the two photon polarisation and the axion, respectively. The time-evolution of the density matrix is given by for the evolution operator of equation (9). Numerical solutions in this formalism proceed equivalently to the Schrödinger equation (i.e. by finding the z-evolution operator), and there is no significant conceptual or calculational difference between the two formalisms. In particular, photon absorption is accounted for in an identical manner in the two formalisms, by replacing the unitary evolution operator by a non-unitary transfer matrix, and the perturbative expansion of U can be applies equally to the two formalisms. For this reason, we will not discuss the density-matrix formulation of this problem any further in this paper.

III. FOURIER TRANSFORM FORMALISM: THE MASSIVE CASE
To leading order in the mixing, the amplitude for an initial photon state, linearly polarised along the xdirection, to transition into an axion is given by where the phase is given by with ∆ γ and ∆ a as in equation (8). The phase Φ may both increase and decrease along the trajectory depending on the relative sizes of ∆ γ (z ) and ∆ a . The key property of Φ that we use in this paper is that it can, in general, be factorised into a generalised spatial integration variable of equation (17)  case) multiplied by an independent, 'conjugate' parameter (say 1/ω). In this section, we focus on the simplest case when ∆ γ can be neglected, and the links to Fourier analysis are most apparent. We note that this factorisation is more subtle in the case of non-relativistic axions, which obey a similar Schrödinger-like equation but for which the factors of 1/ω in (8) should be replaced by the inverse of a spatial wave-vector [44], which depends on the spatial coordinate. In this paper, we focus on the simpler relativistic case.

A. The cosine and sine transforms for massive axions
We first consider relativistic axions with a mass larger than the plasma frequency This is the relevant case for gamma-ray searches for axions [41,49], and we will for brevity refer to it as the case of 'massive axions'.
In this case, we set ∆ γ = 0 and note that ∆ a = −m 2 a /(2ω) is independent of z. The transition amplitude becomes Since we are interested in amplitudes evaluated far away from the transition region, we have extended the integral to infinity, assuming that ∆ x is only non-vanishing in a finite region of space. This amplitude is now a halfsided Fourier transform of ∆ x , with ∆ a being the conjugate 'momentum'. This suggests that we may use Fourier analysis to analyse axion-photon mixing.
Functions defined on the real, positive line can be expressed independently using either sines or cosines. For sufficiently well-behaved functions that decay at infinity, 6 the relevant expansions are given by the Fourier sine and cosine transforms: These transform have the convenient property of being, up to a constant, their own inverses: Moreover, the cosine and sine transforms are clearly real if f (z) is real. The transition amplitude is now straightforwardly expressed through the cosine and sine transforms as: where we have dropped the spatial index on the magnetic field, and where the conjugate Fourier variable is given by which is positive semi-definite. So, in the massive case, transition amplitudes are simple Fourier transforms of the magnetic field profile. Since these transforms are real, the conversion probability is conveniently given by without cross-terms between the cosine and sine transforms. This expression means that the oscillatory pattern that axions induce on astrophysical spectra are directly given by the power spectrum of the relevant magnetic fields. Moreover, for a fixed energy and mass, only a single wavelength of the magnetic field affects the conversion probability: the conversion probability is localised in Fourier space. The mathematical implications of equation (26) can be further elucidated by relating the power spectra to the magnetic field autocorrelation function. To do so, we derive the Wiener-Khintchine theorem, as applied to the cosine transform. We begin by noting that the selfconvolution of the standard, exponential Fourier transform implies the identities: where ∆ o x and ∆ e x respectively denote the odd and even extensions of ∆ x to negative values of the spatial coordinate. We now use that and identify as the auto-correlation function of ∆ x . We now have that This compact equation is one of the main results of this paper. It says that the form of axion-induced modulations in high-energy spectra is encoded in the magnetic auto-correlation function, through a simple transform.
Since the cosine transform is its own inverse (up to a constant), the spectrum of the probability as defined by the cosine transform, i.e. F c (P γx→a ), is directly given by the magnetic autocorrelation at a given spatial length scale. This means that questions about the spectrum of oscillations in the conversion probability map onto sharp questions about the real-space properties of the magnetic field. Through equation (30), we can translate established properties of deterministic autocorrelation functions into general statements about the perturbative conversion probability. For example, the autocorrelation function peaks at zero: c B (0) > c B (L) for L > 0, and decays at large L (at least when the magnetic field is of physical origin, or when the Fourier transform is applicable). This means that the zero mode of P γ→a (η) is always the largest Fourier (cosine) component, and the conversion probability decays as η → ∞.
We close this section by commenting on two additional new relations that follow from this formalism. First, we note that equation (30) implies a new equation for extrema of the conversion probability. We have that P γx→a (η) = 0 whenever This equation can be of phenomenological interest as the distribution of extreme probability points is a measure of how featured the probability curve is.
Second, we note a non-trivial equality for the integrated conversion probability. Plancherel's formulas imply the identities: (32) It follows that the integrated probability over all wavelengths (η) equals the real-space integral over the squared magnetic field:

B. Examples
We now apply the general formalism described in section III A to a series of examples. We find that all types of models that have previously been considered in the literature on axion-photon conversion (i.e. 'cell models', turbulent models defined in Fourier space, and regular models defined from Taylor expansions) can be solved analytically. Indeed, since any magnetic field with a finite extent can be expressed as a Fourier series for which we find the general solution in Example 4, this analysis constructively solves the problem of perturbative axionphoton mixing in the massive case. We note that all figures presented in this section (and similar figures in sections IV and V) correspond to analytic solutions evaluated in the perturbative regime, so agree excellently with the full solution. The simplest example of axion-photon conversion corresponds to a constant magnetic field over a finite range. We take B constant for 0 ≤ z ≤ z max , while being zero elsewhere. The real and imaginary parts of the transition amplitude are given by where η = m 2 a /2ω. Squaring the amplitude, the conversion probability is This result (of course) agrees with the standard singledomain formula to leading order in g aγ where ∆ osc = η 2 + 4∆ 2 x . The alternate way to calculate the transition probability is to compute the magnetic autocorrelation function, and then take its cosine transform, cf. equation (30). The autocorrelation function of a rectangular box is a simple, linear function: Here, we have introduced the step functions and The cosine transform of the magnetic autocorrelation is then given by: Equation (30) then reproduces the correct conversion probability: It is interesting to note how the conversion probability changes with z max : the larger the z max , the more rapid the probability oscillations in η-space. From equation (40) we can infer that this is a general property following from the relation between the conversion probability and the magnetic autocorrelation function: long-ranged autocorrelations map to rapidly oscillating modes in P γ→a (η).

Example 2: General cell models
In an often considered class of 'cell models', the magnetic field is taken to be constant within a series of domains along the particle trajectory. In general cell models, the magnitudes and directions of the magnetic field and the size of each cell are treated as independent variables, that may e.g. be generated through some specified probability distributions (cf. e.g. [16,19,40,43,50]). The benefits of cell models are that they are quick to generate and enable rather efficient numerical evaluation of multi-scale magnetic fields, which is important when marginalising over magnetic field realisations to understand the astrophysical variability of the axion predictions. Due to their close resemblance with the constant domain model, cell models have also been used to gain intuition about general properties of axion-photon conversion by simply summing the probabilities from each cell (as opposed to the amplitudes), neglecting the interference terms [43].
In this section, we determine the conversion probability analytically in general cell models. This allows us to address the questions of the importance of the interference terms, and how the predictions from a a cell-model magnetic field differs from a smoothed version.
The (x-component of the) magnetic field in a general cell model can be defined as with the constant parameters B i and where z i < z i+1 . The real and imaginary parts of the transition amplitude are given by To find the conversion probability, one may of course sum the squares of the real and imaginary parts of the amplitude. However, in this case, it is arguably simpler to directly compute the autocorrelation function and its cosine transform, as we now show.
The magnetic field auto-correlation function is given by, This equation means that in a general cell model, the magnetic autocorrelation function is a continuous and piece-wise linear function of L, with discontinuities in its first derivative. This generalises the well-known result that the autocorrelation function of a step-like signal is a triangle function.   The conversion probability is given by the cosine transform of this function, according to equation (30). Explic-itly, we obtain where d i = z i+1 − z i is the cell length and The conversion probability is obtained as P γx→a = g 2 aγ 2 F c c B , and setting η = m 2 a /(2ω). We identify the first term in equation (45) as an incoherent sum of oscillatory transition probabilities from each cell. For each term, the oscillation frequency in ηspace is simply set by the cell size, d i .
All other terms are due to interference effects, and are also given by simple oscillatory functions (possibly with a prefactor of η). For the interference terms, the frequency of oscillation in η-space is set by the separations between domain boundaries, which need not be adjacent. Thus, interference terms can contribute with a wide range of oscillation frequencies to P γx→a (η), including rapid oscillations. Moreover, unless there is a hierarchy of magnetic field strengths, the interference terms are unsuppressed. Figure 2 shows the effects of the interference terms in equation (45) are particular apparent at low energies (large η), where they generate the fast oscillations in the conversion probability. The magnitude of this discrepancy depends on the number of cells, and here we illustrate just one simple example. However, we believe that an inclusion of the interference effects may further improve interesting analyses that have neglected them, such as [43].
In figure 3 we compare the results of equation (45) for a cell model and the complete numerical solution for the same field smoothed out on the ∼ 1 kpc scale. The corresponding magnetic autocorrelation functions are qualitatively similar, but the auto-correlation of the cell-model is only piece-wise linear, with clear discontinuities in its first derivative. These jagged features of c Bx (L) translate into additional support at large η for P γx→a (η). This is visible also in figure 3d: at high energies (small η) the shape of the conversion probabilities are similar, and mostly differ due to the different magnetic field strengths of the cell model and the smoothed field. At lower energies (large η), both the conversion probabilities are oscillatory, but that of the cell model is more 'featured', and decays more slowly. We conclude that the differences between cell models and smoothed versions of the magnetic field are mostly confined to comparatively low energies. Observations that are only sensitive to axions in the high-energy region where the conversion probability is the largest are likely to be rather insensitive to the differences between cell models and smoothed versions.

Example 3: Single mode oscillatory magnetic fields
A well-known result of [42] is that an oscillating magnetic field can lead to an enhanced axion-photon conversion probability, analogously to magnetic resonance in quantum mechanics, where an oscillating magnetic field provides the additional energy for rapid transition between Zeeman split energy levels. For axion-photon conversion, the basic observation is that an oscillating magnetic field of the form leads to an amplitude involving the terms e i(η+k)R and e i(η−k)R . For η ≈ k and ηR 1, these correspond, respectively, to a rapidly and a slowly oscillating term. Neglecting the former by appealing to the rotating phase approximation, one finds a conversion probability of the form [42] which peaks at η = k. We will now discuss the relation between our more general result of equation (30) and equation (48). Assuming a magnetic field of the form of equation (47), the real and imaginary parts of the transition amplitude are given by where ξ ± = (η ± k)R/2 and η = m 2 a /2ω. The conversion probability is calculated by squaring this amplitude Note that as R increases (keeping k fixed), the interference term becomes negligible and the probability approaches P aγ = g 2 aγ B 2 0 R δ(k − η)/16. Alternatively, one can find the conversion probability by computing the autocorrelation function, (51) and taking its cosine transform. The first term of equation (51) then gives the first two terms of equation (50), including the 'resonant' term for ξ − → 0. This provides a new perspective on the resonance: its characteristic sin 2 ξ/ξ 2 behaviour emerges when the autocorrelation function includes a cosine mode with linearly changing amplitude ∼ L cos(kL)

Example 4: General, 'turbulent' magnetic fields
In astrophysical environments, the components of the magnetic field inevitably involve more than one Fourier mode which leads to interference effects that are absent in the single-mode case. Indeed, 'turbulent' magnetic fields are often modelled as Gaussian random fields, by drawing the amplitudes of a set of modes from some specified power spectrum. 7 We consider the most general magnetic field defined for 0 < z < R as expressed through its Fourier series The wavenumbers are k n = 2πn/R with n = 0, 1, . . ..
The real and imaginary parts of the amplitude for the n-th mode is simply obtained by direct integration, or equivalently by taking the cosine and sine transforms as in (19): where ξ ± n = (η ± k n )R/2 and η = m 2 a /2ω. The conversion probability is then given by The terms proportional to sin 2 (ξ ± n )/(ξ ± n ) 2 , corresponds to an incoherent addition of contributions of the form (48) and would be the only term present in a strict application of the rotating phase   approximation.
The magnetic autocorrelation function of (52) has the general form where the functions C are defined as and Equations (54)-(55) provide the most general solution for relativistic axion-photon conversion of axions with m a ω pl , and equations (55)-(58) provide the corresponding cosine transforms. In figure 4, we show a simple, 7-mode example of the magnetic field, and its solution. We have tested our analytical solution against numerical simulations, and found perfect agreement.

Example 5: Monomial magnetic fields
Sufficiently slowly varying fields are often conveniently Taylor expanded to some finite order in z/R. In this example, we consider the simplest case of a linear dependence on the radius, and in section III B 6 we determine the conversion probability for a general polynomial.
The simplest monomial magnetic field is linear and the transition amplitude is and the conversion probability is

Example 6: 'Regular' magnetic fields
In this section, we consider a broad class of magnetic fields that can be modelled by finite-order polynomials within a finite radius, i.e. as where b(z) can be expressed as a polynomial in z/R. The purpose of this section is to provide an algorithmic prescription for finding the conversion probability of such models analytically, using our Fourier analysis approach.
As an example, we also calculate the conversion probability explicitly in the regular magnetic field model proposed in [41] (however, see our discussion about the astrophysical consistency of this model in section VII B). The transition amplitude involves the sine and cosine transforms of b z R Θ(R − z). It is convenient to define u = z/R and write b(u) = nmax n=0 b n u n . The cosine and sine transforms then involve terms like 1 0 du u n cos(ηu) and whereη = ηR. These integrals respectively, and rather formally, evaluate to instances of the regular and generalised hypergeometric functions. For our purposes, it will simpler to evaluate them by repeated integration by parts, i.e. by using The resulting amplitude then has a simple form where s i and c i denote the polynomials obtained from the sine and cosine transforms, respectively. 8 The conversion probability is simply found by squaring and summing the real and imaginary parts of the amplitude, as given by equation (66). Alternatively, the conversion probability can be found by first evaluating the magnetic field autocorrelation function. Writing = L/R, we have that Each integral in this double sum can be expressed using incomplete Euler beta-functions, or alternatively, by expanding the factors of (u + ) m and doing the integral term-by-term. The latter puts c B in polynomial form: From the autocorrelation function, the conversion probability can be found by taking the cosine transform. For a correlation function of the form (68), this involves hypergeometric functions, however, we again use repeated applications of equation (64) to simplify the expression and determine its general form. We find that where the p i (η) are polynomials. The shape of the conversion probability is then determined by these polynomials and the low-frequency oscillations from sin(η) and cos(η), which typically results in a slow-varying conversion probability.
We now consider a particular example of axion-photon conversion in slowly varying magnetic fields. Recently, reference [41] adapted a model of the magnetic field in a radio bubble from [51] as a model of the cluster magnetic field in the Perseus cluster, suggesting that it provides a limiting case of observationally viable magnetic fields. We will return to the astrophysical issues of this nonstandard, and somewhat controversial, generalisation in   [41,51] and section III B 6 (Example 6). The solid and dashed curves respectively correspond to the φ and θ components of the regular magnetic field, cf. equations (72) and (71). See however section VII B for a discussion of the interpretation of these results.
section VII. Here, we simply show that the methods of this paper can be used to find the perturbative conversion probability of this model analytically. The magnetic field model of [51] was found as a solution to the Grad-Shafranov equation in equilibrium magnetohydrodynamics (MHD). In spherical coordinates, the explicit form of the components of the magnetic field for r ≤ R were determined to where A solves a transcendental equation and is approximately given by A = 5. 76. The function f (r) is given by Here u = r/R ≤ 1 is the re-scaled radial coordinate, and c 1 is a normalisation factor that sets the overall magnitude of the magnetic field. The parameter values adopted in reference [41] in applying this model to the entire Perseus cluster were θ = π/4, c 1 = −0.060 µG, and R = 93 kpc. Since the magnetic field is slowly varying in u, this system is easily solved to high accuracy by approximating the components of B by finite-order Taylor expansions. The main results are summarised in figure  5.
The two paths to the conversion probability are illustrated in figure 5: by taking the cosine and sine transforms of the magnetic field components as in equation (24), we obtain the transition amplitudes for linearly po-larised photons converting into axions, as plotted in figure 5b. Taking the squared norm of the amplitudes give the polarised conversion probabilities, plotted in figure 5d. Alternatively, we directly calculate the magnetic autocorrelation function as shown in figure 5c. The cosine transform of these autocorrelation functions again give the conversion probability.
This example magnetic field model is inconsistent with observations of the Perseus cluster, as we explain in section VII B. However, it is still interesting to interpret it through the Fourier formalism. We have seen that longranged magnetic autocorrelations translate into support for rapid-oscillation modes in P γ→a (η), however, this does not automatically translate into oscillations of the full function P γ→a (η) at those frequencies; Fourier analysis decomposes any function into its oscillatory components, but of course, this does not mean all functions are oscillatory. The bumpy features of the conversion probabilities of figure 5d indicate the broadly preferred scales of the autocorrelation function, c B (L). Heuristically, we expect that featured or oscillatory autocorrelation functions at large L correspond to a featured and oscillatory conversion probability at high frequencies.

IV. THE FOURIER TRANSFORM FORMALISM: THE MASSLESS CASE
The assumption of m 2 a ω 2 pl used in the preceding section simplifies the phase-factor of the transition amplitude to Φ(z ) = m 2 a z 2ω = ηz . The linear appearance of the integration variable z makes the connection to Fourier transforms evident. When the plasma frequency cannot be neglected, the analysis is more subtle. In this section, we focus on the next-simplest case of massless axions. We will see that in this case, a change of variable again makes the simple Fourier formulas applicable, although at a slightly modified form.

A. The transition amplitude and the conversion probability
In this section, we set m a = 0 and consider ω 2 pl (z) = 0. The transition amplitude is again given by but now the phase factor reads where we have defined λ = 1/ω and Since ω 2 pl is positive definite, ϕ is monotonically increasing with z , has the range 0 ≤ ϕ ≤ ∞, and is a welldefined coordinate, replacing z . The measure transforms as dz = dϕ/ω 2 pl and the full amplitude can again be expressed using Fourier cosine and sine transforms: (77) whereĜ denotes the transform of the function We emphasise that B x and ω pl in this equation are functions of ϕ, as defined by (76). Since theĜ are real, the conversion probability is simply given by which is the power spectrum of G(ϕ). Our derivation of the Wiener-Khintchine theorem leading to equation (30) now implies that where c G (ψ) denotes the autocorrelation of G in ϕ-space, i.e.
In sum, also in the massless case there are two possible routes to calculate the conversion probability: either by calculating the cosine and sine transforms of G, then taking the sum of squares, or by computing the autocorrelation function of G and taking its cosine transform. A novel feature of the massless case is that, for the purpose of axion-photon conversion, distances are more naturally measured in terms of the phase ϕ instead of the spatial coordinate z. Moreover, the Fourier cosine and sine transforms are taken of the function G, which includes a factor of 1/ω 2 pl . The appearance of this factor is mathematically intuitive: when ω 2 pl is large, the phase in equation (74) varies quickly and the integral tends to wash out. For smaller values of ω 2 pl , the oscillations are slower and it is easier to build up a non-vanishing amplitude. When changing integration variable from z to ϕ, the phase Φ grows at a constant rate, by definition. The impact of ω 2 pl on the transition amplitude is instead captured by the explicit factor of 1/ω 2 pl in G(ϕ). This factor makes the impact of the plasma frequency on the transition probability apparent. In particular, this form of the amplitude suggests that fluctuations in the electron density along the direction of propagation can be an important source of axion-photon oscillations. We expect to return to this phenomenon in future work.

B. Examples
We have seen that axion-photon transitions for massless axions proceed similarly to the case of massive axions (cf. e.g. equations (26) and (30) to equations (79) and (80)). The examples worked out for massive axions in section III B applies with small modifications to massless axions if the plasma frequency is constant. More generally however, one needs to account for the non-trivial coordinate change from z to ϕ, as we here illustrate by an example.

Example 7: Massless axion in an oscillating magnetic field
Consider an axion with m a = 0 in an environment where the electron density decreases like n e ∼ R/z. The squared plasma frequency is then given by where ω 2 pl (R) = 4πe 2 me n e (R). From equation (76), we have that For this simple form of the electron density, we can invert this expression to find z(ϕ) explicitly: where we have defined the dimensionless variableφ = ϕ/(ω 2 pl (R)R). We consider a damped, oscillatory magnetic field of the form The explicit damping factor mimics a common parametrisation for galaxy cluster magnetic fields (cf. e.g. [52]), but is here mainly motivated by convenience: since n e (z)/n e (R) = ω 2 pl (z)/ω 2 pl (R), the function G becomes wherek = kR.
The real part of the transition amplitude is given by where we have introduced η = λω 2 pl (R)R, and where C and S denote the corresponding Fresnel integrals. The imaginary part of the amplitude is given by The transition probability is, as usual, obtained from the sum of squares of equations (87) and (88). The autocorre-lation function of G(ϕ) is readily found to be whereψ 2 = 2 +ψ − 2ψ 2 . Using equation (80) and that the cosine transform is its own inverse, cf. equation (22), we note that this expression gives the spectrum of oscillations of P γx→a , as understood through the cosine transform. This spectrum would be non-trivial to obtain explicitly from directly transforming the conversion probability, and similarly, in this case it's challenging to analytically take the cosine transform of the autocorrelation function to obtain the conversion probability.

V. THE FOURIER TRANSFORM FORMALISM: THE GENERAL CASE
In general, the axion/photon trajectory may pass through dense regions with ω pl > m a as well as dilute regions with ω pl < m a . The phase Φ will then increase in some regions and decrease in others, which prohibits a one-to-one, global change of coordinate from z to Φ. Moreover, at points of stationary phase, dΦ/dz = 0, the photon and axion are mass degenerate, and can interconvert resonantly. The full amplitude is then the sum of resonant and non-resonant contributions, and in sections V A-V C we discuss how to carefully calculate these contributions in the perturbative formalism.
Our approach is to split the integration domain into sub-regions in which the phase Φ is monotonic, either increasing or decreasing, or approximately stationary. Points of stationary phase will produce a 'resonant' contribution to the amplitude that depends only on the local environment around that point. For each region in which Φ is monotonic, we can identify a suitable coordinate extension and express the contribution to the amplitude as a Fourier transform. The total amplitude is simply given by the sum over all resonant and non-resonant contributions. Figure 1 shows a particular example of the general case (there denoted 'Case II'): the plasma frequency crosses the axion mass at three points, resulting in four contributions to the non-resonant amplitude (labelled by 1 through 4) and three resonant contributions.
Our notation for the general case is as follows: we factorise the phase into the (dimensionful) 'coordinate' ϕ and its Fourier conjugate λ = 1/ω as We denote the points of stationary phase by z * i (for i = 1, . . . , N ), so that ϕ (z * i ) = 0. We also define a small interval around each stationary point as z i < z * i < z u i . As we discuss in section V A, these intervals should be chosen to be small enough so that ϕ is well-approximated by a leading-order Taylor expansion around z * i . For ease of notation, we also define z u 0 = 0 and z N +1 = z. In this notation and without loss of generality, the full amplitude is given by where iA non−res Equation (92) picks up the 'resonant' amplitudes from stationary points, while equation (93) collects the nonresonant integrals. We will now discuss these contributions in turn.
A. Resonant contributions: standard stationary phase approximation and ω/ma enhancement In this section we provide context for the general discussion on resonant conversion in section V B, and point out the subtleties relating to the use of the stationary phase approximation when evaluating resonant amplitudes, that can lead to large inaccuracies at sufficiently high energies.   The resonance condition, m a = ω pl (z * ), implies that the eigenvalues of H 0 are degenerate at the resonance point, and the splitting of the eigenvalues of the total Hamiltonian is only due to the off-diagonal interaction H I . This corresponds to an 'avoided level crossing' and under certain assumptions, transitions occurring at avoided level crossings can be solved non-perturbatively in g aγ . In our notation, if one assumes that m 2 a − ω 2 pl (z) varies linearly with z − z * , that the magnetic field is constant, and that the boundaries of integration can be taken to z → ±∞, then the conversion probability is described by the Landau-Zener formula [53,54]: where γ = ∆ 2 x (z * )/|∆ γ (z * )|. When γ 1, the conversion probability is non-perturbatively large and referred to as 'adiabatic'. In the opposite, 'non-adiabatic' regime, γ 1 and the exponential can be approximated order-by-order in γ, with the linear-order probability given by: This expression is equivalent to what one obtains by applying the stationary phase approximation to a perturbative solution of the Schrödinger-like equation, as we now show. The premise of the method of stationary phase is that integrals of rapidly oscillating functions, such as (91), tend to cancel, and so the dominant contribution is expected from where the phase is approximately constant. Suppose that the amplitude A res γx→a has a point of stationary phase at z * , so that Φ (z * ) = λϕ (z * ) = 0. In the limit of λ → ∞, the amplitude is dominated by the resonant contribution, which can be evaluated by Taylor expanding ϕ(z) to second order and treating the non-exponential part of the integrand as constant: where Using the explicit expressions for the matrix elements in (8), we can write the resonant conversion lengths as where in the last step we have defined the energyindependent conversion length L = 2π |ω pl | , and set ω pl (z * ) = m a . The conversion probability is given by with all quantities evaluated at the point of stationary phase. Clearly, P res γx→a = 2πγ ≈ P LZ γx→a , so that the probability obtained from the perturbative amplitude identical to the Landau-Zener probability in the non-adiabatic limit, cf. equation (95). This is unsurprising, since these expressions are derived under the same assumptions, although taken in different order.
An important aspect of equation (99) is that the conversion length is energy-dependent, and can grow large for highly relativistic axions. As we will now discuss, this signals the breakdown of the stationary phase method at sufficiently high energies, and makes equation (99) inaccurate. In this section, we propose an approximate, simple formula that regulates the probability at high energies, and in section V B, we instead carefully re-derive the amplitude and find an analytic formula that is accurate, but more complicated.
The factor of ω/m a in equation (99) comes from the square-root dependence of the resonance length on the energy, or equivalently the 1/ √ λ in equation (97). The stationary phase approximation is expected to be accurate for λ → ∞, but in the relativistic case, we are interested in small m a λ, where equation (96) is no longer guaranteed to be applicable. Indeed, the first step of equation (96) approximates the non-exponential part of the integrand, in our case ∆ x (z), as a constant over distances large compared to L res . This is reasonable if L res is sufficiently small, but less so if the resonance length extends over astronomical distances, as is the case when ω/m a is sufficiently large. When the magnetic field varies significantly over the resonance length, the amplitude may partially cancel and the stationary phase approximation becomes inaccurate.
A better approximation of the resonant amplitude for relativistic axion-photon conversion can be obtained by replacing L res by a regulated resonance length, L reg . This corresponds to assuming that contributions to the amplitude further away than a distance of L max /2 from the resonance point cancels: where ω max /m a = L max /L. This expression involves the free parameter L max , which can be approximated by the length scale associated with variations in the nonexponential part of the integrand, in our case: L max = ∆ x (z * )/|∆ x (z * )|. We expect equation (100) to be accurate for ω ω max . For ω ≥ ω max , using L reg instead of L res in equation (99) for the conversion probability should lead to an improved estimate, though a more detailed calculation is required for accuracy. Moreover, we expect this simple regularisation to be more accurate than simply discarding contributions where L res becomes large compared to other physical scales of the problem, as done in e.g. [55,56].
Two conceptual points emerge from this analysis. First, with all else the same, resonant axion-photon conversion of relativistic axions is more relevant in environments with slowly varying magnetic fields. That is, the resonance length is not only dependent on the plasma frequency as in the standard analysis, but also the magnetic field properties. Second, the spectral shape of the resonant amplitude is in all cases very simple: it increases with √ ω up to a maximal energy that is determined by the spatial variation of the magnetic field, after which it is expected to plateau. To contextualise these results, we close this section with a brief review of previous work on relativistic axionphoton conversion in astrophysical environments, focussing on the extent to which resonant conversion has been accounted for.
Reference [57] discussed resonant production of axions in cosmological magnetic fields, considering primarily the Cosmic Microwave Background (CMB) as the source. In the adiabatic limit considered in [57], the full expression of (94) need to be used to find the conversion probability, so that the linear dependence of L 2 res on the mode energy (cf. equation (99)) appears in the exponential. The issue of the validity of the Landau-Zener approximation when accounting for the spatial variation of primordial, cosmic magnetic fields was not discussed in [57]. This work was later updated in [55], which used the stringent limits on CMB spectral distortions from COBE/FIRAS to constrain resonant axion-photon conversion, now pressed into the non-adiabatic limit. The enhancement of the resonance length with ω/m a appears in [55], but contributions where L res was larger than the estimated coherence scale of the magnetic field were simply discarded.
At gamma-ray energies, reference [58] provided an overview of astrophysical environments that could lead to substantial conversion probabilities, considering in particular adiabatic, resonant mixing. More recently, reference [59] discussed the mixing of relativistic axions from magnetars with hard X-rays, and considered also the resonant case. These references did not explicitly discuss the dependence of L res on the mode energy.
Axion-photon conversion at X-ray energies in galaxy clusters can be highly efficient, and has been studied by several groups at X-ray and gamma-ray energies. For axion masses that fall in the range between the minimal and maximal plasma frequencies in the cluster, the conversion amplitude will include one or more resonant contributions, in addition to non-resonant contributions. The relevant mass range has been studied in Xray searches for axions [13][14][15][16][17][18][19]36] (but not gamma-ray searches, which have focused on the higher-mass region where non-resonant conversion can induce large irregularities). However, many studies have neglected the resonant contribution, either explicitly or implicitly. In part, this is an immediate consequence of the reliance on cell models, or comparatively coarsely sampled GRF models, for which the mixing equations have been solved numerically. In such models, the resonance condition is typically not satisfied in any patch, and the resonant contribution is effectively omitted. If by chance the resonance condition is nearly satisfied in a patch, but the cell length is large compared to L res , the resonant conversion probability is overestimated, cf. [36] for a recent discussion. This issue has been briefly discussed in [16] using the singledomain formula, where it was argued that unless the coherence length of the magnetic field grows very large near the resonance point, the resonant contribution is negligible. Our discussion below will make this argument much more precise, but will confirm the general conclusion.
Finally, a related issue appears when considering the resonant conversion of non-relativistic axions and photons. The mixing is still governed by the Schrödinger-like equation, but the resonance length is now enhanced by a factor of ∼ 1/ √ v a . The perturbative amplitude can be evaluated using the standard stationary phase approximation, cf. [56,60,61] in the case of dark matter conversion in neutron star magnetospheres. However, the standard stationary phase approximation again suffers analogous theoretical corrections, as recently discussed in [47].

B. Resonant contributions: general treatment
The standard treatment of the resonant amplitude assumes that it suffices to make a quadratic Taylor expansion of the phase, while extending the limits of integration to ±∞. This procedure is motivated by the assumption that contributions from regions far away from the resonance point cancel, and erases information about nonresonant contributions. However, non-resonant axionphoton conversion can be non-negligible or even dominate resonant contributions, and a careful calculation requires both types of contributions to be accounted for. We here show how this can be done, and that the width of the finite region attributed to the resonant contribution is arbitrary if taken sufficiently small, and the non-resonant contributions are always important for an accurate result.
To evaluate the resonant amplitude (92), we treat the non-oscillatory part of the integrand as slowly varying (or constant), and Taylor expand the phase λϕ(z) to second order around z . The integral runs from z = z * − ∆z to z u = z * + ∆z. However, to make the build-up of the resonant contribution more apparent, we replace the upper limit by an intermediate z satisfying z < z ≤ z u . We then have that The error function grows linearly for small argument, and performs damped oscillations as the argument increases, and finally asymptotes to 1, consistently with equation (96). So far, we have kept ∆z arbitrary, however, we expect that accurate results are obtained when ∆z/L 1 so that ∆z/L res 1 for all relevant energies, and the Taylor expansion of the phase can be trusted. In this limit, the resonant amplitude simplifies further: Note that ∆z is assumed small but is still arbitrary, and varying ∆z shuffles around contributions to the amplitude between the resonant and non-resonant parts. This implies that the non-resonant contributions will generically always be important: when insisting on taking ∆z small enough for the Taylor expansion of the phase to be a good approximation, the resonant contributions never dominate the non-resonant contributions.

C. Non-resonant contribution
We now turn to the non-resonant amplitudes of equation (93). The terms are of the form and ϕ is either monotonically increasing or decreasing in the entire interval [z u i−1 , z i ]. To interpret this integral as a Fourier transform, we would like change integration variable from z to ϕ, and extend the range of ϕ to ±∞.
A suitable extension, ϕ i (z ), of ϕ can be defined by having ϕ i coincide with ϕ for z ∈ [z u i−1 , z i ], and extend linearly outside this range: With this definition, ϕ i is continuous and has a continuous first derivative, and an unbounded range for z ∈ (−∞, ∞). Clearly, this is a one-to-one map from z to ϕ i , which provides the change of coordinates that we seek. We furthermore denote the even extension of ∆ x (z ) by ∆ e B (z ), defined for z ∈ (−∞, ∞). Equation (103) can now be written as Using dϕ i = ϕ i dz , this may be written as where for a monotonically increasing phase we have ϕ min = ϕ(z u i ) and ϕ max = ϕ(z i+1 ), and for a decreasing phase we have ϕ min = ϕ(z i+1 ) and ϕ max = ϕ(z u i ).
We've slightly abused the notation to write ∆ e B (z ) = ∆ e B (z (ϕ i )) = ∆ B,i (ϕ i ). Equation (107) shows that I i (λ) is, up to a factor of 2π, an inverse Fourier transform of the function G i (ϕ i ), defined as (108) Given that ϕ i is now merely a coordinate, we can relabel it by φ, dropping the i, and write the total amplitude from (93) concisely as where In equation (109), the circumflex denotes the ordinary, complex Fourier transform, for which we use the convention: We note that the transform of equation (109) is welldefined for any real λ, but only positive values have the physical interpretation as inverse energies: λ = 1/ω. Furthermore, this equation makes it evident that the nonresonant amplitude in general can be understood as (an inverse) Fourier transform of G, when expressed as a function of the wavelength. Clearly, questions about the oscillations in the amplitude, and more specifically its spectrum, map directly to questions of the real-space properties of G. The spectrum of the amplitude is given by its Fourier transform: (113) The coordinate ψ has dual interpretations: it is the frequency of probability-oscillations in wavelength-space, and it's the phase corresponding to a set of positions along the trajectory. The map to position space is oneto-many and given by φ → ϕ i → z . Thus, the Fourier transformed amplitude reads off, and sums, the values of ∆ x /|ϕ | at up to N + 1 locations along the line of sight. Clearly, if ψ falls outside the range of all window functions in G, there is no transition amplitude.
The total conversion probability from the non-resonant contributions is given by the power spectrum of G P non−res γx→a (λ) = Ĝ (−λ) (114) where in the last line we've defined the individual probabilities from each regions, P i = |Ĝ i | 2 , and the interference terms, P ij = Re(Ĝ iĜ * j ). The Wiener-Khintchine theorem can again be applied to express the conversion probability as the Fourier transform of the auto-correlation function of G: where Note that c G is an even function of its argument.

D. Examples
Fully solvable analytical examples are more sparse in the general case than in the simpler cases of very massive or massless axions. Here, we provide one example in which both the resonant and non-resonant contributions can be calculated exactly (to this order in perturbation theory), and which illustrates the application of the Fourier formalism in the general case.

Example 8: Resonant and non-resonant contributions
We consider a model which is simple enough that some of the approximations of the stationary phase method actually hold exactly, and the resulting amplitude can be calculated in three ways: as a purely resonant contribution; as an arbitrary mix of resonant and non-resonant contributions; or as a purely non-resonant contribution. We assume the following plasma frequency profile leading to a resonance m a = ω pl for u = z/z * = 1. The resonance point is assumed to lie within a region with a constant magnetic field of magnitude B 0 , which extends up to u max > 1. The (dimensionful) phase coordinate is then given by . Clearly, this phase is a quadratic function of the spatial coordinate around the resonance point.
We first evaluate the amplitude by using our version of the stationary phase approximation, as derived in section V B. The resonance length is given by where we in the last term have introduced the full phase at the resonance point: Φ * = −λϕ * . The resonant contribution from some small region [z * − ∆z, z * + ∆z] is now (exactly) given by In this simple example where the magnetic field is constant and the phase quadratic in z, we can also calculate the full amplitude (i.e. over [0, z max ]) as a resonant amplitude by using equation (101): In this example, the standard treatment of the stationary phase approximation, which neglects the Error functions and uses equation (96), is a good approximation when Φ * 1, but not otherwise.
The transition amplitude can also be computed using the Fourier techniques developed in section V C. The first step is then to change coordinate from z to ϕ. Inverting equation (118), gives two branches: where we have introduced the notationφ = ϕ/|ϕ * |. Following the procedure outlined in section V C, we divide the integration domain in u into three regions: [0, 1 − δ], [1 − δ, 1 + δ] and [1 + δ, u max ]. The first and last of these can be evaluated with Fourier techniques as non-resonant amplitudes, and the middle contribution is a resonant amplitude that evaluates as in equation (120).
In the first region ϕ drops from 0 to −|ϕ * | + δϕ (while u = u − ), and in the third region the phase increases from −|ϕ * | + δϕ to some ϕ fin = ϕ(u max ) (while u = u + ). To account for the change of measure, we need dϕ/dz: The function f i (ϕ i ) = ∆ x /|ϕ i | are in this simple case identical between the two regions: The full function G(φ) is then given by (125) The non-resonant amplitude is now given by the inverse Fourier transform of G, and the result can be expressed using Fresenel's C and S integrals, incomplete Euler gamma functions, or, as we do here, error functions: To get the full amplitude, we add the resonant contribution from (120) to this equation, which simply cancels the last term of equation (126), giving a total amplitude that agree with the resonant-only calculation of equation (121). Finally, by taking ∆z → 0 so that there is no resonant contribution, the last term of equation (126) goes to zero, and the non-resonant amplitude reproduces the full result of equation (120). We conclude that all three ways of calculating the amplitude are in exact agreement, and the conversion probability P γ→a = |A non−res γ→a + A res γ→a | 2 , is independent of the split between resonant and nonresonant.
The conversion probability can also be calculated directly from the autocorrelation function of the function G(φ), taking ∆z = 0. The autocorrelation function is given by for 0 ≤ψ ≤ 1, and for 1 <ψ ≤ 2. Since the autocorrelation function is an even function defined over the entire real line, this completely defines c G . In figure 7, we show the function G(φ), the real and imaginary parts of the transition amplitude, the autocorrelation function c G , as well as the resulting conversion probability. A notable feature is the discontinuity in G(φ) at φ = 0 where G goes from having support from two functions (G 1 + G 2 ) to just a single non-vanishing function (G 2 ). This is also reflected as a kink in the autocorrelation function of G at ψ = ±1. In figure 8, we show how the conversion probability builds up with distance across the resonance point.

VI. NUMERICAL TESTS OF THE FORMALISM
Expressing axion-photon mixing as Fourier transforms makes it possible to leverage the highly effective numerical techniques of the Fast Fourier Transform (FFT) to determine the predictions. In this section, we show how FFT techniques can be applied to complex examples with massive (m a ω p ) and massless (m a ω p ) axions. 9 We use these numerical examples to test the applicability of the Fourier formalism, and compare it to the traditional method of solving the Schrödinger-like equation directly. We find that the Fourier-based techniques are much faster than direct simulations (as expected), and that, in our examples, the perturbative expansion holds very well over a wide range of interesting energies and for   couplings up to about an order of magnitude larger than the current observational limit. This suggests that numerical Fourier methods are superior to traditional methods when searching for axions using high quality data.

A. The discrete cosine and sine transforms
In sections III-V, we presented a series of examples where the relevant Fourier transforms can be computed analytically. In more complex models for the magnetic field and the plasma density, this may not be possible, and we can instead make use of the discrete cosine and discrete sine transforms (DCT and DST, respectively) [62]. The DCT and DST can be computed quickly using a modified FFT, and implementations of these algorithms are available as standard in numerical libraries such as scipy [63,64] or the GNU scientific library [65]. We here focus on the DCT, but equivalent equations can be written for the DST. There are four different types of the DCT [66], and here we use the 'type-I' definition which simply corresponds to discretising the cosine transform over a finite range, as we now discuss.
The DCT can be obtained by replacing the continuous coordinate (z or ϕ) by an evenly sampled list of points: z n = n ∆z for n = 0, . . . , N , so that z max = N ∆z. The conjugate variable (η or λ) is also discretised with η m = m π/z max where m = 0, . . . , N . Note that the product z n η m = n m π/N . The DCT of type I is given by: which can be understood as the midpoint Riemann sum (divided by ∆z) of the continuous cosine transform. Clearly, the resolution of the conjugate variable is ∆η = π/z max and the maximal value is η N = π/∆z. In the context of axion-photon oscillations, obtaining sufficient energy resolution and coverage of the conversion probability, requires sampling with high enough spatial resolution and out to sufficiently large values of z max . From the examples we have considered, a modestly large value of N = 10 5 gives extremely good resolution in the X-ray band. The value z max is naturally no smaller than the region over which the magnetic field is non-vanishing, but it can be made arbitrarily large by 'zero-padding', i.e. appending the set of f n with zero-valued data points.

B. Numerical implementations applied to the Perseus cluster
For our numerical tests, we consider classes of models that have previously been used to study photon-axion conversion in astrophysical settings, adopting magnetic field profiles and density laws appropriate to the Perseus cluster. We first consider two examples with massive axions (m a ω pl ): a general cell model, and a turbulent, divergence-free model based on Gaussian random fields. As shown in section III, any such massive model can be solved analytically to leading order in perturbation theory, so here we compare the results of the discrete Fourier transform with direct numerical simulations of the Schrödinger-like equation. We then consider an example with a massless axion and a 'beta-model' plasma density, which illustrates how the discrete Fourier transform can be applied to the phase variable ϕ.
In all three examples in this section, we adopt g aγ = 10 −13 GeV −1 and use a domain size of 500 kpc. The DCT is calculated using 2 × 10 4 and 10 5 samples in the massive and massless case, respectively. We will focus on the gamma-ray regime in the massive ALP case and the X-rays in the massless ALP case, but we stress that the formalism can be applied to a much broader class of problems, and over any energy regime.

Cell model, massive axions
We first consider massive axions (m a ω pl ) in a general cell model, cf. section III B 2 for the analytical solution. These types of models are commonly used to model photon-axion conversion in complex astrophysical environments, such as galaxy clusters. Specifically, we consider a realisation of the stochastic 'Model B' of [19], where cells of random size L are generated according to a power-law probability distribution p(L) ∝ L −1.2 between 3.5 and 10 kpc, and the magnetic field direction in each cell is random and isotropic. Such magnetic domains are assumed to extend to 500 kpc from the centre. The total magnetic field strength of each cell is assumed to decline with radius as where the function h(z) is given by a double-β law of the form h(z) = 3.9 × 10 −2 (131) In [19], h(z) corresponds to the electron density, a form originally given by [67], but, in this test, we explicitly set n e = 0 everywhere so as to ensure the massive axion case. This cell-based model is a slightly altered version of the model used by [68], and is based on very long baseline array (VLBA) observations of NGC 1275, the central AGN in the Perseus cluster.
The magnetic field profile is shown in figure 9, together with the associated DCT and DST (giving the real and imaginary parts of the ampllitude), the magnetic autocorrelation function, and the resulting conversion probability. The probability curve obtained from squaring the amplitude or taking the DCT of the magnetic autocorrelation function agrees well with the numerical solution from the Schrödinger-like equation, but is significantly quicker to compute.

Gaussian random field, massive axions
Next, we consider the mixing of photons and massive axions in a divergence-free Gaussian random field (GRF) model. The magnetic field profile is generated using a variation of the method of Tribble [69], which has been used in various guises in the astrophysical and axion literature [39,70,71]. Briefly, the field is generated by drawing random Fourier components of the vector potential, from which the real-space A is obtained through the inverse Fourier transform. The vector potential is then rescaled to account for the decrease of the field strength with radius: specifically we use the scaling of 'Model B' from [19], as in the previous subsection. Finally, a divergence-free field magnetic field is generated from the curl of the rescaled vector potential. To  generate the field, we have to specify the power spectrum; we assume a Kolmogorov power spectrum such that E(k) ∝ k −5/3 , where E(k) is the energy contained in the interval (k, k + dk) with wavenumber k. We use minimum and maximum scale lengths 3.5 kpc and 30 kpc, respectively, and the field is sampled at 1.75 kpc intervals for the solution of the Schrödinger-like equation. As in the previous section, we set n e = 0 as we are considering the massive axion regime.
In figure 10, we show B x of a specific realisation of the magnetic field, the real and imaginary parts of the amplitude (obtained from the DCT and DST), the magnetic autocorrelation function, and the resulting conversion probability. The GRF has significantly more smallscale structure in the magnetic field than the previously studied cell-model, due to the presence of a large range of Fourier components. This leads to finer structure in the DCT/DST, and the conversion probability. The conversion probability curve calculated from the DCT of c Bx agrees very well with the result from the numerical solution of the equation of motion.

Gaussian random field with β-law density, massless axions
We now consider the case of massless axions. A nontrivial aspect of the evaluation of the amplitude is the change of coordinates z → ϕ, which depends on the plasma density. In this example, we consider the famous class of electron densities described by 'β-models': Here n 0 denotes the electron density at the origin, and r c and β are positive constants. β-models and their variants  are ubiquitously used as simple approximations of the gas density in galaxy clusters [72,73]. We consider a source of photons located at the centre of the cluster, and calculate ϕ from equation (76) as: where we have made a change of integration variable from z to u = z/r c . The integral evaluates to the 'ordinary', or Gaussian, hypergeometric function: The inverse function, z(ϕ), is readily obtained numerically from this expression, cf. figure 11. Four our specific example, we take β = 1, n 0 = 0.05 cm −3 , r c = 125 kpc and a total domain size of 4r c = 500kpc, a choice which gives densities and scale lengths comparable to those in typical cool-core clusters.
We assume that the magnetic field is the GRF-example of section VI B 2. The function G(ϕ) is plotted in figure  12a, and should be compared to the magnetic field of figure 10a: in contrast to the magnetic field that decays with radius, G(ϕ) keeps a rather constant magnitude of oscillations over a large potion of its range, and then increases for large values of ϕ. This reflects how the transition probability is more sensitive to the magnetic field in (real-space) regions where the phase Φ varies more slowly, i.e. where the plasma density is suppressed.
The real and imaginary parts of the transition amplitude are plotted in figure 12b, and the magnetic autocorrelation function, c G (ψ), is shown in Fig 12c. The conversion probability is shown in Fig 12d, again compared to a numerical solution of the Schrödinger-like equation. Once more, the agreement is extremely good.

C. Performance Improvement
To compare the numerical efficiency of the traditional Schrödinger approach to our new, numerical Fourier approach, we consider the following problem. Suppose that we would like to calculate the conversion probability for N m axion masses and N g axion-photon couplings in the perturbative regime. We will for simplicity consider massive axions in this example, and assume that the background magnetic field and autocorrelation function are known. Suppose further that to obtain sufficient energy resolution over the specified range, the calculation should include at least N ω beam energies, and that axion-photon trajectory needs to be sampled with at least N z domains for numerical accuracy. We first consider the computational time required in the traditional approach, before comparing with that of the Fourier approach.
To solve this problem in the traditional approach requires calculating the z-evolution operator through matrix multiplication of the N z operators for the uniform domains. This has to be done for each mass, energy and mode energy, so if the calculation time for a single uniform domain is ∆t 1 , then the total time to solve this problem is given by By contrast, in the perturbative approach, g aγ is an overall prefactor which trivially rescales the probability (at negligible computational cost), so the discrete Fourier transform can be performed at a fixed coupling. Similarly the mass only appears inside η = m 2 a /(2ω), so that the probability can be calculated for a fixed mass, and then re-interpreted for other masses (again at negligible numerical cost). Moreover, a single discrete Fourier transform samples all spatial points and generates the probability at all mode energies. The number of points that the discrete Fourier transform runs over, N , must then satisfy both N ≥ N z and N ≥ N ω . If N ≥ N ω > N z , taking the Fourier transform requires 'zero-padding' as described above. If N ≥ N z > N ω , the Fourier transform tends to give a higher energy resolution than the minimal requirement. Given the usual scaling in the DCT algorithm of O(N log 2 N ) [63], we expect the total calculation time to scale as where ∆t 2 is now a single multiplication in the DCT algorithm, and ∆t 2 ∆t 1 . The fractional numerical gain can in this example be estimated as This performance improvement is model-dependent but appreciable for all relevant values of the parameters, and tends to be maximised when N z ≈ N ω , so that the Fourier approach does not lead to 'superfluously' high spatial or energy sampling. We now consider a concrete example where we numerically compare the actual computational time. We take N m = N g = 1, and set g aγ = 10 −13 GeV −1 and m a = 5×10 −9 eV. For the magnetic field, we consider the model of GRF model described in section VI B 2. To illustrate the dependence on N z , we 'truncate' the magnetic field at the maximal radius N z ∆z with ∆z = 1.75 kpc and vary N z from 1 to 501 at intervals of 20. Thus, the parameter N z in this case parametrises the size of the magnetised region. To ensure sufficient energy resolution in the band 1 − 100 GeV, we take N ω = 1104. This requires the number of discrete Fourier transform samples to be N = 4096 for the chosen value of the mass. We note that this example is rather conservative since N ω N z , and the Fourier approach samples the spatial profile more densely than the minimal requirement.
The speed-up factor is shown in Fig 13 as a function of N z . The speed improvement is linear in N z , because ∆t sle ∝ N z and ∆t dct is constant for constant N and N ω . The speed-up factor quickly exceeds 100 even for modest numbers of cells, and this particular test provides quite a conservative estimate of the performance improvement. For broader energy ranges, the speed-up is even more significant since more solutions to the Schrödinger-like equation must be calculated (the factor N ω /N increases). The speed-up is further increased by several orders of magnitude when N m and N g are taken to be sufficiently large to ensure a dense sampling of the axion parameter space.
Overall, our results show that the Fourier formalism is significantly quicker -often by orders of magnitudefor models with more than a few cells, and therefore also for the types of magnetic field models commonly used for astrophysical ALP searches.

D. Perturbativity
A central aspect of our analysis in this paper is perturbation theory, where we have consider the leading-order contributions to the conversion probability (i.e. A ∼ g aγ and P γ→a ∼ g 2 aγ ). As discussed in section II, parametrically, the perturbative expansion should be good when conversion probabilities are small or moderately small, which includes a large fraction of all cases of phenomenological interest. In this section, we elaborate on this general argument through two examples: the simplest case of a constant magnetic field, and a complex example with a turbulent magnetic field. In both cases, we seek to identify the breakdown of the perturbative expansion when the mixing grows large.
First, we consider the single domain of uniform B and length z max discussed in section III B 1. This is one of a very small set of examples in which both the full and the perturbative conversion probabilities can be calculated analytically. The ratio between the leading-order and full conversion probabilities is: The perturbative expansion in g aγ B can be interpreted as an approximation using 4∆ 2 x /η 2 1. Clearly, when ηz max is sufficiently small to allow for a Taylor expansion of the sines, An interesting special case is when 4∆ 2 x /η 2 1 but ηz max is sufficiently large that ∆ 2 x z max /η π. In this case, the oscillations of the conversion probability are very fast, and the LO approximation is out-of-phase with the exact result. This discrepancy is not always observable however as detector resolution or Doppler broadening limits the detectability of high frequency oscillations, and the cycle averaged predictions at LO agree with that of the full model. In figure 14a we show a comparison of the first-order, Fourier-like formula, which can be calculated from the cosine and sine transforms, to the full formula accurate to all orders. The calculation is conducted for a beam energy of 10 keV, with m a = 10 −11 eV, z max = 10 kpc and B = 10 µG. We show the conversion probabilities as a function of g aγ B x z max , together with the absolute fractional error. The relative error is at the few % level or lower until the conversion probability exceeds ∼ 0.01; beyond this point errors can be significant and approach unity or higher -moreover the probability exceeds 1 at very large g aγ .
Second, we consider a more complex cell model with many domains, taking the example show in figure 9a for concreteness (but now calculated using m a = 10 −11 eV). Here, we compare the result obtained from a numerical solution of the Schrödinger-like equation to a numerical DCT of the autocorrelation function of the field, c Bx (L). The conversion probabilities are shown as a function of energy for a few different values of gradually increasing g aγ , in figure 14b. Overall, the outcome is similar to the single-cell case, with good agreement until the amplitude of the conversion probability exceeds a few %, where second-order effects start to imprint themselves on the curve. In fact, the values of g aγ above which the first-order predictions from the DCT deviate substantially from the more detailed calculation are already ruled out in X-ray searches using similar models [19]. The tests presented here are specific and not exhaustive, so caution should be exercised when applying our scheme to any situation whether the conversion probability exceeds a few percent.

A. Magnetic Fields in Clusters
The formalism developed in this paper provides new insights into the dependence of axion-photon conversion probabilities on the structure of the magnetic field along the particle trajectory. Given the special role that galaxy clusters play in these studies, it is appropriate to review our current understanding of the structure and strength of the magnetic field in the ICM.
To state the current paradigm up-front -multiple lines of evidence suggest that the ICM is a turbulent magnetized plasma that is, almost everywhere, in approximate hydrostatic equilibrium in the gravitational potential of the cluster's dark matter halo. In the bulk of the ICM, the magnetic pressure appears to be approximately 1% of the thermal pressure and the magnetic field is tangled on scales of a few-kpc, consistent with characteristic spatial scales in the turbulent power spectrum. There is little evidence for large-scale regular magnetic field structures in the central regions of clusters. In the rest of this subsection, we summarize in brief the key evidence underpinning this paradigm.
Direct observational evidence for magnetic fields in the ICM comes principally from radio observations [74]. Faraday rotation, i.e. the characteristic rotation of the polarization angle as a function of wavelength due to propagation through a magnetized plasma, is clearly seen in the radio band towards the radio-emitting lobes of AGN embedded in the ICM. The corresponding Rotation Measure (RM) provides a direct measure of the line-ofsight magnetic field B z integrated along the path from the source to observer. For a given RM, the inferred typical magnetic field strength depends upon the number of field reversals along the path and hence the field coherence length. In a number of clusters, the RM mapping across large radio lobes reveals the typical coherence length of the field to be a few kpc [75][76][77]. Adopting these in-plane coherence lengths as being typical of the line-of-sight structure allows typical field strengths to be estimated; we find B ∼ 1 − 10µG in the cores of coolcore clusters, corresponding to a ratio of the thermal-tomagnetic pressures of ∼ 100.
Independent support for this picture comes from observations of radio halos and minihalos in the cluster cores and radio relics in cluster outskirts [78]. The radio emission from these diffuse structures is synchrotron radiation from a population of relativistic electrons gyrating in the magnetic field of the ICM (with the relativistic population likely arising from a central AGN for the central minihalos, and shock acceleration for the ra-  (35), accurate to leading order, are compared to the fully accurate analytic formula, as a function of the dimensionless quantity gaγBxzmax. dio relics). Measurements (or upper limits) on the X-ray inverse Compton scattering of the relativistic electron population can be combined with measurements of the synchrotron emissivity to estimate (or set a lower limit on) the magnetic field strength. An example of such a study for the radio halo in the Perseus cluster is provided by [79]; note that the reported detection of the inverse Compton X-rays in this work was later discovered to be due to an error in the Chandra mirror calibration. Thus the reported measurements of the magnetic field strength in [79] should be taken as strict lower bounds, implying B > 3µG in the core of Perseus. The polarization of these structures provides information about the configuration of the magnetic field. Synchrotron radiation in a highly ordered magnetic field can provide extremely high levels of polarization (with polarization fractions exceeding 50%). By contrast, the observed polarization fractions of these structures is observed to be significantly lower (< 10%), especially for the radio minihalos in the cores of relaxed clusters [80]. This is interpreted as the effect of magnetic field tangling on a scale smaller than the resolution of the radio observations, i.e. beam depolarization.
The kpc-scale tangling of the magnetic field implied by the radio observations is intimately linked to turbulence in the ICM plasma. Again, there are multiple lines of evidence suggesting that the ICM in even relaxed clusters is turbulent with typical velocity fluctuations δv ∼ 0.1 − 0.2c s , where c s is the sound speed. The density fluctuations associated with turbulence directly translates into fluctuations of the X-ray surface brightness, and these have become a powerful tool for constraining ICM turbulence. For example, [81] employ the particularly high-quality Chandra X-ray Observatory data for the Virgo and Perseus clusters, showing that both have a turbulent ICM with a Kolmogorov-like spectrum of fluctuations extending to scales well below 10kpc. The most direct detection of turbulence in the ICM of a cluster came from measurements of the Doppler broadening of X-ray emission lines by the Soft X-ray Spectrometer (SXS) on the Hitomi observatory [82] 10 . The (1-d) velocity dispersion across most of the core was found to be 100-150 km s −1 [83], to be compared with the sound speed of c s ≈ 1000 km s −1 and in good agreement with the findings of the surface brightness fluctuation analysis. This suggests that the turbulent energy density is a few percent of the thermal energy density, and so very similar to the magnetic energy density.
From a fluid dynamics perspective, the existence of MHD turbulence in the ICM is entirely expected. While the microphysics of viscosity in ICM-like plasmas is still an area of active work [84], it is clear that it is significantly suppressed below the näive textbook value [85], making the ICM atmosphere a moderate-to-high Reynolds number (Re) system. Dynamics in the ICM atmosphere is driven by a variety of external processes including merging sub-clusters which induces ICM sloshing [86], the orbiting galaxies [87], and jet-activity from the central AGN that drives ICM shocks and inflates ICM cavities [88]. There are also purely internal processes within the ICM that can drive dynamics, including conduction-driven buoyancy instabilities [89] and smallscale kinetic dynamos [90]. Turbulence is a natural outcome of such driving in a high-Re atmosphere.

B. Galaxy cluster magnetic field models in axion searches
Having discussed the evidence that the ICM is a turbulent magnetized plasma, here we briefly review the specific magnetic field models that are used for axionphoton conversion in galaxy clusters. These range from simple cell-models to more sophisticated treatements of turbulent fields. We critically discuss a recent proposal to model galaxy cluster magnetic fields as regular [41] for the purpose of axion-photon conversion, and point out that this model is inconsistent with observations.
• General cell models. The most commonly used class of magnetic field models for axion searches in galaxy clusters is the general cell model, which we discussed (and solved) in section III B 2 for massive axions. An attractive feature of this class of models is its close link to the single-domain model (cf. section III B 1), which is commonly used to build intuition around axion-photon mixing also in complex environments. In a general cell model, the sizes of the domains can be generated stochastically as in section VI B 1, which allows mimicking some properties of smooth, 'multi-scale', turbulent magnetic fields. Moreover, with a general cell-model, it is easy (though not very fast) to numerically solve the Schrödinger-like equation.
However, as an astrophysical model, cell models are extremely simplified, and do not satisfy some of the most basic properties of magnetic fields, such as being continuous and divergence free. Indeed, the non-vanishing divergence of cell models has been argued to provide a route to rule them out observationally through its impact on Faraday rotation measurements [77]. Furthermore, while qualitative agreement has been found between the predictions for axion-photon conversion in cell models and in more sophisticated, turbulent models [39], there is no clear dictionary between these models that would allow an unambiguous translation of assumptions and constraints. See also [31,40,91,92] for a discussion of variants of cell models.
• Gaussian random field models. The second class of magnetic field models commonly used for axion searches are the 'GRF' models that we described in section VI B 2. These are smooth and divergence free, and since the power spectrum is used when generating the vector potential, it is easy to connect these models to other theoretical and observational studies of magnetic field structure. GRF models have been used numerically for axion searches in e.g. [22,24,39,46].
The magnetic fields generated from a GRF are rather featureless and lack large-scale coherent filamentary structures arising from intermittency and cooling instabilities which may relate to cold filaments that are often observed in galaxy cluster cores. Thus, the GRF models are still too simple to fully capture the complex structure of galaxy cluster magnetic fields. Moreover, axionphoton conversion in GRF models has not been well-understood analytically (though our discussion in section III B 4 addresses this issue for massive axions). Finally, GRF models are slightly more involved to simulate than cell models.
• MHD models. Magnetohydrodynamic (MHD) dynamos transfer mechanical energy of the plasma into magnetic energy, and is believed to amplify small seed fields into galaxy cluster magnetic fields of O(µG) strength. Cosmological MHD simulations indicate that such magnetic fields are tangled, can feature turbulence-induced non-Gaussian structures, and extend over several decades in Fourier space (cf. [93] for a review). To date, there has been no study of axion-photon conversion in MHD magnetic fields. We note that our formalism will be very useful in developing an understanding of how axion-photon conversion in MHD magnetic fields differs from GRF models. In particular, coherent structures that are clearly evident in MHD simulations are completely erased when the phases of the (exponential) Fourier modes are randomised [94]. However, as we've discussed in this paper, axion-photon mixing is only sensitive to the magnetic field through the autocorrelation function, which is independent of the phases. This suggests that magnetic field models used for axion searches can be visibly different from more realistic magnetic fields, but still produce the same predictions.
Outflows from AGN at the centre of galaxy clusters provide an additional source of energy and complexity to the intracluster medium (ICM). In particular, jets extending from AGNs can inflate bubbles that rise through the ICM by buoyancy, and can reach large radii before 'pancaking' or dispersing. Such bubbles, or cavities, are visible in high-resolution images of galaxy clusters as regions with suppressed X-ray emission, indicating that the thermal ICM has been displaced from the cavity. For example, in the Perseus cluster, a number of such cavities have been identified in various directions within 100 kpc from the central galaxy, NGC 1275 [79,[95][96][97][98]. In particular, reference [79] used deep Chandra observations to identify a high-abundance ridge, or 'shell', around 93 kpc from the centre of the cluster (coincident with the edge of the radio mini-halo), and interpreted this ridge as a relic of a collapsed radio bubble.
Recently, reference [41] appealed to the evidence for X-ray cavities to propose a galaxy cluster model with slowly varying, 'regular' fields in Perseus, and then used this model to argue that the uncertainties in gamma-ray and X-ray limits on axions had been underestimated, possibly by several orders of magnitude. Here, we critically re-examine the motivation and consistency of the astrophysical model proposed in [41] (see section III B 6 for the analytic, perturbative solution of axion-photon mixing in this magnetic field model).
The central assumptions of the model of [41] are that: 1. The centre of the Perseus cluster consists of an X-ray cavity of radius 93 kpc.
2. The cavity can be modelled using the analytical solution of [51] (that we discussed in section III B 6), and the electron density is given by the deprojection analysis of [67].
However, both these assumptions are problematic: 1. The Perseus cluster is the brightest cluster in the X-ray sky. The azimuthally averaged surface brightness peaks towards the centre of the cluster [67], and so, the central 93 kpc of Perseus is certainly not an X-ray cavity. As mentioned above, observations by a number of groups have identified several localised cavity regions, including the high-abundance ridge at 93 kpc from the centre [79], but these do not encapsulate the cluster or even fall along the line of sight from earth to the centre of the cluster. Thus, the central assumption of [41] is in stark conflict with observations.

2.
A key prediction of the analytical, axisymmetric solution of [51] 11 is that the gas pressure dips inside the cavity, and attains its peak value at the centre and at the boundary. This prediction is incompatible with the radially decreasing gas pressure found from Chandra observations of Perseus in [99]. Moreover, the cavity model of [51] assumes that the electron density inside the cavity is smaller than the surrounding ICM, while [41] implicitly assumes that it is filled by the ICM when using the observational electron density from [67].
We conclude that the regular magnetic field model of [41] is incompatible with observations, which should be taken into account when examining its consequences for axion-photon mixing.
We note in closing that magnetic field modelling is currently the dominant systematic in studies of axionphoton mixing, and further theoretical, numerical and observational studies will be critical to further improving the sensitivity to axions.

C. Outlook
In sections III-V, we discussed the mathematical and conceptual consequences of the formalism that we have developed, and in section VI, we demonstrated how it allows for drastic improvements of numerical simulations of axion-photon mixing. In this short section, we briefly outline how our formalism suggests entirely new analysis methods for axion searches.
The current paradigm for studying axion-photon mixing involves specifying the magnetic field and the plasma density, which are then used to numerically solve the modified Klein-Gordon and Maxwell equations to find the conversion probability. In many astrophysical applications, this process must be repeated many times as one scans over the allowed magnetic field configurations.
The key message of this paper is that the perturbative axion-photon conversion probability is given by the magnetic power spectrum (in some cases rescaled by 1/ω 2 pl ), or equivalently, as a Fourier-type transform of the corresponding magnetic autocorrelation function. The autocorrelation function contains less information than the magnetic field, and visibly distinct magnetic fields can share the same autocorrelation function. This suggests that scanning over magnetic fields is redundant, and that more efficient methods could be developed based directly on the autocorrelation function. Here we consider three possibilities along these lines.
First, when defining an ensemble of magnetic fields (e.g. using one of the classes discussed in section VII B), one implicitly defines an ensemble of magnetic autocorrelation functions. This suggests that instead of generating magnetic field profiles along the axion/photon trajectory, one may directly generate the magnetic autocorrelation function. In simple cases, the distribution of autocorrelation functions may be known analytically. In more complex cases, machine learning techniques based on Gaussian processes are likely to be very suitable for learning the distribution. An autocorrelation function generated in real space need only to be Fourier transformed to yield the conversion probability. In some cases, it may be possible to generate the Fourier coefficients of the autocorrelation function directly, which simply corresponds to writing down the conversion probability, without numerically solving the equations of motion.
Second, whenever high quality radio observations of extended radio halos are available, it may be possible to dispense of explicit magnetic field modelling all together, and move directly from radio observations to predictions for the axion-photon conversion probability at high energies. In the case of Gaussian, isotropic magnetic turbulence, maps of RMs of an extended radio background source can be used to reconstruct the statistics of the magnetic autocorrelation tensor, M ij = δB i (x) δB j (x + r) , and [76]: Here M L , M N and M H respectively denote the longitudinal, normal and helical components of the tensor. We note however that M H is irrelevant for both RMs or axion-photon mixing. The remaining components, M N and M L , can be reconstructed from the spatial distribution of the observed RMs, and determine the statistics of axion-photon conversion along any trajectory through the Faraday screen. In this context, the conversion probabilities could be obtained as Fourier transforms of autocorrelation functions generated from the observationally inferred tensor M ij .
Third, photons 'disappearing' into axions lead to distortions of astrophysical spectra. These distortions will appear as oscillatory residuals in fits that don't take axions into account. Recently, reference [43] suggested that it may be beneficial to search for these imprints by expressing them as functions of their wavelength, and then performing a Fourier transform to hopefully isolate the oscillatory features. Our formalism suggests that this procedure is a map from the observational residuals to an approximation of the magnetic autocorrelation function (up to corrections due to noise and the fitting procedure). For unpolarised fluxes, the relevant conversion probability is the superposition of the polarised conversion probabilities, cf. equation (16), and the relevant correlation function will be the sum of two components of the autocorrelation tensor. This method opens up the possibility of using mathematical, physical and observational properties of the autocorrelation function to improve the sensitivity of axion searches.

VIII. CONCLUSIONS
We have considered the mixing of relativistic axions and photons to leading order in the coupling, g aγ , and asked what properties of the magnetic field are reflected in the observationally interesting conversion probability, P γ→a . Our conclusion is that P γ→a is given by the magnetic power spectrum, or equivalently, by the Fourier transform of the magnetic autocorrelation function. This result has notable conceptual, calculational, numerical and methodological consequences.
Conceptually, our formalism maps questions about the oscillatory spectrum of the conversion probability into questions about real-space magnetic correlations. Moreover, the mature framework of Fourier analysis and auto-correlation theory lead to new, interesting identities for axion-photon conversion.
Calculationally, we have demonstrated in a series of examples that large classes of magnetic fields models, that have previously only been studied numerically, can be solved analytically. This is particularly striking in the case of axions more massive than the plasma frequency (which is of direct relevance for gamma-ray searches for axions), where one of our examples comprises the general solution to axion-photon mixing in any magnetic field.
Numerically, the current paradigm for computing axion photon mixing is to simulate the equations of motion for each magnetic field realisation. Our formalism allows one to circumvent numerical simulation, and the techniques of the Fast Fourier Transform can be leveraged to efficiently determine the predictions of each model.
Methodologically, our formalism suggests new approaches for generating the predictions of the model, tightening the connection with radio observations of rotation measures, and for analysing observational data that may carry spectral imprints from axions.
Moreover, in the light of our new formalism, we have critically reviewed the types of magnetic field models used in axion searches in the literature. We have argued that the commonly used cell models qualitatively agree well with more sophisticated, turbulent GRF models, though with certain differences that can be explained by the zigzaggy nature of the cell-model autocorrelation function. We have also pointed out some inconsistencies in a recently proposed regular model for the magnetic field in the Perseus cluster.
The focus of this paper is the interconversion of photons and relativistic axions, e.g. as those that may be produced from high-energy astrophysical fluxes. Our formalism extends straightforwardly to non-relativistic axions, such as axion dark matter, in the case of a constant plasma density. The general problem of interconversion of dark matter axions and photons is more subtle, and not discussed in this paper. Similarly, we do not discuss the opposite regime of extremely energetic particles in strong magnetic fields, where the QED birefringence contribution cannot be neglected.
Our results draw on the quantum mechanical analogy for axion-photon mixing developed in [42], and in the context of quantum mechanical perturbation theory, it is well-known that the leading-order, infinite-time transition amplitude is given by a Fourier transform of the interaction Hamiltonian. To the best of our understanding, the equivalent Wiener-Khintchine form in which the transition probability is given by the Fourier transform of the autocorrelation of the interaction Hamiltoinian, is less explored in the literature.
Finally, our approach is based on perturbation theory and will be most useful when searching for spectral irregularities that are sufficiently small. The precise breaking point of perturbation theory depends on the magnetic field model, but our analysis indicate that our approach could consistently be used to search for axions with cou-plings of about an order of magnitude larger than the current observational limit, using currently available precision X-ray data. Thus, our first-order formalism can already be used for large portions of the parameter space considered in axion searches. Moreover, our method is even more relevant to future studies with the Athena mission [37] and the gamma-ray CTA observatory [38].