Detecting ALP wiggles at TeV energies

Axions and axion-like-particles (ALPs) are characterised by their two-photon coupling, which entails so-called photon-ALP oscillations as photons propagate through a magnetic field. These oscillations lead to distinctive signatures in the energy spectrum of high-energy photons from astrophysical sources, allowing one to probe the existence of ALPs. In particular, photon-ALP oscillations will induce energy dependent oscillatory features, or “ALP wiggles”, in the photon spectra. We propose to use the discrete power spectrum to search for ALP wiggles and present a model-independent statistical test. By using PKS 2155-304 as an example, we show that the method has the potential to significantly improve the experimental sensitivities for ALP wiggles, and that the ALP wiggles may be detected using the Cherenkov Telescope Array (CTA) for optimistic values of the photon-ALP coupling constant and the magnetic field. Moreover, we discuss how these sensitivities depend on the modelling of the magnetic field. We find that the use of realistic magnetic field models, due to their larger cosmic variance, substantially enhances detection prospects compared to the use of simplified models.


Introduction
Axions are well motivated beyond the standard model particles that can explain a variety of unsolved problems in physics, such as the strong CP problem [1,2] and the nature of dark matter [3][4][5].These particles are mainly characterised by their twophoton coupling g aγ from the interaction term L = 1  4 g aγ aF µν F µν = g aγ aE • B, and by their small mass m a obtained through pion mixing [6,7].The relationship between g aγ and m a is thus fixed as g aγ GeV ∼ 10 −16 m a / µeV up to a O(1) factor [8,9].A more general class of light pseudoscalar particles which share the same two-photon coupling as the axion but have an arbitrary mass m a , is known as axion-like particles (ALPs).Although ALPs do not solve the strong CP problem, they are nevertheless interesting as they, e.g., arise naturally in string theories and other extensions of the standard model [10][11][12].
The majority of ALP searches are based on photon-ALP mixing in a magnetic field (see Refs. [13,14] for two recent reviews): Due to the characteristic two-photon-ALP vertex, a photon/ALP may interact with a virtual photon provided by the magnetic field and convert into an ALP/photon.Currently, the most solid and extensive exclusions at sub-eV masses, g aγ < 6.6 × 10 −11 GeV −1 , are set by the CAST helioscope (m a ≲ eV) [15] by attempting to convert solar ALPs into photons on Earth.A comparable limit is found for m a ≲ keV by studying the lifetime of horizontal branch stars [16,17].The planned "shining light through a wall" experiment ALPS-II [18] and the solar helioscope IAXO [19] are expected to improve upon these limits immensely.Significantly stronger limits around m aγ ∼ 10 −6 eV are obtained for ALP dark matter in haloscope experiments, such as ADMX [20] and the upcoming ABRACADABRA experiment [21], or by conversion near neutron stars [22,23].
The strongest limits (g aγ ≲ 10 −11 -10 −13 GeV −1 ) at low masses (m a ≲ 10 −6 eV) are set by observations of astrophysical photon sources using the signatures that photon-ALP oscillations will imprint on photon spectra: First, photon-ALP oscillations will induce "irregularities".The non-detection of such spectral irregularities has been used to constrain the parameter space using, e.g., gamma-ray observations by HESS [24] and Fermi-LAT [25], observations of the Galactic diffuse gamma-rays by Tibet ASγ and HAWC [26], and using X-ray observations from Chandra [27][28][29].The Cherenkov Telescope Array (CTA) is expected to improve the limits from HESS and Fermi-LAT [30].Second, ALPs that are produced near or in the source can convert back into photons in, e.g., the Galactic magnetic field, thus inducing an additional gamma-or X-ray flux which has been searched for in SN1987A [31], Betelgeuse [32], and super star clusters [33].Moreover, ALPs that are sourced near the polar caps in pulsars and resonantly converted to photons have recently been used to set leading limits [34].Third, photon-ALP oscillations will increase the linear polarisation of photons, which can be used to set limits using e.g.optical photons or X-rays 1 from magnetic white dwarfs and neutron stars [36,37].In Ref. [38], it was shown that the measurement of linear polarisation in white dwarf spectra excludes g aγ ≳ 5.4 × 10 −12 GeV −1 for m a ≲ 3 × 10 −7 eV, which is the strongest existing limit for ALP masses between ∼ 10 −9 and ∼ 10 −6 eV.At masses m a ≲ 10 −11 eV, the best limit is set by the non-detection of spectral irregularities in X-ray data from Chandra [27].Finally, photon-ALP oscillations will effectively increase the mean-free path of photons at TeV energies since ALPs travel practically without any interactions with the extragalactic background light (EBL) [39].This effect has been recently used to set strong limits with HAWC [40].Moreover, this effect is important in combining fit analyses, such as in the recent limit set using FERMI flat radio quasars [41].
All the limits discussed in the previous paragraph are, however, strongly dependent on the treatment of the magnetic fields [42][43][44][45].Therefore, one either needs a reliable description of the magnetic fields, or knowledge of how uncertainties in the magnetic fields affect the results (see e.g. the discussions in Refs.[28,46,47]).This is particularly important for the turbulent component of the magnetic fields, since oversimplified models are often used to describe these fields.
In Ref. [47] we introduced the idea of using the discrete power spectrum to probe photon-ALP oscillations in photon spectra.In this work, we further discuss and exemplify this concept.In particular, we introduce a statistical procedure that has the potential to significantly improve current detection prospects for irregularities induced by photon-ALP oscillations, which we name "ALP wiggles".The statistical method has two main applications: First, it can be used to search for ALP wiggles without specifying the EBL distribution and the magnetic field model.Second, the method is a convenient way to analyse the effect of various magnetic field models on the expected ALP wiggles.We find that this method is more robust than a standard χ 2 comparison with data.
In order to stay as concrete as possible, the examples focus on gamma-rays at TeV energies, of relevance for the upcoming CTA experiment.However, the same considerations and discussions also apply to other energy ranges, relevant for e.g.Fermi, HESS and Chandra.Furthermore, we only include the effect of the extragalactic magnetic field and fix B 0 = 10 −9 G in Sec. 3, while we use B 0 = 5 × 10 −9 G in Sec. 4.Although these are (over-) optimistic choices for a space-filling, primordial magnetic field, the magnetic fields in filaments can easily reach higher values.For example, the turbulent and regular components of the magnetic fields in galaxies and clusters of galaxies are expected to be as large as O(10 −6 ) G. In addition, we fix g aγ = 10 −11 GeV −1 , which at m a ≲ 10 −6 eV is excluded by a factor of ∼ 2, and at m a ≲ 10 −11 eV by more than an order of magnitude.Nevertheless, this choice is appropriate for our purposes: We advocate for a model independent approach to an ALP searches, since current limits based on astrophysical photon-ALP oscillations depend strongly on the modeling of the magnetic fields.Moreover, it suffice to show that the approach is more sensitive than the standard χ 2 search for residuals, which have already been used to set competitive and leading limits [24][25][26][27][28][29].
2 Photon-ALP oscillations in astrophysical magnetic fields

Equation of motion
Physically, one can interpret the photon-ALP oscillation as a mixing between two mass eigenstates, similar to neutrino oscillations.The mixing strength and oscillation length depend on the effective mass of the photon which in turn is determined by the propagation environment (i.e. the surrounding magnetic field, plasma, photon bath, etc.) and photon energy.A photon and an ALP with energy E propagating in z direction can be described by the linearised equation of motion [48], where ϕ = (A ⊥ , A ∥ , a) T is the wave function describing the two photon polarisation and ALP states.The mixing matrix can be written as where ∆ ⊥/∥ = (n ⊥/∥ − 1)E (n being the refractive index of the photon), ∆ a = −m 2 a /(2E) and ∆ aγ = g aγ B T /2.The transverse magnetic field B T is the component of the magnetic field perpendicular to the propagation direction, and the index ⊥ (∥) refers to the direction perpendicular (parallel) to B T .
In this work, we consider only photons in the sensitivity range of CTA (∼ 10 11 -10 14 eV) and low ALP mass (m a ≲ 10 −10 eV).Then the dominant contribution to the photon refractive index is the one from the EBL [47], given by [49] ∆ EBL ≃ ∆ CMB ≃ 0.5 × 10 −42 E. (2.3) All turbulent magnetic fields lead to similar dependencies on the oscillation parameters [47], and the discussions in this paper can therefore be applied to other energy ranges and magnetic field2 .

ELMAG
We simulate the propagation of photons using ELMAG [50,51] which is a Monte Carlo program that simulates electromagnetic cascades of high-energy photons, electrons and positrons created by their interactions with the EBL.We have implemented ALPs into ELMAG [47], thereby allowing for a consistent treatment of cascading and oscillations.This advantage is however at the cost of being significantly more computationally demanding than the alternative Python packages GammaALP [52] and ALPro [28], which are based on transfer matrices.Compared to Ref. [47] we have added the following features ELMAG3 : Gaussian turbulent fields with a broken power-law as power spectrum [see Eq. (2.4)] can be modelled, the magnetic field strength can be distributed as a top-hat function with a given filling factor, and the computation time is significantly reduced.

Magnetic field models
High-energy photons will encounter a variety of turbulent magnetic fields on their path towards Earth, with strengths varying from B ∼ 1 G near jets of AGNe, fields on galactic scales with ∼ µG, within galaxy clusters (∼ 0.1-10 nG) and finally the intergalactic magnetic field, see e.g.Refs.[53,54] for recent reviews.The energy in the turbulent magnetic fields in galaxies and galaxy clusters are believed to be generated at large scales 0.1-10 kpc through, e.g., "mechanical stirring", large-and small-scale dynamos, and compression.The energy is in turn transported to smaller length scales through an energy cascade, leading to a power-law spectrum of the turbulent magnetic field.It is common to assume that the magnetic field either has a Kolmogorov (γ = −5/3) or Kraichnan (γ = −2/3) spectrum.At small k, a Batchelor spectrum (β = 5) is expected, but other spectral indices have been suggested too.
In order to take into account the stochastic nature of the turbulence, we will describe it as a divergence-free Gaussian turbulent field with zero mean and RMSvalue B 2 rms = ⟨B 2 ⟩.Following the approach of Refs.[55,56], we describe the magnetic field as a superposition of n left-and right-circular polarised Fourier modes.The modes will be distributed according to the power-law spectrum between k min and k max .The parameter k 0 determines the break in the power law which is visible in the magnetic field spectra shown in Fig. 1.In the case of astrophysical magnetic fields, L 0 = 2π/k 0 corresponds to the injection scale.The field modes extend down to the dissipation scale L min = 2π/k max which is below any astrophysical scale of interest.In practise, one cuts off therefore the spectrum at a value of L min which is much smaller than the smallest relevant scale of the problem in question 4 .We use k max = 100k 0 , fix k min by the condition B(k min ) = B(k max ) and use 33 modes per decade.The field is normalised such that B 2 rms /2 coincides with the energy density stored in the field.We define the coherence length L c of the turbulent fields as For comparison, we will also consider a simple domain-like field which is often used in the literature due to its simplicity [27,29,57,[57][58][59][60][61].In this approach, the magnetic field is split into patches with a size equal to the coherence length L c .Within each patch, the magnetic field is homogeneous with a randomly chosen direction.This model is unphysical and may lead to a bias in the strength of the ALP signatures deduced [46,47].
As already mentioned, we will focus on the effect of the intergalactic magnetic field.From the non-detection of electromagnetic cascades from blazars, it was concluded that the extragalactic space must be filled with an turbulent magnetic field with a strength of B ≳ 10 −14 G with a large filling factor [62][63][64], while an upper limit of B ≲ 10 −9 G is derived from Faraday rotation measurements [65].The nature and the production of the extragalactic magnetic field remain however unknown: A large range of magnetic field strengths, spectral index at small k (i.e.β) and coherence lengths are made possible by the many conceivable production mechanisms.For example, if produced during inflation, the initial magnetic spectrum will be scale invariant (β = 0).Its coherence length is currently limited by hydrodynamical turbulence decay from below (∼ kpc) and the Hubble radius from above.Meanwhile, the range of allowed magnetic field strengths is slowly closing, and it has been argued that the remaining parameter space can be completely eliminated by the non-detection of magnetic halos from misaligned blazars [66].As a solution, Ref. [66] proposed that the electromagnetic cascades are quenched by plasma instabilities, what, if confirmed, would re-open large parts of the parameter space for intergalactic magnetic fields.
In the remaining of this work, we will focus on the effect of a primordial intergalactic magnetic field with a field strength B ∼ 10 −9 G.This value is chosen to highlight the signatures and the effect of the statistical method introduced in section 3.3.Although this value is arguably over-optimistic for primordial fields, similar field strengths can easily be obtained in filaments between clusters of galaxies, or in Galactic magnetic fields which we for concreteness do not include.Although we consider only TeV photons, all of the discussions and considerations made in this paper can be applied to other energies and astrophysical environments [47], taking into account that the energy dependence of the refractive index scales as E −1 at low energies and as E 1 at high energies.Due to the many uncertainties, and since the expected signal depends strongly on the treatment of the magnetic fields, we will in this work advocate for an experimental approach independent of the modelling of the magnetic fields and the source spectrum.

Parameter space
Photon-ALP oscillations will lead to two important signatures on high-energy photon spectra at E ∼ TeV.First, they will perturb the photon spectrum by energy dependent oscillations with k ∼ ∆ osc [see Eq. (2.6)], even for a turbulent magnetic field [47].Second, the mean free path length of photon will increase since ALPs will propagate without interacting with the EBL.In this work, we focus on the former effect.In this subsection, we will estimate the ALP and magnetic field properties needed to observe ALP wiggles with CTA, and motivate our focus on intergalactic magnetic fields.The conditions discussed here can be deduced graphically from Fig. 3 in Ref. [47].
For a homogeneous magnetic field, the oscillation probability is given by The oscillation length is then defined as L osc = 2π/∆ osc .The oscillatory features-which we name "ALP wiggles"-described by the solution in Eq. (2.6) are present also in turbulent magnetic fields provided that the coherence length is on the same order of magnitude or larger than the oscillation length.At TeV energies, this happens when 2π/∆ CMB ≲ L c , or (2.7) The coherence length of the intergalactic magnetic field is practically unconstrained from above, for concreteness we will use L c ∼ 5 Mpc as default value.Meanwhile, the Galactic magnetic field has a turbulent component which coherence length is usually assumed to be around 20 pc, and the regular component should be comparable to the size of the Galaxy, ∼ 10 kpc.Thus, the turbulent component of the Galactic magnetic field is expected to contribute little to the ALP wiggles at CTA energies.
The ALP wiggles are most prominent around the transition from the strong mixing regime, occurring when ∆ CMB ∼ ∆ aγ or Since CTA is most sensitive in the range between 10 11 and 10 14 eV, one should ideally have 10 11 eV ≲ E crit .This yields (2.9) Furthermore, the ALP mass should be small enough that there exists a strong mixing regime for the given magnetic field strength.This leads to the condition m a ≲ 10 −10 B T /nG.Since the onset of the wiggles is determined by the weakest magnetic field and photon spectra are usually steeply falling, the intergalactic magnetic field may prove to lead to the strongest wiggles.In Eqs.(2.8) and (2.9),only the combination g aγ B T is of importance.This implies that if we change g aγ or B rms , we will change the energy at which the wiggles are most prominent-weaker magnetic fields or lower g aγ implies that we should look at lower energies with a different detector.However, the exact morphology and distribution of magnetic field strengths is unknown.In general, a combination of the Galactic magnetic field (B = O(µG)), magnetic fields around galaxy clusters (B = O(µG)), and the extragalactic magnetic field (B = 10 −9 -10 −14 G) will influence ALP oscillations.Therefore, in the next two sections, we introduce and advocate for a search for ALP wiggles in which the magnetic field does not need to be specified.The example parameters that will be used are chosen such that the wiggles will be prominent in the CTA energy range (g aγ = 10 −11 GeV, B rms = 10 −9 [Sec.3], B rms = 5 × 10 −9 [Sec.4]), and the results should therefore only be considered as a proof of principle.

Statistical tests for ALP wiggles 3.1 The χ 2 test for irregularities
The ALP wiggles induced by photon-ALP oscillations will be perceived as "irregularities" in the photon spectrum.Thus, one can use as a probe the χ 2 test, where f data (E) is the measured binned energy spectrum (with photon-ALP oscillations if they exist), σ data is the experimental uncertainty, N bins is the number of data points, and f (E) is the modelled spectrum (without ALP oscillations) [59].However, even though this method is statistically sound, it can only measure whether the photon spectrum is more irregular than statistically expected.In the simulated examples in this work we 'model' instead the spectrum by fitting the function where β(x) is a fifth order polynomial and b is the spectral index5 to the un-binned spectra by minimising the maximum likelihood estimate (MLE) in order to isolate the effect of the wiggles in a model independent way.

The discrete power spectrum
The photon-ALP oscillations will perturb the photon spectrum by energy dependent oscillations, k ∼ ∆ osc , even for a turbulent magnetic field.At energies above the strong mixing regime, the ALPs with thus lead to wiggles with k ∼ E in the observed photon spectra.Likewise, below the strong mixing regime, the wiggles have the wavenumber k ∼ E −1 .In Ref. [47], we suggested therefore to use the windowed discrete power spectrum, to extract information on the wiggles.The sum in Eq. (3.3) goes over the N detected photon events.Only photons with energies between E min and E max are included, and we use η = E/E min to resemble the expected energy dependence of the wiggles above the strong mixing regime.A similar concept was introduced in Ref. [67].Importantly, one can use the discrete power spectrum to search for ALPs without specifying the magnetic field.However, for a turbulent magnetic field, the ALP signal is a broadened peak whose location and width is a priori unknown.While this makes a detection more challenging, it enables the extraction of information on the magnetic field.Note that the signal strength depend on the choice of E min : It should be chosen close to the transition from the strong mixing regime, which a priori is unknown.This means, on the other hand, that the conditions (2.7) and (2.9) can in principle be used to deduce the ALP parameters from a detected photon-ALP oscillation signal: The combination g aγ B ⊥ will for most astrophysical environments determine the onset of the oscillations ∆ CMB = 2∆ a∥ (see Fig. 3 in Ref. [47]).This means that g aγ B ⊥ can be fixed by finding the value of E min that optimises the observed oscillations.The mass m a can likewise be determined by X-ray measurements.
We consider the test statistic (TS) given by the goodness-of-fit measure compared to an estimated background, Here, G B N (k) and σ B N (k) are the estimated background power spectrum and its 1σ variation (see Sec. 3.3).We choose ∆k = 6 to reduce the contributions from random fluctuations at large k.While Eq. (3.4) shares similarities with the χ 2 statistics, one should emphasise that one expects a longer tail in this test statistics since we are integrating over a range in which there statistically is expected to be random peaks, i.e. the probability that there is a random peak at any k is larger than the probability that there is a peak at a fixed k.The TS can in principle be improved if the shape, position and width is taken into account, using for example machine learning.

Statistical procedure and examples
In this section, the use of the discrete power spectrum to detect ALP wiggles will be exemplified.We will focus on the effect of the magnetic field modelling on the ALP wiggles, thereby illustrating the importance of a proper treatment of the magnetic fields in modelling photon-ALP oscillations.Based on the discussions above, we define the following statistical procedure: 1. Photons are sampled according to the chosen source spectrum using ELMAG.The simulations are stopped when a given number N of photons has reached the detector within the considered energy range.The energy of the simulated photons that reached the detector is used to compute the discrete power spectrum G.
2. The "measured spectrum" is modelled by minimising the maximum-likelihoodestimate (MLE) of the fit function (3.2) to the simulated data.
3. The background power spectrum and its statistical variation is in turn found by drawing N × 10 3 energies using the fitted spectrum as a probability distribution.
In all the scenarios considered in this section, we repeat this procedure for 10 3 realisations of the magnetic field in order to obtain the distribution of TS values.
For concreteness, we will consider N = 10 4 detected photons and in the energy range E ∈ (10 12 , 10 14 ) eV with the injection spectrum dN/dE ∝ E −1.2 .Moreover, we fix g aγ = 10 −11 GeV −1 and use a turbulent magnetic field with B rms = 5 nG.The power spectrum using a Gaussian turbulent field with γ = 5/3, β → ∞ and L c = 5 Mpc (default parameters) is shown in the left pane of Fig. 2. In all plots in this section, the various scenarios will be labelled using the parameters that differs from the default parameters.The results for 50 realisations with photon-ALP oscillations are shown in orange lines, those without photon-ALP oscillations in blue.The averages and the 1σ statistical variance (black lines) were computed using the full set of 10 3 realisations.For comparison, the results using a domain-like field are shown in the right pane of Fig. 2.
The results in Fig. 2 show the power of the statistical procedure: There are clear peaks in the power spectrum including photon-ALP oscillations compared to the case without photon-ALP oscillations 6 .Interestingly, due to the lack of cosmic variance in the simple domain-like field, there is a lack of variance in the photon spectra which represents itself as a clear signal in the discrete power spectrum, even after averaging over many realisations of the magnetic field.This becomes even clearer for larger coherence lengths, as shown in Fig. 3 for L c = 10 Mpc.As such, the use of simplified magnetic field models, such as the domain-like field, may lead to a bias in searches for ALP wiggles and impact the estimated limits on g aγ .However, the larger variance in more realistic magnetic field models-in these examples represented by Gaussian turbulent fields-increases the rate of random encounters of regions of magnetic fields that may enhance the wiggles (see also the discussion in, e.g., Ref. [68]).Thus, a more realistic modelling of the magnetic fields may, in fact, improve the detection prospects by such random encounters.The detection prospects could be further improved choosing more suitable fitting functions.Moreover, a constant windowing function was used.By varying the minimal energy, E min , one may hope to further increase the detection sensitivity For visualisation and to better understand the essence of the method, we plot in Fig. 4 the binned energy spectrum (green errorbars) and the fitted7 spectrum (blue dashed line) for one random realisation of the magnetic field, both with and without photon-ALP oscillations.In addition, the spectrum averaged over all simulations and its standard deviation is shown (orange region).It is clear that photon-ALP oscillations increase the variation in the energy spectrum.The task of the generic fitting function (3.2) is to reduce the effect of unknown features in the source spectrum, such as uncertainties in the modelling of the EBL or unresolved features in the source spectrum.This leads to a caveat of this approach, well visualised in Fig. 4 with the spectrum that yielded the highest TS value in this analysis: The spectrum may be "over-fitted", i.e. part of the signal will be incorporated into the fit function, weaken- ing thereby the signal.This applies especially for the wiggles extending over a larger energy range.Since the true injection spectrum of the source is not known, a detailed modelling of the source would be required in such cases to distinguish between intrinsic and ALP induced features in the energy spectrum.In Fig. 5, we plot the distribution of the TS (3.4), for the default parameters, L c = 1 Mpc, L c = 10 Mpc, and a domain-like field.With the chosen TS, the domainlike field is difficult to distinguish from the non-ALP scenario.The Gaussian turbulent field, however, has a clear tail in the TS distribution, which distinguishes the ALP from the non-ALP scenario.While increasing the coherence length improves the detection prospects, details of the magnetic field like the values of γ and β have only a minor influence on the TS distribution and we therefore do not vary them in the figure.The reason for the weak dependence on these parameters is that the integrated magnetic field distributions, or the filling factor, is independent of the magnetic field spectrum for a Gaussian turbulent field with the same B rms and L c .In order to more clearly quantify the differences, we list in Tab. 1 the probability that a signal is detected with a confidence level of 2σ, denoted as C 95 , and the 99 % quantile for the various magnetic field scenarios considered.Although there are only minor differences in the C 95 value, there are noticeable differences in the tails of the distributions, registered in the 99 % quantiles.
Note that the results of our examples clearly show that the sensitivities depend strongly on whether a Gaussian turbulent or a domain-like field is used: Due to the cosmic variance in more realistic magnetic field models, there is a chance that photons and ALPs propagate through a region of magnetic fields favourable for photon-ALP oscillations, thus enhancing the probability for detection.The same conclusion can be drawn by, e.g., the results in Ref. [43] wherein limits where set using a domain-like magnetic field and using cosmic MHD simulations.Likewise, in Ref. [68], it was found that cosmic MHD simulations have a larger change of such "rare encounters" than a Gaussian turbulent field.Thus, using a more realistic field than the Gaussian magnetic field considered here may improve the sensitivities even further.In this section, we will consider PKS 2155-304 [69] at redshift z = 0.116 as a concrete example.Its photon spectrum can be approximated by [70] with N 0 = 15.4 × 10 −12 cm −2 s −1 MeV −1 , E b = 1136 MeV, α = 1.77 and β = 0.035.
Since CTA [71] is not yet operational and its sensitivities are preliminary, we assume conservatively an energy-independent effective collection area A = 10 5 m 2 .We take, however, into account the energy resolution of the detector by scrambling the detected energies using a normal deviate with an energy dependent half-width given by the preliminary energy resolution of CTA.Furthermore, we consider an energy range E ∈ [10   Histograms of the TS (3.4) obtained using the statistical method described in subsection 3.3.The various colored lines are obtained using different parameters for the magnetic field; the labels indicate the parameter changed compared to the default parameters (see the main text for a description).The results in the no-ALP scenario is plotted as a dashed black line.
independent 8 .The expected number of photons detected by CTA in this energy range from PKS 2155-304 can then be approximated as A being the effective detector area and ∆t the detection time.In Fig. 6, we plot the spectrum obtained from one simulation of PKS 2155-304 for an observation time ∆t = 50 h and a Gaussian turbulent field with L c = 5 Mpc.
To get an idea of the detectability of ALP wiggles from PKS 2155-304, we follow the statistical procedure from section 3.3.The result is shown in Fig. 7 for observation times ∆t = {50, 100, 400} h.As expected, increasing the observation time increases the detection prospects.As a basis for comparison, we consider in Fig. 8 the TS distribution using Eq.(3.1) with the same binning as in Fig. 6.
In order to more clearly quantify the differences between the statistical methods used in Figs.7 and

10
Table 1.The probability that a signal is detected with a confidence level of 2σ, denoted as C 95 , for the various magnetic field scenarios considered with Eq. (3.4) as TS. with a confidence level of 2σ, denoted as C 95 , for the various magnetic field scenarios considered.From these values, we conclude that the use of the discrete power spectrum  leads to better detection prospects compared to a standard irregularity search in the energy spectrum, especially for low observation times.Note, however, that the search for irregularities will depend on the binning, and can thus be improved.On the flip side, the optimal windowing function for the power spectrum will depend on the ALP coupling and the magnetic field strength.There is, in any case, a clear advantage of using the discrete power spectrum compared to a standard search for residual: While a high χ 2 value merely indicates that the data are more irregular than expected, a signal in the discrete power spectrum is a clear indication that the data have wiggles with the same energy dependence as expected for photon-ALP oscillations.In this section, we have considered an optimistic value for the intergalactic magnetic field, B rms = 5 × 10 −9 G, and an excluded coupling strength, g aγ = 10 −11 GeV −1 .It is thus useful to check to what extent the detectability worsens when the magnetic field strength is decreased 9 .Therefore, we plot in Fig. 9 the histograms of the TS obtained for varying magnetic field strength, B rms = {5, 1, 0.5} nG.The corresponding detection probabilities are C 95 = {0.947,0.547, 0.521}.Two effects lead to the quick reduction in C 95 with decreasing magnetic field strength: First, the wiggles are strongest close to the strong mixing regime, and decreasing the magnetic field strength shifts the strong mixing regime to lower energies.From the condition in Eq. (2.9), it it expected that our choice of default parameters leads to the strongest wiggles.This reduction in sensitivity may be partly compensated for by changing the lowest energy considered, for example by considering a different detector.Second, the mixing strength is proportional to the magnetic field strength, which can only be compensated for by increasing the observation time.At lower energies, however, the number of photons is usually significantly higher.

Summary and conclusion
Photon-ALP oscillations will imprint energy-dependent oscillatory features, which we name "ALP wiggles", on photon spectra from distant high-energy sources.We have therefore proposed to use the discrete power spectrum (3.3) to directly probe such wiggles in experimental data.Such a search will be independent of the modelling of magnetic fields and theoretical uncertainties in, e.g., the EBL.This work serves as a first proof of principle, and there is room for improvement: We only considered the simple test statistic (3.4) which only measures the residual of the measured discrete power spectrum compared to the estimated background.Furthermore, in order to stay as concrete as possible, and since the onset of the ALP wiggles at TeV energies is determined by the weakest magnetic field that contribute to the oscillations, we considered only an intergalactic magnetic field.In a complete analysis, one should furthermore consider photon-ALP oscillations in, e.g., the Milky Way and the host galaxy, and in the source itself.As a second step, the discrete power spectrum can be used to extract information about the magnetic field, more specifically, it can be related to the two-point correlation function of the magnetic field [47,72].
We have compared two different treatments of the magnetic field: a Gaussian turbulent field, and a simple and unphysical domain-like field.We found that the increased cosmic variance of the Gaussian turbulent field may significantly improve the detection prospects.Varying the shape of its power spectrum, we did not observe a strong dependence on the resulting axion wiggles, as long as the effective coherence length remained constant.As a concrete example, we considered the detection of ALP wiggles in the energy spectrum from PKS 2155-304 using conservative estimates of the Cherenkov Telescope Array (CTA) sensitivity.Our analysis indicates that ALP wiggles can be detected by CTA for optimistic values of the magnetic field and photon-ALP coupling.Importantly, the method is an improvement compared to a standard search for "irregularities" in photon spectra, which is currently used to set leading limits.The statistical method can furthermore be optimised choosing an appropriate windowing function: Since the extragalactic magnetic field strength currently is only weakly constrained, one cannot know at which energy the first ALP wiggle occur.Moreover, the simple test statistic considered does not take into account the size, shape and location of the peak.

Figure 1 .
Figure 1.Visualisation of different magnetic field spectra that can be modelled in ELMAG.

Figure 2 .
Figure2.The power spectrum with the estimated background subtracted is plotted using a Gaussian turbulent field (left) and a domain-like field (right).The results for 50 realisation of the magnetic field with (orange) and without (blue) photon-ALP oscillations is shown, and the averages and the statistical standard deviations from a sample of 10 3 realisations are shown in black lines.The parameters used in the simulations are discussed in the main text.

Figure 4 .
Figure 4.The simulated data (green errorbars) for one random realisation of the magnetic field are plotted with and without photon-ALP oscillations.The spectra are multiplied by a constant to improve visibility.Furthermore, the function (3.2) fitted to the (un-binned) data is shown as a blue dashed line together with the average obtained from the complete sample of 1000 realisations is shown.In addition, the spectrum from the simulation that yielded the highest TS is shown.

Figure 5 .
Figure 5.Histograms of the TS (3.4) obtained using the statistical method described in subsection 3.3.The various colored lines are obtained using different parameters for the magnetic field; the labels indicate the parameter changed compared to the default parameters (see the main text for a description).The results in the no-ALP scenario is plotted as a dashed black line.

Figure 7 .
Figure 7. Histograms of the TS (3.4) obtained using the statistical method described in subsection 3.3 on the simulated data from PKS 2155-304 for the observation times ∆t = 50 h (blue), 100 h (orange), 200 h (green) and 400 h (red).The corresponding no-ALP cases are shown with dashed lines.

Figure 9 .
Figure 9. Histograms of the TS (3.4) obtained using the statistical method described in subsection 3.3 on the simulated data from PKS 2155-304 for an observation time ∆t = 200 h.The results obtained with magnetic field strengths B rms = 5 nG (blue), B rms = 1 nG (orange) and B rms = 0.5 nG (green) are shown.The no-ALP scenario is shown in dashed line.
12, 10 14 ] eV, for which the energy resolution ∆E/E is approximately energy 8, we list in Tab. 2 the probability that a signal is detected

Table 2 .
The table indicates the probability that ALPs are detected in an observation of PKS 2155-304 using CTA with a confidence level larger than 2σ, denoted as C 95 , using the power spectrum [Eq.(3.4)] and a standard χ 2 search [Eq.(3.1)].The columns corresponds to the detection times used in Fig.7.