Axion-like Particles Implications for High-Energy Astrophysics

We offer a pedagogical introduction to axion-like particles (ALPs) as far as their relevance for high-energy astrophysics is concerned, from a few MeV to 1000 TeV. This review is self-contained, in such a way to be understandable even to non-specialists. Among other things, we discuss two strong hints at a specific ALP that emerge from two very different astrophysical situations. More technical matters are contained in three Appendices.


Introduction
The Standard Model (SM) of strong, weak and electromagnetic interactions based on the gauge group SU(3) C ⊗ SU(2) W ⊗ U(1) Y with three sequential families of quarks and leptons has had wonderful success in explaining all known processes involving elementary particles, and the detection of the Higgs boson with the right properties in 2012 represents its crownig.
However nobody would seriously regard the SM as the ultimate theory of fundamental interactions. Apart from more or less aesthetic reasons, electroweak and strong interactions are not unified, and gravity is not even included. The natural expectation from a final theory would rather be a full unification of all four fundamental interactions at the quantum level. Moreover, the need for an extension of the SM is made compelling by the fact that it cannot account for the observational evidence for non-baryonic dark matter-ultimately responsible for the formation of structures in the Universe-as well as for dark energy, presumably triggering the present accelerated cosmic expansion.
Thus, the SM is presently viewed as the low-energy manifestation of some more fundamental and complete theory of all elementary-particle interactions including gravity. Any specific attempt to accomplish this task is characterized by a set of new particles, along with a specific mass spectrum and their interactions with the standard world. This point will be outlined in detail in Section 2.
Basically, ALPs are very light, pseudo-scalar bosons-denoted by a-which mainly couple to two photons with a strength g aγγ according to the Feynman diagram in Figure 1. Owing to their very low mass m a , they are effectively stable particles (even if they were not, their lifetime would be much longer than the age of the Universe). Additional couplings to fermions and gauge bosons may be present, but they are irrelevant for our forthcoming considerations and will be discarded.
Consider now the Feynman diagram in Figure 1.
A possibility is that one photon is propagating but the other represents an external magnetic field B. In such a situation we have γ → a and a → γ conversions, represented by the Feynman diagram in Figure 2. Note that B could be replaced by an external electric field E, but we will not be interested in this possibility. Because we can combine two diagrams of the latter kind together-as shown in Figure 3this means that γ ↔ a oscillations take place in the presence of an external magnetic field. They are quite similar to flavour oscillations of massive neutrinos, apart from the need of the external magnetic field B in order to compensate for the spin mismatch [30][31][32][33]. Over the last fifteen years or so, ALPs have attracted an ever growing interest, basically for three different reasons.

2.
In another region of the parameter plane (m a , g aγγ )-which can overlap with the previous one-ALPs give rise to very interesting astrophysical effects (for a very incomplete list of references, see . In particular, we shall see that ALPs can be indirectly detected by the new generation of gamma-ray observatories, such as CTA (Cherenkov Telescope Array) [130], HAWC (High-Altitude Water Cherenkov Observatory) [131], GAMMA-400 (High-Altitude Water Cherenkov Observatory) [132], LHAASO (High-Altitude Water Cherenkov Observatory) [133], TAIGA-HiSCORE (Hundred Square km Cosmic Origin Explorer) [134] and HERD (High Energy cosmic-Radiation Detection) [135]. 3.
The last reason is that the region of the parameter plane (m a , g aγγ ) relevant for astrophysical effects can be probed-and ALPs can be directly detected-in the laboratory experiment called shining through the wall within the next few years thanks to the upgrade of ALPS (Any Light Particle Search) II at DESY [136] and by the STAX experiment [137]. Alternatively, these ALPs can be observed by the planned IAXO (International Axion Observatory) observatory [138], as well as with other strategies developed by Avignone and collaborators [139][140][141]. Moreover, if the bulk of the dark matter is made of ALPs they can also be detected by the planned experiment ABRACADABRA (A Broadband/Resonant Approach to Cosmic Axion Detection with an Amplifying B-field Ring Apparatus) [142].
Our aim is to offer a pedagogical and self-contained account of the most important implications of ALPs for high-energy astrophysics.
The paper is structured as follows. Section 2 contains a brief outline of what can be considered as the standard view about the relation of the SM with the ultimate theory.
Section 3 describes the most important properties of ALPs, which will be used in the subsequent discussions, along with most of the bounds on m a and g aγγ .
The aim of Section 4 is to provide all the necessary astrophysical background needed to understand the rest of the paper, which is therefore fully self-contained.
Section 5 describes the propagation of a photon beam in the γ-ray band-emitted by a far-away source with redshift z-in extragalactic space and observed on Earth at energy E 0 . Extragalactic space is supposed to be magnetized, and so γ ↔ a oscillations should take place in the beam as it propagates. Moreover, the extragalactic magnetic field B is modeled as a domain-like structure-which mimics the real physical situation-with the strength of B nearly equal in all domains, the size L dom of all domains being taken equal and set by the B coherence length while the direction of B jumps randomly and abruptly from one domain to the next. Because of the latter fact, this model is called domain-like sharp edge model (DLSHE), and is adequate for beam energies currently detected (up to a few TeV) since the γ ↔ a oscillation length L osc is very much larger than L dom . It is just during propagation that an important effect of the presence of ALPs comes about. What happens is that the very-high-energy (VHE, 100 GeV < E 0 < 100 TeV) beam photons scatter off background infrared/optical/ultraviolet photons, which are nothing but the light emitted by stars during the whole evolution of the Universe, called extragalactic background light (EBL). Owing to the Breit-Wheeler process γ + γ → e + + e − [143], the beam undergoes a frequency-dependent attenuation 1 . However in the presence of γ ↔ a oscillations photons acquire a 'split personality': for some time they behave as a true photon-thereby undergoing EBL absorption-but for some time they behave as ALPs, which are totally unaffected by the EBL and propagate freely. Therefore, now the optical depth τ ALP γ (E 0 , z) is smaller than in conventional physics. However since the corresponding photon survival probability is P ALP γ→γ (E 0 , z) = exp −τ ALP γ (E 0 , z) , even a small decrease in τ ALP γ (E 0 , z) gives rise to a much larger photon survival probability as compared to the case of conventional physics. This is the crux of the argument, first realized in 2007 by De Angelis, Roncadelli and Mansutti [52].
Section 6 addresses the so-called VHE BL LAC spectral anomaly. Basically, even if the EBL absorption is considerably reduced in the presence of ALPs, it nevertheless produces a frequency-dependent dimming of the source. Owing to the Breit-Wheeler process γ + γ → e + + e − , the beam undergoes a frequency-dependent attenuation. Thus, if we want to know the emitted spectrum Φ em E 0 (1 + z) we have to EBL-deabsorb the observed one Φ obs (E 0 , z). Correspondingly, we obtain the emitted spectrum which slightly departs from a single power law, hence for simplicity we perform the best-fit Φ em E 0 (1 + z) ∝ E 0 (1 + z) −Γ em (z) . When we apply this procedure to a sufficiently rich homogeneous sample of VHE sources and perform a statistical analysis of the set of the emitted spectral slopes {Γ em (z)} for all considered sources, in the absence of ALPs we find that the best-fit regression line is a concave parabola in the Γ em − z plane. As a consequence, there is a statistical correlation between Γ em (z) and z. However how can the sources get to know their z so as to tune their Γ em (z) in such a way to reproduce the above statistical correlation? At first sight, one could imagine that this arises from selection effects, but this possibility has been excluded, whence the anomaly in question. So far, only conventional physics has been used. Nevertheless, it has been shown that by putting ALPs into the game with m a = O 10 −10 eV and g aγγ = O 10 −11 GeV −1 such an anomaly disappears altogether! Section 7 discusses a vexing question. A particular kind of active galactic nucleus (AGN)named Flat Spectrum Radio Quasars (FSRQs)-should not emit in the γ-ray band above 30 GeV according to conventional physics, but several FSRQs have been detected up to 400 GeV! It is shown that in the presence of γ ↔ a oscillation not only FSRQs emit up to 400 GeV, but in addition their spectral energy distribution comes out to be in perfect agreement with observations! Basically, what happens is that the above mechanism that increases the photon survival probability in extragalactic space works equally well inside FSRQs.
Section 8 is superficially similar to Section 5, but with a big difference. In 2015 Dobrynina, Kartavtsev and Raffelt [145] realized that at energies E 15 TeV photon dispersion on the CMB (Cosmic Microwave Background) becomes the leading effect, which causes the γ ↔ a oscillation length to get smaller and smaller as E further increases. Therefore, things change drastically whenever L osc L dom , because in this case a whole oscillation-or even several oscillations-probe a whole domain, and if it is described unphysically like in the DLSHE model then the results also come out as unphysical. The simplest way out of this problem is to smooth out the edges in such a way that the change of the B direction becomes continuous across the domain edges. After a short description of the new model, the propagation of a VHE photon beam from a far-away source is described in detail.
Section 9 presents a full scenario, which also includes the γ → a conversions also in the source. Actually, the VHE photon/ALP beam emitted by the considered sources crosses a variety of magnetic field structures in very different astrophysical environments where γ ↔ a oscillations occur: inside the BL Lac jet, within the host galaxy, in extragalactic space, and finally inside the Milky Way. For three specific sources a better agreement with observations is achieved as compared to conventional physics. Section 10 addresses another characteristic effect brought about by photon-ALP interaction, namely the change in the polarization state of photons. Less attention has been paid so far to this effect in the literature. Very recent and interesting results on this topic are discussed.
Finally, in Section 11 we draw our conclusions.

The Standard Lore in Particle Physics
As already stressed, the Standard Model (SM) of particle physics is presently regarded as the low-energy manifestation of some more fundamental theory (FT) characterized by a very large energy scale, much closer to the Planck mass M 10 19 GeV than to the Fermi scale G −1/2 F 250 GeV.

General Framework
In order to be somewhat specific, let us suppose that the FT is a string theory. As is well known, in order for such a theory to be mathematically consistent it has to be formulated in ten dimensions. As a consequence, six dimensions must be compactified. Unfortunately, the number of different compactification patterns is O 10 105 , and so far no criterion has been found to decide which one uniquely leads to the SM at low-energy: many of them appear to be viable. However each one predicts a new physics beyond the SM. Nevertheless, one thing is clear. Generally speaking, every compactification pattern leads to a certain number of ALPs as pseudo-Goldstone bosons, namely Goldstone bosons with a tiny mass due to string non-perturbative instanton effects. Basically, they arise from the topological complexity of the extra-dimensional manifold, and their properties depend on the specific compactification pattern. While for a Calabi-Yau manifold one expects O(100) ALPs, there are manifolds which give rise to the so-called Axiverse, a plenitude af ALPs with one per decade with mass down to 10 −33 eV [16]. Below, we do not commit ourselves with any specific string theory.
After compactification at the very high scale Λ, the resulting four-dimensional field theory is described by a Lagrangian. We collectively denote by φ the SM particles together with possibly new undetected particles with mass smaller than G −1/2 F -the above ALPs are an example-while all particles much heavier than G −1/2 F are collectively represented by Φ. Correspondingly, the Lagrangian in question has the form L FT (φ, Φ), and the generating functional for the corresponding Green's functions reads where J and K are external sources and N is a normalization constant. The resulting low-energy effective theory then emerges by integrating out the heavy particles in Z FT [J, K], so that the low-energy effective Lagrangian L eff (φ) is defined by Evidently, the SM Lagrangian L SM (φ SM ) is contained in L eff (φ), and-in the absence of any new physics below G −1/2 F -it will differ from L eff (φ) only by non-renormalizable terms involving the φ SM particles alone, that are suppressed by inverse powers of Λ.
In any theory with a sufficiently rich gauge structure-which is certainly the case of L FT (φ, Φ)-some further ALPs can arise. Indeed, some global symmetries G invariably show up as an accidental consequence of gauge invariance. Since the Higgs fields which spontaneously break gauge symmetries often carry nontrivial global quantum numbers, it follows that the group G undergoes spontaneous symmetry breaking as well. As a consequence, some Goldstone bosons-which are collectively denoted by a if G is non-abelian-are expected to appear in the physical spectrum and their interactions are described by the low-energy effective Lagrangian, in spite of the fact that G is an invariance group of the FT. Because the instanton-like effect explicitly break G by a tiny amount, additional ALPs show up. Moreover, ALPs can arise simply because the effective low-energy theory does not respect the symmetry which gives rise to the Goldstone bosons. Therefore, by splitting up the set φ into the set of SM particles φ SM plus the ALPs a, the low-energy effective Lagrangian has the structure L eff (φ SM , a) = L SM (φ SM ) + L nonren (φ SM ) + L ren (a) + L ren (φ SM , a) + L nonren (φ SM , a) , (3) where L ren (a) contains the kinetic and mass terms of a, and L ren (φ SM , a) stands for renormalizable soft-breaking terms that can be present whenever G is not an automatic symmetry of the low-energy effective theory (we recall that automatic symmetry means any global symmetry which is implied by the gauge group, such as lepton and baryon numbers in the SM). Clearly, the SM is contained in L SM (φ SM ).
Needless to say, it can well happen that between G −1/2 F and Λ other relevant mass scales {Λ i } i=1,2,3,... exist. In such a situation the above scheme remains true, but then G may be spontaneously broken stepwise at various intermediate scales.
Finally, we recall that pseudo-Goldstone bosons such as ALPs are necessarily pseudoscalar particles [146].

Axion as a Prototype
A characteristic feature of the SM is that non-perturbative effects produce the term ∆L θ = θ g 2 3 G µν aGaµν / 32π 2 in the QCD Lagrangian, where θ is an angle, g 3 and G µν a are the gauge coupling constant and the gauge field strength of SU C (3), respectively, andG µν a ≡ 1 2 µνρσ G aρσ . All values of θ are allowed and theoretically on the same footing, but a nonvanishing θ values produce a P and CP violation in the strong sector of the SM. An additional source of CP violation comes from the chiral transformation needed to bring the quark mass matrix M q into diagonal form, and so the total strong CP violation is parametrized byθ = θ + arg Det M q (for a review, see e.g., [147][148][149][150][151]). Observationally, a nonvanishingθ would show up in an electric dipole moment d n for the neutron. Consistency with the experimental upper bound |d n | < 3 · 10 −26 e cm requires |θ| < 10 −9 [147][148][149][150][151]. Thus, the question arises as to why |θ| is so unexpectedly small. A natural way out of this fine-tuning problem-which is the strong CP problemwas proposed in 1977 by Peccei and Quinn [152,153]. Basically, the idea is to make the SM Lagrangian invariant under an additional global U(1) PQ symmetry in such a way that the ∆L θ term can be rotated away. While this strategy can be successfully implemented, it turns out that the U(1) PQ is spontaneously broken, and so a Goldstone boson is necessarily present in the physical spectrum. Actually, things are slightly more complicated, because U(1) PQ is also explicitly broken by the same non-perturbative effects which give rise to ∆L θ . Therefore, the would-be Goldstone boson becomes a pseudo-Goldstone boson-the original axion [154,155], which was missed by Peccei and Quinn-with nonvanishing mass given by m 0. 6 10 7 GeV f a eV , (4) where f a denotes the scale at which U(1) PQ is spontaneously broken. We stress that Equation (4) has general validity, regardless of the value of f a . Qualitatively, the axion is quite similar to the pion and it has Yukawa couplings to quarks which go like the inverse of f a . Moreover-just like for the pion-a two-photon coupling aγγ of the axion a is generated at one-loop via the triangle graph with internal fermion lines, which is described by the effective Lagrangian µνρσ F ρσ . The two-photon coupling g aγγ entering Equation (5) has dimension of the inverse of an energy and is given by with k a model-dependent parameter of order one [156]. Note that g aγγ ∝ 1/ f a and it turns out to be independent of the mass of the fermions running in the loop. Hence, the axion is characterized by a strict relation between its mass and two-photon coupling m = 0.7 k g aγγ 10 10 GeV eV , In the original proposal [152,153], U(1) PQ is spontaneously broken by two Higgs doublets which also spontaneously break SU W (2) ⊗ U Y (1), so that f a G −1/2 F . Correspondingly, from Equation (4) we get m 24 keV. In addition, the axion is rather strongly coupled to quarks and induces observable nuclear de-excitation effects [157]. In fact, it was soon realized that the original axion was experimentally ruled out [158].
A slight change in perspective led shortly thereafter to the resurrection of the axion strategy. Conflict with experiment arises because the original axion is too strongly coupled and too massive. However, given the fact that both m and all axion couplings go like the inverse of f a the axion becomes weakly coupled and sufficiently light provided that one chooses . This is straightforwardly achieved by performing the spontaneous breakdown of U(1) PQ with a Higgs field Φ which is a singlet under SU W (2) ⊗ U Y (1) [159,160], and with a vacuum expectation value Φ G −1/2 F . We have told the axion story in some detail because the latter condition leads to the conclusion that the U(1) PQ symmetry has nothing to do with the low-energy effective theory to which the axion belongs-namely L nonren (φ SM , a) in Equation (3)-but rather arises within an underlying more fundamental theory (alternative strategies to make the axion observationally viable were developed in [161][162][163]).
Thus, we see that the axion strategy provides a particular realization of the general scenario outlined in Section 2.1, with G = U(1) PQ , f a = Λ i (for some i) and L nonren (φ SM , a) including L aγγ among other terms involving the SM fermions and gauge bosons. This fact also entails that new physics should lurk around the scale f a at which U(1) PQ is spontaneously broken. Incidentally, the same conclusion is reached from the recognition that the Peccei-Quinn symmetry is dramatically unstable against any tiny perturbation-even for f a at the Planck scale-unless it is protected by some discrete gauge symmetry which can only arise in a more fundamental theory [164][165][166][167].
Finally, we remark that the first strategy to detect the axion was suggested in 1984 by Sikivie [168], but nowadays a lot of experiments are looking for it.

Emergence of ALPs
ALPs are very similar to the axion, but important differences exist between them, mainly because the axion arises in a very specific context; in dealing with ALPs the aim is to bring out their properties in a model-independent fashion as much as possible. This attitude has two main consequences.

•
Only photon-ALP interaction terms are taken into account. Nothing prevents ALPs from coupling to other SM particles, but for our purposes they will henceforth be discarded.
Observe that such an ALP coupling to two photons aγγ is just supposed to exist without further worrying about its origin.
• The parameters m a and g aγγ are to be regarded as unrelated for ALPs, and it is merely assumed that m a 1 eV and g −1 As a result, ALPs are described by the Lagrangian where E and B have the same meaning as before. Note that L 0 ALP is included in L ren (a) + L nonren (φ SM , a).
Throughout this Review, we shall suppose that E is the electric field of a propagating photon beam while B is an external magnetic field. Because a couples to E · B, the effective coupling is proportional cos β, where β is the angle between E and B, and we denote by B T the B-component transverse to the beam momentum. Thus, what matters is B T and not B. As a consequence, if the photon beam propagates in an external magnetic field with varying direction, the beam polarization changes. This effect will be addressed in Section 10.
When the strength of either E or B is sufficiently large, also QED vacuum polarization effects should be taken into account, which are described by Heisenberg-Euler-Weisskopf Lagrangian [169,170] where α is the fine-structure constant and m e is the electron mass, so that ALPs are described by the full Lagrangian L ALP = L 0 ALP + L HEW (no confusion between the energy and the electric field will arise since the electric field will never be considered again).
So far, we have been concerned with the emergence of ALPs as pseudo-Goldstone bosons from the physics after compactification. However ALPs can also arise due to the topological complexity of the extra-dimensional manifold, and their properties depend on the specific compactification pattern. An example thereof is the so-called String Axiverse [16]. Finally, it has been shown that ALPs can be one of the decay products of the fundamental (moduli) fields Φ (for a review, see e.g., [19] and the references therein).

General Properties of ALPs
We are presently concerned with a monochromatic, polarized photon beam of energy E and wave vector k propagating in a cold plasma which is both magnetized and ionized. We suppose for the moment that an external homogeneous magnetic field B is present and we denote by n e the electron number density. We employ an orthogonal reference frame with the y-axis along k, while the x and z axes are chosen arbitrarily.

Beam Propagation Equation
It can be shown that in this case the beam propagation equation following from L ALP can be written as [32] where A x (y) and A z (y) denote the photon amplitudes with polarization (electric field) along the xand z-axis, respectively, while a(y) is the amplitude associated with the ALP. It is useful to introduce the basis {|γ x ≡ (1, 0, 0) T , |γ z ≡ (0, 1, 0) T , |a ≡ (0, 0, 1) T }, where |γ x and |γ z represent the two photon linear polarization states along the xand z-axis, respectively, and |a denotes the ALP state. Accordingly, we can rewrite ψ(y) as and the real, symmetric photon-ALP mixing matrix M 0 entering Equation (10) has the form where we have set While the terms appearing in the third row and column of M 0 are dictated by L ALP and have an evident physical meaning, the other ∆-terms require some explanation. They reflect the properties of the medium-which are not included in L ALP -and the off-diagonal ∆ xz = ∆ zx term directly mixes the photon polarization states giving rise to Faraday rotation, while ∆ xx and ∆ zz will be specified later.
As already stated, in the present paper we are interested in the situation in which the photon/ALP energy is much larger than the ALP mass, namely E m a . As a consequence, the short-wavelength approximation will be appropriate and greatly simplifies the problem. As first shown by Raffelt and Stodolsky [32], the beam propagation equation accordingly takes the form which is a Schrödinger-like equation with time replaced with y. We see that a remarkable picture emerges, wherein the beam looks formally like a three-state non-relativistic quantum system. Explicitly, they are the two photon polarization states and the ALP state. The evolution of the pure beam states-whose photons have the same polarizationis then described by the three-dimensional wave function ψ(y), with the y-coordinate replacing time, which obeys the Schödinger-like equation (15) with Hamiltonian Denoting by U 0 (y, y 0 ) the transfer matrix-namely the solution of Equation (15) with initial condition U 0 (y 0 , y 0 ) = 1, the propagation of a generic wave function can be represented as where y 0 represents the initial position. Moreover, we can set where U 0 (y, y 0 ) is the transfer matrix associated with the reduced Schödinger-like equation

A simplified case
Because B is supposed to be homogeneous, we have the freedom to choose the z-axis along B, so that B x = 0. The diagonal ∆-terms receive in principle two different contributions. One comes from QED vacuum polarization described by Lagrangian (9), but since here we suppose for simplicity that B is rather weak this effect is negligible [32]. The other contribution arises from the fact that the beam is supposed to propagate in a cold plasma, where charge screening produces an effective photon mass resulting in the plasma frequency [32] ω pl = 3.69 · 10 −11 n e cm −3 which entails Finally, the ∆ xz , ∆ zx terms account for Faraday rotation, but since we are going to take E in the X-ray or in the γ-ray band Faraday rotation can be discarded. Altogether, the mixing matrix becomes with the superscript (0) recalling the present choice of the coordinate system and We see that A x decouples away while only A z mixes with a.
Application of the discussion reported in Appendix A with M → M (0) 0 yields for the corresponding eigenvalues where we have set As a consequence, the transfer matrix associated with Equation ( where the matrices T 0,1 (0), T 0,2 (0) and T 0,3 (0) are just those defined by Equations (A17)-(A19) as specialized to the present situation. Actually, a simplification is brought about by introducing the photon-ALP mixing angle since then simple trigonometric manipulations allow us to express the above matrices in the simpler form Now, the probability that a photon polarized along the z-axis oscillates into an ALP after a distance y is evidently and in complete analogy with the case of neutrino oscillations [33] it reads which shows that ∆ osc plays the role of oscillation wave number, thereby implying that the oscillation length is L osc = 2π/∆ osc . Owing to Equations (27) and (32) can be rewritten as which shows that the photon-ALP oscillation probability becomes both maximal and independent of E and m a for ∆ osc g aγγ B , and explicitly reads P (0) 0,γ z →a (y) sin 2 g aγγ By 2 .
This is the strong-mixing regime, which-from the comparison of Equations (25) and (34)turns out to be characterized by the condition and so it sets in sufficiently above the energy threshold Observe that in the strong-mixing regime the mass term ∆ aa = −m 2 a /(2E) and the plasma term ∆ pl = −ω 2 pl /(2E) should be omitted in the mixing matrix, just for consistency. Below E L the photon-ALP oscillation probability becomes energy-dependent, oscillates typically over a decade in energy and next monotonically decreases becoming rapidly vanishingly small. The reader should keep this point in mind, since it will be used to put astrophysical bounds on g aγγ .

A More General Case
In view of our subsequent discussion it proves essential to deal with the general case in which B is not aligned with the z-axis but forms a nonvanishing angle ψ with it. Correspondingly, the mixing matrix M 0 presently arises from M (0) 0 through the similarity transformation operated by the rotation matrix in the x-z plane, namely This leads to indeed in agreement with Equation (13) within the considered approximation. Therefore the transfer matrix reads U 0 (y, and its explicit representation turns out to be with T 0,2 (ψ) ≡   sin 2 θ sin 2 ψ sin 2 α sin ψ cos ψ − sin α cos α sin ψ sin 2 α sin ψ cos ψ sin 2 α cos 2 ψ − sin α cos α cos ψ − sin α cos α sin ψ − sin α cos α cos ψ T 0,3 (ψ) ≡   sin 2 ψ cos 2 α sin ψ cos ψ cos 2 α sin α cos α sin ψ sin ψ cos ψ cos 2 α cos 2 ψ cos 2 α sin α cos α cos ψ sin ψ cos α sin α cos ψ sin α cos α sin 2 α   .
As a result, a beam initially containing only photons propagating in an external magnetic field, after a distance of many oscillation lengths L osc will contain ALPs. Moreover, the three states |γ x , |γ z , |a will equilibrate, namely the beam will be composed by 2/3 of photons and of 1/3 of ALPs.

Complications
So far we have considered the implications of Lagrangian (8) alone in order to introduce the reader in a rather gentle way into the present formalism. Nevertheless, we shall meet cases in which also Lagrangian (9) has to be taken into account. Moreover, we shall see that in some instances, other terms must be included into the mixing matrix, one of which is imaginary. As a consequence, the mixing matrix-which will be simply written as M-will not be self-adjoint, and correspondingly the transfer matrix will be denoted by U (y, y 0 ) and will be not unitary anymore. Thus, the beam looks formally like a three-state non-relativistic unstable quantum system.

Unpolarized Beam
So far, our discussion was confined to the case in which the beam is in a polarized state (pure state in the quantum mechanical language). This assumption possesses the advantage of making the resulting equations particularly transparent but it has the drawback that it is too restrictive for our analysis. Indeed, photon polarization cannot be measured in the VHE band with present-day and planned detectors, and so we have to treat the beam as unpolarized. As a consequence, it will be described by a generalized polarization density matrix rather than by a wave function ψ(y). Remarkably, the analogy with non-relativistic quantum mechanics entails that ρ(y) obeys the Von Neumann-like equation associated with Equation (19). Thus, the propagation of a generic ρ(y) is still given by and the probability that a photon/ALP beam initially in the state ρ 1 at position y 0 will be found in the state ρ 2 at position y is provided that ρ 2 is measured, since we are assuming as usual that Trρ 1 = Trρ 2 = 1.

Parameter Bounds
Before proceeding further, it seems important to report the observational bounds on the parameters g aγγ and m a , which have been considered free so far. Moreover, we would like to stress that in the absence of direct couplings to fermions ALP interact neither with matter nor with radiation, in spite of the two-photon coupling (the proof will be provided in Appendix B). * * * ALPs were searched for by the CAST experiment at CERN by using a superconductive magnet (decommissioned from the Large Hadron Collider) pointing towards the Sun. The upper side of the magnet was closed, in order to stop the X-rays produced in the outer part of the Sun, since ALPs are expected to be produced in the Sun core through the Primakoff process γ + ion → a + ion-represented in Figure 4-with an energy also in the X-ray band. However since they do not interact with matter, they travel unimpeded through the Sun and enter the magnet, which has an X-ray photon detector at its bottom. Owing of the strong magnetic field of the magnet, ALPs should convert into X-ray photons and be detected. Unfortunately, no detection has taken place and from this fact the resulting upper bound has been set as g aγγ < 0.66 · 10 −10 GeV −1 for m < 0.02 eV at the 2σ level [107]. Coincidentally, exactly the same bound has been derived from the analysis of some particular stars in globular clusters [92].
Let us next turn to the other bounds on g aγγ and m a . Basically, use is made of the observed absence of the characteristic oscillating behavior of the individual realizations of the beam propagation around E L (as explained in Section 3.2). Several bounds have been derived, among which the most relevant ones for us are the following.
Since 1996, Supernova1987A in the Large Magellanic Cloud has been used to set bounds on ALPs. Basically, the idea is that when the supernova exploded, a burst of ALPs produced by the Primakoff effect considered above was released, and when some of these ALPs entered the Milky Way they should have transformed into γ-rays and been detected by the Solar Maximum Mission satellite. From the failure of detection, upper bounds on m a and g aγγ can be established [42,43]. This issue was reconsidered in more detail in 2015 by Payez et al., getting the improved bound m a 4.4 · 10 −10 eV and g aγγ 5.3 · 10 −12 GeV −1 , without reporting the confidence level [94]. Although we do not trust such a bound since too many uncertainties are still present in our precise knowledge of a protoneutron star and the above analysis is rather superficial, two remarks are in order.
• Payez et al. 2015 say that ALPs are emitted simultaneously with neutrinos, and this is repeated by everybody. Concerning our ALPs which are supposed to interact only with two photons, this is not true. Since they do not interact either with matter nor with radiation (see Appendix B), they escape as soon as they are produced, while neutrinos remain trapped. Thus, this weakens the supernova bound. • Because the value of m a is rather uncertain-basically depending on where the strongmixing regime sets in-we will consider throughout this Review m a = O 10 −10 eV and g aγγ = O 10 −11 GeV −1 . Hence, we can easily take m a > 4.4 · 10 −10 eV no contradiction exists even apart from the previous remark.

Astrophysical Context
In order to make this Rewiev self-contained, we are going to provide all information that will be used in the next Sections.

Blazars
A mechanism for the very-high-energy (VHE) emission from active galactic nuclei (AGNs) that attracted a lot of interest consists of a binary system made of a supermassive black hole (SMBH) and a massive star: the latter provides a sort of reservoir of matter that accretes onto the SMBH. In a sense, the situation is similar to a Type IA supernova, but the explosion is replaced with matter falling into the SMBH. Correspondingly, about 10% of AGNs possess an accretion disk around the SMBH and supports two opposite relativistic jets emanating from the SMBH and perpendicular to an accretion disk, propagating from the central regions out to distances that-depending on the nature of the source-are in the range 1kpc-1Mpc. Ultrarelativistic particles (leptons and/or hadrons) in the gas carried by these jets are accelerated by relativistic shocks and emit non-thermal radiation extending from the radio band up to the VHE band. Aberration caused by the relativistic motion makes the emission strongly anisotropic, mainly in the direction of motion. Hence, in this case the emission is extremely beamed. When one jet happens to point towards us just for chance, the AGN is called a blazar, otherwise a radio galaxy. A schematic picture of a blazar (radio galaxy) is shown in Figure 5. As a matter of fact, the phenomenology of the AGNs is quite involved, and the complication of their classification is an indication that many of their aspects have not yet been understood [172][173][174][175] (a nice updated and concise discussion of this and related topics is contained in [176]). For instance, it is totally unclear why some blazars become flaring-from a few hours to a few days-and next come back to their state of smaller luminosity. As far as our considerations are concerned, we will restrict our attention to blazars and to their VHE photon emission mechanisms. It turns out that VHE blazars are typically hosted by elliptical galaxies in small groups.
According to the current wisdom, two competing VHE non-thermal photon emission mechanisms can work.

•
One is called leptonic mechanism (syncrotron-self Compton). Basically, in the presence of the magnetic fields inside the AGN jet, relativistic electrons emit synchrotron radiation and the produced photons are boosted to much higher energies by inverse Compton scattering off the parent electrons (in some cases also external photons from the disc participate in this process). The resulting emitted spectral energy distribution (SED) νF ν (ν)-F ν (ν) is the specific apparent luminosity-has two humps: the synchrotron one-somewhere from the IR band to the X-ray band-while the inverse Compton one lies in the γ-ray band [177][178][179][180]. • The other mechanism is named hadronic mechanism (proton-proton scattering). As far as the synchrotron emission is concerned the situation is the same as before, but the gamma hump is produced by hadronic collisions. The resulting π 0 immediately decays as π 0 → γ + γ, while the π ± produce neutrinos and antineutrinos [181,182]. Thus, the detection of these neutrinos can discriminate between the two mechanisms. In 2017 the IceCube neutrino telescope has detected one neutrino coming from the flaring blazar TXS 0506+056, thereby demonstrating that the hadronic mechanism does work [183].
As far as blazar photon polarization is concerned, the situation is as follows. In the X-ray band and below, where photons are produced via synchrotron emission, they are partially polarized with a realistic degree of linear polarization (see Section 10.1 for the definition) Π L =0.2-0.4, as argued e.g., in [184]. Instead, in the HE and VHE bands, where photons are likely produced via an inverse Compton process, they are expected to be unpolarized [185].
It turns out that VHE blazars fall into two sharply distinct classes. BL Lacs: They are named after the prototype BL Lacertae discovered in 1929 by Hoffmeister [186], but they were originally believed to be a variable star inside the Milky Way. Only in 1968 was it realized that they are instead a radio source, and in 1974 their redshift was found to be z = 0.07, corresponding to a distance l 300 Mpc. They lack broad optical lines, which entails that the broad line region (BLR) depicted in Figure 5 is absent. Moreover, this fact makes their redshift determination very difficult, which explains why only in the 1970s did BL Lacs starte to be understood. Their jets contain a magnetic field and extend out to about 1 kpc. Flat spectrum radio quasars (FSRQs): They are considerably more massive than BL Lacs, and their jets also contain a magnetic field but they extend out to about 1 Mpc. Their structure is by far more complicated than that of BL Lacs, especially in the outer region, with the presence of radio lobes and hot spots where a magnetic field is also present. Because of the very high density of ultraviolet photons in the BLR-effectively centred on the SMBH with radius (0.1-0.3) pc-and infrared photons emitted by the torus, the VHE photons produced at the jet base undergo the process γ + γ → e + + e − . As a result, the FSRQs should be invisible above (20-30) GeV [187][188][189][190]. However, observations with IACTs have shown that such an expectation is blatantly wrong, since In view of our subsequent analysis, we carefully address the propagation of a monochromatic photon beam emitted by a blazar at redshift z and detected at energy E 0 within the standard ΛCDM cosmological model, so that the emitted energy is E 0 (1 + z) owing to the cosmic expansion. Regardless of the actual physics responsible for photon propagation, two important quantities are the observed and emitted spectra (number fluxes) Φ ≡ dN/(dtdAdE). They are related by where P γ→γ (E 0 , z) is the photon survival probability throughout the whole trip from the source to us. Moreover, the SED is related to the observed spectrum by where F ν (ν, E 0 , z) is the specific apparent luminosity. We suppose hereafter that E 0 lies in the VHE γ-ray band.

Conventional Photon Propagation
Within conventional physics the photon survival probability P CP γ→γ (E 0 , z) is usually parametrized as where τ γ (E 0 , z) is the optical depth, which quantifies the dimming of the source. Note that in general τ γ (E 0 , z) increases with z, since a greater source distance entails a larger probability for a photon to disappear from the beam. Apart from atmospheric effects, one typically has τ γ (E 0 , z) < 1 for z not too large, in which case the Universe is optically thin up to the source. However depending on E 0 and z it can happen that τ γ (E 0 , z) > 1, so that at some point the Universe becomes optically thick along the line of sight to the source. The value z h such that τ γ (E 0 , z h ) = 1 defines the γ-ray horizon for a given E 0 , and it follows from Equation (52) that sources beyond the horizon tend to become progressively invisible as z further increases past z h . Owing to Equations (50) and (52) becomes Whenever dust effects can be neglected, photon depletion arises solely when hard beam photons of energy E scatter off soft background photons of energy permeating the Universe and isotropically distributed-we shall come back later to their nature-and produce e + e − pairs through the Breit-Wheeler γγ → e + e − process, represented by the Feynman diagram in Figure 6 [143]. Needless to say, in order for this process to take place enough energy has to be available in the centre-of-mass frame to create an e + e − pair. Regarding E as an independent variable, the process is kinematically allowed for where ϕ denotes the scattering angle, c is the speed of light and m e is the electron mass. Note that E and change along the beam in proportion of 1 + z. The corresponding Breit-Wheeler cross-section is [195] which depends on E, and ϕ only through the dimensionless parameter and the process is kinematically allowed for β 2 > 0. The cross-section σ γγ (E, , ϕ) reaches its maximum σ max γγ 1.70 · 10 −25 cm 2 for β 0.70. Assuming head-on collisions for definiteness (ϕ = π), it follows that σ γγ (E, , π) gets maximized for the background photon energy where E and correspond to the same redshift.
Within the standard ΛCDM cosmological model τ γ (E 0 , z) arises by first convolving the spectral number density n γ ( (z), z) of background photons at a generic redshift with σ γγ (E(z), (z), ϕ) along the line of sight for fixed values of z, ϕ and (z), and next integrating over all these variables [196][197][198]. Hence, we have where the distance travelled by a photon per unit redshift at redshift z is given by with Hubble constant H 0 70 Km s −1 Mpc −1 , while Ω Λ 0.7 and Ω M 0.3 represent the average cosmic density of matter and dark energy, respectively, in units of the critical density ρ cr 0.97 · 10 −29 g cm −3 .
Once n γ ( (z), z) is known, τ γ (E 0 , z) can be computed exactly, even though in general the integration over (z) in Equation (58) can only be performed numerically.
Finally, in order to get an intuitive insight into the physical situation under consideration it may be useful to discard cosmological effects (which evidently makes sense for z small enough). Accordingly, z is best expressed in terms of the source distance D = cz/H 0 and the optical depth becomes where λ γ (E) is the photon mean free path for γγ → e + e − referring to the present cosmic epoch. As a consequence, Equation (52) becomes and so Equation (53) reduces to Note that we have dropped the subscript 0 for simplicity.

Extragalactic Background Light (EBL)
Blazars detected so far with IACTs-or detectable in the near future-lie in the VHE range 100 GeV < E 0 < 100 TeV, and so from Equation (57) it follows that the resulting dimming is expected to be maximal for a background photon energy in the range 0.009 eV < 0 < 9 eV (corresponding to the frequency range 2.17 · 10 3 GHz < ν 0 < 2.17 · 10 6 GHz and to the wavelength range 0.14 µm < λ 0 < 1.38 · 10 2 µm), extending from the ultraviolet to the far-infrared. This is just the extragalactic background light (EBL). We stress that at variance with the case of the CMB, the EBL has nothing to do with the Big Bang. Rather, it is the radiation produced by all stars in galaxies during the whole history of the Universe and possibly by a first generation of stars formed before galaxies were assembled. Therefore, a lower limit to the EBL level can be derived from integrated galaxy counts [199].
Throughout this paper, we adopt the Franceschini-Rodighiero (FR) EBL model mainly because it supplies a very detailed numerical evaluation of the optical depth based on Equation (58), which will henceforth be denoted by τ γ (E 0 , z) [200], since it is in very good agreement with e.g., model of Dominguez et al. [201] (for a review, see e.g., [202]). Regretfully, the errors affecting τ γ (E 0 , z) are unknown. The dimming of a source at redshift z s due to the EBL as a function of the observed energy has been computed in [203] and shown in Figure 7.

Extragalactic Magnetic Field
Unfortunately, the morphology and strength of the extragalactic magnetic field B ext are totally unknown, and it is not surprising that various very different configurations of B ext have been proposed [204][205][206][207]. Yet, the strength of B ext is constrained to lies in the range 10 −7 nG B ext 1.7 nG on the scale O(1) Mpc [208][209][210] Nevertheless, a very realistic scenario for the extragalactic magnetic field exists has existed for a long time, which has become a classic and relies upon energetic galactic outflows.
What happens is that ionized matter from galaxies gets ejected into extragalactic space. The key-role is the fact that the associated magnetic field is frozen in, and amplified by turbulence, thereby magnetizing the surrounding space. Such a scenario was first proposed in 1968 by Rees and Setti [211] and in 1969 by Hoyle [212] in their investigations of radio sources. A more concrete and refined picture was considered in 1999 by Kronberg, Lesch and Hopp [213]. They proposed that dwarf galaxies are ultimately the source of B ext . Specifically, they start from the clear consensus that supernova-driven galactic winds are a crucial ingredient in the evolution of dwarf galaxies. Next, they show that shortly after a starburst, the kinetic energy supplied by supernovae and stellar winds inflate an expanding superbubble into the surrounding interstellar medium of a dwarf galaxy. Moreover, they demonstrate that the ejected thermal gas and cosmic-rays-significantly magnetized-will become mixed into the surrounding intergalactic matter, which will coexpand with the Universe. This picture leads to B ext = O(1) nG on the scale O(1) Mpc. Actually, in order to appreciate the relevance of dwarf galaxies, it is useful to recall that our Local Group-which is dominated by the Milky Way and Andromeda-contains 38 galaxies, 23 of which are dwarfs. Because the Local Group has nothing special, it follows that dwarf galaxies are roughly 10 times more abundant than bright Hubble type galaxies. The considered result is in agreement with observations of Lymanalpha forest clouds [214]. A similar situation was further investigated in 2001 by Furlanetto and Loeb [215] in connection with quasars outflows, which still predicts B ext = O(1) nG on the scale O(1) Mpc. Moreover, also normal galaxies possess this kind of ionized matter outflows-especially ellipticals and lenticulars-due to the central AGN (see e.g., [216] and the references therein) and supernova explosions. Remarkably, this picture is in agreement with numerical simulations [217]. Uncontroversial evidence of galactic outflows comes from the high metallicity (including strong iron lines) of the intracluster medium of regular galaxy clusters, which are so massive that matter cannot escape.
A classic and simple modeling of such a magnetic field configuration consists of a domainlike network, made of identically domains of size L dom equal to the B ext coherence length, all having roughly the same strength of B ext and with the direction of B ext uniform in each domain but randomly jumping from one domain to the next (for a review, see [204,205]).
Quite remarkably, the fact that in all above scenarios the seeds of B ext are galaxies indeed explains the three main features of the considered model for B ext . (1) Its cell-like morphology arises from its galactic origin. (2) It is quite plausible that B ext has nearly the same strength around each galaxy, and so in all domains. (3) Because galaxies are uncorrelated, it looks natural that the B ext direction changes randomly from the neighborhood of a galaxy to that of another, thereby explaining the random jump in direction of B ext from one domain to the next. Now, we have to address a critical point. Up until 2017, all models describing the propagation of a photon/ALP beam in extragalactic space assumed that the edges between adjacent domains are sharp, namely that the change in the direction of B ext from a domain to the next one is abrupt, which causes its components to be discontinuous. Models of this kind will be referred to as domain-like sharp edge models (DLSHE). Even though they are obviously a mathematical idealizations, they have been successfully used because the domain size L dom was invariably much larger than the oscillation length L osc . Because coherence is maintained only inside a single domain, this means that only a very small fraction of an oscillation probes a single domain, and so the discontinuity at the interface of two adjacent domains is not felt by the oscillations. However, if it happens that instead L osc L dom , a whole oscillation probes a single domain, and this model breaks down. In such a situation it becomes compelling to modify the model by replacing the sharp edges with smooth ones, so as to avoid any abrupt jump in the direction of B ext , which gives rise to a more complicated model called domain-like smooth edge model (DLSME) and developed in [114] and briefly described in Section 8.2.
A totally different approach to the extragalactic magnetic field is based on magnetohydrodynamic cosmological simulations (see e.g., [218,219] and the references therein). The strategy is as follows. An initial condition for a cosmological B ext is chosen arbitrarily during the dark age and its evolution as driven by structure formation is studied. The link with the real world is the condition to reproduce regular cluster magnetic fields, which fixes a posteriori the initial condition of B ext . As a by-product, a prediction of the magnetic field B fil inside filaments in the present Universe emerges. But this cannot be the whole story. Apart from failing to answer the question of the seed of primordial magnetic fields, galactic outflows are missing. This issue has a two-fold relevance: inside galaxy clusters and in extragalactic space. Indeed, galactic outflows are a reality in regular clusters, owing to the strong iron line. in 2009 Xu, et al. [220] claimed that the magnetic field ejected by a central AGN during the cluster formation can be amplified by turbulence during the cluster evolution in such a way to explain the observed cluster magnetic fields. Still in 2009, Donnert, et al. [221] have found that the strength and structure of the magnetic fields observed in clusters are well reproduced for a wide range of the model parameters by galactic outflows. Therefore, the requirement to reproduce regular cluster magnetic fields can totally mislead the expectations based on cosmological hydrodynamic simulations, in particular the present value of the magnetic field inside filaments and the model described in [106]. For this reason, we prefer to stick to the previous scenario.

Galaxy Clusters
According to Abell, Galaxy clusters contain from thirty up to more than thousands galaxies with a total mass in the range 10 14 -10 15 M and represent the largest gravitationally bound structures in the Universe. They are classified as (i) regular clusters, (ii) intermediate clusters, and (iii) irregular clusters. We consider regular clusters here, since most of them are spherical in first approximation, and so they are very easy to describe. In view of our subsequent needs, we focus our attention on the strength and morphology of their magnetic field B clu and electron number density n clu e . Faraday rotation measurements and synchrotron radio emissions tell us that B clu = O(1 − 10) µG [222,223]. The structure of B clu is believed to have a turbulent nature with a Kolmogorov-type turbulence spectrum M(k) ∝ k q , with k the wave number in the interval [k L , k H ] (more about this, in Section 10.2) and index q = −11/3 [224]. The behavior of B clu is modeled by [224,225] where B is the spectral function describing the Kolmogorov-type turbulence of the cluster magnetic field (for more details see e.g., [86]), B clu 0 and n clu e,0 are the central cluster magnetic field strength and the central electron number density, respectively, while η clu is a cluster parameter. The behavior of n clu e is modeled by the model n clu e (y) = n clu e,0 1 + where β clu is a cluster parameter and r core is the cluster core radius. Galaxy clusters are divided into two main categories: cool-core (CC) and non-cool-core (nCC) clusters. CC clusters generally host an AGN, while nCC clusters do not usually contain an active SMBH. Many aspects differentiate CC and nCC galaxy clusters (see e.g., [226]). However, for our studies their central electron number density n clu e,0 represents the real crucial quantity. We take the following average values for the two classes: n clu e,0 = 5 · 10 −2 cm −3 for CC clusters and n clu e,0 = 0.5 · 10 −2 cm −3 for nCC ones [226]. Other models with a larger number of parameters have been proposed in the literature, but they do not influence much our final results about γ ↔ a oscillations inside clusters.
In the cluster central region photons are produced by several processes in different energy ranges: thermal Bremsstrahlung is responsible in the X-ray band [227], while inverse Compton scattering, neutral pion decay are believed to produce photons in the HE range (see e.g., [228][229][230][231]). Photons produced by all these processes turn out to be effectively unpolarized.

Propagation of ALPs in Extragalactic Space-1
Let us consider a far away VHE blazar at redshift z which is presently detected by an IACT. We stress that H.E.S.S., MAGIC and VERITAS are sensitive to photons with energy from about 100 GeV up to a few TeV. As a consequence, we have E 0 m a and so we can apply the formalism developed in Section 3, but a small extension is needed in order to take EBL absorption into account.
Our ultimate task is the computation of the photon survival probability P ALP γ→γ (E 0 , z) in the presence of ALPs. We have seen that in conventional physics photons undergo EBL absorption, which severely depletes the photon beam when z is sufficiently large. Clearly, now-owing to the presence of the extragalactic magnetic field-γ ↔ a oscillations will take place in the beam. This means that during its propagation a photon acquires a 'split personality': for some time it behaves as a true photon-thereby undergoing EBL absorption-but for some time it behaves as an ALP, and so it is unaffected by the EBL and propagates freely. Therefore, the optical depth in the presence of ALPs τ ALP γ is now smaller than in the conventional case. But since the corresponding photon survival probability is recalling (52) we conclude that P ALP γ→γ (E 0 , z) is much larger than P γ→γ (E 0 , z) evaluated in conventional physics: this is the crux of the argument. As a consequence, far-away sources that are too faint to be detected according to conventional physics would become observable.
In order to be definite-and in view of the discussion to be presented in the next Sectionwe choose the values of some parameters in agreement with the subsequent needs.
As far as the extragalactic magnetic field is concerned, we assume a domain-like structure described in Section 4.4 with a DLSHE model, since we shall see that L osc L dom .

Strategy
Thanks to the fact that B is homogeneous in every domain, the beam propagation equation can be solved exactly in every single domain. But due to the nature of the extragalactic magnetic field, the angle of B in each domain with a fixed fiducial direction equal for all domains (which we identify with the z-axis) is a random variable, and so the propagation of the photon/ALP beam becomes a N d -dimensional stochastic process, where N d denotes the total number of magnetic domains crossed by the beam. Moreover, we shall see that the whole photon/ALP beam propagation can be recovered by iterating N d times the propagation over a single magnetic domain, changing each time the value of the random angle. Therefore, we identify the photon survival probability with its value averaged over the N d angles.
Our discussion is framed within the standard ΛCDM cosmological model with Ω M = 0.3 and Ω Λ = 0.7, and so the redshift is the natural parameter to express distances. In particular, Accordingly, the overall structure of the cellular configuration of the extragalactic magnetic field is naturally described by a uniform mesh in redshift space with elementary step ∆z, which is therefore the same for all domains. This mesh can be constructed as follows. We denote by L (n) dom = L (n − 1)∆z, n∆z the proper length along the y-direction of the generic n-th domain, with 1 ≤ n ≤ N d . Note that N d is the maximal integer contained in the number z/∆z, hence N d z/∆z. In order to fix ∆z we consider the domain closest to us, labelled by 1 and-with the help of Equation (66) dom z. Second, the proper length of the n-th domain along the y-direction follows from Equation (66) * * * Manifestly, in order to maximize P ALP γ→γ (E 0 , z) we choose m a in order to be in the strongmixing regime for E 0 100 GeV. Incidentally, when the external magnetic field is homogeneousas it is in fact in each single domain-a look at Lagrangian (8) shows that all results depend on the combination g aγγ B and not on g aγγ and B separately. It is therefore quite convenient to employ the parameter in terms of which Equation (37) can be rewritten as Because we would like to be in the strong-mixing regime almost everywhere within the VHE band, we take E L = O(100) GeV. What about m a ? We should keep in mind that ω pl is unknown, but the upper bound on the mean diffuse extragalactic electron density n e < 2.7 · 10 −7 cm −3 is provided by the WMAP measurement of the baryon density [232], whichthanks to Equation (20)-translates into the upper bound ω pl < 1.92 · 10 −14 eV. Moreover, in order to fix ξ we use the fact that the result to be derived in the next Section requires ξ = 0.5. As a consequence, we get m a = O(10 −10 ) eV.

Propagation over a Single Domain
We have to determine λ (n) γ and the magnetic field strength B (n) in the generic n-th domain. The first goal can be achieved as follows. Because the domain size is so small as compared to the cosmological standards, we can safely drop cosmological evolutionary effects when considering a single domain. Then as far as absorption is concerned what matters is the mean free path λ γ for the reaction γγ → e + e − , and the term i/2λ γ should be inserted into the 11 and 22 entries of the M matrix. In order to evaluate λ γ , we imagine that two hypothetical sources located at both edges of the n-th domain are observed. Therefore, we apply Equation (53) which upon combination imply that the flux change across the domain in question is But owing to Equation (62) mutatis mutandis implies that Equation (72) should have the form and the comparison with Equation (72) ultimately yields where the optical depth is evaluated by means of Equation (58) or more simply taken from [200].
As for the determination of B (n) , we note that because of the high conductivity of the IGM medium the magnetic flux lines can be thought as frozen inside it [204,205]. Therefore, the flux conservation during the cosmic expansion entails that B scales like (1 + z) 2 , so that the magnetic field strength in a domain at redshift z is B(z) = B(z = 0)(1 + z) 2 [204,205]. Hence in the n-th magnetic domain we have Thus, at this stage the mixing matrix M as explicitly written in the n-th domain reads where ψ n is the random angle between B (n) and the fixed fiducial direction along the z-axis (note that indeed M † = M). Observe that since we are in the strong-mixing regime, ω pl and m a can be neglected with respect to the other terms in the mixing matrix. So-apart from ψ nall other matrix elements entering M (n) are known. Finding the transfer matrix corresponding to M (n) is straightforward even if tedious by using the results reported in Appendix A. The result is with where ξ n is just ξ as defined by Equation (68) and evaluated in the n-th domain.

Calculation of the Photon Survival Probability in the Presence of Photon-ALP Oscillations
As we said, our aim is to derive the photon survival probability P ALP γγ (E 0 , z) from the source at redshift z to us in the present context. So far, we have dealt with a single magnetic domain but now we enlarge our view so as to encompass the whole propagation process of the beam from the source to us. This goal is trivially achieved thanks to the analogy with non-relativistic quantum mechanics, according to which-for a fixed arbitrary choice of the angles {ψ n } 1≤n≤N d -the whole transfer matrix describing the propagation of the photon/ALP beam is Moreover, the probability that a photon/ALP beam emitted by a blazar at z in the state ρ 1 will be detected in the state ρ 2 for the above choice of {ψ n } 1≤n≤N d is given by with Trρ 1 = Trρ 2 = 1.
Since the actual values of the angles {ψ n } 1≤n≤N d are unknown, the best that we can do is to evaluate the probability entering Equation (84) as averaged over all possible values of the considered angles, namely indeed in accordance with the strategy outlined above. In practice, this is accomplished by evaluating the r.h.s. of Equation (84)  Because the photon polarization cannot be measured at the considered energies, we have to sum the result over the two final polarization states Moreover, we suppose here that the emitted beam consists 100% of unpolarized photons, so that the initial beam state is described by the density matrix We find in this way the photon survival probability P ALP γ→γ (E 0 , z) A final remark is in order. It is obvious that the beam follows a single realization of the considered stochastic process at once, but since we do not know which one is actually selected the best we can do is to evaluate the average photon survival probability.

VHE BL Lac Spectral Anomaly
After all these preliminary considerations, we are now in a position to discuss an important effect. Actually, we show that VHE astrophysics leads to a first strong hint at ALPs.
Basically, we are going to demonstrate that conventional physics leads to a paradoxical situation concerning the EBL-deabsorbed BL Lac spectra as a function of the source redshift z. But such a situation disappears altogether once the ALPs consistent with the observational bounds enter the game.
As a first step, we have to select a sample of BL Lacs which is suitable for our analysis. They have to meet the following conditions.

1.
We focus our attention on flaring blazars, which show episodic time variability with their luminosity increasing by more than a factor of two, on the time span from a few hours to a few days: the reason is both their enhanced luminosity-which entails in turn their detectability [233,234]-and our desire to consider a homogeneous sample of BL Lacs.

2.
As we shall see, our analysis requires the knowledge of the redshift, the observed spectrum and the energy range wherein every blazar is observed. This information is available only for some of the observed flaring sources.

3.
In order to get rid of evolutionary effects inside blazars we restrict our attention to those with z ≤ 0.6.

4.
It seems a good thing to deal with sources that are as similar as possible. Therefore, we consider only intermediate-frequency peaked (IBL) and high-frequency peaked (HBL) flaring BL Lacs with observed energy E 0 100 GeV.
We are consequently left with a sample S of 39 flaring VHE BL Lacs, which are listed in Appendix C.
In first approximation, all observed spectra of the VHE blazars in S are fitted by a single power-law-neglecting a possible small curvature of some spectra in their lowest energy part-and so they have the form where E 0 is the observed energy, E ref is a common reference energy whileK obs (z) and Γ obs (z) denote the normalization constant and the observed slope, respectively, for a source at redshift z. Actually,K obs (z) is generally defined at different energies for different sources. So, for the sake of comparison among all observed spectra normalization constants we need to perform a rescalingK obs (z) → K obs (z) in the observed spectrum of the considered blazars in such a way that K obs (z) coincides with Φ obs (E 0 , z) at the fiducial energy E 0, * = 300 GeV for every source in S. Accordingly, Equation (89) becomes As already emphasized, these spectra are strongly affected by the EBL, hence if we want to know the shape of the emitted spectra they have to be EBL-deabsorbed. We know that the emitted and observed spectra are related by Equation (53), which we presently rewrite as where CP stands throughout this Section for conventional physics.

Conventional Physics
Let us start by deriving the emitted spectrum of every source in S, starting from each observed one. As a preliminary step, thanks to Equation (53) we rewrite Equation (91) as Because of the presence of the exponential in the r.h.s. of Equation (92), Φ CP em E 0 (1 + z) cannot behave as an exact power law. Yet, it turns out to be close to it. Therefore, we best-fit (BF) Φ CP em E 0 (1 + z) to a single power-law expression over the energy range ∆E 0 (z) where a source is observed, and so E 0 varies inside ∆E 0 (z) (which changes from source to source). Correspondingly, the resulting values of Γ CP em (z) are plotted in Figure 8.
We proceed by performing a statistical analysis of all values of Γ CP em (z) as a function of z, by employing the least square method and try to fit the data with one parameter (horizontal straight line), two parameters (first-order polynomial), and three parameters (second-order polynomial). In order to test the statistical significance of the fits we evaluate the corresponding χ 2 red,CP . The values of the χ 2 red,CP obtained for the three fits are 2.37 (one parameter), 1.49 (two parameters) and 1.46 (three parameters). Thus, data appear to be best-fitted by the second-order polynomial Γ CP em (z) = − 5.33 z 2 − 0.66 z + 2.64 .
The best-fit regression line given by Equation (94) turns out to be a concave parabola shown in Figure 9.  This is the key-point. In order to appreciate the physical consequences of Equation (94) we should keep in mind that Γ CP em (z) is the exponent of the emitted energy entering Φ CP em (E). Hence, in the two extreme cases z = 0 and z = 0.6 we have thereby implying that the hardening of the emitted flux progressively increases with the redshift. More generally, we have found a statistical correlation between the {Γ CP em (z)} and z.
However, this result looks physically absurd. How can the sources get to know their z so as to tune their Γ CP em (z) in such a way to reproduce the above statistical correlation? We call the existence of such a correlation the VHE BL Lac spectral anomaly, which of course concerns flaring BL Lac alone. According to physical intuition, we would have expected a straight horizontal best-fit regression line in the Γ em − z plane.
The most natural explanation would be that such an anomaly arises from selection effects, but it has been demonstrated that this is not the case [118].

ALPs Enter the Game
As an attempt to get rid of the VHE BL Lac spectral anomaly, we put ALPs into play, with parameters consistent with the previously mentioned bounds. Because the presently operating IACTs reach at most e few TeV, the oscillation length is much larger than the magnetic domain size L dom and so the propagation model in extragalactic space considered in Section 5 is fully adequate.
Basically, we go through exactly the same steps described above. That is to say, we rewrite Equation (92) Next, we still best-fit Φ ALP em E 0 (1 + z) to a single power law expression over the energy range ∆E 0 (z) where a source is observed, hence E 0 varies within ∆E 0 (z). Such a best-fitting procedure is performed for every benchmark value of ξ and L dom , namely L dom = 4 Mpc, ξ = 0.1, 0.5, 1, 5, and L dom = 10 Mpc, ξ = 0.1, 0.5, 1, 5. Moreover, we carry out again the above statistical analysis of the values of Γ ALP em (z) for all blazars in S, for any benchmark value of ξ and L dom .
Finally, the statistical significance of each fit can be quantified by computing the corresponding χ 2 red,ALP , whose values are reported in Table 1 for L dom = 4 Mpc, ξ = 0.1, 0.5, 1, 5, and in Table 2 for L dom = 10 Mpc, ξ = 0.1, 0.5, 1, 5. In both Tables the values of χ 2 red,CP are reported for comparison.  The relevance of such a statistical analysis is to single out two preferred situations (corresponding to the minimum of χ 2 red,ALP ): one for L dom = 4 Mpc and the other for L dom = 10 Mpc. In either case, the results are ξ = 0.5 and a straight best-fit regression line which is exactly horizontal. More in detail, for L dom = 4 Mpc we get χ 2 red,ALP = 1.29 and Γ ALP em = 2.54, while for L dom = 10 Mpc we find χ 2 red,ALP = 1.25 and Γ ALP em = 2.60. Manifestly, both cases turn out to be very similar. We plot the values of Γ ALP em (z) in Figure 10 only for the two considered situations.
Because ξ = 0.5 is our preferred value, we are now in a position to make a sharp prediction of the ALP parameters. Correspondingly-owing to Equation (69) with E L = O(100) GeVthe ALP mass must be m = O(10 −10 ) eV, since in Section 5.1 we have seen that ω pl < 1.92 · 10 −14 eV. Moreover, by recalling Equation (68) with ξ = 0.5 and the upper bounds on g aγγ and B quoted in Section 3.6 we get 2.94 · 10 −12 GeV −1 < g aγγ < 0.66 · 10 −10 GeV −1 .
Remarkably, these parameters are consistent with the bounds reported in Section 3.6.
In conclusion, we have indeed succeeded in getting rid of the VHE BL Lac spectral anomaly, since the Γ ALP em are on average independent of z. We stress that it is an automatic consequence of the ALP scenario, and not an ad hoc requirement.
A final remark is in order. It is obvious that by effectively changing the EBL level-this is what the ALP actually does-the best-fit regression line also changes. But that it transforms from a concave parabola into a perfectly straight horizontal line looks almost a miracle!

A New Scenario for Flaring BL Lacs
Besides getting rid of the VHE BL Lac spectral anomaly, the ALP scenario naturally leads to a new view of flaring BL Lacs.
In order to best appreciate this point, it is enlightening to fit the values of Γ CP em (z) by a horizontal straight regression line, at the cost of relaxing the best-fitting requirement. Accordingly, the scatter of the values of Γ CP em (z) for 95% of blazars belonging to S is less than 20% of the mean value set by horizontal straight regression line in Figure 11, namely equal to 0.47. Superficially, the VHE BL Lac spectral anomaly problem would be solved-but in reality it is not-since we correspondingly have χ 2 red,CP = 2.37 which is by far too large. The result obtained in the presence of γ ↔ a oscillations and ξ = 0.5 leads to a similar but much more satisfactory picture. In the first place, we are dealing with a horizontal straight best-fit regression line, and in addition the corresponding χ 2 red,ALP turns out to be considerably smaller. Specifically, the scatter of the values of Γ ALP em (z) for 95% of the considered blazars is now less than 13% about the mean value set by Γ ALP em = 2.54 for L dom = 4 Mpc and less than 13% about the mean value set by Γ ALP em = 2.60 for L dom = 10 Mpc, namely equal to 0.33 for L dom = 4 Mpc and equal to 0.32 for L dom = 10 Mpc. This situation is shown in Figure 12. We argue that the small scatter in the values of Γ ALP em (z) implies that the physical emission mechanism is the same for all flaring blazars, with the small fluctuations in Γ ALP em (z) arising from the difference of their internal quantities: after all, no two identical galaxies have ever been found! On the other hand, the larger scatter in the values of K ALP em (z) as derived in [118]presumably unaffected by photon-ALP oscillations when error bars are taken into account-is naturally traced back to the different environmental state of each flaring source, such as for instance the accretion rate.
A natural question finally arises. How is it possible that the large spread in the {Γ obs (z)} distribution (see Appendix C) arises from the small scatter in the {Γ ALP em (z)} distribution shown in Figure 10? The answer is very simple: most of the scatter in the {Γ obs (z)} distribution arises from the large scatter in the source redshifts.

New Explanation of VHE Emission from FSRQs
As emphasized in Section 4.1, according to conventional physics FSRQs do not emit at energy larger than about 30 GeV. But the IACTs have detected them up to an energy of 400 GeV [235][236][237]. This fact poses a formidable problem to VHE astrophysicists, who have developed contrived models as a way out of this conundrum, but none of them is really satisfactory. The most iconic example of these FSRQs is represented by PKS 1222+216, and we shall focus our attention on it.

Detection of PKS 1222+216 in the VHE Range
The detection of an intense, rapidly varying emission in the energy range 70 GeV-400 GeV from the FSRQ PKS 1222+216 at redshift z = 0.432 represents a challenge for all blazar models. Since the surrounding of the inner jet in FSRQs is rich in optical/ultraviolet photons emitted by the BLR (see Section 4.1), a huge optical depth for γ rays above 20 GeV-30 GeV is expected (see e.g., [187][188][189]). However, photons up to 400 GeV have been observed by MAGIC [236]. In addition, the PKS 1222+216 flux doubles in only about 10 minutes, thereby implying an extreme compactness of the emitting region. These features are very difficult to explain by the standard blazar models.
The only solution within conventional physics to solve both issues-the detection of PKS 1222+216 in the VHE band, and its rapid variation-appears to deal with a two-blob model: a larger one located in the inner region of the source which is responsible for the emission from IR to X-rays, and a smaller compact blob (r ∼ 10 14 cm) accounting for the VHE emitting region detected by MAGIC beyond the BLR-which is therefore far from the central engine-in order to avoid absorption [238][239][240]. Manifestly, this is an ad hoc solution.
Because PKS 1222+216 has also been simultaneously detected by Fermi/LAT in the energy range 0.3 GeV-3 GeV [241], it is compelling to find a realistic SED that fits both the Fermi/LAT and the MAGIC observations, besides to explaining the VHE γ-ray emission.
We want now to inquire whether a similar two-blob model can produce a physically consistent SED, with the key-difference that also the smaller blob is now located close to the center.

Observations and Setup
The relevant physical parameters for PKS 1222+216 are as follows. We assume a disk luminosity L disk 1.5 · 10 46 erg s −1 , a radius of the BLR R BLR 0.23 pc, and standard values for the cloud number density n c 10 10 cm −3 and the temperature T c 10 4 K of the BLR (see e.g., [189]). However, the average electron number density n e , which is relevant for the beam propagation, gets considerably reduced to n e 10 4 cm −3 .
We now estimate the BLR absorption by evaluating the optical depth τ(E) of the beam photons inside the BLR according to conventional physics. By following the same procedure developed in [189], the optical depth reads where E is the energy of a γ ray, x is the distance from the center of the BLR, µ ≡ cos θ where θ is the scattering angle between a γ ray and a soft photon of energy of the BLR, dΩ = − 2πdµ, while n ph ( , Ω, x) (which is computed by means of the standard photo-ionization code CLOUDY as in [242]) is the spectral number density of the BLR radiation field at position x per unit solid angle and σ γγ (E, , µ) is the pair-production cross-section of Equation (55). The resulting τ(E) is represented by the blue long-dashed line in Figure 13, which shows that we cannot expect photons from PKS 1222+216 in the energy range 70 GeV-400 GeV. But MAGIC has detected such photons. The calculated τ(E) is affected by some degree of uncertainty coming from the uncertainty of the input parameters, and in particular the luminosity of the disk L disk . However, since τ ∝ L 1/2 disk (R BLR ∝ L 1/2 disk , see [243] and references therein), the final impact of these uncertainties is moderate. In addition, scattered disk photons [240] show a maximal absorption of τ 0.2 at about 200 GeV, so that their contribution to the total τ(E) can safely be neglected.

An ALP Model for PKS 1222+216
A natural explanation of the VHE observations arises if γ ↔ a oscillations are put into the game, according to the following scenario. They take place within the BLR, ALPs cross this region unimpeded and re-conversion into photons occurs either in the magnetic field of the source or in that of the host galaxy (more about this, later). Thanks to γ ↔ a oscillations, we can stay within the standard blazar emission models. The resulting SED looks quite realistic, and simultaneously fits both the Fermi/LAT and MAGIC spectra.
We use the same conventions of the previous Sections. In particular, four different regions are crossed by the photon/ALP beam: (i) the inner region where B is homogeneous in first approximation, (ii) the large scale jet where B possesses a smooth y-dependence, (iii) the host galaxy where B has a domain-like structure (as we shall see) and (iv) the extragalactic space where B presents again a domain-like structure.
In the inner region the magnetic field strength is so strong that the photon one-loop vacuum polarization effects coming from L HEW in Equation (9) are not negligible, which is a further complication with respect to the simple scenario outlined in Section 3.2. Correspondingly, to the 11 and 22 entries of the mixing matrix M 0 in Equation (13) two new terms must be added, which read and respectively, where B cr 4.41 · 10 13 G is the critical magnetic field. Thus, we can introduce also a high-energy cutoff above which the photon/ALP beam does not propagate within the strong-mixing regime and presents an energy dependent behavior. We will now address γ ↔ a oscillations in the several regions crossed by the beam. We consider very light ALPs as in [53,55,56,62,70]. We take m a = O(10 −10 ) eV as in Section 5.

Photon-ALP Oscillations inside the BLR
We start by evaluating the photon/ALP beam propagation in the inner part of the blazar, which is the region extending from the centre to R BLR 0.23 pc, to be referred to as region 1.
Because the magnetic field profile along the jet decreases starting from the center and possesses a complicated morphology, we prefer to assume its strength and orientation to be constant and equal to their average values from the center to the edge of the BLR, since their precise estimate is very difficult due to the presence of strong shocks and relativistic winds. Thus, we take B 0.2 G (B 2.2 G at the base of the jet [244]) and an angle of 45 • with the beam direction, since γ ↔ a oscillations vanish if B is exactly along the beam, while it is maximal for B transverse to the beam. Whence B T = 0.14 G.
With the previous parameter choice, using the CAST bound [107] and employing Equations (37) and (101) we see that for Fermi/LAT data we are in the strong-mixing regime but for MAGIC observations we are beyond E H and thus in the weak mixing regime. Still, we will observe that γ ↔ a oscillations are relevant well above E H .
We calculate the mean free path inside the BLR as where τ(E) is the optical depth for the γγ → e + e − process reported in Figure 13 and represented by blue long-dashed line. As a consequence, to leading order the various terms entering the mixing matrix M 0 of Equation (13) are ∆ aγ = 1 2 g aγγ B T 1.37 · 10 −23 g aγγ 10 11 GeV eV , It is now possible to evaluate the transfer matrix U 1 (R BLR , 0; E) in this region by means of the procedure developed in Appendix A.

Photon-ALP Oscillations in the Large Scale Jet
We now estimate the γ ↔ a oscillations in the jet beyond R BLR , which we call region 2. In this zone B is believed to possess a toroidal behavior so that B(y) ∝ y −1 [245,246] and B T (y) reads B T (y) 0.14 R BLR y G 3.22 · 10 −2 pc y G .
Region 2 extends up to R * , which is defined as the distance where B T (y) in Equation (107) reaches the value assumed by the strength of turbulent magnetic field in the host elliptical galaxy, whose typical strength is 5 µG (more about this, later). Accordingly, we get R * 6.7 kpc.
By employing Equations (37) and (101) with the previous choice of the parameter values and g aγγ = O(10 −11 ) GeV −1 , we can conclude that the photon/ALP beam propagates in the strong-mixing regime over the whole region 2. Therefore, the plasma, the ALP mass and the QED one-loop terms can be safely neglected. In this region the same is true concerning photon absorption, so that the various terms entering the mixing matrix M 0 in Equation (13) are ∆ aγ = 1 2 g aγγ B T 3.1 · 10 −24 pc y g aγγ 10 11 GeV eV .
In the present situation, the transfer matrix can be obtained by analytically solving the beam propagation equation and reads where B T (R BLR ) is given by Equation (107).

Photon-ALP Oscillations in the Host Galaxy
As anticipated, FSRQs are usually located in elliptical galaxies, whose B is known with great uncertainty. Nevertheless, it is believed that B possesses a turbulent nature and can be modeled with a domain-like structure, with strength 5 µG and domain size equal to 150 pc [247]. Thus, the photon/ALP beam propagates in this region, which we call region 3 from R * up to the radius of the host galaxy R host . The system is in the strong-mixing regime in region 3. By observing that also in this case photon absorption is negligible, the terms entering the mixing matrix M 0 in Equation (13) become simply ∆ aγ = 1 2 g aγγ B T 3.4 · 10 −28 g aγγ 10 11 GeV eV , (112) while the transfer matrix in this region U 3 (R host , R * ; E) can be calculated by means of a strategy very similar to the one developed in Section 5 and in Appendix A. However, it turns out that in practice γ ↔ a oscillations in region 3 are totally negligible.

Photon-ALP Oscillations in Extragalactic Space
Finally, the photon/ALP beam propagates in extragalactic space from R host to us. We call this region 4. It goes without saying that the treatment of γ ↔ a oscillations in this region is identical to the one developed in Section 5, and we have nothing to add. We denote transfer matrix in the extragalactic space by U 4 (D L , R host ; E), where D L denotes the luminosity distance of PKS 1222+216 (corresponding to z = 0.432).

Overall Photon-ALP Oscillations
The whole transfer matrix U tot (D L , 0; E) associated with the propagation of the photon/ALP beam from the inner jet of PKS 1222+216 to us can be obtained by multiplying in the correct order the transfer matrices calculated in the previous Sections. Correspondingly we obtain Now, let us consider the total probability that a photon/ALP beam emitted by the considered FSRQ in the state ρ 1 will be detected in the state ρ 2 for a fixed configuration of the B direction in each domain of the host galaxy and of the extragalactic space. Unless otherwise stated, we suppose that the angles {ϕ m } 1≤m≤N g and {ψ n } 1≤n≤N d of B in the magnetic domains of the host galaxy and if extragalactic space, respectively are held fixed. Then Equation (49) entails with Trρ 1 = Trρ 2 = 1. But since their values are unknown-the situation is formally identical to that considered in Section 5-we are actually dealing with a N g N d -dimensional stochastic process. As already discussed in Section 5, the total photon survival probability is therefore given by where for clarity we have explicitly written all angles that are averaged over.

Results
We now want to clarify the role of the γ ↔ a oscillations for two issues concerning PKS 1222+216.

•
The explanation of why MAGIC data have been observed [236]. • The fit to both low energy observations from Fermi/LAT [241] and to those at high energy from MAGIC [236] with a physically motivated SED.
In order to investigate the first item, we consider the photon survival probability from the center of PKS 1222+216 up to R * by means of the transfer matrix The corresponding photon survival probability in the presence of γ ↔ a oscillations is given by (65), which presently can be written as We have considered three benchmark cases for B and g aγγ .
The curves for τ ALP γ (E) corresponding to the above cases are reported in Figure 13, where the red solid line corresponds to case (1), the green short-dashed line to case (2), the violet dashed-dotted line to case (3), while the blue long-dashed line represents the case of conventional physics (as in [189]).
From Figure 13, we see that γ ↔ a oscillations greatly reduce the optical depth in the optically thick range. In particular, in the case (1)-which represents our best option-the effective optical depth is almost constant in the MAGIC band around τ ALP γ 4, which cor-responds to a photon survival probability of about 2%. Instead, in the optically thin region below ∼30 GeV, the optical depth with γ ↔ a oscillations is larger than in conventional physics because a fraction of emitted photons are converted to ALPs close to the source but cannot be converted back. Thus, we have a solution to the first problem, namely why MAGIC data have been observed [236].
Let us turn our attention to the second issue. In order to fit both Fermi/LAT and MAGIC data at once with a physically motivated SED, it is instructive to define the quantity P ≡ log P ALP γ→γ (R * , 0; 1 GeV) P ALP γ→γ (R * , 0; 300 GeV) , (118) measuring the ratio between P ALP γ→γ at 1 GeV, which is representative of Fermi/LAT observations, and P ALP γ→γ at 300 GeV which accounts for MAGIC ones. In order to get an acceptable shape for the emitted SED at VHE, we need to have P as low as possible. Indeed, a low value of P not only would reduce the discrepancy between Fermi/LAT and MAGIC data allowing us to produce a smooth SED, but would also give rise to a big correction to the MAGIC spectrum. By exploring the parameter space we conclude that case (1) looks like the best one.
We now present our model for the emitted spectrum of PKS 1222+216 in the whole γ-ray band. The emitted spectrum Φ em E 0 (1 + z) is linked to the observed one Φ obs (E 0 ) by where z is the redshift of PKS 1222+216 and P ALP γ→γ (D L , 0; E 0 ) is computed from the center of PKS 1222+216 to us.
As already stressed, in order to explain the behavior of PKS 1222+216 within conventional physics a two-blob model has been proposed. We recall that a larger blob is located inside the source and is responsible for the emission from IR to soft γ-rays, and a smaller compact blob-accounting for the emission detected by MAGIC-is present well outside the BLR in order to avoid absorption [238]. We want now to inquire whether a similar two-blob model can produce a physically consistent SED with the key difference that the smaller blob is now located close to the center.
Within case (1)-which we have found to be the best situation-we consider electrons radiating through synchrotron, and inverse Compton processes (with both the internally produced synchrotron radiation and the external radiation from the BLR). In each blob several parameters must be chosen: size r, magnetic field B, bulk Lorentz factor Γ, electron normalization K, minimum, break and maximum Lorentz factors γ min , γ b , γ max , and the low energy (n 1 ) and the high energy (n 2 ) slope of the smoothed power law electron energy distribution. Concerning the larger blob we adopt the same parameters of the original model [238], while for the compact VHE γ-ray emitting blob we take r = 2.2 · 10 14 cm, B = 0.008 G, Γ = 17.5, K = 6.7 · 10 9 , γ min = 3 · 10 3 , γ b = 1.2 · 10 5 , γ max = 4.9 · 10 5 and electron slopes n 1 = 2.1, n 2 = 3.5. The resulting SED is reported in Figure 14 for the above parameters.
Actually, Figure 14 shows that our model not only explains the detection of PKS 1222+216 by MAGIC but also leads to a physical motivated SED. Moreover, according to case (1)which we restate to be our best option-the VHE luminosity is L γ = 6 · 10 48 erg s −1 . While other choices such as case (2) give rise to a satisfactory SED shape (see [76] for details), the corresponding VHE luminosity L γ = 10 51 erg s −1 looks unrealistic, since it is about 100 times larger than that of the brightest VHE blazars (see e.g., [248]). In addition, it turns out that γ ↔ a oscillations in extragalactic space is not important for our model, since a SED similar to the one reported in Figure 14 emerges even if such oscillations are discarded (see [76] for details).  (119). The gray data points below 10 20 Hz are irrelevant for the present discussion (details can be found in [238]). In addition, the dashed and solid curves show the SED resulting from the considered two blobs which account for the γ-ray emission at high energy and VHE, respectively. (Credit [76]).
In conclusion, without invoking ad hoc models we have solved the problems caused by PKS 1222+216. This fact represents a second strong hint at the existence of an ALP with the properties specified above, since also in this model we have taken m a = O(10 −10 ) eV. Note that our model can be applied to all other VHE FSRQs with analogous results.

Propagation of ALPs in Extragalactic Space-2
So far, we have been dealing with the DLSHE model for the extragalactic magnetic field since at the VHE currently probed by the IACTs the photon-ALP oscillation length L osc is much larger than the size L dom of the magnetic domains.
However, in 2015 Dobrynina, Kartavtsev and Raffelt [145] realized that at even larger energies photon dispersion on the CMB (Cosmic Microwave Background) becomes the leading effect, which implies L osc to decrease as E further increases. Therefore, things completely change whenever L osc L dom , since in this case a full oscillation-or even several oscillationsprobe a whole domain, and if it is described unphysically like in the DLSHE model then the results come out unphysical as well. Manifestly, this would be a disaster for the VHE observatories of the next generations, which will reach energies up to 100 TeV or even larger.
This problem can be solved by smoothing out the edges in order to make the change of the magnetic field B direction continuous across the domain edges, even if it is still random, as already stressed in Section 4.4. Hence, in both cases only a random single realization of the beam propagation process is observable at once. We still suppose that photon-ALP oscillations are present in the beam from a blazar at redshift z, and so the photon survival probability is denoted by P ALP γ→γ E 0 , z; φ(y), θ(y) , where φ(y) and θ(y) are the two angles that fix the direction of B(y) in space at a generic point y along the beam and perpendicularly to it. In order to achieve our goal we have to resort to a domain-like smooth-edges (DLSME) model-mentioned in Section 4.4-wherein the beam propagation equation within a single domain becomes threedimensional and very difficult to solve analytically. But as shown in [114] such an equation becomes effectively two-dimensional. Moreover, according the above two models [213,215] the strength of B should vary rather little in different domains, hence we average it over many domains and attribute in first approximation the resulting value to each domain, denoting it for simplicity again by B. Finally, we consistently we take the transverse magnetic field component The two-dimensional beam propagation equation has been solved exactly and analytically [114]. It turns out that such a solution is indistinguishable from the numerical solution of the above three-dimensional exact equation (more about this in [114]). Physically, this amounts to the the whole physics of the problem being confined inside the planes Π(y) perpendicular to the beam rather than being spread out throughout the full three-dimensional space. As a consequence, P ALP γ→γ E 0 , z; φ(y), θ(y) → P ALP γ→γ E 0 , z; φ(y) , where φ(y) is the angle between B T (y) and a fixed fiducial z-direction equal in all domains (namely in all planes Π(y)).

Preliminary Remarks
Broadly speaking, what we said in Section 3 remains unchanged, apart from two facts. One is that the mixing matrix depends on y also in a single domain, and its explicit form is The meaning of the terms of M(E, y) is as follows. The contribution from photon dispersion on the CMB is ∆ CMB (E) = 0.522 · 10 −42 E [145], the contribution from the EBL absorption is ∆ abs (E) = i/ 2λ γ (E) where λ γ (E) denotes the corresponding photon mean free path inside a single domain (more about this, later), the contribution from the plasma frequency of the ionized intergalactic medium is ∆ pl (E) = − ω 2 pl /(2E) while the remaining terms are ∆ aγ = g aγγ B T /2 and ∆ aa (E) = − m 2 a /(2E), just like in Section 3.2. The other fact is that we now have three regimes, separated by the low-energy threshold which is equal to Equation (37), and the high-energy threshold E H ≡ 1.92 · 10 42 g aγγ B T 3.74 · 10 2 ξ GeV .
Specifically we have: • E < E L -This is the low-energy weak-mixing regime, wherein the terms ∝ E −1 dominate. Correspondingly, we have and However, since we will not consider this case throughout the paper.
• E L < E < E H -This is the intermediate-energy or strong-mixing regime in which the E = constant term dominates. Accordingly, we obtain L osc 2π g aγγ B T 2.05 · 10 2 ξ −1 Mpc , (125) P γ→a (L dom ) sin 2 g aγγ B T L dom 2 sin 2 1.54 · 10 −2 ξ L dom Mpc . (126) Clearly, L osc and P γ→a (L dom ) are independent of all the energy-dependent terms, and P γ→a (L dom ) becomes maximal: observe that m a enters E L and nowhere else. • E > E H -This is the high-energy weak-mixing regime, which is in a sense a sort of reversed low-energy weak-mixing regime, where however the term 0.522 · 10 −42 E dominates over g aγγ B T . Correspondingly, we get Manifestly, L osc (E) decreases with increasing E and P γ→a (E, L dom ) exhibits oscillations in E: this means that the individual realizations of the beam propagation are also oscillating functions of E. Moreover-since P γ→a (E, L dom ) ∝ E −2 -as E increases the photon-ALP oscillations become unobservable at some point.

Domain-like Smooth-Edges (DLSME) Model
As we said, we are going to apply the DLSME model to a monochromatic photon/ALP beam of energy E emitted by a far-away blazar, propagating through extragalactic space and reaching us [115].
Therefore we briefly summarize this model (see [114] for a full account). We suppose that there are N d domains between the blazar and us, and we number them so that domain 1 is the one closest to the blazar while domain N d is the one closest to us. Momentarily, we take all domains with the same length. We denote by {y D,n } 0≤n≤N d the set of coordinates which defines the beginning (y D,n−1 ) and the end (y D,n ) of the n-th domain (1 ≤ n ≤ N d ) towards the blazar.
Because of our ignorance about the strength of B in every domain and since it is supposed to vary rather little in different domains, we average B over many domains, and next we attribute the resulting value to each of them, so that B-but not B-will henceforth be regarded as constant in first approximation.
As we said the problem is actually a two-dimensional one, since what matters is only B T (y). Therefore, we denote by {φ n } 1≤n≤N d the set of angles that B T (y) forms with the fixed fiducial z-direction in the middle of every domain. Thanks to the previous assumptions also B T -but not B T -can be taken as constant in all domains.
Given the fact that B T (y) changes randomly from one domain to the next, in order for B T (y) to be continuous all along the beam it is compelling that it has equal values on both sides of every edge, e.g., the one between the n-th and the (n + 1)-th domain. Thus, the emerging picture is that B T (y) is homogeneous in the central part, but as the distance from the edge with the (n + 1)-th domain decreases we assume that B T (y) linearly changes thereby becoming equal to B T (y) on the same edge but in the (n + 1)-th domain. Accordingly, the continuity of the components of B T (y) along the whole beam is ensured.
A schematic view of this construction is shown in Figure 15. In practice, it is useful to define the two quantities y 0,n and y 1,n as where σ ∈ [0, 1] is the smoothing parameter. The interval [y 0,n , y 1,n ] is the region where the angle φ(y) changes smoothly from the value φ 0,n ≡ φ n in the n-th domain to the value φ 0,n+1 ≡ φ n+1 in the (n + 1)-th domain. Manifestly, for σ = 0 we have y 0,n = y 1,n , and we recover the DLSHE model. On the other hand, for σ = 1 then y 0,n becomes the midpoint of the n-th domain, and likewise y 1,n becomes the midpoint of the (n + 1)-th domain: now the smoothing is maximal, because we never have a constant value of φ in any domain. The general case is intermediaterepresented by a value of 0 < σ < 1-so that in the central part of a domain the angle is constant (φ 0,n ) and then it linearly joins the value of the constant angle in the next domain (φ 1,n ). Hence, in a generic interval [y 1,n−1 , According to our conventions, the blazar redshift is z ≡ z 0 , the points y D,n−1 and y D,n defining the n-th domain have redshift z n−1 and z n (z n < z n−1 ), respectively, and we set z n ≡ z n−1 + z n /2 for the average redshift of the n-th domain. Similarly, the emitted beam has energy E 0 , whereas the beam at points y D,n−1 and y D,n has energy E n−1 and E n (E n−1 > E n ), respectively. Finally, we define the average energy in the n-th domain as E n ≡ E n−1 + E n /2, and the observer has energy E N d . As usual, E n = (1 + z n )E N d (E 0 = (1 + z 0 )E N d ). We stress that at variance with the previous conventions, the observed beam energy is E N d .

Propagation over a Single Domain
We proceed along the same lines of Section 5.2, namely we have to account for the EBL absorption and to determine the magnetic field strength B T,n in the generic n-th domain of size L dom,n .
The only novelty is that instead of taking the length of all domains strictly equal, we allow for a small spread. Thus, we assume a probability distribution for the {L dom,n } 1≤n≤N . Owing to the properties of B at redshift z = 0, we take for the probability density in question the power law ∝ L −1.2 dom inside the range 0.2 Mpc-10 Mpc, which means that L dom = 2 Mpc, indeed allowed by present bounds [209]. Needless to say, such a choice is largely arbitrary and the corresponding histogram is shown in Figure 16. In order to accomplish this task, we just employ the discussion reported in Section 5.2. Accordingly, we find with τ γ again given by the FR model, and where B T,N d (y) is the strength of B T (y) in the local Universe, namely in the domain closest to the observer (z = 0). So, the above mixing matrix M(E, y) in a single n-th domain is fully determined, and now has E → E n = E N d 1 + z n and all terms replaced with those evaluated within the n-th domain. It can next be inserted into the reduced Schrödinger-like Equation (19).

Solution of the Beam Propagation Equation
Our present job is two-fold. In the first place, we have to solve such a reduced Schrödinger equation, which amounts to find its transfer matrix in a single domain. This is the hard part of the game, which is actually the main achievement of [114]. Next, the overall transfer matrix emerges by properly multiplying the transfer matrices pertaining to all domains, which is a trivial implication of quantum mechanics.
The first one amounts to solving the beam propagation equation inside a single n-th domain. The solution reads U n E n ; z n , z n−1 ; φ 1,n , φ 1,n−1 = (134) = U var,n E n ; z n , z n−1 ; φ 1,n , φ 0,n U const,n E n ; z n , z n−1 ; φ 0,n , for an arbitrary choice of the angle φ 0,n . Unfortunately, the explicit forms of the two transfer matrices in Equation (134) is much too cumbersome to be reported here, and the reader can found them in [114] (see its Equations (54) and (91) with E → E n and the appropriate conversions in order to go over from physical space to redshift space). The second point consists in the evaluation of the whole transfer matrix from the blazar to us, namely along a single arbitrary realization of the whole beam propagation process. Starting from Equation (134) the desired equation presently has the form U var,n E n ; z n , z n−1 ; φ 1,n , φ 0,n U const,n E n ; z n , z n−1 ; φ 0,n .
Note that this product must be ordered in such a way that the transfer matrices with smaller n must be closer to the source.

Results
We are finally in a position to derive the desired result, namely the photon survival probability from the blazar to us along an arbitrary realization of the whole beam propagation process. This goal is again trivially achieved thanks to the analogy with non-relativistic quantum mechanics, namely by employ again Equation (84).
As before, the photon polarization cannot be measured at the considered VHE, hence we have to start with the unpolarized beam state and sum the result over the two final polarization states. So, for the reader's convenience we revert to the same, common notation used in Section 5.2, namely E N d → E 0 , z N d → 0, z 0 → z. Accordingly, we have Below, the photon survival probability P ALP γ→γ,unp E 0 ; ρ x , ρ z ; z, ρ unp ; {φ n } 1≤n≤N is plotted versus the observed energy E 0 for 7 simulated blazars at z = 0.02, 0.05, 0.1, 0.2, 0.5, 1, 2, assuming for each one our benchmark values ξ = 0.5, 1, 2, 5. In order to simplify notations, we will denote P ALP γ→γ,unp E 0 ; ρ x , ρ z ; z, ρ unp ; {φ n } 1≤n≤N as P ALP γ→γ (E 0 , z). Thousand random realizations of the propagation process have been considered for each choice of z, ξ, E 0 . In all figures a random distribution of the domain length L dom has been taken: a power law distribution function has been chosen with exponent α = −1.2 and domain length in the interval between the minimal value L min dom = 0.2 Mpc and the maximal value L max dom = 10 Mpc. The resulting average domain length is L dom = 2 Mpc. Our results are shown in Figures 17-23. We take the smoothing parameter σ = 0.2 for the transition between two adjacent magnetic domains. The dotteddashed black line corresponds to conventional physics, the solid light-gray line to the median of all realizations of the propagation process and the solid yellow line to a single realization with a random distribution of both the domain lengths and of the orientation angles of the magnetic field inside the domains. The filled area represents the envelope of the results on the percentile of all the possible realizations of the propagation process at 68 % (dark blue), 90 % (blue) and 99 % (light blue), respectively. In the upper-left panel we have taken ξ = 0.5, in the upper-right panel ξ = 1, in the lower-left panel ξ = 2 and in the lower-right panel ξ = 5 [115]. A rather similar approach has been developed in 2017 by Kartavtsev, Raffelt and Vogel, where however only the average photon survival probability is considered [102].

A full Scenario
Up until this point we have especially addressed two specific topics which provide two strong hints at the existence of an ALP with m a = O(10 −10 ) eV and g aγγ = O(10 −11 ) GeV −1 .
In addition, we have considered the propagation of a photon/ALP beam in extragalactic space both for VHE energies currently observed by the IACTs (Section 5) and for energies to be measured by the next generation of VHE observatories (Section 8).
A partial scenario-complementary to the one discussed in Section 5-has been put forward in 2008 and consists in the conversion γ → a inside a blazar and the reconversion a → γ in the Milky Way, neglecting any possible γ ↔ a oscillation in extragalactic space [56].
Our present aim is to discuss a full scenario wherein a VHE photon/ALP beam is described from its origin inside a BL Lac to its detection on Earth. An early attempt towards this goal was done in 2009 [62], but since then much progress has been done. So, our analysis will be performed in the light of the most up-to-date astrophysical information and for energies up to above 50 TeV [116].
We are going to consider three sources.
• Markarian 501 at z = 0.034. • The extreme BL Lac 1ES 0229+200 at z = 0.1396. • A simulated source like BL Lac 1ES 0229+200 but at z = 0.6 Recall that BL Lacs have been observed up to z 1. Manifestly, the emitted VHE photon/ALP beam from the BL Lacs in question crosses a variety of magnetic field structures in very different astrophysical environments: inside the BL Lac jet, within the host galaxy, in extragalactic space, and finally inside the Milky Way. Accordingly, we shall have to evaluate the transfer matrix in each of these structures.
Our strategy is to assume a realistic emitted spectrum for the considered three BL Lacs, and derive their observed spectrum up to above 50 TeV.

Propagation in the BL Lac Jet
We denote by R VHE the region where the VHE photons originate inside the BL Lac jet, with y VHE denoting its distance from the central supermassive black hole (SMBH). So, our first step is to evaluate the transfer matrix in the jet region R jet between y VHE and the end of the jet y jet , which we denote as U R jet (E; y jet , y VHE ).
The region R VHE is rather far from the central SMBH, and the jet axis is supposed to coincide with the direction y (as usual). In order to evaluate the photon/ALP beam propagation inside the jet we must know three quantities: (1) the distance y VHE from the central SMBH, (2) the transverse magnetic field profile B T,R jet (y) from y VHE to y jet , (3) the electron density profile n e,R jet (y) from y VHE to y jet . The Synchrotron Self Compton (SSC) diagnostics as applied to the SED of BL Lacs [249] allows us to derive realistic values for these quantities. Inside R VHE we find 0.1 G B T,R VHE 1 G and for definiteness we choose B T,R VHE = 0.5 G. Moreover, we get n e,R VHE 5 · 10 4 cm −3 , leading in turn to a plasma frequency of ω pl 8.25 · 10 −9 eV, thanks to Equation (20). Although there is no direct way to infer a precise value of y VHE , this quantity can be estimated from the size of R VHE -which is assumed as measure of the jet cross-section-thus finding 10 16 cm y VHE 10 17 cm. For definiteness, we shall take y VHE 3 · 10 16 cm. Once produced, VHE photons propagate unimpeded out to y jet 1 kpc where they leave the jet, entering the host galaxy. Within R jet , what is relevant is the toroidal part of the magnetic field which is transverse to the jet axis [245,250,251]. Its profile is Concerning the electron density profile, due to the conical shape of the jet our expectation is n e,R jet (y) = n e,R VHE y VHE y 2 .
The knowledge of the above quantities allows us to compute the entire propagation process of the photon/ALP beam within the jet, namely U R jet (E; y jet , y VHE ).
It should be kept in mind that in R jet we consider the photon/ALP beam in a frame co-moving with the jet, so that we must apply the transformation E → γE to the beam in order to go to a fixed frame-as it will be performed in the next regions-with γ being the Lorentz factor. We take γ = 15.

Propagation in the Host Galaxy
All the three considered blazars are hosted by elliptical galaxies, which we denote by R host . We have already addressed the propagation of a VHE beam in these galaxies in Section 7.6, finding that even if the beam is in the strong-mixing regime the effect of γ ↔ a oscillations is totally negligible. Therefore, denoting by y in,host ≡ y jet and by y out,host the points on the y axis where the beam enters and exits from the host galaxy, respectively, we have U R host (E; y out,host , y in,host ) = 1.

Propagation in the Extragalactic Space
We let R ext be the region where the photon/ALP beam propagates in the extragalactic space, i.e., from y out,host up to the border of the Milky Way y MW . Now the beam behaviour in R ext is affected by the morphology and strength of the extragalactic magnetic field B ext . We have already repeatedly considered this issue in great detail, and for the present purposes in Section 8.

Propagation in the Milky Way
We denote by R MW the region where the photon/ALP beam propagates inside the Milky Way, i.e., from y MW up to the Earth, whose position is denoted by y ⊕ .
By closely following the strategy described in [74], we compute U R MW (E; y ⊕ , y MW ). In order to take into account the structured behaviour of the Galactic magnetic field B MW we adopt the recent Jansson and Farrar model [252,253], which includes a disk and a halo component, both parallel to the Galactic plane, and a poloidal 'X-shaped' component at the galactic center. Its most updated version is described in [254], where newer polarized synchrotron data and use of different models of the cosmic ray and thermal electron distribution are employed.
The alternative model of the Galactic magnetic field existing in the literature is the one in [255]. However this model is based mainly on data along the Galactic plane so that the Galactic halo component of B MW is not accurately determined. For this reason we prefer to use the Jansson and Farrar model. In any case, we have tested the robustness of our results by employing also this model and-even if with some little modifications-our results are qualitatively unchanged.
While the Jansson and Farrar model allows also for a random and a striated component of the field, it turns out that only the regular component is relevant in the present context, since the γ ↔ a oscillation length is much larger than the coherence length of the turbulent field.
Inside the Milky Way disk the electron number density is n e 1.1 · 10 −2 cm −3 , resulting in a plasma frequency ω pl 3.9 · 10 −12 eV owing to Equation (20): this emerges from a new model for the distribution of the free electrons in the Galaxy [256]. Moreover, the Galaxy is modeled by an extended thick disk accounting for the so-called warm interstellar medium, a thin disk standing for the Galactic molecular ring, spiral arms (inferred from a new fit to Galactic HII regions), a Galactic Center disk and seven local features counting the Gum Nebula, the Galactic Loop I and the Local Bubble. The model includes also an offset of the Sun from the Galactic plane and a warp of the outer Galactic disk. The Galactic model parameters are obtained from the fit to 189 pulsars with independently determined distances and DMs.
Accordingly, we compute U R MW (E; y ⊕ , y MW ) for an arbitrary direction of the line of sight to a given blazar.

Overall Photon Survival Probability
Because all the transfer matrices in each region are now known, the total transfer matrix U (E; y ⊕ , y VHE ) describing the propagation of the photon/ALP beam from the VHE photon production region in the BL Lac jet up to the Earth reads where of course we have y in,host ≡ y jet and z. Since photon polarization cannot be measured in the VHE gamma-ray band, we have to treat the beam as unpolarized. Therefore, we must use the generalized polarization density matrix ρ(y) = (A x (y), A z (y), a(y)) T ⊗ (A x (y), A z (y), a(y)). As a consequence, the overall photon survival probability becomes and Below-merely for notational convenience-we shall replace P ALP γ→γ E; y ⊕ , ρ x , ρ z ; y VHE , ρ unp simply by P ALP γ→γ E, z . In order to give the reader a feeling of what happens in the various regions crossed by the photon/ALP beam, in Figure 24 we plot how the oscillation length L osc varies as a function of the energy E in the jet, in the extragalactic space and in the Milky Way. As the upper panel of Figure 24 shows, the behaviour of L osc versus E is strongly affected by the value of B T,R jet (y): as expected (see also [114]), as B T,R jet (y) decreases-when the distance from the emission region increases-the maximal value of L osc increases and the energy where the QED vacuum polarization effect becomes important increases as well. Instead, in the central panel of Figure 24 what happens in extragalactic space is that L osc starts to decrease because of the effect of the photon dispersion on the CMB, which becomes more and more important as E increases (for more details see [114]). Finally, in the lower panel of Figure 24 we see that in the Milky Way L osc is almost constant with E since the QED vacuum polarization effect and photon dispersion on the CMB become subdominant as compared to the photon-ALP mixing in almost all the considered energy range.

Blazar Spectra
Starting from the intrinsic spectra, we are now in a position to employ the overall photon survival probability in order to derive the observed spectra of the considered three blazars, and from them to infer the corresponding SED νF ν in the presence of γ ↔ a oscillations all the way from inside the blazar to us. We can thus compare our findings with the results from conventional physics.
The observable physical quantity is the blazar spectrum pertaining to a single random realization of the photon/ALP propagation process. Nevertheless, it is enlightening to contemplate several realizations at once and to compute some of their statistical properties-the median and the area containing the 68%, 90% and 99% of the total number of realizations-in order to check the stability of the result against the distribution of the angles of the extragalactic magnetic field inside each domain with respect to a fixed fiducial direction z equal for all domain. We recall that these angles are independent random variables.
For the three blazars in question, we model their intrinsic spectrum with a power law exponentially truncated at a fixed cut-off energy E cut as where Φ 0 is a normalization constant accounting for the blazar luminosity, E 0 is a reference energy and Γ is the spectral index.
• Markarian 501-This source is a high-frequency peaked blazar (HBL) at redshift z = 0.034. We use the observational data points from HEGRA [257] in a condition where Markarian 501 was observed in a high emission state, which allows us to have a very good quality spectrum up to ∼30 TeV. This fact is important for testing our model, since at such high energies it starts to make predictions which depart from conventional physics. In Figure 25 we report its observed SED when conventional physics alone is considered, and when γ ↔ a oscillations are at work. In order to obtain the SED we take E cut = 10 TeV, E 0 = 1 TeV and Γ = 1.8 in Equation (145). • 1ES 0229+200-This is a BL Lac at redshift z = 0.1396. This is the prototype of the so-called 'extreme HBL' (EHBL) [258,259], which exhibit a rather hard VHE observed spectrum up to at least 10 TeV. This fact is particularly interesting since the observed data points at such high energies allow to distinguish between the models based on conventional physics and those containing γ ↔ a oscillations. Future observations with the CTA that can eventually reach energies up to 100 TeV can provide a definitive answer. In Figure 26 we plot its observed SED both when only conventional physics is taken into account and in the case in which also γ ↔ a oscillations are present. The SED is obtained by taking in Equation (145) E cut = 30 TeV in the case of conventional physics, and E cut = 10 TeV when also γ ↔ a oscillations are considered, while in either case we choose E 0 = 1 TeV and Γ = 1.4. Note that Γ is in agreement with the one derived for the Fermi/LAT spectrum in the recent analysis of [259]. • Extreme BL Lac at z = 0.6-BL Lacs have been observed also at redshift z ≥ 0.6, and so we assume the existence of an EHBL at redshift z = 0.6. For this blazar we take a SED similar to the one of 1ES 0229+200, namely E cut = 30 TeV, E 0 = 1 TeV and k = 1.4 in Equation (145) for both cases (conventional physics alone, presence of γ ↔ a oscillations). We consider two possibilities: (1) such BL Lac is observed in the sky along the direction of the galactic pole: in Figure 27 we plot its observed SED for both cases of presence/absence of photon-ALP interaction; (2) in Figure 28 we exhibit the corresponding observed SED for the same BL Lac instead observed in the sky along the direction of the galactic plane for both cases of presence/absence of photon-ALP interaction.   Figure 25 but for 1ES 0229+200. The dark gray squares are the spectrum detected by Fermi/LAT [260] while the light gray squares are the spectrum observed by HESS [261]. (Credit [116]).

Results
Our results about the SED of the above-considered BL Lacs are exhibited in Figures 25 and 28. Generally speaking, the γ ↔ a oscillations give rise to a harder observed spectrum for all three sources as compared to the outcome of conventional physics. We stress that this fact becomes increasingly evident as E or z (or both) get larger.
Our findings strongly suggest that γ ↔ a oscillations inside the magnetic field of the BL Lac jet play a key role in starting the propagation in extragalactic space with a sizable amount of already produced ALPs, whose relevance depends both on E and on z. This is a rather subtle point and deserves a clear explanation. Superficially, one might expect P ALP γ→γ (E, z) to increase with g aγγ , according to the physical intuition. This is certainly true as long as the EBL does not play an important role, namely for E and z low enough. Needless to say, γ → a conversions in the BL Lac and a → γ back-conversions in the Milky Way help increasing P ALP γ→γ (E, z), but not that much. Suppose instead that both g aγγ and z are fairly large but that E is not, so that photon dispersion on the CMB can be discarded. Accordingly, the conversion probability increases so that inside each single magnetic domain many γ → a and a → γ conversions take place. But since z is supposed to be fairly large the EBL level is high, so that most of the photons get absorbed. Such a behaviour is very clearly shown in Figures 27 and 28 around E 3 TeV. As the energy increases, photon dispersion on the CMB becomes dominant, which causes a much smaller number of γ → a and a → γ conversions to occur in extragalactic space. By and large, most of the ALPs produced in the BL Lac survive until they enter the Galaxy, whose strong magnetic field allows them to convert to photons. Whence the peak in Figures 27 and 28 around E = (10-30) TeV. All figures show that-as E progressively increases beyond 70 TeV-the area covered by the realizations of the photon/ALP propagation process gradually reduces. The reason is that the EBL absorption becomes so high at those energies that almost all the photons in each extragalactic magnetic field domain are absorbed and only the ones reconverted from ALPs inside the Galaxy are observed (as previously mentioned). Therefore, the parameter space of the model-B ext orientation angles, domain lengths L dom -gets reduced, and this fact reduces the available area that can be covered by the realizations of the propagation process.
Observe that in all figures we draw the CTA sensitivity curve for the South site and 50 h of observation. Because the sensitivity curves are relay on conservative criteria [262,263] we expect that the theoretical spectral features-look at the peak in Figures 27 and 28 around E ∼ 20 TeV-which are close to the sensitivity curve should anyhow be detectable by the CTA. * * * For the reader's convenience, we would like to briefly summarize our results. We have investigated the propagation of a photon/ALP beam originating well inside a BL Lac jet and traveling in the jet magnetic field, in the host galaxy magnetic field, in the extragalactic magnetic field, and in the Milky Way magnetic field up to us. We see from Markarian 501 (see Figure 25) that conventional physics does not fit the highest energy point of the SED while the model including γ ↔ a oscillations naturally matches the data. For 1ES 0229+200 (see Figure 26) the model including γ ↔ a oscillations fits the data remarkably well, especially the highest energy data points of the SED. As far as the simulated 1ES 0229+200 at z = 0.6 is concerned, the situation is striking: only γ ↔ a oscillations predict the peak around E = (10-30) TeV of the SED, while conventional physics prediction is many order of magnitude below.
As it is evident from Figures 27 and 28-as the redshift increases-at high energies the difference between the results from conventional physics alone, and the model including γ ↔ a oscillations becomes more and more dramatic. This is especially true when sizable γ → a conversions take place inside a blazar, since then most of the emitted ALPs can become photons only inside the Milky Way magnetic field. In particular, for very distant BL Lacs we predict a peak in the energy spectra at E = (10-30) TeV as it is evident from Figures 27 and 28 for a BL Lac at z = 0.6. In addition, the energy oscillations in the observed spectrum-clearly recognizable in the figures-are a clear-cut feature of our scenario, which can be observed provided that the detector has enough energy resolution: they arise from the photon dispersion on the CMB.
A competitive scenario capable to reduce the optical depth is the Lorentz invariance violation (LIV) which could predict a somewhat similar peak in the BL Lac spectra above ∼20 TeV [264,265]. But the two scenarios can be distinguished since the LIV does not predict any spectral energy oscillatory behavior [117].
At this point some remarks look compelling.
• The jet parameters (y VHE , B T,R VHE ) are affected by uncertainties, and the amount of produced ALPs in this region clearly reflects this fact. Nevertheless, we have checked that the final spectra qualitatively possess the above-mentioned features regardless of the choice of the jet parameters, provided of course that they are realistic. • Even if we consider very low values of the extragalactic magnetic field-namely B ext 10 −9 G-the considered model predicts the above-mentioned features even if partially reduced, in particular concerning the amplitude of the energy oscillations. However, the peak in the spectra at E = (10-30) TeV remains unaffected at high redshift. • The electromagnetic cascade proposed to mimic γ ↔ a oscillation effects in blazar spectra [266] can work only for B ext O(10 −15 ) G, which is indeed quite close to the B ext lower limits [208][209][210]. Still, for B ext O(10 −15 ) G the charged particles produced in the cascade are deflected by B ext and the resulting additional photon flux turns out to be very likely irrelevant (for more details, see e.g., [267]). This argument also applies to the possible additional e + , e − pairs produced in the process γ + γ → e + + e − . • For E 100 TeV the infrared radiation from dust present inside the Milky Way could play a moderate role in absorbing photons [268]. But this effect is irrelevant for us and can be safely discarded. Basically, the resulting absorption is substantial only inside the Galactic plane and a few degrees above and below it, hence only ALPs converted to photons in the Galactic plane close to the outer border of the Milky Way disk fully undergo such an effect. Actually, two points should be be stressed. (1) It goes without saying that when the line of sight to the blazar lies outside the galactic plane the considered effect is totally irrelevant.
(2) Even for the photon/ALP beam entering the Milky Way along the Galactic plane the γ ↔ a oscillations reduce photon absorption, thereby considerably weakening dust absorption.

Polarization Effects
As stressed in the Introduction and in Section 2.3, not only the photon-ALP interaction produces γ ↔ a oscillations in the presence of an external magnetic field-resulting in several consequences for astrophysical spectra (transparency modification, flux excess, spectrum irregularities)-but it also gives rise to the change of the polarization state of photons. Less attention has been paid to the latter effect in the literature so far. Yet, it gives rise to effects which are potentially detectable from current and planned satellite missions.

ALP Effects on Photon Polarization
In order to describe the consequences of photon-ALP interaction on the final photon polarization we employ the generalized density matrix ρ(y) associated with the photon-ALP system defined by Equation (46). It can be specialized to describe pure photon states in the x and z direction, which can be expressed by respectively, the ALP state reading and unpolarized photons represented by Instead, the polarization density matrix characterizing partially polarized photons shows an intermediate functional expression between Equations (146) and (148).
Once the transfer matrix of the photon-ALP system U is known, the final polarization density matrix ρ can be computed by means of Equation (48) when the initial polarization density matrix ρ 0 is specified. It is useful to express the photonic part of ρ in terms of the Stokes parameters as [269] ρ γ = 1 2 The photon degree of linear polarization Π L is defined as [270] Π L ≡ (Q 2 while the polarization angle χ reads [270] It is simple algebra to express Equations (150) and (151) in terms of the elements of the ρ matrix ρ ij with i, j = 1, 2 as and respectively. The photon degree of linear polarization at emission Π L,0 is linked to the photon conversion and survival probabilities P ALP γ→a and P ALP γ→γ in the presence of the photon-ALP interaction. In particular, some theorems stated and proven in [127] show what follows.

2.
Under the previous conditions, Π L,0 can be viewed as the measure of the overlap between the values assumed by P ALP γ→a and P ALP γ→γ . If photons are emitted unpolarized (Π L,0 = 0), then P ALP γ→a and P ALP γ→γ have no overlap apart from the value 1/2, at most. From item 2 we envisage that photon-ALP interaction can be used to measure emitted photon degree of linear polarization when photon absorption is absent. In order for this strategy to be implemented, the photon-ALP system must be in the weak mixing regime: since P ALP γ→a and P ALP γ→γ possess an oscillatory behavior in energy, their values oscillate within the bounds ensured by the above theorems and Π L,0 can be inferred. In particular, from an observed astrophysical spectrum Φ obs we can extract P ALP γ→γ as where Φ em is the emitted spectrum, which is supposed either known or derivable from Φ obs (for more details see [127]). Moreover, we have P ALP γ→a = 1 − P ALP γ→γ since there is no photon absorption. Now, Π L,0 is simply given by the measure of the interval where P ALP γ→γ and P ALP γ→a overlap, as Figure 29 shows for a typical shape of P ALP γ→γ and P ALP γ→a . Such a procedure represents the only known possibility to measure the initial polarization of photons emitted by astrophysical sources, since all other methods can only detect the final Π L . In addition, only flux observations are needed.

ALP-Induced Polarization Effects in Galaxy Clusters
We now address the impact of photon-ALP interaction on photon polarization for galaxy clusters. We combine the previous results obtained concerning the photon/ALP beam propagation in the different astrophysical backgrounds (galaxy cluster, extragalactic space, Milky Way).
As discussed in Section 4.5, photons produced inside nCC regular clusters are emitted unpolarized in each energy band, so that they have Π L,0 = 0. The photon/ALP beam propagates inside the cluster, in the extragalactic space and in the Milky Way. While crossing these magnetized media, photon-ALP interaction produces a modification in the final Π L and χ. Concerning the parameters of the photon-ALP system, in order to be specific we take g aγγ = 0.5 · 10 −11 GeV −1 and m a = 10 −10 eV. Since a diffuse photon emission in the galaxy cluster is considered, we assume a typical nCC regular cluster at a redshift z = 0.03. The cluster magnetic field profile B clu and the electron number density profile n clu e are given by Equations (63) and (64) of Section 4.5, respectively. We consider the following parameters entering Equations (63) and (64) and regarded as typical [224,225]: B clu 0 = 15 µG, k L = 0.2 kpc −1 , k H = 3 kpc −1 , n clu e,0 = 0.5 · 10 −2 cm −3 , η clu = 0.75, β clu = 2/3, r core = 100 kpc, and a cluster radius of 1 Mpc. By evaluating the photon/ALP beam propagation from the cluster central region up to the cluster radius, we compute the transfer matrix of the photon-ALP system inside the cluster U clu . The propagation of the photon/ALP beam in extragalactic space directly follows from Section 8. Thus, by taking an extragalactic magnetic field strength B ext = 1 nG with coherence length L ext dom in the range (0.2-10) Mpc and average L ext dom = 2 Mpc, we compute the transfer matrix U ext in this region. We recall that γ ↔ a oscillations in the Milky Way have been investigated in Section 9.4, which we closely follow. We consider the cluster located in the direction of the galactic pole, where the Milky Way magnetic field strength is minimal, so as to be conservative about γ ↔ a oscillations. Hence, we are in a position to evaluate the transfer matrix of the photon-ALP system in the Milky Way U MW .
By combining the previous transfer matrices in the correct order, we obtain the whole transfer matrix U tot associated with the propagation of the photon/ALP beam from the cluster core to us U tot = U MW U ext U clu , (155) while the final photon survival probability with γ ↔ a oscillations P ALP γ→γ is given by with ρ in ≡ ρ unpol reading from Equation (148) and ρ x and ρ z from Equation (146). Then, the corresponding final degree of linear polarization Π L emerges from Equation (152).
In the left panel of Figure 30 we report P ALP γ→γ in the MeV energy band for a particular realization of the photon/ALP beam propagation process, in the central panel of Figure 30 we plot the corresponding Π L , while in the right panel of Figure 30 we present the probability density function f P of Π L associated with several realizations at the benchmark energy E 0 = 10 MeV. From Figure 30 we observe that the photon/ALP beam propagates in the weak mixing regime, P ALP γ→γ and Π L show an oscillatory energy behavior and Π L > 0. Since Π L = 0 is never the most probable result, as the right panel of Figure 30 shows, we can conclude that a signal of Π L > 0 can be detected by the planned missions like COSI [271], e-ASTROGAM [272,273] and AMEGO [274]. Similar results have been obtained by considering Coma [129]. Note that P ALP γ→γ in Figure 30 satisfies theorems enunciated and demonstrated in [127] and recalled above.
We have also performed the same procedure for the photon/ALP beam propagation in the VHE range (see [128]). For energies where photon absorption due to the EBL is strong, we have found a feature: photons are fully polarized. The reason is as follows. Since all photons propagating in the extragalactic space are absorbed, only the ALPs reconverted back to photons inside the Milky Way can be detected. Since the Milky Way magnetic field component responsible for photon-ALP conversion is the regular part, photons are fully polarized. In case of a detection of photons completely polarized, this fact would represent a proof for ALP existence. However, observatories cannot measure Π L up to so high energies, yet [275].

ALP-Induced Polarization Effects in Blazars
In Section 4.1 we have described the blazar properties, which are important for the photon-ALP system. Concerning polarization we have seen that in the X-ray band photons are expected to be partially polarized, while in the MeV and VHE range they are expected to be emitted unpolarized with Π L,0 = 0. In particular, we consider BL Lacs and we take the same parameters of Section 9.1. Thus, we can calculate the transfer matrix of the photon/ALP beam in the blazar jet U jet . Concerning the photon/ALP beam propagation inside the host galaxy we closely follow what we have described in Section 7.6, so that we get the transfer matrix U host . The transfer matrices of the photon-ALP system in the galaxy cluster U clu , in the extragalactic space U ext and in the Milky Way U MW exactly read from the previous Section with the same parameters apart from n clu e,0 = 5 · 10 −2 cm −3 since we are considering now a CC galaxy cluster. By evaluating the latter three transfer matrices and combining them with U jet and U host in the correct order we can calculate the total transfer matrix U tot associated to the propagation of the photon/ALP beam from the jet base up to us while the final photon survival probability in the presence of photon-ALP interaction P ALP γ→γ reads from Equation (156) and the corresponding final degree of linear polarization Π L is given by Equation (152). The left panel of Figure 31 shows P ALP γ→γ in the MeV energy range for a peculiar realization of the photon/ALP beam propagation process, the central panel of Figure 31 exhibits the corresponding Π L , while the right panel of Figure 31 reports the probability density function f Π of Π L associated to several realizations at the benchmark energy E 0 = 10 MeV.  Figure 31 shows that the photon/ALP beam is in the weak mixing regime and P ALP γ→γ presents a pseudo-oscillatory behavior with respect to the energy. Correspondingly, we observe in the central panel of Figure 31 that Π L > 0 with a high variation. From the right panel of Figure 31, we infer that Π L = 0 is never the most probable result. Therefore, we can conclude that a signal of Π L > 0 can be observed by observatories such as COSI [271], e-ASTROGAM [272,273] and AMEGO [274]. The present situation is similar but even better with respect to photon production in the cluster central zone (see the previous Section). It is possible to verify that P ALP γ→γ in Figure 31 satisfies theorems enunciated and demonstrated in [127] and recalled above.
Furhtermore, in the present case of photon production at the blazar jet base, we have calculated the photon/ALP beam propagation in the VHE range (see [128]). We obtain the same results reported in the previous Section about the production of fully polarized photons at VHE, when photon-ALP oscillations are present.

Conclusions and Outlook
In the present Review we have tried to summarize the most important implications of ALPs for High Energy Astrophysics (in a broad sense). Two new strong hints at an ALP with m a = O 10 −10 eV and g aγγ = O 10 −11 GeV −1 have emerged.
On the other hand a direct detection can be achieved by means of the experiment called shining through the wall within the next few years, either using a laser at optical frequency such as in the ALP II [136] experiment at DESY or by employing a laser at radiofrequency such as in the planned STAX experiment [137]. Alternatively opportunities are the planned IAXO observatory [138] or the strategies developed by Avignone and collaborators [139][140][141]. Finally, if most of the dark matter is made of the considered ALPs, they can be discovered by the planned experiment ABRACADABRA [142].
Only time will tell!
where the coefficients s, t, u and v are supposed to be complex numbers. We start by diagonalizing M. Its eigenvalues are and it is straightforward to check that the corresponding eigenvectors can be taken to be Correspondingly, any solution of Equation (A1) can be represented in the form ψ(y) = c 1 X 1 e iλ 1 (y−y 0 ) + c 2 X 2 e iλ 2 (y−y 0 ) + c 3 X 3 e iλ 3 (y−y 0 ) , where c 1 , c 2 , c 3 and y 0 are arbitrary constants. As a consequence, the solution with initial condition emerges from Equation (A10) for It is a simple exercize to recast the considered solution into the form ψ(y) = U (y, y 0 ; 0) ψ(y 0 ) (A15) with U (y, y 0 ; 0) = e iλ 1 (y−y 0 ) T 1 (0) + e iλ 2 (y−y 0 ) T 2 (0) + e iλ 3 (y−y 0 ) T 3 (0) , where we have set from which it follows that the desired transfer matrix is just U (y, y 0 ; 0) as given by Equation (A16).

Appendix B
The key-point of the ALP scenario-already stressed in Section 3.6-is that ALPs do neither interact with the EBL-in spite of the fact that they couple to two photons-nor with the ionized intergalactic medium. The proof is as follows ALPs might interact with the EBL only through two processes: a + γ → a + γ and a + γ → f + f , where f denotes a generic charged fermion.
Consider first the process aγ → aγ represented by the Feynman diagram in Figure A1 in the s-channel. A simple estimate gives σ(aγ → aγ) ∼ s g 4 aγ , and enforcing the CAST bound g aγγ < 0.66 · 10 −10 GeV −1 we find σ(aγ → aγ) E γ E ALP /GeV 2 10 −68 cm 2 , which shows that this process is negligibly small for any reasonable choice of E γ and E ALP . Let us next turn our attention to the process aγ → f f represented by the Feynman diagram in Figures A2 in the s-channel. Accordingly we have σ(aγ → f f ) ∼ α g 2 aγγ , which-thanks to the CAST bound-yields σ(aγ → f f ) 10 −50 cm 2 . Finally, we address the process a f → γ f represented by the Feynman diagram in Figures A2 in the t-channel. Manifestly, we have again σ(a f → γ f ) ∼ α g 2 aγγ and so σ(a f → γ f ) 10 −50 cm 2 . Therefore, for all practical purposes ALPs neither interact with photons nor with any fermion.    Other processes discussed in [144] are totally irrelevant for the energy range considered in this paper.