Axion-Like Particle Searches with IACTs

The growing interest in axion-like particles (ALPs) stems from the fact that they provide successful theoretical explanations of physics phenomena, from the anomaly of the CP-symmetry conservation in strong interactions to the observation of an unexpectedly large TeV photon flux from astrophysical sources, at distances where the strong absorption by the intergalactic medium should make the signal very dim. In this latter condition, which is the focus of this review, a possible explanation is that TeV photons convert to ALPs in the presence of strong and/or extended magnetic fields, such as those in the core of galaxy clusters or around compact objects, or even those in the intergalactic space. This mixing affects the observed ${\gamma}$-ray spectrum of distant sources, either by signal recovery or the production of irregularities in the spectrum, called"wiggles", according to the specific microscopic realization of the ALP and the ambient magnetic field at the source, and in the Milky Way, where ALPs may be converted back to ${\gamma}$ rays. ALPs are also proposed as candidate particles for the Dark Matter. Imaging Atmospheric Cherenkov telescopes (IACTs) have the potential to detect the imprint of ALPs in the TeV spectrum from several classes of sources. In this contribution, we present the ALP case and review the past decade of searches for ALPs with this class of instruments.


Axion and Axion-Like-Particles
The presentation of a new fundamental particle called 'axion' traces back to the late 1970s, when Peccei and Quinn [1] introduced it as a possible solution to the otherwiseunexplained missing CP-simmetry violation in strong interactions. The term 'axion' was first used by Weinberg [2], who classified it as a "light, long-lived, pseudoscalar boson" together with Wilczek [3]. Since then, axions were subject of strong scrutiny, from both theory and observation; however, half a century later, they remain one of the most compelling solutions to this so-called strong CP problem.
Although there is nothing in the theory forbidding it, and, therefore, it is expected, a violation of the Charge × Parity (CP) symmetry in Quantum Chromo-Dynamics (QCD) was never experimentally observed. The term of the Lagrangian corresponding to CP violation can be written as where θ QCD is a phase parameter of QCD, G is the gluon field strength tensor, a indicates trace summation over the SU(3) colors and g 2 is the QCD coupling constant. θ QCD = 0 in case of no CP violation.
For example, the electrical dipole moment of the neutron, d n , which shows a dependence on the θ QCD angle, and is, therefore, sensitive to the CP violation term, is experimentally bound [4] to be |d n | ≤ 1 × 10 −26 e cm, which translates into θ QCD < 10 −10 , revealing a fine-tuning problem.
The Peccei-Quinn (PQ) mechanism solves the Strong CP Problem by introducing a new global symmetry, known as U(1) PQ symmetry, which makes the CP-violating term (Equation (  1)) in the QCD Lagrangian negligible. Axions are, therefore, pseudo-Nambu-Goldstone bosons associated with the breaking of the U(1) PQ symmetry [2,3].
In the PQ formalism, the axion is a particle of mass m a and decay constant f a , related to the decay amplitude, i.e., to the coupling. In the original model of the axion proposed by Peccei and Quinn [1], Weinberg [2], Wilczek [3], the axion decay constant f a is of the order of the electroweak scale (∼246 GeV), and the mass of the axion m a is inversely proportional to this. Its mass was, therefore, expected to be rather large, i.e., of the order of 100 keV m a 6 × 10 −6 eV 10 12 GeV f a .
Using experimental limits based on the stellar evolution and rare particle decays, this first model was ruled out. Soon after, two new models, abbreviated as KSVZ [5,6] and DFSZ [7,8], emerged. They had in common the fact that the energy scale of the symmetry breaking was instead proposed to be large, i.e., close to the "Grand Unification scale", with the energy of 10 15 GeV. This translated into a very light axion, with mass m a 10 −9 eV. These axions would be very weakly coupled, hence the name currently used to dub them: "invisible axions". Taking into account their mass and coupling, these axions have eluded several experiments to date. Furthermore, a similar but strictly massless pseudoscalar Goldstone particle was also considered, and named arion [9,10].
At present, after many unsuccessful searches for axions (see Figure 1 for a collection of limits), the axion model was extended to a wider group of particles, called Axion-Like Particles (ALPs), in which the decay constant is no longer coupled with the axion mass, in contrast with the original axion (Equation (2)) [11]. ALPs are also often found in SM extensions, motivated by string theory. A real "treasure" for the experimental detectability of ALPs is the term representing the axion coupling to photons through the two-photon vertex, shown in Figure 2. The mentioned term is L a γγ = − g aγγ 4 F µνF µν a = g aγγ E · Ba, where g aγγ is the photon-ALP coupling, F µν the strength tensor of the electromagnetic field, F µν its dual, a is the axion field with mass m a , E is the electric field of a beam photon, and B is the external 1 magnetic field. This effect, explained as the photon-ALP conversion, occurs in magnetic fields and is the basis of many experiments in the search for ALPs. More recently, axions were also proposed as viable Dark Matter (DM) particle candidates. The reason for this relies upon their small mass, combined with a possibly large decay constant f a 10 12 GeV. Since they are connected to spontaneous symmetry breaking, they could have been produced in the early Universe via "misalignment" mechanisms. As such, they could represent a substantial fraction of DM. Arias et al. [11] report that, in order to explain the current amount of DM with ALPs, the axion coupling, dependent on the mass of axions, has to be g aγγ < 10 −12 m a 1 neV Together with hidden photons, axions pose as viable candidates for DM, and are named Very Weakly Interacting Slim Particles (WISPs).

Experimental Searches for ALPs
A wide class of axion searches are performed with special helioscopes, i.e., instruments pointing at the Sun, such as the well-known CERN Axion Solar Telescope (CAST) [12]. Axion helioscopes search for axions produced in the interior of the Sun by the conversion of plasma photons in the Coulomb field of charged particles, the so-called Primakoff process. By creating a strong magnetic field in the instrument and placing an X-ray detector at the far end, these detectors aim to reveal the reconversion of axions into X-ray photons [13]. CAST uses a dipole magnet with a strength of ≈9 T and length L = 9.26 m. The latest constraint on the coupling of photons to axions obtained with CAST [12] is g aγγ < 6.6 × 10 −10 GeV −1 . Progress in this detection technique is expected from the new-generation axion helioscope International Axion Observatory (IAXO) [14]. Methods to constrain solar axions can be obtained using the M'ossbauer [15] and axioelectric effects [16], among others.
Alternative methods are pursued in so-called 'light-shining-through-the-wall (LSW)' experiments, in which photons from a strong laser beam are searched beyond a wall that can be crossed over by ALPs but not photons, asas wdone with The Optical Search for QED Vacuum Bifringence (OSQAR) [17] at CERN. There are also experiments based on the expected axion-induced birefringence of the vacuum, such as ALPS [18]. As proposed by Sikivie [13], another observable phenomenon could be the conversion of axions to photons in a resonant cavity. This study laid the theoretical ground for modern experiments such as Axion Dark Matter eXperiment (ADMX) [19].
Astrophysics searches for axion and ALPs use cosmic magnetic fields and ample photon fluxes present in the cosmos. Clusters of galaxies, for instance, have magnetic fields at their cores that are orders of magnitude larger than the average intracluster medium [20][21][22][23][24]. Magnetic fields in active galactic nuclei or pulsars could also be considered as a possible "medium" for the conversion of photons in axions or ALPs. In the following section, we will focus on astrophysics experiments in the gamma-ray range.

Phenomenology of the Mixing between Gamma-Rays and ALP and Propagation in the Astrophysical Environment
The existence of axions and ALPs can be probed by their imprints on the spectra of astrophysical sources. This is due to the fact that, in the presence of magnetic fields, ALPs couple with photons. Therefore, TeV gamma rays travelling over cosmological distances can oscillate to photons due to the interaction with magnetic fields, and/or convert to ALPs in strong magnetic fields and, as such, cross astrophysical distances until they possibly encounter another strong magnetic field, such as that of the Milky Way, in which they can convert back into observable gamma rays. All these conversion/reconversion processes are governed by a probability term for the mixing P γγ , which depends on the actual ALP mass and coupling, as well as the magnetic field characteristics.

ALP Propagation
In order to understand the phenomenon of conversion, it is necessary to compute the term P γγ . The Lagrangian of the photon-ALP system can be written as where the first term relates to the photon-ALP coupling L aγγ term discussed in Equation (3), followed by terms relating to the effective Euler-Heisenberg Lagrangian L EH for corrections of QED loops in photon propagators due to an external magnetic field [25], where the last term L a describes the kinetic and mass term of the axionic field. To model the propagation, we consider the motion of the ALP in the x 3 direction in a cold and ionized plasma. Generally, for polarized photons and relativistic ALPs, the equations of motion can be written as: where M is the photon-ALP mixing matrix. A 1 (x 3 ) and A 2 (x 3 ) represent the photon linear polarization amplitudes along the x 1 and x 3 axis, respectively, and a(x 3 ) is the axion field strength [26]. The solution to this equation is the transfer function T (x 3 , 0; E) using the condition T (0, 0; E) = 1.
In case a homogeneous magnetic field transverse to the propagation direction (laying in x 2 direction) of the photon beam is assumed, then the photon-ALP mixing matrix M can be simplified into where the elements in this matrix are written considering the plasma condition, the QED vacuum birefringence effect, the axion field, and the photon-ALP mixing. They can be written as In the above equations, α is the fine structure constant and ω pl is the plasma frequency, connected to the ambient thermal electron density, and a critical magnetic field term B CR ∼ 4.4 × 10 13 G is defined. The term ∆ aγ represents the photon-ALP mixing and depends on the strength of the interaction g aγγ and the intensity of the transverse magnetic field B ⊥ . Generally, the magnetic field B does not have to be in the x 2 direction, but at an angle ψ from it. In this case, the equations of motion are solved with a transfer function T (x 3 , 0,

Probability of ALP-Gamma Conversion
With the transfer function, we can compute the probability of the conversion of a gamma ray to an ALP in an external magnetic field. The simplest description of the magnetic field is that of a single domain. In this case, the probability of the photon-ALP mixing can be written as [25] P γ→a = (∆ aγ d) 2 sin 2 (∆ osc d/2) (∆ osc d/2) 2 = sin 2 (2θ) sin 2 ∆ osc d 2 , where θ is the rotation angle θ = 1/2 arcsin(2∆ aγ /∆ osc ), d is the size of the domain and ∆ osc is the oscillation wave number, . This term is often written in terms of a critical energy E crit defined as where ω pl,neV is the plasma frequency in units of neV, B µG is magnetic field in microgauss and g 11 = g aγγ /10 −11 GeV −1 . The critical energy is computed such that, around and above this value, the probability of conversion P γ→a in Equation (8)

Gamma-Ray Survival Probability
We are now in the position to compute the gamma-ray survival probability, that is, the fraction of photon that did not convert to ALP.
To compute it, the exact morphology of the magnetic field should be considered and the hypothesis of having just one single magnetic field domain with a fixed orientation is not plausible. A common approach is to divide it into N different domains. By doing this, the transfer matrix can be reformulated see [28] properly, thus providing the total photon survival probability P γγ When we write Equation (8) following the previously introduced substitutions, we can obtain P γ→a = sin 2 (2θ) sin 2 g aγγ Bd As one can see from Equation (11), P γ→a is dependent on the product of domain length d and magnetic field B. Because of this, it is essential to have a well-defined magnetic field model to account for the oscillations in the spectra of astrophysical objects caused by the photon-ALP mixing.

Astrophysical Magnetic Field and Photon Survival
P γγ depends on the strength of the axion coupling to the photon, the intensity, and the coherence scale of the magnetic field in the medium in which the photon/ALP beam is propagating. While the first term is governed by the microscopic nature of ALP, there are several magnetic field realizations in the universe. Therefore, one needs to consider all different magnetic field environments in the path, from the source to the detector: photon-ALP mixing can be assumed in the magnetic field at the source, in the local environment around the source, in the intergalactic magnetic field and, finally, in the galactic magnetic field [29]. Depending on the observed source, different combinations of magnetic fields can be considered, and as reported by Sikivie [13], there are clearly several ways this problem can be approached. For example, one of the most studied cases is that of the Active Galactic Nuclei (AGN) located in the cores of galaxy clusters. Here, once generated, gamma rays from the AGN would encounter the strong magnetic fields of the cluster core and have a sizeable chance of being converted to ALPs. Such an ALP could travel unimpeded along the intergalactic distances, whose magnetic field is extremely low, thus allowing only a moderate photon reconversion. Finally, the ALP, when entering the Milky Way (MW) magnetic field, could (or could not) be reconverted back to gamma rays.
These are, therefore, several kinds of imprint in the original gamma-ray spectrum. In the first case, if an ample fraction of photons is converted at the source into ALP that do not later convert back in the MW, a signal depletion would be observed. In the second case: if an ample conversion happens in the source but then a back conversion happens in the MW, then one could also observe an ampler signal than expected; for example, if the ALP travelled regions of space that are opaque to gamma-rays (for example, regions with strong particles or radiation fields). One should mention that the above signatures would be observed on top of the well-known gamma-ray extinction due to the interaction with the Extragalactic Background Light (EBL) [27,[30][31][32][33] which strongly limits the observation of TeV emission above redshift z ∼ 1. The propagation of VHE photons is affected by pair production processes with the EBL. Depending on the photon energy, they interact with the extragalactic background photons (EBL) or the cosmic microwave background (CMB), producing an electron-positron pair (γ → e + + e − ). The flux attenuation caused by these processes is dominant for photon energies around E γ ≈ 500 GeV and E γ ≈ 10 6 GeV, respectively [30]. In that way, the greater part of photons is absorbed and evades detection: the universe becomes opaque to VHE gamma rays. The above-mentioned cases of ALP signatures are possible in a regime above the critical energy E crit of Equation (9), where the photon-ALP mixing is maximum. A third case is possible at around E crit . In this regime, the oscillatory behaviour in Equation (11) would create 'wiggles' in the spectrum, in correspondence with the probability term. These wiggles would be hardly misinterpreted as being of astrophysical origin and would, therefore, constitute a clear detection. Such a case is extensively discussed by, e.g., Sánchez-Conde et al. [32], de Angelis et al. [33], Hooper and Serpico [34], de Angelis et al. [35].

A Concrete Example of the Photon Survival Probability
As a showcase, in Figure 3, we report the P γγ calculated for m a = 100 neV and g aγγ = 1 × 10 −11 GeV −1 , assuming conversion in the Perseus galaxy cluster magnetic field and in the Galactic magnetic field. The reason for neglect of the intergalactic magnetic field in this case is its strength being restricted to B ∼ (0.1 − 1) neV, and is still not confirmed. In order to probe the photon-ALP conversions in a magnetic field of this strength, one needs to access critical energy E crit 500 TeV or probe significantly low ALPs masses m a < 10 −10 eV [27]. On the other hand, there are works considering the photon-ALP mixing only in the host galaxy cluster magnetic field and the intergalactic magnetic field [36], while some include all three of the mentioned magnetic field environments [26]. For the magnetic field of galaxy clusters, there are usually not well-established values, with the bounds being between 10 −15 G and 10 −9 G, so their strengths are modeled assuming turbulence and using the Kolmogorov power spectrum. Regarding the magnetic field of the Milky Way, a few models are used most often [37][38][39]. Most of these models are based on the Galactic Synchrotron Emission maps and the extragalactic rotation measurements, modelling a disk field and an extended halo field. In one recent work [40], it is shown that, in ALPs searches using observations of BL Lacertae (often shortened as BL Lac, a well known blazar), there is a sizable jet-mixing effect, meaning that the modelling of the BL Lac jet magnetic field is needed. It is shown that the changes in the parameters of the jet model can cause changes in the photon-ALP mixing in a way that it will enlarge the part of ALP parameter space available for study. It is also shown that, in case of the sources embedded in strong cluster magnetic fields of dense environments, this effect is not relevant, so the constrains set by Abramowski et al. [36], Ajello et al. [41] are still valid. In the future, photon-ALP mixing in the blazar jet might become relevant and, with the new generation of Cherenkov telescopes ( [42], e.g., the Cherenkov Telescope Array (CTA)) the detection of more blazars at a higher redshift is expected. In conclusion, very detailed magnetic field models are needed to address the photon-ALP mixing in a more accurate way.

VHE γ-ray Detection and Analysis Techniques
While there have been several early attempts to detect gamma rays at the ground starting, from the 1950s [43], ground-based gamma-ray astronomy officially started with the detection, in 1989, of the Crab Nebula by the Whipple telescope, which has been operating since 1986 [44]. Increasingly new TeV emitters populated the gamma-ray sky, considering one of the last unexplored windows in the electromagnetic radiation from the cosmos. Whipple belongs to the Imaging atmospheric Cherenkov telescopes (IACTs) class. IACTs are suitable for the detection of VHE γ rays, highly energetic photons which can be produced in the environments of astrophysical objects such as Active Galactic Nuclei (AGNs), supernovae, binary stars, pulsars, etc., as the result of highly accelerated (TeV-PeV) cosmic rays such as electrons and protons. The sensitivity of IACTs is in the range ∼50 GeV-50 TeV. At present, Whipple is decommissioned, and there are currently three major operating IACT arrays: the High Energy Stereoscopic System (H.E.S.S.) [45], the Major Atmospheric Gamma-ray Imaging Cherenkov Telescopes (MAGIC) [46] and the Very Energetic Radiation Imaging Telescope Array System (VERITAS) [47]. IACTs measure the energy and direction of γ-rays indirectly: when the γ-ray penetrates the atmosphere, it interacts with the present nuclei and produces a shower of particles. Charged particles belonging to the shower travel faster than the speed of light in the medium atmosphere, consequently producing Cherenkov light. The faint Cherenkov light is collected by the mirror dishes of the telescopes and reflected into a camera positioned in front of the mirror dish. The energy threshold of IACTs is inversely proportional to the signal-to-noise ratio, so it is convenient to maximize the mirror area and throughput of the optical system to minimize the threshold. The shape of the image shower is described by the so-called Hillas parameters [48]. The flux in γ-rays is calculated using MonteCarlo simulations trained with OFF data, taking the collection area of the telescopes and the effective time of the observations into account. The IACT observations are usually performed in the so-called wobble-mode, to allow for the subtraction of the background during the observations [49]. The analysis of the data for the existing IACTs differs at the high level of analysis, when different methods to correct (unfold) the energy spectrum are used in the respective collaborations. The unfolding methods can be based on different algorithms, in order to assign to the γ-rays a true energy, and to calculate the intrinsic spectrum of a source.
In particular, each array of IACTs possess a different configuration and asset so the instrument response function, used to obtain the final spectra, is different. The principles of detection for IACTs are explained in detail in Section 2.2 of [50]. At present, the collaborations are converging towards common software analysis tools, such as ctools 2 and gammapy 3 . Despite a build-up of successes from the early Crab detection, the technique became really mature in the first decade of this century, when not only were an increasing number of targets acquired, but the results also reached a level of precision and significance never achieved before. As an example, in Reference [51] MAGIC reports the spectrum of the Crab Nebula over three orders of magnitude in energy and four orders of magnitude in intensity, able to detect the source in less than 1 min. Along with this ramp-up of performance, the attention moved from purely astrophysical interests to more fundamental questions, such as the possibility of observing the signature of ALPs in gamma-ray spectra. The first decade of the 21st century brought interest in the imprints and modifications that the conversion of photons to ALPs and vice versa could leave on the spectra of astrophysical objects [28,[32][33][34]52,53].

Astrophysical Targets for ALPs Searches with IACTs
In the attempt to maximize the ALP signatures, it is possible to select the best target of observation. These are astrophysical emitters, where both ample, high-energy gammaray photons fluxes are produced, and where the gamma-ray radiation encounters extended regions with significantly intense magnetic fields, which extend over much larger distances than their coherence length [36]. These conditions guarantee that the probability of interaction is maximal (see Equation (8)). Recently, Abdalla et al. [42] quantified the importance of the intensity of the magnetic field and the source brightness, showing that, for example, a factor of 2.5 more intense magnetic field could result in factor 10 stronger constraints on the ALP coupling ( [42] Figure 7). In the gamma-ray TeV sky, sources often display a flaring state, as opposed to a baseline emission state. If possible, flaring states are then preferred to search for ALP. The best candidates for observation are, therefore, Active Galactic Nuclei (AGNs), where particle acceleration and subsequent gamma-ray emission are found in the region around the central supermassive black holes (SMBHs). AGNs are the largest population of TeV targets. An optimal situation is the AGNs being located in the central core of galaxy clusters, especially in a cool core one, in which extended and intense magnetic fields permeate the region around the central galaxies. In this condition, the magnetic field is not only more intense (tens of µG) with respect to that in the intergalactic space, but also more easily experimentally quantifiable. One of the best examples of this is the AGN NGC 1275 at the center of the Perseus Galaxy Cluster, presented above. Another class of objects of interest for ALP searches is compact objects, namely, pulsars and neutron stars, which are also present in binary systems. Here, the magnetic field is more localized, but significantly more intense. We will come back to discussion of source-specific information later in the text. iI order to make a prediction of the ALP-photon interaction pattern, one has to define both the microscopic nature of the ALP (mass and cross-section) as well as the magnetic field. For the former, one has to build a model of the interaction, as done, for example, in the aforementioned, open-source gammaALPs code, and scan the available parameter space. This is, at present, mostly done with grid sampling. For the magnetic field, since the knowledge is not accurate, the procedure used is normally the computation of several random realizations and attempt of a marginalization procedure of the likelihood over this nuisance parameter, as we will show later.
Other targets have been explored for ALPs searches. In case of a supernova explosion, ALPs would be emitted via the Primakoff process and could be observed with γ rays after a possible re-conversion in the magnetic field of the Milky Way. Following the observation of the supernova SN1987A, constraints due to the non-observation of γ rays, coincidental with the neutrino observations, were set [54,55], but affected by the strong uncertainties. Due to this, Payez et al. [56] revisited these papers using a more detailed analysis. Additionally, neutron stars are another possible candidate for ALP searches. Considering the radiative decays of axions produced by nucleon-nucleon bremsstrahlung in neutron stars, [57,58], Berenji et al. [59] have set constraints on the axion mass m a using the Fermi-LAT data of four neutron stars. This phenomenon was investigated in previous works with X-ray [60] and γ-ray data from a supernova [61].

Critical Energy and Parameter Space for γ-ray Studies
This interest in ALP searches in the γ-ray range was firstly encouraged by the unexplained observation of a change in light polarization in a vacuum filled with a magnetic field detected by the Polarization of the Vacuum with Laser (PVLAS) experiment, Zavattini et al. [62], that offered an explanation based on the existence of a light axion. The results of the PVLAS experiment were in tension with the astrophysical limits. In order to reconcile the signal obtained with PVLAS, authors theorized an ALP with mass m a = 1.3 meV and coupling g aγγ = 3 × 10 −6 Gev −1 . Following on this interpretation, Mirizzi et al. [52] included photon-ALP conversion in the magnetic field of our galaxy and, taking the mentioned parameters into account, present the possible distortions in the photon spectra above the energies E γ ≥ 10 TeV. A few months later, de Angelis et al. [63] and Hooper and Serpico [34] extended this approach. Taking into account the possibility of the photon-ALP conversion in and around the gammaray source, as strong astrophysical accelerators, they showed that the critical energy in Equation (9) falls directly in the gamma-ray range. The photon-ALP conversion then depends on the condition g aγγ B s/2 ≥ 1 (12) where B is the magnetic field component aligned with the photon polarization vector and s is the size of the magnetic field domain. If the photon-axion conversion happens at the source, the product B s in Equation (12) is directly connected to the Hillas criterion [64] for the maximum possible acceleration energy of cosmic rays, and taking into account that cosmic rays with energies up to a few times 10 20 eV have been observed, it follows that sources with B G s pc ≥ 0.3 should exist [34]. Hooper and Serpico [34] showed that IACTs such as H.E.S.S., MAGIC and VERITAS could have probed the range of masses of m a = (10 −9 − 10 −3 ) eV with sensitivities stronger than CAST, as shown in Figure 4. The best candidates for observation were identified with AGNs located in the cores of galaxy clusters. One can now compare Figure 4 with Figure 1 to see how Hooper and Serpico [34] were right in their predictions. The first works by de Angelis et al. [30,63] are based on the unexpected transparency of the universe: EBL observations at the time showed higher transparency at higher redshifts than expected [65,66]. Following the idea that, if converted to ALPs, photons could travel through the extragalactic space without interaction with the EBL or CMB photons, be converted back to photons in the Galactic magnetic field and be detected as such, the photon-ALP conversion could reduce the opacity of the universe to VHE gamma rays, as discussed above. In order to explain the possible detection of TeV photons from a source located at z = 0.44, which was not expected by conventional physics of photon propagation at the time, Sánchez-Conde et al. [32] laid out a similar model. They built a model combining both the mixing near or in the source and mixing in the intergalactic space, stressing the importance of observations, both in the lower and highest energies in order to better constrain the intrinsic spectra of the sources, the EBL attenuation and explore the morphology of the considered magnetic fields. The photon flux attenuation was investigated by varying and combining the photon energy, magnetic field intensity, source redshift and ALPs parameters, showing that these effects could be observed in the spectra of AGNs at the higher energies, E γ ≥ 1 TeV, especially if combined with the Fermi-LAT energy regime [32]. After MAGIC detected the surprising rapidly varying emission from the flat spectrum radio quasar (FSRQ) PKS 1222+216 [67], Tavecchio et al. [68] performed a combined ALPs study using the MAGIC and Fermi-LAT data. The aim of [68] was to present the emission model, including the photon-ALP oscillations mechanism, and explain the mentioned detection. The results showed an agreement with the previously introduced De Angelis, Roncadelli and Mansutti (DARMA) scenario that includes photon-ALP oscillations triggered by large-scale magnetic fields to effectively reduce the EBL attenuation at the energies above 100 GeV [33,35]. These results showed the possibility of explaining such emissions with photon-ALPs oscillations by applying them to the other detected FSRQs.
The challenge related to the detection of spectral features induced by ALPs in the gammaray spectra is due to the the number of statistical and systematics fluctuations that shape the spectrum, even in the case of no ALP effect. First of all, the intrinsic spectrum is shaped by the absorption by the EBL, as discussed above. Such an effect is not-negligible for targets farther than z ∼ 0.1, but many models have been created based on EBL observations in the UV-infrared. Therefore, it is possible to correct the spectra for EBL absorption at different redshifts. The effects of LIV on the flux in photons could also compete with ALPs conversion, but the power of a given source to constrain LIV increases with its distance, its variability in time and the hardness of its energy spectrum, so not all the considered targets are also good targets for studying LIV. The energy reconstruction is generally performed with IACTs at about 10-20% precision, depending on the energy. Finally, the data are affected by a variety of systematics due to the instrument itself (e.g., telescope mirror reflectivity) as well as external factors (atmospheric optical depth). While the former are estimated more accurately, less accurate results were obtained for the latter. Observing irregularities in the spectrum, such as those caused by the ALPs-see Figure 3-is, therefore, challenging.

H.E.S.S. Results with PKS 2155-304
After the first predictions of Hooper and Serpico [34], one of the first attempts to constraint ALP with gamma rays was made by H.E.S.S., using the data from the BL Lac object PKS 2155-304 [36]. In this work, a search for irregularities induced by the photon-ALP mixing in the spectrum was performed, and schematically shown in Figure 5. The problem is in searching for ALP-induced spectral patterns on top of a spectrum generated by the main astrophysical processes at the source. Normally, these generate rather smooth and featureless spectra, such as power-laws, with or without a cutoff or log-parabolic shape. Abramowski et al. [36] assumed a power-law function as a local spectral model, justified by the processes explaining the acceleration and radiation in the extreme astrophysical sources, such as BL Lacs [69]. For an estimation of the irregularities, Wouters and Brun [70] proposed a reduced χ 2 test with the null hypothesis build without the ALP (φ w/oALP ( θ)): where d is the number of degrees of freedom, k runs over the N bins, and φ w/oALP ( θ) is a global fit without ALPs with spectral parameters θ. This method relies on the accuracy of the assumed shape of the spectrum, and is, therefore, subject to possible bias, but can be used in the case when the global fit represents a good estimate on the spectrum [71]. Expanding on this, Abramowski et al. [36] searched for irregularities avoiding a global fit and using only a spectral shape over three adjacent points in the energy spectrum (a triplet i): where (φ i − φ i 2 is the residual of the middle bin in the triplet, φ i the measured flux,φ i the flux in the median bin expected from the power-law fit to the side bins, C i covariance matrix for the triplet and d T i = ∂φ i . Although both methods showed consistent results, Abramowski et al. [36] evaluated that, due to its independence of the global spectral model assumption, the sum of residuals over three adjacent spectral bins is preferred for this kind of analysis. This estimator is calculated for each set of ALPs parameters and 1000 spectra are simulated in order to take the randomness of both the intergalactic magnetic field and the galaxy cluster magnetic field into account. The distribution of values of the spectral irregularity estimator for both the observed spectrum and spectra folded with photon-ALP oscillations for different ALPs parameters are compared, and exclusions of the ALPs parameter space were obtained at 95 % confidence level. The results ( Figure 5, right) yielded constraints on the photon-ALP coupling value g aγγ < 2.1 × 10 −11 GeV −1 for masses of the ALPs m a in the range  neV [36].

Studies on Spectral Irregularities of NGC 1275
The IACT results were completed at lower energies, making use of the Fermi-LAT instrument data. Ajello et al. [41] analyzed 6 years of NGC 1275 data, collected with Fermi-LAT, using the Pass 8 event analysis, and produced ALP predictions by including the photon-ALP conversion in the intracluster magnetic field and in the galactic magnetic field of the Milky Way. A fit of the time-averaged spectrum of NGC 1275 and ALPs models was made, and a likelihood analysis was performed. In Figure 6, one can see the likelihood of one of the event types, together with the best spectral fit with and without ALPs. To evaluate the ALPs hypothesis, Ajello et al. [41] exploited a likelihood ratio test statistics (TS). In the procedure, a time-averaged spectrum is modelled by a smooth function, and likelihood is extracted for each reconstructed energy bin k , L(µ k , θ|D k ), where µ k is the expected number of photons in the photon-ALP conversion scenario, θ are the nuisance parameters of the fit, and D k is the observed photon count. For each set of ALPs parameters and magnetic field, the joint likelihood of all reconstructed energy bins k is maximized and the best-fit parameters are determined. Among the different turbulent magnetic field realizations, simulated by accounting for its randomness, the one corresponding to the 0.95 quantile of the likelihood distribution is chosen. The likelihood ratio test is performed as where the null hypothesis is the no-ALP scenario (including the EBL attenuation) with expected photon count µ 0 and nuisance parametersθ, and the alternative hypothesis of ALP, shows an expected photon count µ 95 and nuisance parametersθ [41]. Aside from the degeneracy of the photon-ALP conversion in coupling and magnetic fields, and non-linearly scaled irregularities considering the ALPs parameters, in comparison with the ALP hypothesis, the null-hypothesis is independent of the realisations of the magnetic field. Considering this, the null distribution needs to be derived from Monte Carlo simulations [41]. The exclusion threshold value, above which the set of ALPs parameters can be excluded with the 95% confidence level statistics, is also calculated from Monte Carlo simulations. The result of this research was the exclusion of the ALP coupling values in the range 0.5 × 10 −11 GeV −1 ≤ g aγγ ≤ ×10 −11 GeV −1 for ALPs masses 0.5 neV ≤ m a ≤ 5 neV and g aγγ ≥ 1 × 10 −11 GeV −1 for 5 neV ≤ m a ≤ 10 neV [41], as seen in Figure 6.

Combined Fermi-LAT and H.E.S.S. Observations of PKS 2155-304
Another study using the Fermi-LAT data from the PKS 2155-304 was carried out by Zhang et al. [73]. The used data were taken from Fermi-LAT observations in the energy range of 100 MeV-500 GeV. Photon-ALP oscillations in the inter-cluster magnetic field and the galactic magnetic field of the Milky Way are included. For different sets of couplings in the range of 10 −12 GeV −1 ≤ g aγγ ≤ ×10 −10 GeV −1 and mass of ALPs 10 −1 neV ≤ m a ≤ 10 2 neV and 800 different realizations of the inter-cluster magnetic field, a binned likelihood analysis similar to [41] was performed. The best fits with and without ALPs were compared to the observed spectrum, and the result is shown in Figure 7. A joint likelihood was calculated; parameter space regions were excluded with 99.9% confidence level and compared with the previous results from H.E.S.S. [36] and with the Fermi-LAT observations of NGC 1275 [41] in Figure 7.

H.E.S.S. Study with Galactic Sources
More recently, H.E.S.S. data of galactic TeV γ-ray sources were used to search for ALP oscillation effects [74]. Ten sources, mainly supernova remnants and pulsar wind nebulae studied by H.E.S.S., were utilized. By using sources in the galactic plane, one can probe the ALPs parameter space with higher ALP mass, m a > 10 −7 neV. This is due to the strength of the galactic magnetic field, an important factor for the photon-ALP oscillation, as seen in Equation (9). The ALP model was obtained by multiplying a spectral fit without ALPs with the P γγ for a certain parameter set (m a , g aγγ ), and including the instrument energy resolution. As above, for each set of parameters (m a , g aγγ ) the ALP model was fitted to the observed spectrum, and a χ 2 value was calculated and compared to the best fit over the whole parameter space. The best parameters were deduced from the calculation of the χ 2 ; however, as the photon-ALP conversion is degenerate in the coupling and magnetic field, and that the induced irregularities are not linearly scaled with the ALPs parameters, a threshold value for excluding the ALPs parameters was derived using the Monte Carlo simulations. In, e.g., in [74] the threshold value is calculated from Monte Carlo simulations and compared to the difference in χ 2 values for each set of ALPs parameters and the best fit over the whole parameter space. Since a scan of the whole parameter space is not feasible [41], it is assumed that the overall shape probability distribution of the alternative hypothesis (with ALPs) can be approximated with the null distribution (no ALPs). It has been shown that such an approach yields conservative limits [41]. The results of Liang et al. [74] were consistent with other limits, but were uniquely sensitive towards the higher mass range. This showed that using galactic observations of TeV sources can improve and further constrain the high-mass part of the ALPs parameter space.
Other studies using Fermi-LAT observations combined with IACTs results have been carried out, using the MAGIC [75] and the H.E.S.S. [76] data. In [75], both signatures induced by the photon-ALP oscillations and step-like flux suppression at the energies E γ > E crit in the spectrum of NCG1275 were investigated. As can be seen, the irregularity estimator in Equation (13) is the reduced-χ 2 . For its general applicability in testing fits to the observed data, and simplicity of calculation, the χ 2 test has been used in several works [74,75]. For each set of the considered ALPs and each magnetic field realization, and photon survival probability is calculated and multiplied by the best fit of the time-averaged spectrum, not including the ALPs effects. χ 2 values for each of these fits are calculated. Testing of the ALPs hypothesis is performed using the ∆χ 2 , defined as ∆χ 2 = χ 2 wALP − χ 2 w/oALP . Based on the distribution of these values for each set of ALPs parameters, the exclusion region is evaluated under specific criteria and ALPs parameters are excluded. Malyshev et al. [75] considered 1000 different random realizations of the cluster magnetic field (modelled as in [41]) for a range of ALP parameters, coupling 10 −14 GeV −1 ≤ g aγγ ≤ 10 −9 GeV −1 and mass of ALPs 10 −2 neV ≤ m a ≤ 10 2 neV. By combining observations of Fermi-LAT, the MAGIC energy range available is extended, and by observing both the patterns of the spectrum, a higher sensitivity to the photon-ALP coupling values is reached, dropping down to g aγγ ∼ 10 −12 GeV −1 . The result was the exclusion of the broader part of ALPs parameter space, compared to the previous analysis of the Fermi-LAT data alone. The excluded region also included the part of the ALPs parameter space which can be assigned to the possible ALP Dm. This showed the potential of combining data obtained by different instruments for the purpose of increasing part of the available ALPs parameter space and increasing the sensitivity. Following previous interest in the effects ALPs oscillations could have on the BL Lac spectra [77,78], recent works investigated the same using the simulations for the upcoming experiments and showed that BL Lac could be used for future studies of the ALPs oscillations [79].

Supernova Remnants
Expanding their previous work, [80], Xia et al. [81] performed a search for spectral irregularities in three galactic supernova remnants, combining GeV data from Fermi-LAT and TeV data from IACTs (H.E.S.S., MAGIC and VERITAS). The broadband spectra were fitted with models with and without photon-ALP conversion in the galactic magnetic field. The ALP hypothesis was tested using the χ 2 analysis and the combined limits were again shown to be inconsistent with limits already set by CAST. The authors speculated that a possible reason for this result could be the uncertain connection between the Fermi-LAT spectrum and the IACT observations, which are not easily calibrated in energy, and also the systematic uncertainties of the instruments that were not taken into account [81]. This approach is likely to be revisited once CTA start taking data.

Studies Obtained Comparing Data from Different Blazars
In [76], the Fermi-LAT and H.E.S.S. data of two BL-Lacs are analyzed. Two different EBL models are also probed. The ALP model included mixing of the inter-cluster magnetic field modeled as a Gaussian turbulent field with zero mean and variance σ B , as in [82], and the Galactic magnetic field [37]. The ALPs hypothesis was evaluated in a similar way, as in [41], using a likelihood ratio test. The results showed the improvement of the fit when ALP models are included and set constraints on the ALPs parameter space consistent with the previously obtained ones.
In [83] the highest energy spectra of AGN studied by Fermi-LAT and IACTs are compared, showing that the inclusion of proton-ALPs oscillation effects improves the agreement of the standard AGN model with the data. Recently, the analysis of Ajello et al. [41] was revisited by Cheng et al. [71] using a different analysis method, calculating the irregularity of the spectrum of NGC 1275. Aiming to measure the irregularity of the spectrum, an estimator needs to be chosen. Looking back to the article by H.E.S.S [36], one could decide to use the estimator from Equation (14). A possible problem arises in the case of a large number of energy bins (∼100) (as in [71]), since the ALPs signatures might become wider than the bin size, making this kind of estimator insensitive to such alternations. Using the energy windows, instead of the spectral points triplets, and following the assumption of a power-law model in those energy windows, Cheng et al. [71] proposed an alternative version of the estimator, where i and j represent the energy window and bin, respectively, while φ pl is the flux assumed by the power-law spectral fit in each energy window, and φ and σ are the measured values of the flux and uncertainty, respectively [71]. Each of the simulated ALPs models were fitted assuming a baseline log parabola. From the assumption that the observed irregularity can be explained by the photon conversion connected to a given set of ALPs parameters, exclusion limits were set. This study included mixing in the intracluster magnetic field and in the galactic magnetic field of the Milky Way. Excluded couplings are g aγγ > 3 × 10 −12 GeV −1 for masses of the ALPs m a ∼ 1 neV at a 95 % confidence level. The results of this search show the possibility of further improving the constraints by combining NGC 1275 observations with observations of another source PKS 2155-304 [71].

ALP-Photon Back Conversion in the Galactic Magnetic Field
As investigated by Long et al. [29], new observations of VHE γ-ray sources could lead to the detection of the flux enhancement due to the ALP-photon back-conversion in the Galactic magnetic field. This enhancement is expected at energies E crit ∼ 100 TeV [29] and could be detected by the Large High-Altitude Air Shower Observatory (LHAASO) [84], CTA, and by the planned Southern Wide-field Gamma-ray Observatory (SWGO) [85]. Long et al. [29] analyzed HE and VHE γ-ray data from three promising AGNs and the spectra were extrapolated to the energies E ∼ 100 TeV. Further on, the assumed intrinsic spectra were folded with the P γγ , assuming the photon-ALP conversions in the source magnetic field and the back-conversion in the magnetic field of the Milky Way. These spectra were compared to the ones obtained only by including the EBL and CMB attenuation. The results showed that, in the respective energy range (above E ∼ 100 TeV), predicted flux enhancement is above one order of magnitude and higher than the sensitivity of the instrument, which will allow for the constraints to be set on the ALPs parameter space [29]. It is also emphasized that, to set stringent constraints, a better estimation of the intrinsic spectra, magnetic fields and EBL attenuation needs to be established, all of which are expected with the upcoming experiments in HE and VHE γ-ray astronomy.
The next generation of IACTs, CTA, is expected to lower the uncertainties in the spectra of astrophysical sources such as active galaxies and BL Lacs and increase the sensitivity to photon-ALP oscillations [42]. In that way, they may surpass the current constraints by broadening the part of the ALPs parameter space available for γ-ray studies, and excluding it even further. A comparable performance is expected from future laboratory experiments: the axion helioscope IAXO [14] and Any Light Particle Search II (ALPSII) [86].

Outlook: The Cherenkov Telescope Array
In the previous section, we discussed the current constraints set by IACTs on ALPs. As noted, to date, several studies using H.E.S.S. and MAGIC data have been performed, but there is still room for improvement with the new upcoming generation of experiments. In particular, CTA is expected to probe the energies up to E γ ≈ 300 TeV, which directly improves the possibility of studying ALPs manifestations. Abdalla et al. [42] created simulations of the observation of the radio galaxy NGC 1275. The magnetic field of the Perseus cluster is modeled following Jansson and Farrar [37], with morphology modeled as a random field with Gaussian turbulence. The conservative value of the central magnetic-field strength was set to 10 µG, along with the other parameters listed in [42]. Using three different sets of ALPs parameters with 100 different magnetic field realizations, P γγ was calculated using the GAMMAALPs code 4 . GAMMAALPs solves the equations of motion for the photon-ALP system using the transfer matrix method. Considering other effects that could impact the photon flux, GAMMAALPs includes the EBL absorption, dissipation in QED, and CMB effects. Observations in both the quiescent and the flaring state are included in a ≈300 h exposure. The authors included the systematic uncertainties of the instrument, and fits are performed both with and without ALPs effects. As the energy binning has a great importance for observing wiggles in the spectrum, three different sets of parameters are used, and fits for each of them are performed by maximizing the likelihood and summing over 40 energy bins. For each set, 100 different magnetic field realizations are computed and likelihood values corresponding to quantile Q = 0.95 of the distributions are chosen. To obtain the confidence intervals of 95% and 99%, Monte Carlo simulations are used.
The results showed that, in contrast to the quiescent state, the flaring state of the source provides a stronger exclusion of the ALPs parameter space, reaching a level where ALPs could constitute the entirety of DM. A probable reason for this is a strong background cut on the quiescent data, which causes the exclusion of the low energy bins from the analysis. On the other side, flaring state observations extend to lower energies. As concluded by the authors, this shows the great importance of observing the high activity states of this and other sources that have yet to be studied. Changes in the magnetic field parameters are also tested and the projected exclusion parameters are presented in Figure 8.
It is important to note that the constraints on the ALPs parameter space are sensitive to changes in the assumed parameters' values in the model of the magnetic field of the Perseus galaxy cluster. Moreover, it is found that finer energy binning gives stronger constraints, as the analysis becomes more sensitive to small and fast oscillations in the spectrum, caused by the photon-ALP oscillations. The projected limits obtained in [42] can be seen in Figure 8. Compared to future laboratory experiments (e.g. ALPSII [86]), CTA exclusions of the ALPs parameter space will be dominated by the systematic uncertainties of the model [42]; CTA is expected to have a similar sensitivity to the planned IAXO [14] and ALPS II experiments [86].

Summary and Conclusions
ALPs are one of the most promising candidates to solvethe strong CP problem, and also a viable solution to the long-lived mystery of astrophysics, the existence of DM. In this review, searches for ALPs are presented, focusing on VHE γ-ray astronomy and IACTs. Current constraints set with IACTs [36] are still viable, and complementing constraints are set by other experiments and instruments, such as axion helioscopes, see, e.g., [12,[17][18][19]41]. Even though a great number of searches have been performed, ALPs parameter space still leaves room for future developments. Probably the most interesting part of yet-unexplored ALPs parameter space is accounting for ALPs which are able to explain and constitute most or all of the DM in the universe. This region is anticipated in the future γ-ray experiments, such as CTA [42] and SWGO [85], or in LHAASO [29], that will be able to explore higher energies of up to about 100 TeV and exclude masses of ALPs of m a ∼ 200 neV [42]. With these and other upcoming laboratory axion experiments, constraints on the ALPs parameter space, or even a possible detection of the ALPs, are increasingly anticipated.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: