Magnetohydrodynamics predicts heavy-tailed distributions of axion-photon conversion

The interconversion of axionlike particles (ALPs) and photons in magnetised astrophysical environments provides a promising route to search for ALPs. The strongest limits to date on light ALPs use galaxy clusters as ALP-photon converters. However, such studies traditionally rely on simple models of the cluster magnetic fields, with the state-of-the-art being Gaussian random fields (GRFs). We present the first systematic study of ALP-photon conversion in more realistic, turbulent fields from dedicated magnetohydrodynamic (MHD) simulations, which we compare with GRF models. For GRFs, we analytically derive the distribution of conversion ratios at fixed energy and find that it follows an exponential law. We find that the MHD models agree with the exponential law for typical, small-amplitude mixings but exhibit distinctly heavy tails for rare and large mixings. We explain how non-Gaussian features, e.g.~coherent structures and local spikes in the MHD magnetic field, are responsible for the heavy tail. Our results suggest that limits placed on ALPs using GRFs are robust.


I. INTRODUCTION
Axion-like particles (ALPs) are ubiquitous in extensions of the Standard Model of particle physics [1][2][3][4][5][6] and provides an increasingly popular candidate for dark matter [7][8][9].ALPs and photons can interconvert in background magnetic fields through the interaction [10,11] where a is the ALP field, g aγ is the ALP-photon coupling, F µν is the electromagnetic tensor, and F µν = 1 2 ϵ µνρσ F ρσ is its dual.Laboratory searches based on this interaction include light-shining-through-the-wall experiments, 'helioscopes' that are sensitive to solar ALPs, and dark matter 'haloscopes', see [12][13][14] for reviews.Moreover, complementary astrophysical searches for ALP-photon mixing have proven exceptionally powerful and have led to some of the strongest limits on g aγ , as we now discuss.
The strength of the ALP-photon mixing grows with the magnitude of the magnetic field and the size of the magnetised region.Galaxy clusters are both magnetised (with µG field strengths) and large (spanning hundreds of kiloparsecs), and are known to be efficient ALP-photon converters: a significant fraction of high-energy photons travelling through a cluster could emerge as ALPs, and vice versa.It follows that the high-energy photon spectra of a bright Active Galactic Nucleus (AGN) located within or behind a galaxy cluster would be distorted by ALPphoton conversions, which makes it possible to use X-ray and gamma-ray satellites to search for ALPs.For example, precision X-ray spectrometry of the cluster-hosted quasar H1821+643 limits the amplitude of possible ALPinduced distortions to the ≲ 2.5% level, which implies strong constraints on ALP theories [15].A comparable precision has been reached for the AGN NGC1275 at the centre of the Perseus cluster [16], and the several similar studies have been performed in the X-ray and gamma-ray ranges  A critical step in translating the absence of spectral distortions into limits on ALPs is the modelling of the cluster magnetic field.For some clusters, Faraday rotation measure studies constrain key properties of the magnetic field, but the detailed spectral shape of the ALP-photon conversion ratio is given by the full spatial autocorrelation function of the magnetic field [50], which is not observationally accessible.The magnetic field must therefore be treated as a nuisance parameter to be marginalised over in the statistical analysis.Traditionally, rather simple models for the cluster magnetic field have been used in the literature, e.g. by taking the field to be constant within a series of cells along the direction of propagation (i.e.'cell-models'), or by using divergence-free Gaussian random fields (GRFs).Recent studies have found the predictions to be rather robust to the choice of model [40].Still, it is important to note that these simple magnetic field models differ significantly from more realistic models found by solving the dynamical magnetohydrodynamic (MHD) equations arXiv:2208.04333v2[hep-ph] 16 Nov 2023 of motion.In particular, MHD simulations of turbulent plasmas exhibit coherent magnetic structures on moderately large scales, which have no counterpart for GRF or cell-model fields.This prompts the question: can the predictions of the simpler models be trusted, given their lack of coherent structure?
This paper presents the first systematic study of ALPphoton conversion in MHD models of cluster magnetic fields (see also Ref. [51] for a similar study with ENZO code).We carefully analyse the predicted distributions of ALP-photon conversion ratios and compare them to those of similar GRF models.Importantly, we find that the typical predictions are in excellent agreement between the MHD and GRF models, but that the distributions differ significantly for rare fluctuations.While the fluctuations to large conversion ratios are exponentially suppressed in GRF models, the MHD models predict distinctly heavy tails.We establish the origin of the heavy tails to be coherent structures with large amplitude magnetic fields, which reflect the non-Gaussianity of the MHD field.By contrast, we find the helicity of the magnetic field to be unimportant for the conversion ratio.
In Sec.II we summarise the formalism of ALP-photon conversions in the perturbative regime.In Sec.III we discuss how MHD simulations model the environment encountered in galaxy clusters.In the literature, these turbulent magnetic fields are modeled with simplified GRF models, discussed in Sec.IV.In Sec.V we compare the results obtained for ALP-photon conversions in these different magentic field models and conclude.

II. ALP-PHOTON OSCILLATIONS
The linearised, classical field equations of electromagnetism and relativistic ALPs is given by a Schrödingerlike equation [11]: where z denotes the spatial coordinate along the direction of propagation.The components of the three-level 'state vector' Ψ(z) = (A x , A y , a) T consists of the transverse components of the vector potential, i.e. 'photon', and the ALP field.In the presence of a background magnetic field, B(x), the ALP-photon interaction of Eq. ( 1) induces off-diagonal elements that mix the ALP with the photons, where . We have neglected Faraday rotation and the QED birefringence effect, which are both negligible at X-ray energies.
Due to the mixing by ∆ aγi , a photon beam of energy ω that is linearly polarised along the i-direction loses a fraction of P aγi (ω) of its flux into ALPs, where P aγi denotes the 'quantum mechanical' conversion probability calculated from Eq. (2).Throughout this paper, we refer to P aγ as the conversion ratio.In the bulk of this paper, we focus on the most instructive, 'massive' regime where m a ≫ ω pl , but the general conclusions hold for arbitrary ALP masses (see App. C).This is not surprising, given that the mixing of ALPs and photons in an environment with constant plasma frequency is formally identical to the massive regime, and that the mixing for varying plasma frequencies results in well-understood, and typically mild, modifications to the calculations.
Throughout this paper, we will be interested in the small mixing regime where this problem can be solved perturbatively in g aγ [11], resulting in the leading order transition amplitude where and where, without loss of generality, the trajectory goes through the origin of the coordinate system.We discuss the conversion ratio for unpolarised fluxes in App.E. Clearly, as L → ∞, the perturbative amplitude is the one-dimensional Fourier transform of the relevant component of the magnetic field along the z-direction, evaluated at the conjugate 'momentum' η a .It follows that the conversion ratio is given by the one-dimensional power spectrum of the magnetic field, and can also be expressed as the Fourier transform of the real-space magnetic autocorrelation function [50].The Fourier representation makes it possible to efficiently evaluate the conversion ratio using celebrated numerical methods such as the Fast Fourier Transform [15,40,50].Moreover, Eq. ( 6) immediately suggests a subtle conceptual question.Extended magnetic structures rely on phase coherence, and disappear when the phases of the Fourier components are randomised. 1 The conversion ratio is only a function of the norm of Bi , so is P aγi independent of coherent structures?The answer to this question is no, as we will now explicitly demonstrate by comparing the predictions of magnetic field models from three-dimensional MHD simulations to GRF models.

III. MHD MODELS OF THE MAGNETIC FIELD
The intra-cluster medium (ICM) is a dilute, magnetised plasma with tangled structures on few-kiloparsec scales, as is evident from radio and X-ray observations.The ICM is near local hydrostatic equilibrium and is characterised by a very high magnetic Reynolds number (∼ 10 30 ) and significant MHD turbulence.The magnetic field is thought to be generated through gas motions by a dynamo process, driven by a variety of internal and external processes such as magneto-thermal instabilities [53][54][55], jet-activity from the central AGN resulting in ICM shocks and cavities [56], turbulent wakes of individual galaxies [57][58][59], and sporadic cluster mergers [60,61].
To analyse ALP-photon conversion in dynamically generated magnetic fields, we perform state-of-the-art MHD simulations using the Pencil Code [62] in a box of size L 3 = (200 kpc) 3 with 512 3 mesh points.The turbulence is assumed to be driven through some external volume forcing, which we model as random sinusoidal waves that are δ-correlated in time, i.e. the forcing function changes at each time step.The wave vectors are taken from a shell of finite thickness and radius k f , which we chose to be close to the smallest wave number of the computational domain k 1 ≡ 2π/L.This results in turbulence at a moderate magnetic Reynolds number that is as large as possible for the given numerical resolution.The ratio of viscosity to magnetic diffusivity is 20.We ignore the density stratification and just consider an isothermal gas with constant sound speed.The simulated magnetic fields do not decrease with radius, as is expected in galaxy clusters, and should not be interpreted as fully realistic models.Rather, the MHD models exhibit turbulence and structure, as are expected in real galaxy clusters, and allow us to test the robustness of the ALP theory predictions.
We initialise the simulations with a weak seed magnetic field.After about 50 turnover times (u rms k f t = 10 where u rms is the rms velocity), the magnetic field begins to grow exponentially.During this phase, the magnetic field is highly non-Gaussian, but the field strength is still weak.To assess the consequences of such a highly non-Gaussian field, we consider a scaled version of this magnetic field, referred to as Run K, because the dynamo is kinematic, i.e. unaffected by magnetic feedback.
When the magnetic energy density reaches values comparable to the kinetic energy density, the Lorentz force begins to affect the turbulence and leads to a saturation of the dynamo.We refer to this state as Run S. The magnetic field is then still non-Gaussian, but the kurtosis is smaller than during the kinematic stage.The density also becomes more strongly affected by the magnetic field.In both runs, the turbulence is isotropic on large length scales to a good approximation (see App. F).
To study the properties of ALP-photon conversion in the MHD magnetic fields, we consider the ensemble of trajectories defined by straight lines in the z-direction through the simulation volume.The corresponding ensemble of conversion ratios is then defined by solving Eq. ( 2) along each trajectory for a range of energies.In practice, this procedure is drastically simplified by the perturbative formalism, as the relevant predictions can be efficiently extracted from the one-dimensional power spectrum of the magnetic field, cf.Eq. ( 6).This way, Runs K and S result in two distinct statistical distributions of conversion ratios, which we now compare with the predictions of GRF models.

IV. GRF MODELS OF THE MAGNETIC FIELD
GRFs provide a convenient mathematical framework for tailoring smooth magnetic fields with an arbitrary power spectrum.GRF models have been used to study the properties of cluster magnetic fields based on Faraday rotations [63,64], and have also been used in ALP searches in the X-ray and gamma-ray energy ranges [15,28,30,40,42,65,66].In this section, we analytically compute the distribution of ALP-photon conversion ratios in GRF magnetic field models.
We consider photons that are linearly polarised along the i-direction and propagate along an ensemble of rays in the z-direction.The perturbative conversion ratio is given by Eq. ( 6), and our goal in this section is to determine the statistical distribution, f P (ηa) , of P aγi at fixed η a .
In order to isolate the intrinsic differences between models from MHD and GRF, we compare magnetic field models that have the same three-dimensional power spectrum.We denote the three-dimensional Fourier transform of the magnetic field as where the index a runs over all three spatial coordinates.For a magnetic field that is statistically homogeneous and isotropic, as our MHD and GRF magnetic field models, the two-point correlation function is given by [63] where M N and H respectively are the normal and helical autocorrelation functions, following the conventions of [67].The three-dimensional power spectrum is given by the trace of the autocorrelation tensor: The autocorrelation of Bi is then given by where The delta-function prefactor in Eq. ( 9) arises in the L → ∞ limit of the expressions sin For finite but large L and η = η ′ , the factor δ(η a − η ′ a ) is regulated to L/(2π).Clearly, helicity is unimportant for the conversion ratio of linearly polarised photons (see App. E).
The one-dimensional power spectrum P 1D fully determines the statistical properties of the Gaussian magnetic field, and thus also those of the conversion ratio.Following the derivation in, e.g.Ref. [39], we determine the probability density function (PDF), f P (ηa) , of the random variable P aγ (η a ) in the ensemble defined by the Gaussian magnetic fields (see App.A for GRFs and App.B for the non-Gaussian case).The resulting PDF takes a very simple, exponential form: This equation provides an ideal starting point to compare the predictions of GRFs to those of MHD simulations: we numerically extract the three-dimensional magnetic power spectrum from the MHD runs, and use this in Eqs.(10) and (11) to determine the semi-analytical GRF prediction for f Paγ at fixed η a .We then compare this prediction with the empirical distribution of conversion ratios calculated from Eq. ( 6) for the ensemble of trajectories through the MHD runs.

V. RESULTS AND CONCLUSIONS
Figure 1 shows a cross-section of the magnetic field for Run S in the left panel, the extracted power spectrum (P 3D ) and the helicity (H) in the right panel, and a realisation of a GRF magnetic field with the same P 3D in the central panel.
Clearly, fields are visually very different due to the appearance of coherent structures in the MHD simulation.For each run we used 5122 straight lines as photon paths to calculate the conversion ratio using the perturbative formalism, cf.Eq. ( 6), which we checked to be in agreement with a complete numerical solution (see Ref. [40,68]) at the ∼ 0.1% level in the relevant regime. 2  Each row of Fig. 2 shows, for a fixed η a , the distribution (PDF) of P aγ (η a ) in Run K (left column) and Run S (right column) together with the analytical, GRF prediction from Eq. ( 11) (black solid line).The red and blue points correspond to different linear polarisations of the photon flux (i.e.i = x, y).The agreement between the red and blue histograms for typical, small values of the probability is a reflection of the approximate statistical isotropy of the MHD magnetic field.
All models agree well for the prediction of the smallamplitude conversion ratios that are the most common, and the MHD runs reproduce the exponential decay of the GRF field.However, small-amplitude features in observational photon spectra can be buried by systematic uncertainties and measurement errors.For sufficiently large conversion ratios, the predictions from MHD runs differ significantly from the GRF predictions: for rare fluctuations, the MHD distributions exhibit a spectral break followed by a distinctly heavy tail.Consequently, large conversion ratios are more frequent in MHD models than in GRFs.We note that the average conversion probability depends only on P 3D and is, therefore, the same for MHD and GRF models.This explains the additional, fractionally-small difference, visible in Fig. 2, between the GRF and MHD conversion ratios at intermediate values of 4P aγ /g 2 aγ .We note that the appearance of the heavy tails does not originate from poor statistics, as is evident from the small error bars of the data points in Fig. 2.
Unsurprisingly, the highly non-Gaussian magnetic field from the kinematic phase of Run K results in the largest deviations from the GRF predictions.In this case, a strong deviation from the exponential law occurs in 1-5% of the realisations for each of the sampled η a .The distribution from the saturated phase of Run S deviates from the GRF predictions more strongly at low η a (i.e.high energy).
The independence of P aγ on the phases of B, cf.Eq. ( 6), ensures that the mean value, ⟨P aγ ⟩, is fixed by the power spectrum and is therefore the same for the MHD and GRF models.However, the heavy tails are reflected in enhanced higher moments of P aγ .
Quantitatively for a random variable µ with mean μ and variance σ 2 , we indicate the the skewness as S = ⟨(µ− μ) 3 ⟩/σ 3 and the kurtosis as K = ⟨(µ− μ) 4 ⟩/σ 4 .These quantities are shown in each panel of Fig. 1 corresponding to the conversion ratio PDFs.In most cases, both S and K are significantly larger than the values for an exponential distribution of the GRFs, S = 2 and K = 9.
There are two distinct non-Gaussian features that can contribute to the physical origin of the heavy tails of f P (ηa) : the presence of extended coherent structures, and larger amplitude peaks (spikes) of the magnetic field compared to a GRF.In the context of TeV-scale ALPphoton conversion in magnetic fields from cosmological simulations, Ref. [51] suggested that coherent structures (alone) can explain heavy-tailed distributions.We show that both effects contribute, see App.D. However, in our MHD simulations, peaks of the magnetic field tend to be located in coherent structures, and we expect the effects to be correlated.
In many cases, the energy-dependent conversion ratio for a single sightline is of direct observational interest.It goes beyond the scope of this paper to discuss the effect of MHD magnetic fields on the limits placed on ALPs as this would require significantly extended modelling (e.g. of the, on average, radially decreasing plasma frequency and magnetic field strength in a galaxy cluster) and data analysis (using real data, e.g. from precision Xray observations, cf.[16]).However, in Appendix G, we construct sets of mock X-ray observations and ask: what is the probability that a random line-of-sight taken from an MHD model is 'sufficiently different' from an analogous GRF realisation?In more detail, we consider a hypothetical point-like source in a galaxy cluster and compare the energy-dependent residuals obtained from MHD and GRF models after ALP-photon interconversion.The (mock) observed photon spectra are generated by multiplying a power-law spectrum, representing the primary X-ray flux, by photon survival probabilities calculated, respectively, with the MHD or GRF models.The mock spectra are fitted by a power-law and the statistical distribution of the residuals from the GRF models is determined.We then use the Anderson-Darling test statistic [69] to quantify the differences between the MHD and GRF models.For 5-8% of the MHD spectra, we can reject the null hypothesis that the residuals come from the GRF distribution at 95% confidence level; for the vast majority of realisations, the Anderson-Darling statistic does not suffice to distinguish MHD models from the GRF distribution.This conclusion holds even for experiments with high energy resolution (∼ O(eV)), like Chandra diffraction grating observations or the future Athena mission [25,70,71].It may be possible to improve the sensitivity to MHD effects by designing dedicated statistical tests, however, we expect GRF modelling to suffice for most applications.
In conclusion, we have presented the first systematic study of ALP-photon conversions in magnetic fields obtained by turbulent MHD simulations and compared the predictions to those of simpler GRF models with the same power spectrum.We showed that the typical predictions agree between these models, even though the non-Gaussian fields from MHD simulations result in heavy-tailed distributions of conversion ratios, corresponding to larger-than-expected oscillations in photon spectra from sources embedded in galaxy clusters.This effect is typically small for observables sensitive only to single sightlines through a cluster, in which case the GRF modelling is expected to be sufficient.
An important future direction is to study MHD models of realistic cluster magnetic fields (including the radial decrease of the field strength), and to test the robustness of published limits.Moreover, it would be interesting to further develop our cluster MHD models by including, e.g.anisotropic viscosity [72] and apply these to searches for ALPs in observational X-ray data.

Data availability.
The source code used for the simulations of this study, the Pencil Code, is freely available from Ref. [62].The simulation setups and the corresponding data along with animations for Runs K and S are freely available from Ref. [73].In this appendix, we briefly review how GRFs with a given power spectrum can be explicitly generated, and how the probability distribution of the conversion ratio is derived.Since the conversion ratio is determined by the onedimensional Fourier transform, we consider a one-dimensional lattice of 2N + 1 points with spacing dz = L/(2N + 1) along the direction of propagation.At each lattice point, we generate independent, random numbers, W h , from a normal distribution with mean zero and unit variance The discrete Fourier transform of this white noise is given by To construct a GRF with a given one-dimensional power spectrum, P 1D (η), we define the Fourier transform of this field as where η j = 2π|j|/(N dz) and P1D (η j ) = (L/2π)P 1D (η j ).The components in Eq. (A3) directly determine the ALPphoton conversion ratio and the statistical distribution of P aγ (η j ) (at fixed η j ) only depends on the distribution of Bj (at fixed j).The moments of the conversion ratio are given by Given the moments, the characteristic function (the expectation value of e itPaγ (ηj ) ) is easy to evaluate, and the probability distribution of the perturbative conversion ratio can be extracted by a simple Fourier transform, Appendix B: Non-Gaussian random field In this appendix, we construct an analytically tractable example of a non-Gaussian magnetic field which we then use to demonstrate how non-Gaussianity leads to heavy tails for the ALP-photon conversion ratio.Non-Gaussian fields can be constructed in many ways.Here, we follow the approach of Ref. [74] and construct a non-Gaussian random field through a non-linear operation on a GRF.We denote the Gaussian field by Wj , defined as in Eq. (A2), and take the non-Gaussian field to be given by where, for j = 1, ...N , The j = 0 mode is defined as f ϵ ( W0 ) = W0 and the remaining, negative modes are fixed by the reality condition for the magnetic field in real space B * j = B−j .We now consider j > 0 unless explicitly stated otherwise.The factor of 1 + ϵ 2 2 in Eq. (B2) ensures that the average magnetic field does not depend on the parameter controlling the strength of the non-Gaussianity, ϵ.The real and imaginary parts of the non-Gaussian field are explicitly given by Re( Bj ) = P1D (η j ) Re( Wj ) + ϵ Im( Wj ) , Im( Wj ) .

(B3)
It is possible to derive the PDF of the non-Gaussian field Bj for an arbitrary power spectrum of W .In general, if n random variables ξ i follow a distribution f ξ (ξ 1 , ...ξ n ), then the new variables defined as ψ j = g j (ξ 1 , ..., ξ n ) have the following PDF From Eq. (B2), the inverse transformation is , and the inverse Jacobian is then the PDF of Bj is Clearly, the Gaussian distribution is recovered as ϵ → 0. Given the PDF f Bj , we now calculate the probability distribution for the conversion ratio, P aγ , by imposing the constraint that the probability is proportional to where p = 4P aγ /g 2 aγ P1D (η j ) and b 1 , b 2 respectively denote the real and imaginary parts of Bj .The resulting PDFs for the non-Gaussian case with ϵ = 1 and the Gaussian case with ϵ = 0 are shown in Fig. 3.The non-Gaussian magnetic field produces a characteristic heavy tail, with increased support for high conversion ratios.The heavy tail is reflected in the higher moments of the distribution of P aγ (η j ): The mean of the conversion ratio is unchanged as the power spectrum of the magnetic field is, by construction, the same for the Gaussian and non-Gaussian magnetic fields.However, the higher-order moment all increase for ϵ ̸ = 0 due to non-Gaussianities, which is indicative for the heavy-tailed probability distribution.

Appendix C: Heavy tails for lighter axions
In the bulk of this paper, we have considered the case of m a > ω pl in which numerical simulations can be matched by detailed analytical results.In this appendix, we extend this discussion to the case of arbitrary ω pl (z)/m a and show that, also in this case, MHD magnetic fields lead to heavy-tailed probability distributions.
The ALP-photon transition amplitude is, to leading order in perturbation theory, given by [50,75] where and λ = 1/ω.The function ϕ(z) is not, in general, linear in z, and is even non-monotonous if there are 'resonance points' (i.e.level-crossings) where ω pl − m a changes sign along the trajectory.We note that the amplitude can still be expressed as a sum of Fourier transforms (of B i /|ϕ ′ |) and efficient numerical evaluation is possible through methods like the Fast Fourier Transform, though we will not use these methods here.The new conceptual issue of this general case is that regions where |ϕ ′ | is small are expected to contribute more to the amplitude than similar regions with large ϕ ′ .Naively, one might expect that resonance points give large contributions to the amplitude, and that the amplitude would be rather insensitive to the magnetic properties far away from the resonance points.If so, non-Gaussian peaks of the magnetic field would only contribute significantly if located at a resonance point, which could lead to very different properties of the distribution of the conversion ratio.We now show that this intuition is not correct: resonant contributions are typically small compared to the cumulative, non-resonant contributions, and the heavy tails of the non-Gaussian magnetic field persist also in this general case.
From our MHD runs, we first extract the plasma frequency as a function of position.The average plasma frequency value for the two Runs is ωpl = 6.59 × 10 −12 eV, and fluctuations around the mean are small, indeed much smaller than the expected range of the plasma frequency in a galaxy cluster.We therefore re-scale the fluctuations to an amplitude of ∼ 3ω pl around the mean value and choose to consider ALP masses that have one or more resonance points within the simulation box.We then calculate the conversion ratio as described in the main text of this paper.
Figure 4 shows the conversion ratio as a function of z along two trajectories within Run S. Each trajectory has two resonance points where ω pl (z)/m a = 1.Clearly, around the resonance points, there is no significant increase in the conversion ratio which, as in the massive case, is determined by the properties of the magnetic field along the entire trajectory.This result is consistent with the general discussion in [50].
Figure 5 shows the numerical distribution of the conversion ratio, P aγ , in the MHD runs with m a = ωpl .The solid lines indicate the analytical predictions for the same ALP mass but considering the massive case (discussed in the main text) where the plasma frequency is negligible.The comparison of the analytical result (solid line) and the numerical simulations (histograms) is therefore not rigorous, but the appearance of heavy tails at large conversion ratios is apparent in both figures.We checked that similar results are obtained in case of very light ALPs, m a ≪ ω pl .This indicates that the conclusions of this work hold over the entire ALP parameter space.where we defined More generally, for GRFs the helicity always enters the correlation functions through an even power of f H , and hence gives non-negative contributions to higher moments, and can contribute to heavy tails of the probability distribution of P aγ (cf. the non-Gaussian example of appendix B. For non-Gaussian fields, the situation is model dependent as there are additional contributions involving powers of f H and correlators of an odd number of magnetic fields.However, as we will now argue, in practice helicity tends to be negligible for ALP-photon conversion also for MHD fields.In general, the magnetic fields generated by small-scale dynamo with non-helical forcing are non-helical.To provide a quantitative comparison for the role of helicity of the magnetic field in P aγ , we consider a scenario where the helicity spectrum (H) is proportional to the power spectrum (M N = P 3D /2) and the ratio 2H/P 3D is either 1 i.e. the fully helical magnetic field case or 0.1.Figure 7 shows the ratio f 2 H /3P 2 1D for these two cases.To estimate this ratio we considered a 3D power spectrum given by [76] P 3D (η a , k ⊥ ) = 1 where k p = (2π/10)kpc −1 is the wavenumber corresponding to the peak of the spectrum.From Fig. 7, it is concluded that the impact of the helicity on the photon to ALP conversion ratio is negligible unless the magnetic field is fully helical and η a > k p .In applications it means that it is more important for low ALP energy, where typically the ALP-photon conversions are suppressed.

Appendix F: Magnetohydrodynamic simulations
In this appendix, we review details of the properties of the MHD simulations.Magnetic fields in galaxy clusters are thought to be generated through gas motions by a dynamo process.The gas motions are turbulent, but their driving is not well understood.Magneto-thermal instabilities have been studied in this context [53,55], and even the driving by turbulent wakes of individual galaxies has been discussed [57,58].Turbulence could also be driven by sporadic cluster mergers [60,61], but the turbulence would be slowly decaying between such events.
The dynamo process itself is a generic phenomenon that is mainly characterized by a large magnetic Prandtl number [77].This means that the magnetic diffusivity η is much smaller than the kinematic viscosity ν.This is numerically difficult to handle.In addition, the magnetic Reynolds number is huge (of the order of 10 30 ).This means that the turbulent magnetic cascade extends down to very small length scales.As a compromise, since the energy contained in the smallest length scales becomes very weak, we consider here turbulence at a moderate magnetic Reynolds number that is as large as possible for a given numerical resolution.The turbulence is assumed to be driven through some external volume forcing as is commonly done in numerical simulations of homogeneous turbulence.We also ignore the density stratification and just consider an isothermal gas with constant sound speed.As advertised above, we consider a dynamo-generated magnetic field B driven by forced turbulence.We solve the induction equation for the magnetic vector potential A with B = ∇ × A in the Weyl gauge, i.e.
Here, we have adopted Gaussian Heaviside units where J = ∇ × B is the current density, and u is the velocity which obeys the momentum equation, where are the components of the traceless rate-of-strain tensor S, f is the forcing function, and ϱ is the density, which obeys the continuity equation, written here as We solve these equations in a periodic domain of size L 3 with 512 3 mesh points using the Pencil Code [62].The forcing function consists of random sinusoidal waves that are δ-correlated in time, i.e. the forcing function changes at each time step.The wave vectors are taken from a shell of finite thickness and radius k f , which we chose to be close to the smallest wave number of the computational domain k 1 ≡ 2π/L.The strength of the forcing function is arranged such that the turbulent Mach number Ma = u rms /c s , where u rms is the rms velocity, is around unity or less.We define the kinetic and magnetic Reynolds numbers as Re = u rms /νk f and Re M = u rms /ηk f , respectively, and the Lundquist number as Lu = v rms A /ηk f , where v rms demonstrate how this works for a mock X-ray observation.We consider hypothetical X-ray observations of a bright, power-law source in a galaxy cluster, observed by instruments with energy resolutions comparable to current and future experiments.We proceed as follows: • We consider a source with intrinsic power-law spectrum with index −2, and modulate these by the energydependent conversion ratios calculated in the MHD and GRF model.The observed photon flux is tabulated as a function of the energy for a set of energy-spacings, where a finer binning roughly corresponds to an X-ray experiment with higher energy resolution.The final results are insensitive to the value of the spectral index.
• The resulting spectrum is fitted with a power law, and the residuals are extracted.
• For the GRF model, we determine the distribution of the residuals by scanning 2,500 realisations.
• For the residuals from each MHD realisation, we test the null hypothesis that the residuals come from the GRF distribution.
The resulting cumulative distribution function (CDF) of the residuals for a GRF is shown in Fig. 9 (dashed lines) and compared with single realisations of each of the two MHD models considered (solid lines).Intuitively, one may expect more negative residuals for the MHD case compared to the GRF model since the heavy tails give larger conversion ratios and deeper 'dips' in the observed spectrum.This effect is indeed visible in Fig. 9, where the empirical CDFs of the MHD realisations grow faster than the GRF CDF.We use the Anderson-Darling test to determine the significance of the heavy tails for single-sightline spectra.The Anderson-Darling test is a non-parametric statistical test that can be used to compare a sample with a theoretical model, to determine whether the sample is generated by the theoretical underlying distribution.This test is more sensitive to differences in the distributions' tails than other distributional equality tests, such as the Kolmogorov-Smirnov test.This makes it particularly useful for detecting differences in extreme values.Its non-parametric nature and sensitivity to differences in the tails of the distributions make it a valuable tool in this context.We determine the distribution of the test statistic [69] where n is the sample size and Y is an ordered list of residuals.This step allows us to determine the 95% confidence level for the variable A 2 .These results are shown in Tab.II for different values of the energy resolution and the two magnetic field models considered.For illustrative purpose, we show the distribution of the statistical test variable A 2 in Fig. 8, comparing the GRF case (red histograms) with run K (blue histogram, left panel) and run S (blue histogram, right panel).Slightly larger values of A 2 are preferred by the MHD models compared to the GRF case.Note that we consider N div bins, determining an equally spaced grid in terms of Log 10 ω, in the 0.1 − 10 keV range in which we calculate the simulated photon signal.The energy resolution ∆E (first column of Tab.II) is the smallest energy interval that we simulate.
The critical values in Tab.II are used to determine the similarity between a MHD realization and the GRF prediction.If the test statistic is outside the interval determined by the critical values, the null hypothesis that the   III.Probability that a line-of-sight of an MHD model is significantly different (as defined in the text) from the GRF case.The first column refers to the energy resolution that we use to simulate the detected signal, the second and third ones to the probability for the two MHD models considered.Here the ALP parameters are ma = 10 −12 eV and gaγ = 5×10 −13 GeV −1 .
the MHD residuals come from the same distribution of the GRF ones, is rejected.This comparison is performed over 2500 lines-of-sight for the MHD models.These lines-of-sight are uniformly distributed in the innermost 190 kpc of the simulation box and connect two opposite faces of the box.Finally, we count the number of times the hypothesis that the residuals are extracted from the same distribution of the GRF is rejected at the 95% confidence level.Therefore, we can estimate the probability that a P aγ vs ω realization of an MHD run is significantly different from the GRF case.These results are shown in Tab.III.This analysis allows us to conclude that, even with high energy resolution, only in the 5 − 8% of the cases the observed line-of-sight differs significantly from the predictions of a GRF model.This suggests that GRF models of turbulent cluster magnetic fields suffice to reliably determine the predictions for X-ray searches for ALPs for observables probing a single sightline.

FIG. 1 .
FIG. 1. Cross-section of the magnetic field for (a) Run S and (b) a GRF model with same the power spectrum.The red regions highlight locations with |B| > 3Brms.(c) Power and helicity spectra for Run S.

FIG. 2 .
FIG. 2. Histogram of the conversion ratio extracted from Runs K (left column) and S (right column) for an initial photon flux that is linearly polarised along the x (blue) and y (red) directions together with the analytical GRF prediction of Eqs.(10)-(11) (black line).Here ηa = 7.8 × 10 −3 kpc −1 (top row), ηa = 0.07 kpc −1 (middle row), ηa = 0.195 kpc −1 (bottom row).Skewness S and kurtosis K of the distributions are shown with the corresponding colour.For comparison, S = 2 and K = 9 for an exponential distribution.Approximate locations and cumulative probabilities of the heavy tails are respectively indicated by dashed vertical lines and percentages.The error bars are obtained by assuming Poissonian fluctuations.

ACKNOWLEDGMENTS
We would like to thank Riccardo Z. Ferreira for discussions.The work of PC, DM and EM is supported by the European Research Council under Grant No. 742104 and by the Swedish Research Council (VR) under grants 2018-03641 and 2019-02337.The work of RS and AB is supported by VR grant 2019-04234.Nordita is sponsored by Nordforsk.We acknowledge the allocation of computing resources provided by the Swedish National Allocations Committee at the Center for Parallel Computers at the Royal Institute of Technology in Stockholm and Linköping.
Appendix A: Gaussian random field (GRF)

FIG. 3 .
FIG. 3. Probability distribution of the conversion ratio for fixed values of gaγ and ηj in the case of a Gaussian field (blue line) and the non-Gaussian model (red line).The analytical predictions of Eqs.(A7) and (B8) are in excellent agreement with the numerical simulation for 5 × 10 4 realisations, ηj = 0.117 kpc −1 , and a magnetic field extended over 100 kpc.

FIG. 4 .FIG. 5 .
FIG. 4. Conversion ratio (solid lines) as function of the distance along the z axis in Run S for two trajectories represented by the two colours.The ALP mass is ma = 6.59 × 10 −12 eV, the energy ω = 10 keV and the dashed lines show ω pl (z)/ma.The grey line indicates the resonance condition, ma = ω pl .

FIG. 7 . 2 1D(
FIG. 7. The ratio of f 2H (contribution of magnetic helicity spectrum to the ALP to photon conversion ratio square) to 3 P 2 1D

A
FIG. 8. Comparison of the distribution for the Anderson-Darling variable A 2 for GRF models (red histograms) with the same power spectrum of run K (blue histogram, left panel) and run S (blue histogram, right panel) models.The vertical black dashed lines denote the 95% CL interval for the test variable.Here we used ∆E = 4.7 eV, ma = 10 −12 eV and gaγ = 5 × 10 −13 GeV −1 .

FIG. 9 .
FIG.9.Comparison of the theoretical CDFs (dashed lines) extracted from GRF models with the same power spectrum of run K (black lines) and run S (red lines) models.This prediction is compared with two CDFs extracted from the two MHD models considered (solid lines).Here we used ∆E = 4.7 eV, ma = 10 −12 eV and gaγ = 5 × 10 −13 GeV −1 .

TABLE I .
Parameters of each magnetic field model used.

TABLE II .
Critical intervals for the stochastic variable A 2 for different energy resolutions.The first column refers to the energy resolution that we use to simulate the detected signal.The 95% of the GRF realizations give a value of A 2 in the range in the second and third columns of the table, referring to the two different MHD models.Here the ALP parameters are ma = 10 −12 eV and gaγ = 5 × 10 −13 GeV −1 .