Turbulent axion-photon conversions in the Milky-Way

The Milky-Way magnetic field can trigger conversions between photons and axion-like particles (ALPs), leading to peculiar features on the observable photon spectra. Previous studies considered only the regular component of the magnetic field. However, observations consistently show the existence of an additional turbulent component, with a similar strength and correlated on a scale of a few 10$\,$pc. We investigate the impact of the turbulent magnetic field on the ALP-photon conversions, characterizing the effects numerically and analytically. We show that the turbulent magnetic field can change the conversion probability by up to a factor of two and may lead to observable irregularities in the observable photon spectra from different astrophysical sources.

The Milky-Way magnetic field can trigger conversions between photons and axion-like particles (ALPs), leading to peculiar features on the observable photon spectra. Previous studies considered only the regular component of the magnetic field. However, observations consistently show the existence of an additional turbulent component, with a similar strength and correlated on a scale of a few 10 pc. We investigate the impact of the turbulent magnetic field on the ALP-photon conversions, characterizing the effects numerically and analytically. We show that the turbulent magnetic field can change the conversion probability by up to a factor of two and may lead to observable irregularities in the observable photon spectra from different astrophysical sources.

I. INTRODUCTION
Conversions between photons and axion-like particles (ALPs) in cosmic magnetic fields are the subject of intense investigations (see, e.g., [1][2][3][4][5][6][7][8][9][10]). If ALPs exist, the analysis of the gamma-ray spectra from Galactic and extragalactic sources may reveal valuable information on their coupling with photons g aγ and mass m a . Several specific cases have been studied so far. Conversions of photons into ALPs in cosmic magnetic fields of Galactic or extra-galactic origin may imprint peculiar deformations on the spectra of very-high-energy photons from faraway sources, such as blazars, active galactic nuclei, pulsars and galaxy clusters (see, e.g. [7][8][9][11][12][13]), and alter the polarization of X-rays in a measurable way [14,15]. Dimming effects on the photons from Supernovae Ia induced by photons oscillating into ALPs and viceversa in the intergalactic magnetic field have also been characterized [16][17][18]. Finally, ALPs produced in the hot core of massive stars such as red supergiants [19,20] and corecollapse supernovae [21][22][23], and converted in the Milky-Way magnetic field, would produce a copious photon flux.
Because of the above considerations, the Milky-Way magnetic field is recognized as a valuable tool to understand the properties of ALPs, triggering ALP-photon conversions for both Galactic and extra-galactic sources (see, e.g. [3,5]). The Milky-Way magnetic field models adopted in these studies (see, e.g. [7,9]) make use of state-of-the-art simulations, based on the combined * Electronic address: pierluca.carenza@ba.infn.it † Electronic address: carmelo.evoli@gssi.it ‡ Electronic address: MGiannotti@barry.edu § Electronic address: alessandro.mirizzi@ba.infn.it ¶ Electronic address: daniele.montanino@le.infn.it analyses of Faraday rotation measurements of extragalactic sources and of polarized galactic diffuse radio emission [24,25]. Typically these simulations resolve the large-scale fields (also called regular or coherent) indicating a magnetic field with strength B ∼ O(1) µG that is smooth on Galactic scales, usually assumed to follow the spiral arms [26]. However, radio observations of synchrotron emission and polarization reveal also a smallscale component of the magnetic field connected to the turbulent Inter Stellar Medium (see, e.g., [27]). This component is known as the random or turbulent magnetic field, and has an amplitude comparable to that of the regular magnetic field, but a much shorter correlation length, l corr ∼ O(10 − 100) pc [28][29][30][31]. The impact of the turbulent component has been neglected in previous investigations, on the basis that the ALP oscillation length would necessarily be much larger than the correlation length of the turbulent field (see, however, [3]) and so the effects of the random field would be averaged to zero over the oscillation length. However, though it is true that in the most interesting cases l osc is much larger than the correlation length l corr of the random field, this condition does not justify neglecting the turbulent component of the magnetic field. On the contrary, we will show that in certain cases and depending on the ALP parameters, the turbulent component may contribute as much as the regular component to the oscillation probability. In general, as we shall see, the presence of a sizable turbulent component in the Galactic B-field leads to peculiar irregularities in the observable photon spectra from different astrophysical sources.
The plan of our work is as follows. In Section II, we present the equations of motion of the ALP-photon mixing in the Galactic magnetic field accounting for the regular and the turbulent component. We discuss our numerical and analytical approach to characterize the conversions in turbulent B-fields and we present different representative examples. In Section III, we show the spectral features induced by the conversion in the turbulent B-field observable in photon spectra from different sources. Finally, in Section IV we discuss our results and draw our conclusions. There follow two Appendixes, where we give details on our analytical recipe to calculate the ALP-photon conversion probabilities in a regular plus turbulent magnetic field configuration in a perturbative (Appendix A) and non perturbative (Appendix B) regime.

A. Equations of motion
The Lagrangian describing the ALP-photon system is The first term in Eq. (1) is the QED Lagrangian for photons is the electromagnetic field tensor, F µν = 1 2 µνρσ F ρσ is its dual, α is the fine-structure constant and m e is the electron mass. The second term on the r.h.s. of Eq. (2) is the Euler-Heisenberg-Weisskopf (HEW) effective Lagrangian [32], which accounts for the one-loop corrections to classical electrodynamics. The Lagrangian for the non-interacting ALP field a in Eq. (1) is where m a is the ALP mass. Finally, the ALP-to-photon interaction in Eq. (1) is represented by the following Lagrangian term [32] where g aγ is the ALP-photon coupling constant. We shall consider a monochromatic photon/ALP beam of energy E propagating along the z-direction in the presence of an external magnetic field B. The linearized equations of motion for the photon/ALP system are given by [32] −i d dz where A x (z) and A y (z) are the two photon linear polarization amplitudes along the x and y axis, respectively, a(z) denotes the ALP amplitude and M represents the ALP-to-photon mixing matrix.
The mixing matrix M takes a simpler form if we consider the case of a photon beam propagating in a single magnetic domain with an homogeneous field inside. We denote by B T the transverse magnetic field, i.e. its component in the plane normal to the beam direction. We can choose the y-axis along B T so that B x vanishes. Under these simplifying assumptions, the mixing matrix can be written as [2] whose elements are [32] ∆ ≡ ∆ pl + 3.5 ∆ QED , with where n e is the electron density in the medium, ω 2 pl = 4παn e /m e is the associated plasma frequency, B cr ≡ m 2 e /e 4.41 × 10 13 G is the critical magnetic field and e denotes the electron charge. Concerning the ALP-tophoton coupling g aγ we took a representative value below the CAST bound for ultralight ALPs from solar searches, g aγ < 6.6 × 10 −11 GeV −1 [33], while the strength of the magnetic field B T and the electron density n e are given by typical conditions in the Milky-Way [34]. For the range of energies we are considering in this work, ∆ QED is negligible and will be neglected hereafter.

B. Conversions in the regular magnetic field
The ALP-photon oscillation probability takes a particularly simple form in the case of a uniform magnetic field. This is a good approximation for the regular magnetic field component, since our ultimate goal is to show the effects of adding a turbulent component over it. Our qualitative results are not modified when this assumption is lifted.
With this assumption, the ALP-photon conversion probability can be calculated analytically. Considering a photon initially polarized along the y axis, the probability of conversion into an ALP after a distance r reads [32] where the mixing angle θ is given by and the oscillation wave number is In order to gain more physical intuition, it is convenient to rewrite Eq. (13) as where the sinc function is defined as sinc x ≡ x −1 sin x. The probability in Eq. (16) has different limiting regimes, depending on the relative importance of ∆ aγ and |∆ a − ∆ | in Eq. (13). Specifically, when ∆ aγ |∆ a − ∆ |, the probability scales as sin 2 (∆ aγ r/2). However, since ∆ aγ r 1.5×10 −2 g 11 (B/µG)(r/kpc) 1 for representative ALP-photon couplings, g 11 ≡ g aγ /(10 11 GeV −1 ) ∼ 1, galactic distances, and typical values of the galactic magnetic field, for all practical purposes the oscillation probability in this case reduces to On the other hand, if ∆ aγ |∆ a − ∆ | we can conveniently rewrite the conversion probability as Eq. (18) shows that for Φ 1, the dependence on the energy in the oscillation probability disappears, so that we recover again Eq. (17). Furthermore, in these conditions the probability grows quadratically with the distance, cancelling the dependence on the distance of the ALP flux. Of course, these results depend on the assumption that the magnetic field is correlated on a scale larger than r.
As obvious from Eq. (19), for a massless ALP the condition Φ 1 is satisfied for distances r 20 kpc at energies of E ∼ 10 keV. A finite ALP mass may spoil this condition. The exact mass threshold for this to happen depends on the distance r, and on the ALP energy. In particular, high energy ALPs may be correlated on a larger distance and for larger masses. On the other hand, for large ALP masses Φ 1 and the oscillation probability acquires a dependence on the energy and decreases rapidly with the axion mass, P aγ ∝ m −4 a , as evident in Figure 1, where we show the behavior of the conversion probability P aγ for a source at distance r = 1 kpc for a value of the regular magnetic field B reg = 3 µG in function of the ALP mass m a . In the Figure, we fixed the axion coupling to g aγ = 10 −11 GeV −1 and take three representative ALP energies: E = 50 keV (black curve), representative of the typical energy of ALPs emitted from a red supergiant star during its late evolutionary stages [20]; E = 100 MeV (red curve), representative of the energies of ALPs emitted from supernova explo- The black curve is the average probability for a turbulent field with Brms = 1 µG and a correlation length lcorr = 10 pc, while the grey-band represents the "1σ" dispersion around the average. The blue line is the conversion probability for one possible random realization of the magnetic field. Note that lcorr = l cell /2, where l cell is the cell length.

C. Conversions in the turbulent magnetic field
A realistic description of the photon-ALP oscillation in the Galactic magnetic field requires a model for the magnetic field with both regular and turbulent components. The most commonly adopted model for the regular component is the Jansson-Farrar model [24] (see, e.g., Refs. [7,9,35]), though the Pshirkov model [25] is a common alternative choice.
The small-scale turbulent field is less well known. It is usually assumed to follow a power law with a certain outer scale, where energy is injected, which then cascades down to smaller turbulent scales until energy dissipates at the dissipation scale [36]. The strength of this random magnetic field component can be estimated from rotation measure fluctuations, combined with an estimate of the thermal electron density, giving a ratio of random to regular magnetic field components of 1 [37,38]. This value is compatible with the predictions of large-scale magnetic field models that include the turbulent magnetic field as a free parameter [39].
Concerning the correlation length, Ref. [29] found a large l corr ∼ 100 pc, using structure functions of rotation measures and averaging over large parts of the sky. Ref. [30] confirmed this large outer scale for interarm regions in the Galactic plane; however, they found a much smaller outer scale ∼ 10 pc in the spiral arms. Also, arrival anisotropies in TeV cosmic-ray nuclei can be best explained by a turbulent magnetic interstellar medium on a maximum scale of about ∼ 1 pc [31]. More recently, an upper limit of ∼ 20 pc for the outer scale of the magnetic interstellar turbulence toward the Fan region has been found in Ref. [28] by measuring with LOFAR the fluctuations in the diffuse synchrotron emission. We now consider different representative situations showing the impact of the turbulent component of the Galactic magnetic field on the ALP-photon conversion. For definiteness, we take a source at a distance r = 1 kpc. There are several Galactic ALP sources, such as red supergiants, within a radius of O(1) kpc [40]. ALPs produced in these sources would propagate in a magnetized medium where they can convert into photons. In our examples below, we take the regular component of the magnetic field to be uniform on such a short length scales. In any case, scales of 1 kpc are not resolved in large scales simulations. On the other hand, for the turbulent component we assume a correlation length l corr = 10 pc. The choice of a very short correlation length is conservative since, as we shall see below, increasing l corr enhances the effect of the turbulent field on the conversion probability [Cf. Eq. (20)]. In the description of the turbulent Galactic magnetic field we adopt a simplified approach, assuming it along a given line-of-sight as a network of domains with a fixed correlation length l corr = 10 pc. 1 Finally, we fix B reg = 3 µG and assume a turbulent component with a root-mean-square (rms) amplitude B rms = 1 µG [24]. As a representative situation, we consider ALPs with coupling g aγ = 10 −11 GeV −1 , a value slightly below the current experimental and astrophysical bounds, and energy E = 50 keV. We will vary the ALP mass in order to probe different cases of ALP oscillation length.
In the conditions described above, an ALP mass above about 1 neV corresponds to a correlation length l osc = 2π/∆ osc 4 pc l corr . The associated conversion probability, in this case, is negligibly small, P aγ 10 −9 . Therefore, we will ignore this case hereafter and limit ourselves to situations in which l osc > l corr . Within this hypothesis, we can show that when the conversion probability is small, the effects induced by the turbulent magnetic field can be evaluated analytically. As shown in Appendix A, the average probability over all possible field configurations in function of the distance r from the source is remarkably simple, where P (0) aγ (r) is the probability assuming only a regular magnetic field. We remark that Eq. (20) is valid also when the regular magnetic field is not homogeneous, as long as P aγ 1. This is the typical situation for conversions in the Milky Way, assuming allowed values for the ALP-photon coupling and a realistic magnetic field. When P aγ ∼ O(1) the previous approach is no longer valid and in general we have to resort to a Monte Carlo simulation. However, in the hypothesis of δ-correlated perturbations it is possible to modify the evolution equations in order to calculate arbitrary moments of the variable P aγ . Although this situation is outside the scope of this work, for completeness we treat this case in Appendix B.
A full derivation of Eq. (20) is provided in Appendix A. However, it is easy to recognize the validity of this result also using very simple arguments. Assuming l corr l osc , it is possible to express the average conversion probability on a given magnetic domain as P (0) aγ,turb (g aγ B rms l cell /2) 2 = (g aγ B rms l corr ) 2 , where the size of a single magnetic domain is given by l cell = 2l corr (see Appendix A). Since the photons travel N = r/l corr magnetic configurations, the incoherent total conversion probability on the turbulent configuration is given by P aγ,turb N P (0) aγ,turb = g 2 aγ B 2 rms l corr r (see also Ref. [3]). We notice that since the turbulent field has a random configuration, P aγ is a random variable itself. In the hypothesis that the turbulent magnetic field has a gaussian distribution, x,y /2B 2 rms ), we can calculate the probability density function F (P aγ ) in the limit in which the diagonal terms in the Hamiltonian in Eq. (6) are small (see Appendix A). If this condition is not fulfilled, it is very difficult to obtain a closed form for the probability distribution function (p.d.f.) but it is still possible to get handy formulas for the calculation of the momenta of P aγ . In particular one can obtain the second momentum of P aγ [Eq. (A11)] and thus the standard deviation of the distribution, σ(r) = [ P 2 aγ (r) − P aγ (r) 2 ] 1/2 . In Appendix A we compare our analytical results with a numerical simulation to show the agreement between the two approaches.
In Figs ios. In particular, in Figure 2 we consider the case m a = 10 −10 eV, corresponding to l osc 400 pc. The red curve is the conversion probability for a purely regular B-field. The black curve is the average probability in presence of regular plus turbulent fields, calculated using Eq. (20). The gray band represents the "1σ" dispersion around the average represented by the grey-band. Here, by "1σ" we mean the standard deviation around the mean value. This cannot immediately translate into a definite confidence level since the p.d.f. is not in general a gaussian. We see that the 1σ dispersion in the probability induced by the turbulent component is seizable, producing variations up to ∼50% with respect to the average. For the purpose of illustration, with the blue line we show the conversion probability for one representative random realization of the magnetic field. The realization of the magnetic field used to draw this representative case is plotted in Fig. 3. We observe that the behavior of P aγ (r) for a single realization is in general irregular, depending on the configuration of the field. In particular, the chosen realization for r 0.5 kpc gives a conversion probability smaller than the purely regular field. In general, some realizations of the turbulent field could result in a conversion probability smaller than the one obtained with a purely regular field. However, the average probability is always larger compared to a regular field.
In Figure 4, we consider an ALP with mass m a = 10 −11 eV, which implies l osc 35 kpc. This case corresponds to an energy-independent oscillation probability. As expected, the conversion probability induced by the regular component of the magnetic field scales as r 2 . The effect of the turbulent component is the spread of the P aγ around the mean. Although on average the effect is very small, for many realizations the difference with the regular case can exceed up to 10% within 1σ. Since in this case ∆ pl , ∆ a ∆ aγ we can use the analytic expression for the p.  Finally, in order to isolate the role of the turbulent field, in Fig. 5 we consider the same case of Fig. 4, but assuming only the turbulent component of the magnetic field. We show the average and the 1σ band for a pure turbulent field with B rms = 1 µG. In this case, as expected from Eq. (20), we see that P aγ (r) grows linearly (at least until P aγ 1). Also in the last two cases in blue we draw a representative conversion probability curves, using the same realizations drawn in Fig. 3.

III. EFFECT ON OBSERVABLE PHOTON SPECTRA
In this Section, we assess the impact of the turbulent component of the Galactic magnetic field in three physically relevant examples. Specifically, we consider i) the case of ALPs produced in the red supergiant Betelgeuse, which can oscillate into photons originating a hard X-ray flux; ii) ALPs produced in a galactic supernova (SN), which can originate a gamma-ray flux; and iii) ALPs produced by conversions from gamma-rays from Galactic pulsars. All these cases have already been considered in the literature. However, the impact of the turbulent magnetic field has always been ignored.
The signatures of the random component of the magnetic field can be revealed through energy-dependent irregularities imprinted on the photon energy spectra. As we shall see, these effects can be more or less pronounced, depending on the ALP parameters, particularly the ALP mass. In general, we expect the impact of the turbulent field to be maximal at the threshold between the energy-dependent and the energy-independent regime [see Eq. (18)], where the conversion probability reaches its maximum before saturating (see Fig. 1).
In all of the examples in this section, we will assume a turbulent field with a strength B rms = 1 µG and correlation length l corr = 10 pc. For simplicity we will consider a single realization of the turbulent magnetic field, different for the three physical cases, avoiding to characterize a  distribution obtained generating different configurations of the field. In all the examples the regular magnetic field is parametrized according to the Jansson and Farrar model [24] along different lines-of-sight.
To start, we consider Betelgeuse, a red supergiant star of about 20 M , at a distance of 200 pc. Betelgeuse's exact evolutionary phase is unknown, except for the fact that it has already exhausted the H in its core. Here, we assume that the star is in the core He-burning stage, which is the most likely scenario since the following stages are considerably shorter. For the numerical model, we refer to [20]. ALPs can be thermally produced in Betelgeuse through the Primakoff process. Their production rate can be expressed as a quasi-thermal spectrum where C is a normalization factor, E 0 represents the average energy, and β is the spectrum index. The values of C, E 0 and β depend on various structural parameters characterizing the core of the star, such as temperature, density and chemical composition. Here, we adopt the parameters in Table I, corresponding to Model 0 in Table S1 of [20] After being produced, the ALPs leave the star unimpeded and may convert into photons in the Galactic magnetic field. From the average ALP energy in Table I it is clear that the resulting photons are expected to have a hard X-ray spectrum. Such a spectrum was searched, with negative results, in a recent dedicated NuSTAR observation [20], leading to the bound g aγ < 2 × 10 −11 GeV −1 for m a < 5 × 10 −11 eV. The analysis, however, ignored completely the random component of the magnetic field.
In the upper panel of Figure 6, we show the ALPphoton conversion probability vs. the energy, for g aγ = 10 −11 GeV −1 and m a = 10 −10 eV. The black curve includes only the regular magnetic field. In this case, we notice an almost periodic behavior for E 30 keV, while at higher energies the energy-independent plateau is reached. The effects of the turbulent component are shown with the red curve in the same figure. In this case, we notice not only a larger conversion probability, expected since the total magnetic field strength is (in average) enhanced. We also see that the periodic pattern observed in the black curve is destroyed and replaced with an irregular behavior.
In the lower panel of Figure 6, we show the observable photon number flux, obtained by convolving the original ALP flux with the P aγ , We also consider a gaussian energy resolution of 1 keV, comparable with the performance of NuSTAR [42]. We observe that the irregular behavior is smeared out by the resolution effect. We find peculiar features in the energy range [10:40] keV, appearing as peculiar bumps and dips. Thus, in principle, the case of just regular field is distinguishable from the case of regular plus turbulent fields. As a second example, we consider the ALP flux produced by Primakoff process in supernovae. The ALP production rate is calculated as in Ref. [43], using a SN model with a 18 M progenitor, simulated in spherical symmetry with the AGILE-BOLTZTRAN code [44,45]. We consider the rate at t pb = 1 s after the bounce. Even in this case, the ALP production rate has the quasi-thermal spectrum given in Eq. (21). The fitting parameters corresponding to the model we are considering are also given in Table I.
Because of the conversion in the Galactic magnetic field, the ALP flux generates a gamma ray flux, with typical energies of O(100) MeV. The lack of a gamma-FIG. 7: Upper panel: ALP-photon conversion probability for a Galactic supernova at r = 10 kpc in function of the energy for gaγ = 10 −12 GeV −1 and ma = 5 × 10 −10 eV for regular (black curve) and regular plus turbulent (red curve) magnetic field. Lower panel: Photon energy spectrum after ALP conversions in regular (black curve) and regular plus turbulent (red curve) magnetic field. A gaussian energy resolution of 10% is assumed. ray signal in the Gamma-Ray Spectrometer (GRS) on the Solar Maximum Mission (SMM) in coincidence with the observation of the neutrinos emitted from SN 1987A therefore provided a bound on ALPs coupling to photons [21,22]. Specifically, for m a = 4 × 10 −10 eV the most recent analysis finds g aγ = 5.3 × 10 −12 GeV −1 [23]. A future ALP burst from a Galactic SN would allow to probe a considerable larger parameter space through a Fermi-LAT observation, if the explosion occurs in its field of view [46].
Here we want to assess the impact of the turbulent magnetic field of these potential observations. For definitiveness, we assume a Galactic SN at r = 10 kpc. We consider a specific line of sight in the Milky Way (b, l) = (75.22 • , 0.11 • ). In the upper panel of Figure 7 we show the ALP-photon conversion probability in function of the energy for g aγ = 10 −12 GeV −1 and m a = 5×10 −10 eV for regular (black curve) and regular plus turbulent (red curve) magnetic field. We observe peculiar wiggles in the energy range E ∈ [1 : 10 2 ] MeV. Larger irregularities are observed in the presence of the turbulent field. In the lower panel we present the corresponding photon spectra assuming an energy resolution of 10%, as expected for Fermi-LAT at those energies [47]. We realize that in the presence of only regular field the energy-FIG. 8: Upper panel: Photon-ALP conversion probability for the Galactic pulsar PSR J2021+3651 at distance r = 6 kpc for gaγ = 3.5 × 10 −12 GeV −1 and ma = 4.4 × 10 −9 eV for regular (black curve) and regular plus turbulent (red curve) magnetic field. Lower panel: Photon energy spectrum after ALP conversions in regular (black curve) and regular plus turbulent (red curve) magnetic field. A gaussian energy resolution of 10% is assumed. dependent modulation is almost washed out by the effect of the resolution, while in the presence of the turbulent component bumpy features below E ∼ 100 MeV would be clearly visible. We notice that this energy range would be below Fermi-LAT sensitivity, but is in the reach of future gamma-ray experiments like eASTROGAM [48] and AMEGO [49].
Finally we consider the case of Galactic pulsars. This case is particularly significant since it was recently claimed by two different groups [9,50] that the unexpected spectral modulation observed by Fermi-LAT in gamma-rays from Galactic pulsars and supernova remnants could be due to conversion of photons into ultralight ALPs in the Milky-Way magnetic field. The bestfit ALP parameters, g aγ = 2.4 × 10 −10 GeV −1 and m a ∼ 4 × 10 −9 eV, are in tension with the CAST bound on solar ALPs. The tension can be lifted in ALP models with environmental dependent mass/coupling [51]. Here, we examine the effects of the previously neglected turbulent magnetic field on these results. For definitiveness, we consider the pulsar PSR J2021+3651 at distance r = 6 kpc from us and with Galactic coordinates (l, b) = (75.22, 0.11) studied in Ref. [9]. The original pho-ton spectrum is given by where N 0 = 0.15×10 −9 MeV −1 cm −2 s −1 , E 0 = 0.8 GeV, Γ = 1.59, E cut = 3.2 GeV. In the upper panel of Figure 8, we show the photon-ALP conversion probability for PSR J2021+3651 for g aγ = 3.5 × 10 −12 GeV −1 and m a = 4.4 × 10 −9 eV corresponding to the best-fit for this source [9] for regular (black curve) and regular plus turbulent (red curve) magnetic field. We recognize that energy modulations are visible in the range E ∈ [10 2 : 10 3 ] MeV, and strongly enhanced when a turbulent component is also present. In the lower panel, we show the photon spectrum after photon-ALP conversions, i.e.
A gaussian energy resolution of 10% is assumed. It is very clear how the modulations, already visible in the case of regular field, are magnified by the presence of the turbulent component. This result would suggest that the analysis of [9] might deserve a revisitation with the inclusion of the turbulent component of the Galactic Bfield.

IV. CONCLUSIONS
The magnetic field in the Milky Way has been recognized as a very powerful tool to probe the low mass ALP parameter space, since it leads peculiar ALP-induced signatures in observable photon spectra from different (extra)galactic sources. The current literature has adopted state-of-the-art models to characterize the morphology of the regular component of Galactic magnetic field. However, the small-scale turbulent component has been neglected so far, in spite of having a strength comparable to the regular one. In our work, we have investigated this important aspect providing numerical and analytical recipes to characterize the impact of the turbulent component of the Galactic magnetic field on the ALPphoton conversions. Here, we summarize and discuss our main results. First, we have shown that on average the effect of the turbulent magnetic field on photon-ALP conversions grows linearly with the distance from the source at least until the conversion probability remains small. This result is summarized in Eq. (20). Furthermore, we have shown that the effect of the turbulent component are especially relevant around the transition energy between the energy-dependent and the energy-independent regime, where the conversion probability reaches its maximum before saturating. Intriguingly, this transition energy selects a specific range of the ALP mass for which we expect energy-dependent irregularities to be imprinted on the photon energy spectra. We have shown that these irregularities could be observable in the photon spectra associated to different galactic sources, such as red supergiants, supernovae and pulsar, on an energy range between 100 keV and 100 GeV. If such features were to be detected in a future observation, that would represent a direct signature of the turbulent component of the Galactic B-fields. Moreover, it would allow to constraint the ALP mass range, which is typically a very difficult task. Our study might be relevant also in other contexts, where one expects a combination of regular and turbulent magnetic fields. An example is the case of Galaxy Clusters, as recently discussed in Ref. [52].
In conclusion, our results show that neglecting the random component of the Galactic magnetic field is not always justified, even in cases when its correlation length is much smaller than the ALP oscillation length. Moreover, we have shown how the random component leads to recognizable features in the photon spectrum, which could reveal information about the ALP parameters and the magnetic field itself. All these results confirm once more the high physics potential in gamma-ray observations to constrain the ALP parameter space.
We can decompose the transverse B T field in two components, a regular one and a turbulent one: B T = B reg (z) + B(z), with B(z) = 0 where ... denotes an average on all possible field configurations. We suppose that the two transverse components of the turbulent field are totally uncorrelated, B x (z 1 ) B y (z 2 ) = 0, while we can define the two-ponits correlation function Reasonably, the correlation functions for the x and y components are identical. The correlation function is in general a function that goes rapidly to zero for |z 1 − z 2 | l corr . We can conveniently define the correlation length l corr as where B 2 rms = C B (0) is the variance of the turbulent magnetic field.
The simplest model of turbulent field is the "cell model" in which the space is divided into cells with dimension ∼ l cell . In each cell B has a constant random value with zero mean and variance B rms and uncorrelated to the value of the field in adjacent cells. In this case is easy to show that C B (ξ) = B 2 rms · (1 − |ξ|/l cell ) for |ξ| ≤ l cell and zero otherwise. In fact, if |ξ| ≥ l cell the points z and z + ξ always fall in different cells and thus their correlation must vanish. Conversely, when |ξ| < l cell , assuming that the point z falls in a given cell, the correlation is proportional to the probability that the point z + ξ falls in the same cell. Is easy to realize that this probability is just the ratio (l cell − |ξ|)/l cell .
Using the definition in Eq. (A1) we see that the correlation length is half of the cell length, l corr = l cell /2. Although this model is very handy for practical purposes, it is unrealistic since cannot satisfy the condition ∇ · B = 0 on the boundary of cells.
A more realistic model is a Kolmogorov-like power-law spectrum whose three-dimensional Fourier transform of the correlation function is given bŷ The one dimensional correlation function can be obtained as in [53] [Eqs. (A. [9][10][11][12][13][14][15][16]]. For the sake of illustration we tale the regular component constant and parallel to the y-axis without loss of generality, although the results can be easily extended to the case of non constant B reg . We notice also that the conversion probability P aγ is generally very small in the Galaxy due to the relatively short distances travelled by the photon. Within this approximation we can solve Eq. (5) perturbatively.
We start from the solution of Eq. (5) for A x,y (z) with ∆ x = 1 2 g aγ B x , ∆ y = 1 2 g aγ (B reg + B y ), and k = ∆ pl +∆ QED −∆ a . For an initial ALP state we can neglect the back-conversion of photons into ALPs, a(z) a(0) = 1. The ALP to photon conversion probability can be approximatively written as with and T x,y (z) = Averaging on all possible configurations of the turbulent field we have with where we have used the Cauchy formula for nested integrals. For z l corr the upper limit in the integral can be approximated to infinity and thus ∆P (z) has a simple linear form, ∆P (z) A + Bz.
A remarkable simplification can be obtained when the correlation length is much smaller than the oscillation wavelength l corr k −1 . In this case we can approximate the correlation function to a δ function, C B (ξ) 2B 2 rms l corr δ(ξ), and thus valid when ∆P (z) 1 (see also [3] for a similar derivation. The second momentum of the distribution can be calculated by squaring Eq. (A4) and taking the average. Further simplifications can be obtained in the hypothesis that the turbulent component has a gaussian distribution. In this case it can be shown that the n-points correlators are vanishing for n odd and are the sum of all permutations of product of 2-point correlators for n even, e.g.
. For gaussian δ-correlations, after straightforward calculations we have where ∆P is given by Eq. (A9). From this relation we can obtain the standard deviation of the distribution, σ(z) = P 2 aγ (z) − P aγ (z) 2 . In Figure 9 we compare the results obtained through a Monte Carlo simulation with 10 realizations of the magnetic field with the analytic formulas in Eqs. (A7), (A9) and (A11). For the Monte Carlo simulation we consider a cell-like structure with l cell = 20 pc. The regular field is chosen B reg = 3 µG while for each cell the transverse magnetic field has a random direction and a random strength with gaussian distribution and a r.m.s B rms = 1 µG 2 , like in Figure 2. Since the oscillation length for the regular field is l osc 400 pc ( l cell ), we can safely consider the turbulent component as δcorrelated. In the right plot we calculate the average probability and 1σ dispersion using Eqs. (A7), (A9) and (A11). The grey band represents the ±1σ spread around 2 A practical way to obtain two random generated components with gaussian distribution with 0 mean and B 2 rms variance is to use the Box-Muller theorem: B T,x = Brms √ −2 log U cos(2πV ), B T,y = Brms √ −2 log U sin(2πV ) with U , V are random numbers with uniform distribution in the interval ]0; 1[ [54]. the average. We remark that this band does not represent a definite confidence level since the distribution is in general not gaussian, as we will see soon. We notice the good agreement between the numerical simulation and the analytic formulas, despite the limited number of realizations.
In the limit k → 0 (i.e., near the resonance point or if ∆ a , ∆ pl ∆ aγ ) the calculation of T m x,y are straightforward, for example: where (. . .)!! is the double factorial and we made use of the property of the 2j-points correlator for a gaussian field. This allows us to calculate all the moments of p x,y = |A x,y | 2 . For p y we have We can build the moment-generating function for the distribution of p y as where L denotes the Laplace transform, and We can prove the following identity (A16) In fact, for j = 0 is easily verified Eq. (A16) can be proven by induction by deriving both members two times respect to the variable t. FIG. 11: Same as in Figure 10 but for a pure turbulent field (Breg = 0).
For p x the moment-generating function can be obtained in the same way, but with P (0) aγ = 0. Since p x and p y are two independent variables, the probability distribution of P aγ = p x + p y is the convolution of f Px and f Py . Consequently, the generating function for P aγ is just the product of the two generating functions for p x and p y . Using the Bromwich integral to calculate the inverse L-transform we have where is chosen in order to have the singularity s = −1/∆P on the left of the integration path. With the change of variable u = 1 + ∆P s we have with > 0. The integral can be calculated closing the path on a half-circle C R in the negative real half-plane and considering that the integral on C R tends to 0 for R → ∞ due to the Jordan Lemma. The integrand has an essential singularity in u = 0. The function can be expanded in Laurent series about u = 0 The residue of the function is the coefficient of where I 0 is the hyperbolic Bessel function of the first kind of order 0. Making use of the residue theorem we can conclude that In Figs. 10 and 11 we compare the P aγ distributions obtained both with a Monte Carlo simulation and with the Eq. (A24) for a regular field B reg = 1 µG and for a pure turbulent field. The parameters used for the calculation are shown in the caption. In the last case the distribution reduces to a pure exponential, F (P aγ ) ∝ e −Paγ /∆P . Finally, we notice that the previous arguments can be easily extended to the case in which B reg is no longer constant. In this case P (0) aγ is just the oscillation probability obtained integrating Eq. (5) for a pure regular field (we omit the proof for simplicity). where M 0 ( M) denotes the Hamiltonian containing the regular (turbulent) part of the field. However, for δcorrelated gaussian perturbations is it possible to modify Eq. (B1) to calculate the average of the density probability. In fact, let us write M as where ∆ x,y = 1 2 g aγ B x,y and Q a,ij = δ ia δ ja . With this position, it can be shown that the average matrix density ρ(z) satisfies the Redfield equation [55,56]  (B3) with β = 1 4 g 2 aγ B 2 rms l corr . To solve this equation is convenient to transform it an a "Schrödinger-like" form. Writing the matrix ρ ij as a 9-component vector R I = ρ ij , with I = 3i + j − 3 and M 0,IJ = M 0,ik δ lj − M 0,lj δ ik and same for Q x,y , we can rewrite the Redfield equation as which has a simple formal solution (in the case of a constant B reg field) Handy subroutines for calculating exponentials of real or complex matrices can be found in the Expokit package [57]. The P aγ conversion probability can be calculated from the density matrix ρ as P aγ = ρ 11 + ρ 22 with the initial condition ρ ij (0) = δ i3 δ j3 .
Higher order moments of P aγ can be calculated by generalizing Eq. (B4) to tensorial products of the matrix ρ. For example, for the second momentum we consider ρ (2) = ρ ⊗ ρ. Using the 81-component vector R (2) = R ⊗ R, that is R (2) I = R I R J , I = 9I + J − 9 and M (2) 0,IJ = M 0,IJ δ lKL + M 0,KL δ IJ and so on, repeating the same argument used to obtain Eq. (B4), it can be shown that the Redfield equation for the tensor product have the same form of Eq. (B4) and thus the same formal solution of Eq. (B5) (we omit here the proof). The variance of P aγ can be calculated as P 2 aγ = (ρ 11 + ρ 22 ) 2 = ρ 11,11 + ρ 22,22 + 2ρ (2) 11,22 . Although the same argument can be used to calculate any moment of P aγ the complexity of Eq. (B4) grows exponentially.
In Figure 12 we show the ALP-to-photon conversion probability for a regular field B reg = 1 µG (red curve) and with a turbulent field B rms = 10 µG and a correlation length l cell = 100 pc (black curve) together with the ±1σ band calculated with the help of Eq. (B4). The value of the turbulent field and the correlation length as well as the z range chosen for this plot are unrealistic for the Milky-Way and are intended only for illustrative purposes. Notice that the distribution of P aγ is neither gaussian nor symmetric in general, so that the gray band should to be intended just as a qualitative range without a defined confidence level (and in fact sometimes the band comes out of the range [0, 1]). We notice the typical effect induced by stochastic term in the Hamiltonian, i.e., the flavor composition tends to be equally distributed among all the degree of freedom, and thus P aγ → 2/3 for z → ∞. This phenomenon is well known for example in neutrino oscillations in presence of a dissipative term (see, e.g., [58]).