{'fraction_non_alphanumeric': 0.028409326210295276, 'fraction_numerical': 0.013749119973495672, 'mean_word_length': 4.411204481792717, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Anisotropy test of the axion-like particle Universe opacity effect: a case for the Cherenkov Telescope Array

The Universe opacity to gamma rays is still an open question, in particular anomalies may have been observed. Assuming that such anomalies find their origin in conventional physics like intrinsic source spectra or the density of the extragalactic background light, they would be evenly distributed over the sky. If they exist, axion-like particles (ALPs) would have a potential effect on the opacity of the Universe to gamma rays, possibly related to the anomalies in the spectral indices of distant gamma-ray sources. In the scenario where ALPs from distant sources convert back to photons in the Galactic magnetic field, their effect on the opacity is expected to depend on the position of the sources. In that case the anomaly is expected to exhibit peculiar correlations on the sky. We propose a method to test the origin of the opacity anomaly, based on angular correlations of spectral softening anomalies. Such a diagnosis requires a wide-field survey of high-energy gamma-ray sources over a broad range of energy. The future Cherenkov Telescope Array (CTA) is perfectly suited to perform such a study. It is shown that while the current sample of sources is not large enough to base conclusions on, with this method CTA will be sensitive to ALP couplings to gamma rays of the order of 3e-11 GeV^-1 for ALP masses below 1e-8 eV.


Introduction
Axion-like particles (ALPs) are hypothetical particles that appear in extensions of the standard model of particle physics [1][2][3]. They are light pseudo-scalar bosons with a two-photon coupling. This coupling leads to a mixing of photon and ALP states when propagation in an external magnetic field is considered [4]. This mixing is used to perform searches for ALPs in various ways, considering ALPs thermally produced in the Sun [5], ALPs produced with Lasers [6], or ALP dark matter [7]. In this article we focus on a possible astrophysical signature of ALPs in the field of γ-ray astronomy. Astrophysical environments do offer magnetic fields over large volumes, which can be in some cases ideal conditions for observing the γ-ray/ALP mixing. In particular, the existence of ALPs could modify the opacity of the Universe to γ rays [8][9][10][11][12][13][14]. This fact raised interest in the last decade after the observation of distant highenergy photon sources having hard spectral indices [15,16]. γ rays are expected to interact with the extragalactic background light (EBL) which fills the Universe in the optical and infrared bands [17,18]. Above some threshold, they are absorbed as they produce e ± pairs. In the Earth frame, this threshold energy for γ rays falls between 100 GeV and a few TeV for sources at redshifts between 0.1 and 0.5. This energy-dependent effect implies a softening of the spectral indices of such distant sources. The observation of hard spectrum sources and distant sources by ground-based imaging atmospheric Cherenkov telescopes (IACTs) such as HEGRA [19], HESS [20], MAGIC [21] and VERITAS [22] is considered in some studies as a possible hint for new physics. One conventional way to solve this puzzle is to revise the EBL models, mainly by lowering the photon density down to values that are close to lower limits set by galaxy counts [23]. Recently, the observation of the spectral features induced by pair production led to measurements of the EBL level at different redshifts by HESS [24] and Fermi [25], showing a possible consistent picture. The solution to the opacity puzzle is still controversial though, as for instance the authors of [26] claim for the persistence of the anomaly even with new EBL models. One possible explanation of the opacity anomaly is provided by the ALP hypothesis [9,14]. In that case the γ rays traveling disguised as ALPs for a fraction of their path, they are less efficient to produce e ± pairs. If the ALPs convert back to photons before they reach the Earth, the high-energy γ rays in excess can give rise to the anomaly. The mixing between photons and ALPs can occur in any magnetic field that exist on the line of sight: around the source itself, in the surrounding structure (e.g. the galaxy cluster), in the intergalactic medium, or in the Milky Way. Still for their tendency to reduce the opacity to γ rays, ALPs have been suggested to be at the origin of a possible mechanism to explain the very high-energy emission of some class of blazars (namely the Flat-Spectrum Radio Quasars, FSRQs) [27,28], although some conventional scenarios might explain it too [29]. For both EBL absorption and absorption in FSRQs, the ALP parameters are invoked in the same region of the parameter space. Concerning the extragalactic opacity, as argued in [10,30] the most natural scenario might be that the γ rays convert into ALPs close to the source, and convert back to γ rays in the Milky-Way magnetic field. This is the scenario we envisage in this paper. Here, we note that if this specific scenario was to explain some opacity anomaly, then the magnitude of the anomaly would depend on the value of the galactic magnetic field (GMF) along specific lines of sight. Conversely, if the anomaly is related to the sources themselves, or to the EBL, there is no reason why it would correlate with the GMF, it would rather be evenly distributed on the sky. In the following, we use the spectral break between the high-energy and the very high-energy bands to quantify the anomaly. Although the sky pattern of the GMF integrated values may depend on the considered GMF model, it turns out that the angular dependence of the anomaly autocorrelation on the sky shows similar trends when ALPs are involved, whatever the model. In every GMF model, if ALP-induced, the anomaly is correlated on small scales on the sky and anti-correlated on large scales. This is the test we propose to perform with the Cherenkov Telescope Array (CTA), which is the only foreseen instrument capable of detecting a sufficient number of sources at both HE and VHE for this test to be conclusive. Indeed as we shall see, the current sample of VHE extragalactic sources is too sparse to be constraining. If positive in the future, i.e. if CTA shows that the anomaly is excessively correlated on small scales and anti-correlated on large scales, such a signal could be interpreted as a strong indication for the existence of ALPs.
The paper is organized as follows. First we recall the basic description of the absorption mechanism on the EBL and the way ALPs can lead to an excess of transparency. In Sec. 3, details are given concerning the scenario we propose to test, in terms of ALP parameters and physical regions where the mixing occurs. Then in Sec. 4, the formal derivation of the test variable is presented together with simulation results showing the expected ALP signal. Sec. 5 presents the results of the proposed test applied to currently available data. In Sec. 6, the method for producing expected samples of extragalactic sources with CTA is described. Then the CTA sensitivity to ALPs using angular correlations is derived in Sec. 7. The results are eventually discussed in Sec. 8.
2 Universe opacity to γ rays and ALPs

Conventional derivation of the attenuation
Gamma rays can interact with cosmic photon backgrounds and produce pairs of electrons and positrons of mass m e . At 1 TeV, the pair production threshold is m 2 e /1 TeV ∼ 0.26 eV, meaning that γ rays are absorbed by pair creation in collisions with photons in the infrared band [17]. This light is the EBL, it originates from stars and interstellar dust emission, and fills the Universe [31]. Because of this pair-creation process, the optical depth τ for γ rays is not zero and their flux is attenuated according to where φ s is the intrinsic source flux at redshift z s and φ obs is the observed flux. The optical depth depends on the cosmological parameters through the distance-redshift relation (z), the EBL density at redshift z and energy n EBL (z, ) and the angle-averaged pair production cross section σ γγ as The cosmological parameters used in (z) are taken from the latest results from Planck for a spatially-flat ΛCDM universe with H 0 = 67.3 km/s/Mpc, Ω m = 0.315 and Ω Λ = 0.685 [32]. The pair-production process results in a cut off in source spectra, similar to an exponential cut off with an additional wiggle. This wiggle is due to the double-bump structure of the EBL spectral energy density in the optical [33,34]. The position of the cut off in the energy spectrum depends on the redshift of the source. An example of such an absorption feature is shown later in Fig. 1 in the case of a source at a redshift of z = 0.1 and an EBL model from [35]. The EBL density and its variation with redshift is very difficult to measure directly because of foreground contamination from the Milky Way, and is subject to large uncertainties. Estimations based on integrated galaxy counts [36][37][38] give lower limits on the density of the EBL. Recently the H.E.S.S. collaboration used γ-ray absorption to perform a measurement of the EBL density [24]. A similar study has been performed by the Fermi-LAT collaboration at higher redshift [25]. Without these measurements, or if one wishes to infer more precise effects, it is necessary to rely on models for the evolution of the EBL density with redshift (for a review, see [18]). As previously mentioned, the observation of sources at redshifts up to ∼0.5 with IACTs challenged this interpretation of the extragalactic absorption [15,16,39]. Although the recent measurements by H.E.S.S. (see also [40]) seem not in conflict with the most recent EBL models, some studies still show evidence for a pair-production anomaly [14,26]. In [26] for instance, the authors claim that data points corresponding to large optical depths are in tension with EBL models. It corresponds to the fact that the VHE part of the extragalactic source spectra are too hard compared to the HE parts. A diagnostic tool for this effect is the difference of spectral indices (absolute values) between HE and VHE ∆Γ = Γ VHE − Γ HE . ∆Γ is therefore a positive number that is too small in case of an anomaly. The conventions regarding the energy bands and the use of this parameter are given in Sec. 4. Because the pair-production process depends on the energy, it could be envisaged that this effect is due to a lower EBL density, or -as τ ∝ c/H 0 -to a Hubble constant that is greater than usually assumed [41,42]. In those two cases, the anomaly is expected to be the same wherever the source is located on the sky. Indeed the Hubble rate is isotropic at the few percent level [43]. Observed from the Earth, the EBL is isotropic at the 15% level on angular scales smaller than 2 • [41,44]. A dipole of the order of 1% is expected [45]. Anyway, an anisotropy of the transparency anomaly would be caused by a variation of the integrated value of the EBL density over the whole line of sight. This is quite unlikely to be sizable as the EBL is expected to be homogeneous over the cosmological distances under consideration here. In any case then, the variation of the anomaly on the sky should not be larger than what is plausible from intrinsic fluctuations from source to source. Another interesting possibility is that ALPs mix with photons in such a way that the Universe happens to be more transparent. As we shall see, in that case specific anomaly patterns can be expected.

How ALPs modify the transparency to γ ray
ALPs are expected to interact with photons through the term where g γa is the ALP-photon coupling strength (in GeV −1 ), F is the electromagnetic tensor, F its dual, and a is the ALP field. It can be shown that in the presence of an external magnetic field, the ALP state mix with the two polarization states of the photon, such that the propagation eigenstates are a mix of ALP and photon [4]. In that case, the propagation of the photon/ALP system can be described by the following equation where A 1 and A 2 are the photon polarization amplitudes, and M is the mixing matrix: where: accounts for the effective mass of the photon, (2.6) accounts for the ALP mass, (2.7) is the γ/ALP mixing term, and (2.8) accounts for the EBL absorption. (2.9) The photon can acquire an effective mass in astrophysical plasmas, such that m 2 γ = 4πα n e /m e , with α being the fine structure constant and n e the electron density. m a is the ALP mass and B t is the value of the strength of the magnetic field projected on the transverse-to-propagation plane. φ is the angle defined by the direction of B t and the first photon polarization component. The i∆ abs term is imaginary, so that the matrix is not Hermitian and the propagation process does not preserve unitarity. The sign of this term in M depends on the sign convention used for the derivative in Eq. 2.4. In the expression 2.9, τ is the optical depth corresponding to Eq. 2.1 and s is the size of the physical region over which the attenuation is computed.
To derive the relevant orders of magnitudes, one can consider a single magnetic domain. In that case, the problem is reduced to a 2-dimensional mixing with a mixing angle θ such that (2.10) In that case the probability for a photon polarized along the transverse component of the magnetic field to convert into ALP after a distance is If θ is small, the mixing is not efficient and photons propagate normally. If θ is close to π/4 though, the mixing is strong and the propagation of photons is affected. Defining the critical energy above which the mixing is strong as [46] the oscillation length can be written (2.13) Current constraints on g γa for low-mass ALPs are set by the CAST experiment [47] and are of the order of 10 −10 GeV −1 . For such couplings, and magnetic fields of a few µG, typical of clusters of galaxies [48], the critical energy lies at the TeV level for ALPs of mass m a ∼ 10 −8 eV. For energies above E c , the strong mixing regime is attained. In the following, only ALP masses such that E c is lower than the observation energy bands are considered. This defines the upper bound of the reachable ALP mass range. In real astrophysical systems, the magnetic field is not homogeneous as considered in the previous example, but rather turbulent. In that case the propagation medium can be described as a collection of patches in which the magnetic field is oriented differently and have different values. The orientations and values of the magnetic field in the patches are related through the turbulence power spectrum, that can be considered for example as a Kolmogorov spectrum. The description of the beam can be done by the means of the density operator [11,49] (2.14) In one domain, the mixing matrix M k includes random values for the parameters φ and B t of Eq. 2.5. After k domains, the composition of the beam is evaluated recursively as where s is the size of the last domain. In the following the initial state is an unpolarized photon beam, described by ρ 0 = diag( 1 /2, 1 /2, 0). The exact configuration of the magnetic field on the line of sight is unknown, so only the statistical properties of the overall transfer function can be predicted. In the case where λ(E) s in each domain, and for a large number N of same-size domains, it can be shown [50] that the average overall transfer function is where P 0 is the transition probability for an unpolarized photon beam in one domain. At the transition between weak and strong mixing regimes, irregularities are expected to appear in the energy spectra of the sources [51]. This effect has been used to set constraints on low-mass ALPs, using H.E.S.S. data [52], and Chandra data [53]. At this point it appears clearly that ALPs modify the propagation of photons in astrophysical media. For ALPs to affect the transparency of the Universe, magnetic fields such that the strong mixing regime is attained must be present on the line of sight. Those include magnetic fields in the source itself, in relation to the acceleration process of the high-energy particles. In case the source lies in a galaxy cluster, the cluster magnetic field can play a role. On larger scales the intergalactic magnetic field (IGMF) can be considered, and at the end of the photon journey, the galactic magnetic field (GMF) is relevant. The case of conversions in the IGMF has been widely studied, e.g. in [11], the authors compute explicitly the mean transfer function and its variance. They show in particular that in general the variance of the transparency effect is large enough to include the conventional non-ALP case. It means that it can happen that the Universe is less transparent with ALPs than without. This is the case for instance if photons are converted to ALPs and do not convert back to photons before detection. It results that the effect can hardly be used to get constraints when IGMF only is considered. In addition, for the transparency to be significantly affected, the value of the IGMF that has to be assumed is of the order of 1 nG. That value is close to the current constraint [54,55] and should be considered as very optimistic. A safer scenario do not involve IGMF but rather magnetic fields close to the source and the GMF. This was proposed in [10] and then studied in more details in [30]. In that case, γ rays originate from the source and are converted into ALPs either in the source itself or in the surrounding galaxy cluster. Then, the ALP part of the beam travels unchanged from the source to the Milky Way, whereas the photons get absorbed on the EBL. When entering the Milky Way, the ALPs convert back into photons in the GMF. This can happen at high energies, where most photons are absorbed on the way, producing a potentially significant boost of the flux at high energy. Figure 1 gives an illustration of the overall transfer function for photons in such a scenario. For that figure, a source of VHE photons in a galaxy cluster located at z = 0.1 is considered. The first conversions are assumed to occur in the galaxy cluster. In that generic case, a typical galaxy cluster magnetic field is considered, with a RMS value of 1 µG and a Kolmogorov power spectrum on scales from 10 kpc to 1 kpc. ALPs then convert back to photons in the GMF. On Fig. 1, the red-dashed line represents the transfer function in the case of EBL absorption only. Other curves are drawn for ALPs of mass m a = 1 neV and coupling g γa = 5×10 −11 GeV −1 . The blue-dashed curve is the mean value of the transfer function in the case of ALPs. One can see in particular that the opacity at high energy tends to be constant, in contrast with the conventional EBL case. The blue-shaded area is the 1-σ envelope for this transfer function, corresponding to different possibilities of the magnetic field configuration on that specific line of sight. The solid black line is the transfer function in one randomly picked realization for the magnetic field. For that single observation, the irregularities appear around E c , and above 1 TeV the transparency becomes larger than in the conventional case. As shown in Fig. 1, the variance of the ALP scenario includes the conventional opacity prediction. It means that it is in principle impossible to derive constraints. Thus a statistical method based on the observation of a sample of sources is required. In the next section, the scenario is further motivated, and more details are given concerning the different conversion regions and the assumptions made in the subsequent analysis.

The tested scenario
The bulk of extragalactic very high-energy γ-ray sources is composed of blazars. These sources are also bright emitters of synchrotron radiation from radio to X rays, witnessing the presence of strong magnetic fields. They are therefore prime targets in the search for ALPs. In this study, γ/ALP oscillations in the source are assumed to be in the strong mixing regime. This assumption is verified if where B is the magnetic field strength, and L is the size of the conversion region [56]. For a value of the coupling strength allowed by current constraints, g γa = 5 × 10 −11 GeV −1 , this condition translates to B × L 60 µG × kpc. As shown in [56] this condition is similar to the Hillas criterion for the acceleration of ultra high energy cosmic rays, implying that the best candidates for the acceleration of UHECRs are also relevant targets to search for ALPs. These sources are typically active galaxy nuclei (AGN). Three different regions containing magnetic fields can be identified in the source. First, the magnetic field in the jets of the AGN, which are believed to host the processes responsible for the γ-ray emission, can be probed using multifrequency radio observations [57,58]. From these observations, a magnetic field strength of a few mG on a scale of 10 pc is typically estimated in the jets [59,60] of AGN. Second, the radio lobes at the extremities of the jets have magnetic fields of a few µG on scales of a few kpc [61]. Third, Fanaroff-Riley I radio galaxies, which are thought to be the parent population of BL-Lac objects, are frequently embedded in a galaxy cluster [62]. Typical magnetic fields in galaxy clusters are at the µG level on scales of a hundred of kpc [48]. Blazars have their jets oriented close to the line of sight of the observer, so that the photons which are produced in these sources will cross all three domains. In each domain, the condition of strong mixing can be verified and globally, if one domain is skipped (the source may not be in a galaxy cluster for instance), the assumption of strong mixing in the source is still valid. In this study, the composition of the beam after propagation in the source is simulated following a random realization of a magnetic field of strength 1µG over a size of 500 kpc and turbulent on a scale of 10 kpc. As long as the strong mixing regime is attained, this simulation correctly describes the statistical distribution of the content of the beam, for any of the above mentioned conversion regions. Some statistical properties of this distribution, the mean and the variance, are formally derived in [11]. The beam of γ rays and ALPs is propagated through the intergalactic medium. The EBL is modeled as in [35]. γ rays and ALPs may also mix toghether in the intergalactic magnetic field (IGMF). The nature and strength of the IGMF still remains unknown so that its modeling is hazardous. Current upper limits on the strength of the IGMF are at the level of 1 nG [54,55], but for a significant mixing to occur over cosmological distances of Gpc scales, a strength of at least 0.1 nG is required following Eq. 3.2 and still assuming g γa = 5 × 10 −11 GeV −1 . This means that a high value of the IGMF strength, close to the current upper limits, is needed to significantly mix photons and axions. Some studies even claim that the IGMF should be lower than 10 −14 G to be in agreement with the high energy spectra of blazars [63]. In the following, a zero IGMF is assumed, and the influence of a strong IGMF on the results is be discussed in Sec. 8.
When entering the Milky Way, at energies corresponding to large τ , the beam is mainly composed of ALPs since the photons have been absorbed on the EBL. During the propagation in the GMF, the ALPs can be back-converted into photons before reaching the Earth. The typical strength of the GMF is a few µG over more than 10 kpc, so that significant γ-ALP oscillations can take place. Here only the regular component of the GMF is of interest. The turbulent component is on average of the same order of magnitude than the regular one but is turbulent on small scales [64,65]. As shown in [66], because of the large number of domains, the conversion does not happen efficiently for small turbulence scales, so that in the case of the GMF, the component turbulent on scales lower than 100 pc can be ignored (see also [30]). In this study, three models for the regular component of the GMF are considered. The model from [67] is based on Faraday rotation measurements from pulsars and extragalactic sources, while the other two models [65,68] also take into account the diffuse polarized synchrotron emission from plasma of relativistic electrons. For each direction on the sky, the conversion probability of ALPs to photons over the line of sight is computed using the full mixing matrix. A sky map of this probability is shown on Fig. 2 for the three models of the GMF. In the regions where this probability is the highest, the process is predicted to be more efficient and in these regions the Universe should appear to be the most transparent, i.e. the anomaly is expected to be there maximal.
In this scenario, the structures of the GMF are correlated with the measured hardening of the TeV spectra. If these structures were known to a high accuracy, it would be possible to use the cross-correlation between the measured spectral hardening and the transverse magnetic field integrated on the same line of sight as a signature of this anomalous transparency scenario. However, as shown on Fig. 2, the structures of the GMF for the three models considered here are not similar. More generally, no consensus is reached on the nature of the different components of the GMF and the resulting back-conversion map is model dependent.
One common feature of these structures is that they show patterns that appear on similarly large angular scales. As a result, the hardening measured from two neighbor directions would appear correlated, whatever the GMF model. Conversely, two sources with large relative angular distance would appear anti-correlated. A model-independent and more robust signature of the scenario is thus provided by the angular autocorrelation of the measured spectral hardening from various sources.

The anisotropy test
The spectral hardening is defined as ∆Γ = Γ VHE − Γ HE where Γ HE and Γ VHE are respectively the spectral indices measured in the high energy (HE, 1 GeV < E < 100 GeV) band and the very-high energy (VHE, E > 100 GeV) band. This observable has been introduced in [40] to study the effect of the EBL absorption in the conventional no-ALP case. The binned spatial autocorrelation of a sample of N sources is computed in the following manner. In each bin of angular distance θ i the autocorrelation C i is computed from the set of the N i couples of sources that are separated by an angular distance falling in the considered bin: where (∆Γ j 1 and ∆Γ j 2 ) are the spectral hardening values for the two sources of the j-th couple, and σ ∆Γ 1 σ ∆Γ 2 are the rms of the ∆Γ distributions. The average spectral hardening over the whole sample is: The theoretical autocorrelation expected from the scenario described in Sec. 3 is shown on Fig. 3 for the three considered GMF models and g γa = 5 × 10 −11 GeV −1 . It is computed by simulating a large number (5000) of sources randomly positioned on the sky and splining the autocorrelation calculated using Eq. 4.1 and 4.2. The level of autocorrelation depends on the GMF model and on g γa but the trends are similar for all models. At small angular distances, the spectral hardening values are correlated because the GMF structures are ordered on small scales. One consequence is that the spectral hardening appears to be slightly anti-correlated at the largest angular scales since the GMF structures will more likely be opposed when looking at scales much larger than their coherence scales. Note that the three models cross zero for quite similar angular distances, between 60 • and 80 • . This is a common feature of the three considered GMF models. The current sample of extragalactic VHE gamma-ray emitters, that will be described in the next section, is far from being in the limit of an infinite number of sources, so binned distributions are used. Moreover, the measurements of the spectral indices in the HE and VHE suffer from some uncertainties. These uncertainties result in a dispersion of the measured spectral hardening and in the end in a dispersion of the autocorrelation. The related uncertainty propagated on C i is then: It has been checked that this procedure correctly implements the errors. In particular, fits of a constant performed over many realizations without ALP signal show that the statistics follow a correct χ 2 distribution.

Results with current data
To date, 37 extragalactic sources are detected by Fermi-LAT and IACTs with measured spectral indices in both bands and with known redshift. Most of those are BL-Lac objects with the exception of 4 radiogalaxies and 3 FSRQs. The list of the sources considered in this study is shown in Tab. 1. The spatial distribution of these sources is shown on the sky map of Fig. 4 in galactic coordinates. The colors stand for the measured break between the spectral indices in the VHE and the HE band. The conversion probability of ALPs to photons in the GMF model of [65] is also shown. The spatial autocorrelation pattern for this set of sources is computed using Fermi-LAT measurements of the spectral index from [105] in the HE band and measurements by IACTs in the VHE band. The measured autocorrelation is shown in Fig. 5. The measurements are compatible with the no-autocorrelation hypothesis with a probability of 0.80. It has been checked that the results do not depend on the choice of the binning. The measurements seem to disfavor the autocorrelation model with an ALP coupling of 5 × 10 −11 GeV −1 and the GMF model of [65] , also shown in Fig. 5. However, the current list of detected sources is still rather limited. This implies a limited sampling of the sky, and an intrinsic randomness in the observation of the ALP effect. It could be for instance that the sources are all located in areas of similar back-conversion probability, leading to an observed transparency effect apparently independent of the position of the sources. This would translate into an autocorrelation pattern compatible with zero. To account for the trial factors for the look-elsewhere effect, simulations are performed where the source positions are scrambled. This allows to estimate the uncertainty that arises from their configuration, for which only one realization is observed. The 1-σ uncertainty on the model from the look-elsewhere effect is shown in   5 as blue boxes. When this effect is taken into account, the observed dataset does not allow to disentangle the ALP transparency model from the no-ALP case. It is found that the probability of the ALP model is 0.7 taking into account the look-elsewhere effect. Given the data points, the probability of the no-autocorrelation hypothesis is 0.8. This does not allow to confirm or reject the ALP hypothesis. Current data cannot exclude the ALP case because the number of detected sources is too small. In the future, the CTA project will enable the discovery of numerous extragalactic sources, thus increasing the sensitivity to the autocorrelation observable.
Γ IACT Γ LAT $ P a→γ # Figure 4. Map of probability of conversion from ALPs to photons in the galactic magnetic field model of [65]. Overlaid is the spatial distribution of sources detected at VHE. The point colors indicate the break between spectral indices at VHE and HE.

Expected sample of extragalactic sources with CTA
CTA is the next-generation array of Cherenkov telescopes. It will consist of two large arrays of tens of telescopes of different sizes, the configuration of which are not yet fully defined [106]. Two CTA sites are foreseen to allow a survey of the whole sky, one in the southern hemisphere and one in the northern one. In the following we assume an even coverage of the whole sky, which can be achieved with a long enough operation of CTA in survey mode [107]. The detectability of AGN with CTA has been thoroughly discussed, for a recent study, see [108]. In [108], it is estimated that at least 500 AGN will be detected by CTA above 30 GeV in ten years in the array B configuration (see Fig. 7 of [108]), which is optimized for low energies, [106]. In the following, it is assumed that 500 extragalactic sources with known redshifts are detected by CTA. In Sec. 8, the consideration of other array configurations is briefly discussed. Simulations of this sample of sources are performed to estimate the sensitivity to the angular autocorrelation observable that could be achieved with CTA. The cumulative distribution function (CDF) for redshifts of the blazars (BL-Lac objects and FSRQs) that should be detected by CTA in the array B configuration is shown in Fig. 7 of [108]. This distribution is obtained by simulating spectral energy distributions of the sources based on a blazar γ-ray luminosity function in agreement with Fermi-LAT results and with conservative assumptions regarding the standard blazar sequence. This prior CDF of the redshift is used to randomly generate the redshifts of the simulated source sample.
The extragalactic sky will mostly be observed by CTA in survey mode. The distribution of the detection significance of the sources in the sample should then closely match the distribution found in [105] for the sample of sources detected by Fermi-LAT which is also operated in survey mode. In this distribution, shown in Fig. 23 of [105], the bulk of detected sources has a significance close to 5 σ, the threshold for a detection. The distribution then rapidly decreases as brighter sources are scarcer. To each source in the simulated CTA sample, a significance is randomly associated with a prior following the Fermi-LAT in the 2 year catalog (2FGL) distribution of source significance. This detection significance is then converted into a number of events in the source region (ON region) using the Eq. 17 of [109] and an estimation of the expected background that would be measured in sample regions around the ON region (OFF regions). The estimation of the background is given in the Monte-Carlo study for CTA of [110]. The number of excess events associated to the source, that are used for the spectrum reconstruction, is eventually obtained by subtracting the number of events in the ON region by the number of events in the OFF regions normalized by the respective region areas.
The spectrum of the excess events is built following a power law with absorption on the EBL and transparency effect due to ALPs. The shape is the following: The first term inside the brackets corresponds to the flux of photons converted into ALPs in the source that are back-converted to photons in the Milky Way before detection. The second term is the flux of photons that are not converted in the source and get absorbed in the intergalactic medium. The slope of the power law Γ is randomly generated following the distributions of spectral indices measured by Fermi-LAT in the 2FGL and shown in Fig. 21 of [105]. The relevance of spectral breaks intrinsic to the source is discussed in Sec. 8. The excess events are distributed following this spectral shape and convolved with the effective area expected for CTA given in [110]. The energy distribution of the excess events is then binned in energy bins of relative size of 10% of the energy. This width is typical of the energy resolution that can be achieved with CTA [110]. The simulated spectrum is eventually obtained by dividing the histogram by the effective area. A re-binning procedure ensures a significance of at least 2 σ per point. A mock spectrum, as it would be observed by CTA, is built through this procedure for each source. Figure 7 shows two examples of would-be measured spectra for sources detected a the 50-σ and 5-σ level, and redshifts of z = 0.1 and z = 0.7 respectively. The simulated sample of sources is then used to estimate the sensitivity of CTA to the spatial autocorrelation of spectral hardening values. P a→γ # Γ VHE Γ HE $ Figure 6. Map of probability of conversion from ALPs to γ rays in the galactic magnetic field model of [65]. Overlaid is the spatial distribution of a simulated sample of 500 sources that could be detected by CTA, assuming g γa = 5 × 10 −11 GeV −1 . The point colors indicate the break between spectral indices at VHE and HE.  Figure 7. Examples of two sources spectra as they would be observed by CTA, for sources located at z = 0.1 and z = 0.7, with significance of detections of 50 σ and 5 σ for the left and right panels respectively. Simulated observations correspond to the ALP scenario with g γa = 5 × 10 −11 GeV −1 .

CTA sensitivity to ALPs using angular correlations
By covering a large energy range, CTA will be able to measure simultaneously the spectral index in the HE band, not affected by EBL absorption, and the spectral index in the VHE band. Here it is assumed that the redshift is known for all considered sources. For each source, the energy separating the two domains is given by the optical depth of absorption on the EBL without ALP-induced transparency effects. The spectral index in the HE band is measured on energy bins for which the optical depth is lower than 0.5. This low value ensures that the spectral index measured in this energy band is not affected by EBL absorption. Conversely, for the VHE band, energy bins with optical depth greater than 1 are used, so that the measured slope will be effectively affected by the absorption on the EBL, and possibly modified by ALP effects. The corresponding energy bands for the two examples of Fig. 7 are shown in the corresponding plots. A sky map of the spectral hardening values measured in one set of 500 simulated CTA source observations is shown in Fig. 6 for the GMF model of [65]. In the example of Fig. 6, ALP effects are considered and simulations are performed assuming g γa = 5 × 10 −11 GeV −1 . A correlation between the spectral hardening and the ALP back-conversion in the Milky Way probability is visible through the color codes. The autocorrelation measured for the simulation of this set of sources is shown in Fig. 8. As in Fig. 5, the blue boxes stand for the theoretical uncertainty from the look-elsewhere effect. Due to the much larger number of sources as compared to the currently detected sample, the uncertainty on the points is reduced, as well as the uncertainty arising from the look-elsewhere effect. When combining these two uncertainties, the probability of the no-autocorrelation hypothesis for the single displayed realization is 8 × 10 −20 , this corresponds to a 8.4-σ deviation. This sample of simulated sources clearly exclude the no-autocorrelation hypothesis, showing that CTA will be sensitive to the autocorrelation observable.  Figure 8. Autocorrelation measured on a simulated sample of 500 sources that could be detected by CTA with g γa = 5 × 10 −11 GeV −1 and the GMF of [65]. The plain black curve shows the ALP prediction for the same coupling strength and GMF model. The blue boxes stand for the uncertainty on the model from the look-elsewhere effect (see text for details). Figures 6 and 8 are shown for a coupling constant of g γa = 5 × 10 −11 GeV −1 . When decreasing g γa , the measured autocorrelation is more compatible with the no-ALP hypothesis. The sensitivity of CTA is computed by estimating the lowest coupling strength that could be excluded assuming that no significant autocorrelation is observed. To estimate this sensitivity, a large number of sample of sources (1000) are simulated with g γa = 0, assuming that no signal is observed. From the corresponding autocorrelation distributions, predicted shapes for the autocorrelation are fitted with increasing g γa until the coupling is excluded at the 95% C.L. For each generated sample, the exclusion set on g γa is different because of the different source positions in the realizations. The averaged exclusion over all the generated sample at the 95% C.L. is g γa < 2.92 × 10 −11 GeV −1 . The variance of the level of exclusion over the whole set of realizations is 2.9 × 10 −12 GeV −1 .
As explained in Sec. 3, when considering µG level magnetic fields, like the magnetic fields in the source and in the Milky Way, the coupling between ALPs and photons is efficient above tens of GeV for ALP masses lower than 10 neV. The corresponding range of ALP parameters accessible with CTA is shown in Fig. 9. On this figure, other constraints from independent analyses and experiments are also shown. For higher masses, the γ-ALP conversion is not in the strong regime and spectral irregularities occur. Current exclusions might be extended by looking at the regularity of the spectrum of bright TeV sources with CTA, which is beyond the scope of the present study. Also shown are various constraints in the same mass range (see Sec. 1 for details). The range of ALP parameters that could explained the opacity anomaly is also shown as "Transparency hint" [14].

Discussion
One possible limitation of the test is that the redshifts of all the sources may not be known. In the 2FGL, two third of the sources have known redshift, which means that one can expect one third of the prospected CTA sample to not have a measurement of the redshift. Knowledge of the redshift of the source is not crucial to the test and it can be worked around if no measurements are available. The redshift of the source is used in the analysis here to define the energy ranges of the HE and VHE domain. For example in Sec. 5, the energy intervals are defined from whether Fermi-LAT or an IACT are used to determine the spectral index, the redshift of the sources are not used at that point. In Sec. 6, the ranges are defined as a function of the optical depth, for which the redshift of the source is required. If the redshift of the source is not known, some other definitions for the HE and VHE ranges can then be used. With CTA, the analysis could be performed with fixed energy intervals. The consequence might be to broaden the ∆Γ distributions, but as long as no spatial correlation are introduced, the autocorrelation observable is still efficient. One can even think of other methods to define these intervals, for instance, the numerical derivative calculated on the spectral points could be used to determine when there is a significant curvature, that would mark the transition between both domains.
It could be that the presence of intrinsic curvature of the source spectra, that would be seen as intrinsic cut offs due to the finite energy resolution, change the results. Such an effect would add to the difference in spectral indices, in a way unrelated to the transparency scenario and could eventually blur the observable. To study the impact of the intrinsic curvature of spectra on the test, the probability density function of intrinsic spectral break between the HE and VHE domains, given by [40] is used. This probability density function is computed by simulating blazar spectral energy distributions (SEDs) following a synchrotron self Compton model with one zone. This model successfully reproduces the SED of BL-Lac objects, which will form the bulk of the extragalactic sources detected by CTA. A break between the HE and VHE spectral indices of each source is randomly added, following this probability density function. Figure 10 shows the autocorrelation prediction computed when spectral breaks are randomly added. The results are compared to the prediction when no break is considered, for g γa = 5 × 10 −11 GeV −1 and the GMF of [65]. The effect is still sizable in spite of the addition of intrinsic spectral break, which is therefore not a limitation to the test. As for the knowledge of the redshifts, the fact that the autocorrelation method still works here is understandable because intrinsic spectral breaks change the distribution of the ∆Γ in just the same manner anywhere on the sky, without strongly modifying the correlations.  Figure 10. Autocorrelation ALP predictions assuming g γa = 5 × 10 −11 GeV −1 and the GMF of [65] with and without intrinsic spectral breaks (see text for details).
In Sec. 7, it is assumed that the strength of the IGMF is low enough so that it could not efficiently mix photons and ALPs. However, in the case of a strong IGMF, close to the upper bound of 1 nG for large turbulence scales of order of 1 Mpc, the mixing is strong enough to influence the test. To study the effect of a strong IGMF, the autocorrelation model is computed for different values of the IGMF strength, assuming a coherence length of 1 Mpc. The result, for g γa = 5 × 10 −11 GeV −1 and the GMF of [65] is shown in Fig. 11. For a large IGMF strength of 1 nG, the autocorrelation pattern is significantly reduced. If the lookelsewhere effect induced by the finite sample of 500 sources expected for CTA is included (see Fig. 8), the autocorrelation is compatible with zero. In that case, the test cannot conclude. For lower values of the IGMF, the sensitivity is retrieved. In Fig. 9, the region of the parameter space that could explain the anomalous spectra claimed in [26] is shown (from [14]). The test proposed in this study covers a significant part of this parameter space, another part being already excluded by constraints from spectral irregularities at TeV energies [52]. Independent studies are required to cover all of this parameter space region. Note that the expected CTA sample of 500 sources used in this study will only be fully complete at the end of operation of CTA. Prospected instruments such as IAXO [111] will also be able to probe this parameter space and perhaps then seal the issue before the end of the full life time of CTA. Nevertheless, the study on the blazar population expected to be detected by CTA [108] shows that in only two years of operation, CTA should detect at least 140 extragalactic sources with known redshifts, thus already more than the triple of the current sample, enabling the use of the test. The CTA simulations used in the previous sections are performed with the array B configuration. In that case the effective area is larger at low energy, and this provides the best prospect in terms of absolute number of detected sources. To estimate the importance of the array configuration, the analysis has been redone completely using the array I. In that case, compared to the B configuration, the effective area is larger at higher energies, at the expense of a larger threshold. As a consequence, less AGN (370) are expected to be detected in that configuration (see Fig. 7 of [108]).
For the sensitivity to the ALP transparency signal through the autocorrelation, the lower number of sources is counterbalanced by a better determination of ∆Γ, in particular thanks to the larger number of collected photons at high energy. Eventually these two antagonist effects compensate, as the sensitivity on the coupling is found to be 3.0 × 10 −11 GeV −1 in the case of the array I, compared to 2.9 × 10 −11 GeV −1 in the case of the array B. Finally, note that when CTA will be running, archival Fermi-LAT data will be available and could be used to improve the analysis, in particular by further lowering the energy threshold and improve the Γ H E measurements.

Conclusions
A new observable related to the ALP-induced anomalous transparency of the universe is proposed. It is based on the measurement of the differences of spectral indices between the HE and the VHE (∆Γ). If the scenario of transparency induced by ALPs is indeed at work, the effect is dependent of the integrated GMF on the line of sight, which offers a possibility to test it through the ∆Γ autocorrelation. The current sample of extragalactic sources detected at very high energies is not large enough for this test to be conclusive. It is shown that CTA, with its expected large sample of 500 extragalactic sources, will be sensitive to the autocorrelation observable. The test is robust with respect to the knowledge of the redshifts, intrinsic spectral breaks and requires that the IGMF is weaker than 1 nG. The detection by CTA of an anomaly autocorrelation signal would be a strong indication for the existence of ALPs. A final test would be to cross-correlate the anomaly pattern with specific GMF models. The proposed observable is complementary to the search for irregularities in the source spectra, and the search for an opacity anomaly based on the magnitude and statistics of the spectral hardening.