Axion-like Particles Implications for High-Energy Astrophysics

Abstract

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.

1. Introduction

The Standard Model (SM) of strong, weak and electromagnetic interactions based on the gauge group $\mathrm{SU}{\left(3\right)}_C\otimes \mathrm{SU}{\left(2\right)}_W\otimes \mathrm{U}{\left(1\right)}_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.

Although it is presently impossible to tell which new proposal—out of so many ones—has any chance to successfully describe Nature, it seems remarkable that several attempts along very different directions such as four-dimensional supersymmetric models [1,2,3,4,5], multidimensional Kaluza-Klein theories [6,7] and especially M theory—which encompasses superstring and superbrane theories—predict the existence of axion-like particles (ALPs) [8,9] (for a very incomplete list of references, see [10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29]).

Basically, ALPs are very light, pseudo-scalar bosons—denoted by a—which mainly couple to two photons with a strength $g_{a\gamma \gamma }$ 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 $\mathbf{B}$. In such a situation we have $\gamma \to a$ and $a\to \gamma$ conversions, represented by the Feynman diagram in Figure 2. Note that $\mathbf{B}$ could be replaced by an external electric field $\mathbf{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 3—this means that $\gamma \leftrightarrow 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 $\mathbf{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.

• In a suitable region of the parameter plane $(m_a,g_{a\gamma \gamma })$ ALPs turn out to be very good candidates for cold dark matter [34,35,36,37,38].
• In another region of the parameter plane $(m_a,g_{a\gamma \gamma })$—which can overlap with the previous one—ALPs give rise to very interesting astrophysical effects (for a very incomplete list of references, see [39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112,113,114,115,116,117,118,119,120,121,122,123,124,125,126,127,128,129]. 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].
• The last reason is that the region of the parameter plane $(m_a,g_{a\gamma \gamma })$ 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\gamma \gamma }$.

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 $\gamma$-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 $\gamma \leftrightarrow a$ oscillations should take place in the beam as it propagates. Moreover, the extragalactic magnetic field $\mathbf{B}$ is modeled as a domain-like structure—which mimics the real physical situation—with the strength of $\mathbf{B}$ nearly equal in all domains, the size $L_{\mathrm{dom}}$ of all domains being taken equal and set by the $\mathbf{B}$ coherence length while the direction of $\mathbf{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 $\gamma \leftrightarrow a$ oscillation length $L_{\mathrm{osc}}$ is very much larger than $L_{\mathrm{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\mathrm{GeV}<E_0<100\mathrm{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 $\gamma +\gamma \to e^{+}+e^{-}$ [143], the beam undergoes a frequency-dependent attenuation1. However in the presence of $\gamma \leftrightarrow 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 ${\tau }_{\gamma }^{\mathrm{ALP}}(E_0,z)$ is smaller than in conventional physics. However since the corresponding photon survival probability is $P_{\gamma \to \gamma }^{\mathrm{ALP}}(E_0,z)=\exp \left(-{\tau }_{\gamma }^{\mathrm{ALP}}(E_0,z)\right)$, even a small decrease in ${\tau }_{\gamma }^{\mathrm{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 $\gamma +\gamma \to e^{+}+e^{-}$, the beam undergoes a frequency-dependent attenuation. Thus, if we want to know the emitted spectrum ${\Phi }_{\mathrm{em}}\left(E_0(1+z)\right)$ we have to EBL-deabsorb the observed one ${\Phi }_{\mathrm{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 ${\Phi }_{\mathrm{em}}\left(E_0(1+z)\right)\propto {\left(E_0(1+z)\right)}^{-{\Gamma }_{\mathrm{em}}\left(z\right)}$. 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 $\left\{{\Gamma }_{\mathrm{em}}\left(z\right)\right\}$ for all considered sources, in the absence of ALPs we find that the best-fit regression line is a concave parabola in the ${\Gamma }_{\mathrm{em}}-z$ plane. As a consequence, there is a statistical correlation between ${\Gamma }_{\mathrm{em}}\left(z\right)$ and z. However how can the sources get to know their z so as to tune their ${\Gamma }_{\mathrm{em}}\left(z\right)$ 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=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$ and $g_{a\gamma \gamma }=\mathcal{O}\left({10}^{-11}\right){\mathrm{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 $\gamma$-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 $\gamma \leftrightarrow 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\gtrsim 15\mathrm{TeV}$ photon dispersion on the CMB (Cosmic Microwave Background) becomes the leading effect, which causes the $\gamma \leftrightarrow a$ oscillation length to get smaller and smaller as E further increases. Therefore, things change drastically whenever $L_{\mathrm{osc}}\lesssim L_{\mathrm{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 $\mathbf{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 $\gamma \to 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 $\gamma \leftrightarrow 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.

2. 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\simeq {10}^{19}\mathrm{GeV}$ than to the Fermi scale $G_F^{-1/2}\simeq 250\mathrm{GeV}$.

2.1. 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 $\mathcal{O}\left({10}^{105}\right)$, 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 $\mathcal{O}\left(100\right)$ 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}\mathrm{eV}$ [16]. Below, we do not commit ourselves with any specific string theory.

After compactification at the very high scale $\Lambda$, the resulting four-dimensional field theory is described by a Lagrangian. We collectively denote by $\varphi$ the SM particles together with possibly new undetected particles with mass smaller than $G_F^{-1/2}$—the above ALPs are an example—while all particles much heavier than $G_F^{-1/2}$ are collectively represented by $\Phi$. Correspondingly, the Lagrangian in question has the form ${\mathcal{L}}_{\mathrm{FT}}(\varphi ,\Phi )$, and the generating functional for the corresponding Green’s functions reads

$$Z_{\mathrm{FT}}[J,K]=N\int \mathcal{D}\varphi \int \mathcal{D}\Phi \exp \left(i\int d^{4}x\left[{\mathcal{L}}_{\mathrm{FT}}(\varphi ,\Phi )+\varphi J+\Phi K\right]\right),$$

(1)

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_{\mathrm{FT}}[J,K]$, so that the low-energy effective Lagrangian ${\mathcal{L}}_{\mathrm{eff}}\left(\varphi \right)$ is defined by

$$\exp \left(i\int d^{4}x{\mathcal{L}}_{\mathrm{eff}}\left(\varphi \right)\right)\equiv \int \mathcal{D}\Phi \exp \left(i\int d^{4}x{\mathcal{L}}_{\mathrm{FT}}(\varphi ,\Phi )\right).$$

(2)

Evidently, the SM Lagrangian ${\mathcal{L}}_{\mathrm{SM}}\left({\varphi }_{\mathrm{SM}}\right)$ is contained in ${\mathcal{L}}_{\mathrm{eff}}\left(\varphi \right)$, and—in the absence of any new physics below $G_F^{-1/2}$—it will differ from ${\mathcal{L}}_{\mathrm{eff}}\left(\varphi \right)$ only by non-renormalizable terms involving the ${\varphi }_{\mathrm{SM}}$ particles alone, that are suppressed by inverse powers of $\Lambda$.

In any theory with a sufficiently rich gauge structure—which is certainly the case of ${\mathcal{L}}_{\mathrm{FT}}(\varphi ,\Phi )$—some further ALPs can arise. Indeed, some global symmetries $\mathcal{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 $\mathcal{G}$ undergoes spontaneous symmetry breaking as well. As a consequence, some Goldstone bosons—which are collectively denoted by a if $\mathcal{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 $\mathcal{G}$ is an invariance group of the FT. Because the instanton-like effect explicitly break $\mathcal{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 $\varphi$ into the set of SM particles ${\varphi }_{\mathrm{SM}}$ plus the ALPs a, the low-energy effective Lagrangian has the structure

$${\mathcal{L}}_{\mathrm{eff}}({\varphi }_{\mathrm{SM}},a)={\mathcal{L}}_{\mathrm{SM}}\left({\varphi }_{\mathrm{SM}}\right)+{\mathcal{L}}_{\mathrm{nonren}}\left({\varphi }_{\mathrm{SM}}\right)+{\mathcal{L}}_{\mathrm{ren}}\left(a\right)+{\mathcal{L}}_{\mathrm{ren}}({\varphi }_{\mathrm{SM}},a)+{\mathcal{L}}_{\mathrm{nonren}}({\varphi }_{\mathrm{SM}},a),$$

(3)

where ${\mathcal{L}}_{\mathrm{ren}}\left(a\right)$ contains the kinetic and mass terms of a, and ${\mathcal{L}}_{\mathrm{ren}}({\varphi }_{\mathrm{SM}},a)$ stands for renormalizable soft-breaking terms that can be present whenever $\mathcal{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 ${\mathcal{L}}_{\mathrm{SM}}\left({\varphi }_{\mathrm{SM}}\right)$.

Needless to say, it can well happen that between $G_F^{-1/2}$ and $\Lambda$ other relevant mass scales ${\left\{{\Lambda }_i\right\}}_{i=1,2,3,\dots }$ exist. In such a situation the above scheme remains true, but then $\mathcal{G}$ may be spontaneously broken stepwise at various intermediate scales.

Finally, we recall that pseudo-Goldstone bosons such as ALPs are necessarily pseudo-scalar particles [146].

2.2. Axion as a Prototype

A characteristic feature of the SM is that non-perturbative effects produce the term $\Delta {\mathcal{L}}_{\theta }=\theta g_3^2{G}_a^{\mu \nu }{\tilde{G}}_{a\mu \nu }/\left(32{\pi }^2\right)$ in the QCD Lagrangian, where $\theta$ is an angle, $g_3$ and $G_a^{\mu \nu }$ are the gauge coupling constant and the gauge field strength of ${\mathrm{SU}}_C\left(3\right)$, respectively, and ${\tilde{G}}_a^{\mu \nu }\equiv \frac{1}{2}{ϵ}^{\mu \nu \rho \sigma }{G}_{a\rho \sigma }$. All values of $\theta$ are allowed and theoretically on the same footing, but a nonvanishing $\theta$ 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 ${\mathcal{M}}_q$ into diagonal form, and so the total strong CP violation is parametrized by $\overline{\theta }=\theta +\arg \mathrm{Det}{\mathcal{M}}_q$ (for a review, see e.g., [147,148,149,150,151]). Observationally, a nonvanishing $\overline{\theta }$ would show up in an electric dipole moment $d_n$ for the neutron. Consistency with the experimental upper bound $|d_n|<3·{10}^{-26}\mathrm{e}\mathrm{cm}$ requires $|\overline{\theta }|<{10}^{-9}$ [147,148,149,150,151]. Thus, the question arises as to why $|\overline{\theta }|$ is so unexpectedly small. A natural way out of this fine-tuning problem—which is the strong CP problem— was proposed in 1977 by Peccei and Quinn [152,153]. Basically, the idea is to make the SM Lagrangian invariant under an additional global $\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$ symmetry in such a way that the $\Delta {\mathcal{L}}_{\theta }$ term can be rotated away. While this strategy can be successfully implemented, it turns out that the $\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$ is spontaneously broken, and so a Goldstone boson is necessarily present in the physical spectrum. Actually, things are slightly more complicated, because $\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$ is also explicitly broken by the same non-perturbative effects which give rise to $\Delta {\mathcal{L}}_{\theta }$. 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\simeq 0.6\left(\frac{{10}^7\mathrm{GeV}}{f_a}\right)\mathrm{eV},$$

(4)

where $f_a$ denotes the scale at which $\mathrm{U}{\left(1\right)}_{\mathrm{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\gamma \gamma$ of the axion a is generated at one-loop via the triangle graph with internal fermion lines, which is described by the effective Lagrangian

$${\mathcal{L}}_{a\gamma \gamma }=-\frac{1}{4}{g}_{a\gamma \gamma }{F}^{\mu \nu }{\tilde{F}}_{\mu \nu }a=g_{a\gamma \gamma }\mathbf{E}·\mathbf{B}a,$$

(5)

where $F^{\mu \nu }\equiv (\mathbf{E},\mathbf{B})\equiv {\partial }^{\mu }{A}^{\nu }-{\partial }^{\nu }{A}^{\mu }$ is the usual electromagnetic field strength and ${\tilde{F}}^{\mu \nu }\equiv \frac{1}{2}{ϵ}^{\mu \nu \rho \sigma }{F}_{\rho \sigma }$. The two-photon coupling $g_{a\gamma \gamma }$ entering Equation (5) has dimension of the inverse of an energy and is given by

$$g_{a\gamma \gamma }\simeq 0.8·{10}^{-10}{k}^{-1}\left(\frac{{10}^7\mathrm{GeV}}{f_a}\right){\mathrm{GeV}}^{-1},$$

(6)

with k a model-dependent parameter of order one [156]. Note that $g_{a\gamma \gamma }\propto 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.7k\left(g_{a\gamma \gamma }{10}^{10}\mathrm{GeV}\right)\mathrm{eV},$$

(7)

In the original proposal [152,153], $\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$ is spontaneously broken by two Higgs doublets which also spontaneously break ${\mathrm{SU}}_W\left(2\right)\otimes {\mathrm{U}}_Y\left(1\right)$, so that $f_a\simeq G_F^{-1/2}$. Correspondingly, from Equation (4) we get $m\simeq 24\mathrm{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 $f_a\gg G_F^{-1/2}$. This is straightforwardly achieved by performing the spontaneous breakdown of $\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$ with a Higgs field $\Phi$ which is a singlet under ${\mathrm{SU}}_W\left(2\right)\otimes {\mathrm{U}}_Y\left(1\right)$ [159,160], and with a vacuum expectation value $\langle \Phi \rangle \gg G_F^{-1/2}$.

We have told the axion story in some detail because the latter condition leads to the conclusion that the $\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$ symmetry has nothing to do with the low-energy effective theory to which the axion belongs—namely ${\mathcal{L}}_{\mathrm{nonren}}({\varphi }_{\mathrm{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 $\mathcal{G}=\mathrm{U}{\left(1\right)}_{\mathrm{PQ}}$, $f_a={\Lambda }_i$ (for some i) and ${\mathcal{L}}_{\mathrm{nonren}}({\varphi }_{\mathrm{SM}},a)$ including ${\mathcal{L}}_{a\gamma \gamma }$ 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 $\mathrm{U}{\left(1\right)}_{\mathrm{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.

2.3. 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\gamma \gamma$ is just supposed to exist without further worrying about its origin.
• The parameters $m_a$ and $g_{a\gamma \gamma }$ are to be regarded as unrelated for ALPs, and it is merely assumed that $m_a\ll 1\mathrm{eV}$ and $g_{a\gamma \gamma }^{-1}\gg G_F^{-1/2}$.

As a result, ALPs are described by the Lagrangian

$${\mathcal{L}}_{\mathrm{ALP}}^0=\frac{1}{2}{\partial }^{\mu }a{\partial }_{\mu }a-\frac{1}{2}{m}_a^2{a}^2-\frac{1}{4}{g}_{a\gamma \gamma }{F}_{\mu \nu }{\tilde{F}}^{\mu \nu }a=\frac{1}{2}{\partial }^{\mu }a{\partial }_{\mu }a-\frac{1}{2}{m}_a^2{a}^2+g_{a\gamma \gamma }\mathbf{E}·\mathbf{B}a,$$

(8)

where $\mathbf{E}$ and $\mathbf{B}$ have the same meaning as before. Note that ${\mathcal{L}}_{\mathrm{ALP}}^0$ is included in ${\mathcal{L}}_{\mathrm{ren}}\left(a\right)+{\mathcal{L}}_{\mathrm{nonren}}({\varphi }_{\mathrm{SM}},a)$.

Throughout this Review, we shall suppose that $\mathbf{E}$ is the electric field of a propagating photon beam while $\mathbf{B}$ is an external magnetic field. Because a couples to $\mathbf{E}·\mathbf{B}$, the effective coupling is proportional $\cos \beta$, where $\beta$ is the angle between $\mathbf{E}$ and $\mathbf{B}$, and we denote by ${\mathbf{B}}_T$ the $\mathbf{B}$-component transverse to the beam momentum. Thus, what matters is ${\mathbf{B}}_T$ and not $\mathbf{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 $\mathbf{E}$ or $\mathbf{B}$ is sufficiently large, also QED vacuum polarization effects should be taken into account, which are described by Heisenberg-Euler-Weisskopf Lagrangian [169,170]

$${\mathcal{L}}_{\mathrm{HEW}}=\frac{2{\alpha }^2}{45m_e^4}\left[{\left({\mathbf{E}}^2-{\mathbf{B}}^2\right)}^2+7{\left(\mathbf{E}·\mathbf{B}\right)}^2\right],$$

(9)

where $\alpha$ is the fine-structure constant and $m_e$ is the electron mass, so that ALPs are described by the full Lagrangian ${\mathcal{L}}_{\mathrm{ALP}}={\mathcal{L}}_{\mathrm{ALP}}^0+{\mathcal{L}}_{\mathrm{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 $\Phi$ (for a review, see e.g., [171] and the references therein).

3. General Properties of ALPs

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

3.1. Beam Propagation Equation

It can be shown that in this case the beam propagation equation following from ${\mathcal{L}}_{\mathrm{ALP}}$ can be written as [32]

$$\left(\frac{d^2}{d{y}^2}+E^2+2E{\mathcal{M}}_0\right)\psi \left(y\right)=0$$

(10)

with

$$\psi \left(y\right)\equiv \left(\begin{array}{c}{A}_x\left(y\right)\\ A_z\left(y\right)\\ a\left(y\right)\end{array}\right),$$

(11)

where $A_x\left(y\right)$ and $A_z\left(y\right)$ denote the photon amplitudes with polarization (electric field) along the x- and z-axis, respectively, while $a\left(y\right)$ is the amplitude associated with the ALP. It is useful to introduce the basis $\left\{\right|{\gamma }_x\rangle \equiv {(1,0,0)}^T,|{\gamma }_z\rangle \equiv {(0,1,0)}^T,|a\rangle \equiv {(0,0,1)}^T\}$, where $|{\gamma }_x\rangle$ and $|{\gamma }_z\rangle$ represent the two photon linear polarization states along the x- and z-axis, respectively, and $|a\rangle$ denotes the ALP state. Accordingly, we can rewrite $\psi \left(y\right)$ as

$$\psi \left(y\right)=A_x\left(y\right)|{\gamma }_x\rangle +A_z\left(y\right)|{\gamma }_z\rangle +a\left(y\right)|a\rangle ,$$

(12)

and the real, symmetric photon-ALP mixing matrix ${\mathcal{M}}_0$ entering Equation (10) has the form

$${\mathcal{M}}_0=\left(\begin{array}{ccc}{\Delta }_{xx}& {\Delta }_{xz}& {\Delta }_{a\gamma }^x\\ {\Delta }_{zx}& {\Delta }_{zz}& {\Delta }_{a\gamma }^z\\ {\Delta }_{a\gamma }^x& {\Delta }_{a\gamma }^z& {\Delta }_{aa}\end{array}\right),$$

(13)

where we have set

$${\Delta }_{a\gamma }^x\equiv \frac{g_{a\gamma \gamma }{B}_x}{2},{\Delta }_{a\gamma }^z\equiv \frac{g_{a\gamma \gamma }{B}_z}{2},{\Delta }_{aa}\equiv -\frac{m_a^2}{2E}.$$

(14)

While the terms appearing in the third row and column of ${\mathcal{M}}_0$ are dictated by ${\mathcal{L}}_{\mathrm{ALP}}$ and have an evident physical meaning, the other $\Delta$-terms require some explanation. They reflect the properties of the medium—which are not included in ${\mathcal{L}}_{\mathrm{ALP}}$—and the off-diagonal ${\Delta }_{xz}={\Delta }_{zx}$ term directly mixes the photon polarization states giving rise to Faraday rotation, while ${\Delta }_{xx}$ and ${\Delta }_{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\gg 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

$$\left(i\frac{d}{dy}+E+{\mathcal{M}}_0\right)\psi \left(y\right)=0,$$

(15)

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 polarization—is then described by the three-dimensional wave function $\psi \left(y\right)$, with the y-coordinate replacing time, which obeys the Schödinger-like equation (15) with Hamiltonian

$$H_0\equiv -\left(E+{\mathcal{M}}_0\right).$$

(16)

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

$$\psi \left(y\right)=U_0(y,y_0)\psi \left(y_0\right),$$

(17)

where $y_0$ represents the initial position. Moreover, we can set

$$U_0(E;y,y_0)=e^{iE(y-y_0)}{\mathcal{U}}_0(y,y_0),$$

(18)

where ${\mathcal{U}}_0(y,y_0)$ is the transfer matrix associated with the reduced Schödinger-like equation

$$\left(i\frac{d}{dy}+{\mathcal{M}}_0\right)\psi \left(y\right)=0.$$

(19)

3.2. A Simplified Case

Because $\mathbf{B}$ is supposed to be homogeneous, we have the freedom to choose the z-axis along $\mathbf{B}$, so that $B_x=0$. The diagonal $\Delta$-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 $\mathbf{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]

$${\omega }_{\mathrm{pl}}=3.69·{10}^{-11}{\left(\frac{n_e}{{\mathrm{cm}}^{-3}}\right)}^{1/2},$$

(20)

which entails

$${\Delta }_{\mathrm{pl}}=-\frac{{\omega }_{\mathrm{pl}}^2}{2E}.$$

(21)

Finally, the ${\Delta }_{xz}$, ${\Delta }_{zx}$ terms account for Faraday rotation, but since we are going to take E in the X-ray or in the $\gamma$-ray band Faraday rotation can be discarded. Altogether, the mixing matrix becomes

$${\mathcal{M}}_0^{\left(0\right)}=\left(\begin{array}{ccc}{\Delta }_{\mathrm{pl}}& 0& 0\\ 0& {\Delta }_{\mathrm{pl}}& {\Delta }_{a\gamma }\\ 0& {\Delta }_{a\gamma }& {\Delta }_{aa}\end{array}\right),$$

(22)

with the superscript $\left(0\right)$ recalling the present choice of the coordinate system and

$${\Delta }_{a\gamma }\equiv \frac{g_{a\gamma \gamma }B}{2}.$$

(23)

We see that $A_x$ decouples away while only $A_z$ mixes with a.

Application of the discussion reported in Appendix A with $\mathcal{M}\to {\mathcal{M}}_0^{\left(0\right)}$ yields for the corresponding eigenvalues

$${\lambda }_{0,1}={\Delta }_{\mathrm{pl}},{\lambda }_{0,2}=\frac{1}{2}\left({\Delta }_{\mathrm{pl}}+{\Delta }_{aa}-{\Delta }_{\mathrm{osc}}\right),\frac{1}{2}\left({\Delta }_{\mathrm{pl}}+{\Delta }_{aa}+{\Delta }_{\mathrm{osc}}\right),$$

(24)

where we have set

$${\Delta }_{\mathrm{osc}}\equiv {\left[{\left({\Delta }_{\mathrm{pl}}-{\Delta }_{aa}\right)}^2+4{\left({\Delta }_{a\gamma }\right)}^2\right]}^{1/2}={\left[{\left(\frac{m_a^2-{\omega }_{\mathrm{pl}}^2}{2E}\right)}^2+{\left(g_{a\gamma \gamma }B\right)}^2\right]}^{1/2}.$$

(25)

As a consequence, the transfer matrix associated with Equation (19) with mixing matrix ${\mathcal{M}}_0^{\left(0\right)}$ can be written with the help of Equation (A16) as

$${\mathcal{U}}_0(y,y_0;0)=e^{i{\lambda }_1(y-y_0)}{T}_{0,1}\left(0\right)+e^{i{\lambda }_2(y-y_0)}{T}_{0,2}\left(0\right)+e^{i{\lambda }_3(y-y_0)}{T}_{0,3}\left(0\right),$$

(26)

where the matrices $T_{0,1}\left(0\right)$, $T_{0,2}\left(0\right)$ and $T_{0,3}\left(0\right)$ 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

$$\alpha =\frac{1}{2}\mathrm{arctg}\left(\frac{2{\Delta }_{a\gamma }}{{\Delta }_{\mathrm{pl}}-{\Delta }_{aa}}\right)=\frac{1}{2}\mathrm{arctg}\left[\left(g_{a\gamma \gamma }B\right)\left(\frac{2E}{m_a^2-{\omega }_{\mathrm{pl}}^2}\right)\right],$$

(27)

since then simple trigonometric manipulations allow us to express the above matrices in the simpler form

$$T_{0,1}\left(0\right)\equiv \left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),$$

(28)

$$T_{0,2}\left(0\right)\equiv \left(\begin{array}{ccc}0& 0& 0\\ 0& {sin}^2\alpha & -sin\alpha cos\alpha \\ 0& -sin\alpha cos\alpha & {cos}^2\alpha \end{array}\right),$$

(29)

$$T_{0,3}\left(0\right)\equiv \left(\begin{array}{ccc}0& 0& 0\\ 0& {cos}^2\alpha & sin\alpha cos\alpha \\ 0& sin\alpha cos\alpha & {sin}^2\alpha \end{array}\right).$$

(30)

Now, the probability that a photon polarized along the z-axis oscillates into an ALP after a distance y is evidently

$$P_{0,{\gamma }_z\to a}^{\left(0\right)}\left(y\right)={\left|\langle a\left|{\mathcal{U}}_0(y,0;0)\right|{\gamma }_z\rangle \right|}^2$$

(31)

and in complete analogy with the case of neutrino oscillations [33] it reads

$$P_{0,{\gamma }_z\to a}^{\left(0\right)}\left(y\right)={\sin }^2\left(2\alpha \right){\sin }^2\left(\frac{{\Delta }_{\mathrm{osc}}y}{2}\right),$$

(32)

which shows that ${\Delta }_{\mathrm{osc}}$ plays the role of oscillation wave number, thereby implying that the oscillation length is $L_{\mathrm{osc}}=2\pi /{\Delta }_{\mathrm{osc}}$. Owing to Equations (27) and (32) can be rewritten as

$$P_{0,{\gamma }_z\to a}^{\left(0\right)}\left(y\right)={\left(\frac{g_{a\gamma \gamma }B}{{\Delta }_{\mathrm{osc}}}\right)}^2{\sin }^2\left(\frac{{\Delta }_{\mathrm{osc}}y}{2}\right),$$

(33)

which shows that the photon-ALP oscillation probability becomes both maximal and independent of E and $m_a$ for

$${\Delta }_{\mathrm{osc}}\simeq g_{a\gamma \gamma }B,$$

(34)

and explicitly reads

$$P_{0,{\gamma }_z\to a}^{\left(0\right)}\left(y\right)\simeq {\sin }^2\left(\frac{g_{a\gamma \gamma }By}{2}\right).$$

(35)

This is the strong-mixing regime, which—from the comparison of Equations (25) and (34)—turns out to be characterized by the condition

$$\frac{|m_a^2-{\omega }_{\mathrm{pl}}^2|}{2E}\ll g_{a\gamma \gamma }B,$$

(36)

and so it sets in sufficiently above the energy threshold

$$E_L\equiv \frac{|m_a^2-{\omega }_{\mathrm{pl}}^2|}{2g_{a\gamma \gamma }B}.$$

(37)

Observe that in the strong-mixing regime the mass term ${\Delta }_{aa}=-m_a^2/\left(2E\right)$ and the plasma term ${\Delta }_{\mathrm{pl}}=-{\omega }_{\mathrm{pl}}^2/\left(2E\right)$ 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\gamma \gamma }$.

3.3. A More General Case

In view of our subsequent discussion it proves essential to deal with the general case in which $\mathbf{B}$ is not aligned with the z-axis but forms a nonvanishing angle $\psi$ with it. Correspondingly, the mixing matrix ${\mathcal{M}}_0$ presently arises from ${\mathcal{M}}_0^{\left(0\right)}$ through the similarity transformation

$${\mathcal{M}}_0=V^{†}\left(\psi \right){\mathcal{M}}_0^{\left(0\right)}V\left(\psi \right)$$

(38)

operated by the rotation matrix in the x–z plane, namely

$$V\left(\psi \right)=\left(\begin{array}{ccc}cos\psi & -sin\psi & 0\\ sin\psi & cos\psi & 0\\ 0& 0& 1\end{array}\right).$$

(39)

This leads to

$${\mathcal{M}}_0=\left(\begin{array}{ccc}{\Delta }_{\mathrm{pl}}& 0& {\Delta }_{a\gamma }sin\psi \\ 0& {\Delta }_{\mathrm{pl}}& {\Delta }_{a\gamma }cos\psi \\ {\Delta }_{a\gamma }sin\psi & {\Delta }_{a\gamma }cos\psi & {\Delta }_{aa}\end{array}\right),$$

(40)

indeed in agreement with Equation (13) within the considered approximation. Therefore the transfer matrix reads

$${\mathcal{U}}_0(y,y_0;\psi )=V^{†}\left(\psi \right){\mathcal{U}}_0(y,y_0;0)V\left(\psi \right)$$

(41)

and its explicit representation turns out to be

$${\mathcal{U}}_0(y,y_0;\psi )=e^{i{\lambda }_1(y-y_0)}{T}_{0,1}\left(\psi \right)+e^{i{\lambda }_2(y-y_0)}{T}_{0,2}\left(\psi \right)+e^{i{\lambda }_3(y-y_0)}{T}_{0,3}\left(\psi \right),$$

(42)

with

$$T_{0,1}\left(\psi \right)\equiv \left(\begin{array}{ccc}{cos}^2\psi & -sin\psi cos\psi & 0\\ -sin\psi cos\psi & {sin}^2\psi & 0\\ 0& 0& 0\end{array}\right),$$

(43)

$$T_{0,2}\left(\psi \right)\equiv \left(\begin{array}{ccc}{sin}^2\theta {sin}^2\psi & {sin}^2\alpha sin\psi cos\psi & -sin\alpha cos\alpha sin\psi \\ {sin}^2\alpha sin\psi cos\psi & {sin}^2\alpha {cos}^2\psi & -sin\alpha cos\alpha cos\psi \\ -sin\alpha cos\alpha sin\psi & -sin\alpha cos\alpha cos\psi & {cos}^2\alpha \end{array}\right),$$

(44)

$$T_{0,3}\left(\psi \right)\equiv \left(\begin{array}{ccc}{sin}^2\psi {cos}^2\alpha & sin\psi cos\psi {cos}^2\alpha & sin\alpha cos\alpha sin\psi \\ sin\psi cos\psi {cos}^2\alpha & {cos}^2\psi {cos}^2\alpha & sin\alpha cos\alpha cos\psi \\ sin\psi cos\alpha sin\alpha & cos\psi sin\alpha cos\alpha & {sin}^2\alpha \end{array}\right).$$

(45)

As a result, a beam initially containing only photons propagating in an external magnetic field, after a distance of many oscillation lengths $L_{\mathrm{osc}}$ will contain ALPs. Moreover, the three states $|{\gamma }_x\rangle$, $|{\gamma }_z\rangle$, $|a\rangle$ will equilibrate, namely the beam will be composed by $2/3$ of photons and of $1/3$ of ALPs.

3.4. 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 $\mathcal{M}$—will not be self-adjoint, and correspondingly the transfer matrix will be denoted by $\mathcal{U}(y,y_0)$ and will be not unitary anymore. Thus, the beam looks formally like a three-state non-relativistic unstable quantum system.

3.5. 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

$$\rho \left(y\right)=\left(\begin{array}{c}{A}_x\left(y\right)\\ A_z\left(y\right)\\ a\left(y\right)\end{array}\right)\otimes {\left(\begin{array}{c}{A}_x\left(y\right)A_z\left(y\right)a\left(y\right)\end{array}\right)}^{*}$$

(46)

rather than by a wave function $\psi \left(y\right)$. Remarkably, the analogy with non-relativistic quantum mechanics entails that $\rho \left(y\right)$ obeys the Von Neumann-like equation

$$i\frac{d\rho }{dy}=\rho {\mathcal{M}}^{†}-\mathcal{M}\rho$$

(47)

associated with Equation (19). Thus, the propagation of a generic $\rho \left(y\right)$ is still given by

$$\rho \left(y\right)=\mathcal{U}(y,y_0)\rho \left(y_0\right){\mathcal{U}}^{†}(y,y_0)$$

(48)

and the probability that a photon/ALP beam initially in the state ${\rho }_1$ at position $y_0$ will be found in the state ${\rho }_2$ at position y is

$$P_{{\rho }_1\to {\rho }_2}\left(y\right)=\mathrm{Tr}\left({\rho }_2\mathcal{U}(y,y_0){\rho }_1{\mathcal{U}}^{†}(y,y_0)\right),$$

(49)

provided that ${\rho }_2$ is measured, since we are assuming as usual that $\mathrm{Tr}{\rho }_1=\mathrm{Tr}{\rho }_2=1$.

3.6. Parameter Bounds

Before proceeding further, it seems important to report the observational bounds on the parameters $g_{a\gamma \gamma }$ 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 $\gamma +\mathrm{ion}\to a+\mathrm{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\gamma \gamma }<0.66·{10}^{-10}{\mathrm{GeV}}^{-1}$ for $m<0.02\mathrm{eV}$ at the $2\sigma$ 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\gamma \gamma }$ 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.

• $g_{a\gamma \gamma }<2.1·{10}^{-11}{\mathrm{GeV}}^{-1}$ for $15·{10}^{-9}\mathrm{eV}<m_a<60·{10}^{-9}\mathrm{eV}$ from PKS 2155-304 [83].
• $g_{a\gamma \gamma }<5·{10}^{-12}{\mathrm{GeV}}^{-1}$ for $5·{10}^{-10}\mathrm{eV}<m_a<5·{10}^{-9}\mathrm{eV}$ from NGC 1275 [97].
• $g_{a\gamma \gamma }<2.6·{10}^{-12}{\mathrm{GeV}}^{-1}$ for $m_a<{10}^{-13}\mathrm{eV}$ from M87 [103].
• $g_{a\gamma \gamma }<{10}^{-11}{\mathrm{GeV}}^{-1}$ for $0.6·{10}^{-9}\mathrm{eV}<m_a<4·{10}^{-9}\mathrm{eV}$ from PKS 2155-304 [111].
• $g_{a\gamma \gamma }<(2·{10}^{-11}$–$6·{10}^{-11}){\mathrm{GeV}}^{-1}$ for $5·{10}^{-10}\mathrm{eV}<m_a<5·{10}^{-7}\mathrm{eV}$ from Mrk 421 [120].
• $g_{a\gamma \gamma }<3·{10}^{-11}{\mathrm{GeV}}^{-1}$ for $1·{10}^{-8}\mathrm{eV}<m_a<2·{10}^{-7}\mathrm{eV}$ from Mrk 421 [121].
• $g_{a\gamma \gamma }<(6–8)·{10}^{-13}{\mathrm{GeV}}^{-1}$ for $m_a<1·{10}^{-12}\mathrm{eV}$ from NGC 1275 (center of Perseus cluster) [123].

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 $\gamma$-rays and been detected by the Solar Maximum Mission satellite. From the failure of detection, upper bounds on $m_a$ and $g_{a\gamma \gamma }$ can be established [42,43]. This issue was reconsidered in more detail in 2015 by Payez et al., getting the improved bound $m_a\lesssim 4.4·{10}^{-10}\mathrm{eV}$ and $g_{a\gamma \gamma }\lesssim 5.3·{10}^{-12}{\mathrm{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 strong-mixing regime sets in—we will consider throughout this Review $m_a=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$ and $g_{a\gamma \gamma }=\mathcal{O}\left({10}^{-11}\right){\mathrm{GeV}}^{-1}$. Hence, we can easily take $m_a>4.4·{10}^{-10}\mathrm{eV}$ no contradiction exists even apart from the previous remark.

4. 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.

4.1. 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–1$\mathrm{Mpc}$. Ultra-relativistic 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) $\nu F_{\nu }\left(\nu \right)$—$F_{\nu }\left(\nu \right)$ 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 $\gamma$-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 ${\pi }^0$ immediately decays as ${\pi }^0\to \gamma +\gamma$, while the ${\pi }^{\pm }$ 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) ${\Pi }_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\simeq 300\mathrm{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)\mathrm{pc}$—and infrared photons emitted by the torus, the VHE photons produced at the jet base undergo the process $\gamma +\gamma \to e^{+}+e^{-}$. As a result, the FSRQs should be invisible above $(20–30)\mathrm{GeV}$ [187,188,189,190]. However, observations with IACTs have shown that such an expectation is blatantly wrong, since they have been detected up to 400 GeV. This fact poses an open problem.

So far, about 85 VHE blazars have been detected by the Imaging Atmospheric Cherenkov Telescopes (IACTs) H.E.S.S. (High Energy Stereoscopic System) [191], MAGIC (Major Atmospheric Gamma Imaging Cherenkov) [192] and VERITAS (Very Energetic Radiation Imaging Telescope Array System) [193] with redshift up to $z\simeq 1$ [194].

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 $\Lambda$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) $\Phi \equiv dN/\left(dtdAdE\right)$. They are related by

$${\Phi }_{\mathrm{obs}}(E_0,z)=P_{\gamma \to \gamma }(E_0,z){\Phi }_{\mathrm{em}}\left(E_0(1+z)\right),$$

(50)

where $P_{\gamma \to \gamma }(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

$$\nu F_{\nu }(\nu ,E_0,z)=E_0^2{\Phi }_{\mathrm{obs}}(E_0,z),$$

(51)

where $F_{\nu }(\nu ,E_0,z)$ is the specific apparent luminosity. We suppose hereafter that $E_0$ lies in the VHE $\gamma$-ray band.

4.2. Conventional Photon Propagation

Within conventional physics the photon survival probability $P_{\gamma \to \gamma }^{\mathrm{CP}}(E_0,z)$ is usually parametrized as

$$P_{\gamma \to \gamma }^{\mathrm{CP}}(E_0,z)=e^{-{\tau }_{\gamma }(E_0,z)},$$

(52)

where ${\tau }_{\gamma }(E_0,z)$ is the optical depth, which quantifies the dimming of the source. Note that in general ${\tau }_{\gamma }(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 ${\tau }_{\gamma }(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 ${\tau }_{\gamma }(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 ${\tau }_{\gamma }(E_0,z_h)=1$ defines the $\gamma$-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

$${\Phi }_{\mathrm{obs}}(E_0,z)=e^{-{\tau }_{\gamma }(E_0,z)}{\Phi }_{\mathrm{em}}\left(E_0(1+z)\right).$$

(53)

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 $\gamma \gamma \to 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

$$ϵ>{ϵ}_{\mathrm{thr}}(E,\phi )\equiv \frac{2m_e^2{c}^4}{E\left(1-cos\phi \right)},$$

(54)

where $\phi$ 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]

$${\sigma }_{\gamma \gamma }(E,ϵ,\phi )\simeq 1.25·{10}^{-25}\left(1-{\beta }^2\right)\left[2\beta \left({\beta }^2-2\right)+\left(3-{\beta }^4\right)\ln \left(\frac{1+\beta }{1-\beta }\right)\right]{\mathrm{cm}}^2,$$

(55)

which depends on E, $ϵ$ and $\phi$ only through the dimensionless parameter

$$\beta (E,ϵ,\phi )\equiv {\left[1-\frac{2m_e^2{c}^4}{Eϵ\left(1-cos\phi \right)}\right]}^{1/2},$$

(56)

and the process is kinematically allowed for ${\beta }^2>0$. The cross-section ${\sigma }_{\gamma \gamma }(E,ϵ,\phi )$ reaches its maximum ${\sigma }_{\gamma \gamma }^{\max }\simeq 1.70·{10}^{-25}{\mathrm{cm}}^2$ for $\beta \simeq 0.70$. Assuming head-on collisions for definiteness ($\phi =\pi$), it follows that ${\sigma }_{\gamma \gamma }(E,ϵ,\pi )$ gets maximized for the background photon energy

$$ϵ\left(E\right)\simeq \left(\frac{900\mathrm{GeV}}{E}\right)\mathrm{eV},$$

(57)

where E and $ϵ$ correspond to the same redshift.

Within the standard $\Lambda$CDM cosmological model ${\tau }_{\gamma }(E_0,z)$ arises by first convolving the spectral number density $n_{\gamma }(ϵ\left(z\right),z)$ of background photons at a generic redshift with ${\sigma }_{\gamma \gamma }(E\left(z\right),ϵ\left(z\right),\phi )$ along the line of sight for fixed values of z, $\phi$ and $ϵ\left(z\right)$, and next integrating over all these variables [196,197,198]. Hence, we have

$$\begin{array}{c}\hfill {\tau }_{\gamma }(E_0,z)={\int }_0^z\mathrm{d}z\frac{\mathrm{d}l\left(z\right)}{\mathrm{d}z}{\int }_{-1}^1\mathrm{d}\left(cos\phi \right)\frac{1-cos\phi }{2}\times \\ \hfill \times {\int }_{{ϵ}_{\mathrm{thr}}(E\left(z\right),\phi )}^{\infty }\mathrm{d}ϵ\left(z\right)n_{\gamma }(ϵ\left(z\right),z){\sigma }_{\gamma \gamma }\left(E\left(z\right),ϵ\left(z\right),\phi \right),\end{array}$$

(58)

where the distance travelled by a photon per unit redshift at redshift z is given by

$$\frac{dl\left(z\right)}{dz}=\frac{c}{H_0}\frac{1}{\left(1+z\right){\left[{\Omega }_{\Lambda }+{\Omega }_M{\left(1+z\right)}^3\right]}^{1/2}},$$

(59)

with Hubble constant $H_0\simeq 70\mathrm{Km}{\mathrm{s}}^{-1}{\mathrm{Mpc}}^{-1}$, while ${\Omega }_{\Lambda }\simeq 0.7$ and ${\Omega }_M\simeq 0.3$ represent the average cosmic density of matter and dark energy, respectively, in units of the critical density ${\rho }_{\mathrm{cr}}\simeq 0.97·{10}^{-29}\mathrm{g}{\mathrm{cm}}^{-3}$.

Once $n_{\gamma }(ϵ\left(z\right),z)$ is known, ${\tau }_{\gamma }(E_0,z)$ can be computed exactly, even though in general the integration over $ϵ\left(z\right)$ 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

$${\tau }_{\gamma }(E,D)=\frac{D}{{\lambda }_{\gamma }\left(E\right)},$$

(60)

where ${\lambda }_{\gamma }\left(E\right)$ is the photon mean free path for $\gamma \gamma \to e^{+}{e}^{-}$ referring to the present cosmic epoch. As a consequence, Equation (52) becomes

$$P_{\gamma \to \gamma }^{\mathrm{CP}}(E,D)=e^{-D/{\lambda }_{\gamma }\left(E\right)},$$

(61)

and so Equation (53) reduces to

$${\Phi }_{\mathrm{obs}}(E,D)=e^{-D/{\lambda }_{\gamma }\left(E\right)}{\Phi }_{\mathrm{em}}\left(E\right).$$

(62)

Note that we have dropped the subscript 0 for simplicity.

4.3. Extragalactic Background Light (EBL)

Blazars detected so far with IACTs—or detectable in the near future—lie in the VHE range $100\mathrm{GeV}<E_0<100\mathrm{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\mathrm{eV}<{ϵ}_0<9\mathrm{eV}$ (corresponding to the frequency range $2.17·{10}^3\mathrm{GHz}<{\nu }_0<2.17·{10}^6\mathrm{GHz}$ and to the wavelength range $0.14\mathsf{\mu }\mathrm{m}<{\lambda }_0<1.38·{10}^2\mathsf{\mu }\mathrm{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 ${\tau }_{\gamma }(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 ${\tau }_{\gamma }(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.

4.4. Extragalactic Magnetic Field

Unfortunately, the morphology and strength of the extragalactic magnetic field ${\mathbf{B}}_{\mathrm{ext}}$ are totally unknown, and it is not surprising that various very different configurations of ${\mathbf{B}}_{\mathrm{ext}}$ have been proposed [204,205,206,207]. Yet, the strength of ${\mathbf{B}}_{\mathrm{ext}}$ is constrained to lies in the range ${10}^{-7}\mathrm{nG}\lesssim {\mathrm{B}}_{\mathrm{ext}}\lesssim 1.7\mathrm{nG}$ on the scale $\mathcal{O}\left(1\right)\mathrm{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 ${\mathbf{B}}_{\mathrm{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_{\mathrm{ext}}=\mathcal{O}\left(1\right)\mathrm{nG}$ on the scale $\mathcal{O}\left(1\right)\mathrm{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 Lyman-alpha forest clouds [214]. A similar situation was further investigated in 2001 by Furlanetto and Loeb [215] in connection with quasars outflows, which still predicts ${\mathrm{B}}_{\mathrm{ext}}=\mathcal{O}\left(1\right)\mathrm{nG}$ on the scale $\mathcal{O}\left(1\right)\mathrm{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 domain-like network, made of identically domains of size $L_{\mathrm{dom}}$ equal to the ${\mathbf{B}}_{\mathrm{ext}}$ coherence length, all having roughly the same strength of ${\mathbf{B}}_{\mathrm{ext}}$ and with the direction of ${\mathbf{B}}_{\mathrm{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 ${\mathbf{B}}_{\mathrm{ext}}$ are galaxies indeed explains the three main features of the considered model for ${\mathbf{B}}_{\mathrm{ext}}$. (1) Its cell-like morphology arises from its galactic origin. (2) It is quite plausible that ${\mathbf{B}}_{\mathrm{ext}}$ has nearly the same strength around each galaxy, and so in all domains. (3) Because galaxies are uncorrelated, it looks natural that the ${\mathbf{B}}_{\mathrm{ext}}$ direction changes randomly from the neighborhood of a galaxy to that of another, thereby explaining the random jump in direction of ${\mathbf{B}}_{\mathrm{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 ${\mathbf{B}}_{\mathrm{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_{\mathrm{dom}}$ was invariably much larger than the oscillation length $L_{\mathrm{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_{\mathrm{osc}}\lesssim L_{\mathrm{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 ${\mathbf{B}}_{\mathrm{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 ${\mathbf{B}}_{\mathrm{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 ${\mathbf{B}}_{\mathrm{ext}}$. As a by-product, a prediction of the magnetic field ${\mathbf{B}}_{\mathrm{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.

4.5. Galaxy Clusters

According to Abell, Galaxy clusters contain from thirty up to more than thousands galaxies with a total mass in the range $\left({10}^{14}–{10}^{15}\right)M_{\odot }$ 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 ${\mathbf{B}}^{\mathrm{clu}}$ and electron number density $n_e^{\mathrm{clu}}$.

Faraday rotation measurements and synchrotron radio emissions tell us that $B^{\mathrm{clu}}=\mathcal{O}(1–10)\mathsf{\mu }\mathrm{G}$ [222,223]. The structure of ${\mathbf{B}}^{\mathrm{clu}}$ is believed to have a turbulent nature with a Kolmogorov-type turbulence spectrum $M\left(k\right)\propto 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 ${\mathbf{B}}^{\mathrm{clu}}$ is modeled by [224,225]

$$B^{\mathrm{clu}}\left(y\right)=\mathcal{B}\left(B_0^{\mathrm{clu}},k,q,y\right){\left(\frac{n_e^{\mathrm{clu}}\left(y\right)}{n_{e,0}^{\mathrm{clu}}}\right)}^{{\eta }_{\mathrm{clu}}},$$

(63)

where $\mathcal{B}$ is the spectral function describing the Kolmogorov-type turbulence of the cluster magnetic field (for more details see e.g., [86]), $B_0^{\mathrm{clu}}$ and $n_{e,0}^{\mathrm{clu}}$ are the central cluster magnetic field strength and the central electron number density, respectively, while ${\eta }_{\mathrm{clu}}$ is a cluster parameter. The behavior of $n_e^{\mathrm{clu}}$ is modeled by the model

$$n_e^{\mathrm{clu}}\left(y\right)=n_{e,0}^{\mathrm{clu}}{\left(1+\frac{y^2}{r_{\mathrm{core}}^2}\right)}^{-\frac{3}{2}{\beta }_{\mathrm{clu}}},$$

(64)

where ${\beta }_{\mathrm{clu}}$ is a cluster parameter and $r_{\mathrm{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_{e,0}^{\mathrm{clu}}$ represents the real crucial quantity. We take the following average values for the two classes: $n_{e,0}^{\mathrm{clu}}=5·{10}^{-2}{\mathrm{cm}}^{-3}$ for CC clusters and $n_{e,0}^{\mathrm{clu}}=0.5·{10}^{-2}{\mathrm{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 $\gamma \leftrightarrow 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.

5. 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\mathrm{GeV}$ up to a few TeV. As a consequence, we have $E_0\gg 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_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z\right)$ 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—$\gamma \leftrightarrow 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 ${\tau }_{\gamma }^{\mathrm{ALP}}$ is now smaller than in the conventional case. But since the corresponding photon survival probability is

$$P_{\gamma \to \gamma }^{\mathrm{ALP}}(E_0,z)=e^{-{\tau }_{\gamma }^{\mathrm{ALP}}(E_0,z)},$$

(65)

recalling (52) we conclude that $P_{\gamma \to \gamma }^{\mathrm{ALP}}(E_0,z)$ is much larger than $P_{\gamma \to \gamma }(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 section—we 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_{\mathrm{osc}}\gg L_{\mathrm{dom}}$.

5.1. Strategy

Thanks to the fact that $\mathbf{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 $\Lambda$CDM cosmological model with ${\Omega }_M=0.3$ and ${\Omega }_{\Lambda }=0.7$, and so the redshift is the natural parameter to express distances. In particular, the proper length $L_{\mathrm{dom}}(z_a,z_b)$ extending over the redshift interval $[z_a,z_b]$ is

$$\begin{array}{cc}& {\displaystyle L(z_a,z_b)\simeq 4.29·{10}^3{\int }_{z_a}^{z_b}\frac{dz}{(1+z){[0.7+0.3{(1+z)}^3]}^{1/2}}\mathrm{Mpc}\simeq }\\ & {\displaystyle \simeq 2.96·{10}^3\ln \left(\frac{1+1.45z_b}{1+1.45z_a}\right)\mathrm{Mpc}.}\end{array}$$

(66)

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 $\Delta z$, which is therefore the same for all domains. This mesh can be constructed as follows. We denote by $L_{\mathrm{dom}}^{\left(n\right)}=L\left((n-1)\Delta z,n\Delta z\right)$ the proper length along the y-direction of the generic n-th domain, with $1\le n\le N_d$. Note that $N_d$ is the maximal integer contained in the number $z/\Delta z$, hence $N_d\simeq z/\Delta z$. In order to fix $\Delta z$ we consider the domain closest to us, labelled by 1 and—with the help of Equation (66)—we write its proper length as $\left(L_{\mathrm{dom}}^{\left(1\right)}/5\mathrm{Mpc}\right)5\mathrm{Mpc}=L(0,\Delta z)=2.96·{10}^3\ln (1+1.45\Delta z)\mathrm{Mpc}$, from which we get $\Delta z\simeq 1.17·{10}^{-3}\left(L_{\mathrm{dom}}^{\left(1\right)}/5\mathrm{Mpc}\right)$. So, once $L_{\mathrm{dom}}^{\left(1\right)}$ is chosen in agreement with such a prescription, the size of all magnetic domains in redshift space is fixed. At this point, two further quantities can be determined. First, $N_d\simeq z/\Delta z\simeq 0.85·{10}^3\left(5\mathrm{Mpc}/L_{\mathrm{dom}}^{\left(1\right)}\right)z$. Second, the proper length of the n-th domain along the y-direction follows from Equation (66) with $z_a\to (n-1)\Delta z,z_b\to n\Delta z$. Whence

$$L_{\mathrm{dom}}^{\left(n\right)}\simeq 2.96·{10}^3\ln \left(1+\frac{1.45\Delta z}{1+1.45(n-1)\Delta z}\right)\mathrm{Mpc}.$$

(67)

* * *

Manifestly, in order to maximize $P_{\gamma \to \gamma }^{\mathrm{ALP}}(E_0,z)$ we choose $m_a$ in order to be in the strong-mixing regime for $E_0\gtrsim 100\mathrm{GeV}$. Incidentally, when the external magnetic field is homogeneous—as it is in fact in each single domain—a look at Lagrangian (8) shows that all results depend on the combination $g_{a\gamma \gamma }B$ and not on $g_{a\gamma \gamma }$ and B separately. It is therefore quite convenient to employ the parameter

$$\xi \equiv \left(\frac{B}{\mathrm{nG}}\right)\left(g_{a\gamma \gamma }{10}^{11}\mathrm{GeV}\right),$$

(68)

in terms of which Equation (37) can be rewritten as

$$E_L=\frac{25.64}{\xi }\left|{\left(\frac{m_a}{{10}^{-10}\mathrm{eV}}\right)}^2-{\left(\frac{{\omega }_{\mathrm{pl}}}{{10}^{-10}\mathrm{eV}}\right)}^2\right|\mathrm{GeV}.$$

(69)

Because we would like to be in the strong-mixing regime almost everywhere within the VHE band, we take $E_L=\mathcal{O}\left(100\right)\mathrm{GeV}$. What about $m_a$? We should keep in mind that ${\omega }_{\mathrm{pl}}$ is unknown, but the upper bound on the mean diffuse extragalactic electron density $n_e<2.7·{10}^{-7}{\mathrm{cm}}^{-3}$ is provided by the WMAP measurement of the baryon density [232], which—thanks to Equation (20)—translates into the upper bound ${\omega }_{\mathrm{pl}}<1.92·{10}^{-14}\mathrm{eV}$. Moreover, in order to fix $\xi$ we use the fact that the result to be derived in the next section requires $\xi =0.5$. As a consequence, we get $m_a=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$.

5.2. Propagation over a Single Domain

We have to determine ${\lambda }_{\gamma }^{\left(n\right)}$ and the magnetic field strength $B^{\left(n\right)}$ 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 ${\lambda }_{\gamma }$ for the reaction $\gamma \gamma \to e^{+}{e}^{-}$, and the term $i/2{\lambda }_{\gamma }$ should be inserted into the 11 and 22 entries of the $\mathcal{M}$ matrix. In order to evaluate ${\lambda }_{\gamma }$, we imagine that two hypothetical sources located at both edges of the n-th domain are observed. Therefore, we apply Equation (53) to both sources. With the notational simplifications ${\Phi }_{\mathrm{obs}}(E_0,z)\to \Phi \left(E_0\right)$ and ${\Phi }_{\mathrm{em}}\left(E_0(1+z)\right)\to \Phi \left(E_0(1+z)\right)$, we have

$$\Phi \left(E_0\right)=e^{-{\tau }_{\gamma }\left(E_0,(n-1)\Delta z\right)}\Phi \left(E_0[1+(n-1)\Delta z]\right),$$

(70)

$$\Phi \left(E_0\right)=e^{-{\tau }_{\gamma }(E_0,n\Delta z)}\Phi \left(E_0(1+n\Delta z)\right),$$

(71)

which upon combination imply that the flux change across the domain in question is

$$\Phi \left(E_0[1+(n-1)\Delta z]\right)=e^{-\left[{\tau }_{\gamma }(E_0,n\Delta z)-{\tau }_{\gamma }\left(E_0,(n-1)\Delta z\right)\right]}\Phi \left(E_0(1+n\Delta z)\right).$$

(72)

But owing to Equation (62) mutatis mutandis implies that Equation (72) should have the form

$$\Phi \left(E_0[1+(n-1)\Delta z]\right)=e^{-L_{\mathrm{dom}}^{\left(n\right)}/{\lambda }_{\gamma }^{\left(n\right)}\left(E_0\right)}\Phi \left(E_0,(1+n\Delta z)\right),$$

(73)

and the comparison with Equation (72) ultimately yields

$${\lambda }_{\gamma }^{\left(n\right)}\left(E_0\right)=\frac{L_{\mathrm{dom}}^{\left(n\right)}}{{\tau }_{\gamma }(E_0,n\Delta z)-{\tau }_{\gamma }\left(E_0,(n-1)\Delta z\right)},$$

(74)

where the optical depth is evaluated by means of Equation (58) or more simply taken from [200].

As for the determination of $B^{\left(n\right)}$, 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\left(z\right)=B(z=0){(1+z)}^2$ [204,205]. Hence in the n-th magnetic domain we have $B^{\left(n\right)}=B^{\left(1\right)}{\left(1+(n-1)\Delta z\right)}^2$.

Thus, at this stage the mixing matrix $\mathcal{M}$ as explicitly written in the n-th domain reads

$${\mathcal{M}}^{\left(n\right)}=\left(\begin{array}{ccc}i/2{\lambda }_{\gamma }^{\left(n\right)}& 0& B^{\left(n\right)}\sin {\psi }_n{g}_{a\gamma \gamma }/2\\ 0& i/2{\lambda }_{\gamma }^{\left(n\right)}& B^{\left(n\right)}\cos {\psi }_n{g}_{a\gamma \gamma }/2\\ B^{\left(n\right)}\sin {\psi }_n{g}_{a\gamma \gamma }/2& B^{\left(n\right)}\cos {\psi }_n{g}_{a\gamma \gamma }/2& 0\end{array}\right),$$

(75)

where ${\psi }_n$ is the random angle between ${\mathbf{B}}^{\left(n\right)}$ and the fixed fiducial direction along the z-axis (note that indeed ${\mathcal{M}}^{†}\ne \mathcal{M}$). Observe that since we are in the strong-mixing regime, ${\omega }_{\mathrm{pl}}$ and $m_a$ can be neglected with respect to the other terms in the mixing matrix. So—apart from ${\psi }_n$—all other matrix elements entering ${\mathcal{M}}^{\left(n\right)}$ are known. Finding the transfer matrix corresponding to ${\mathcal{M}}^{\left(n\right)}$ is straightforward even if tedious by using the results reported in Appendix A. The result is

$$\begin{array}{cc}& {\displaystyle {\mathcal{U}}_n(E_n,{\psi }_n)=}\\ & {\displaystyle =e^{i{E}_n{L}_{\mathrm{dom}}^{\left(n\right)}}\left[e^{i\left({\lambda }_1^{\left(n\right)}{L}_{\mathrm{dom}}^{\left(n\right)}\right)}{T}_1\left({\psi }_n\right)+e^{i\left({\lambda }_2^{\left(n\right)}{L}_{\mathrm{dom}}^{\left(n\right)}\right)}{T}_2\left({\psi }_n\right)+e^{i\left({\lambda }_3^{\left(n\right)}{L}_{\mathrm{dom}}^{\left(n\right)}\right)}{T}_3\left({\psi }_n\right)\right]}\end{array}$$

(76)

with

$$\begin{array}{cc}& {\displaystyle T_1\left({\psi }_n\right)\equiv }\\ & {\displaystyle \equiv \left(\begin{array}{ccc}{cos}^2{\psi }_n& -sin{\psi }_{n}cos{\psi }_n& 0\\ -sin{\psi }_{n}cos{\psi }_n& {sin}^2{\psi }_n& 0\\ 0& 0& 0\end{array}\right),}\end{array}$$

(77)

$$\begin{array}{cc}& {\displaystyle T_2\left({\psi }_n\right)\equiv }\\ & {\displaystyle \equiv \left(\begin{array}{ccc}\frac{-1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}{sin}^2{\psi }_n& \frac{-1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}sin{\psi }_{n}cos{\psi }_n& \frac{i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}sin{\psi }_n\\ \frac{-1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}sin{\psi }_{n}cos{\psi }_n& \frac{-1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}{cos}^2{\psi }_n& \frac{i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}cos{\psi }_n\\ \frac{i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}sin{\psi }_n& \frac{i{\delta }_n}{\sqrt{1-4{{\delta }_n}^2}}cos{\psi }_n& \frac{1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}\end{array}\right),}\end{array}$$

(78)

$$\begin{array}{cc}& {\displaystyle T_3\left({\psi }_n\right)\equiv }\\ & {\displaystyle \equiv \left(\begin{array}{ccc}\frac{1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}{sin}^2{\psi }_n& \frac{1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}sin{\psi }_{n}cos{\psi }_n& \frac{-i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}sin{\psi }_n\\ \frac{1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}sin{\psi }_{n}cos{\psi }_n& \frac{1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}{cos}^2{\psi }_n& \frac{-i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}cos{\psi }_n\\ \frac{-i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}sin{\psi }_n& \frac{-i{\delta }_n}{\sqrt{1-4{\delta }_n^2}}cos{\psi }_n& \frac{-1+\sqrt{1-4{\delta }_n^2}}{2\sqrt{1-4{\delta }_n^2}}\end{array}\right),}\end{array}$$

(79)

where we have set

$${\lambda }_1^{\left(n\right)}\equiv \frac{i}{2{\lambda }_{\gamma }^{\left(n\right)}\left(E_0\right)},{\lambda }_2^{\left(n\right)}\equiv \frac{i}{4{\lambda }_{\gamma }^{\left(n\right)}}\left(1-\sqrt{1-4{\delta }_n^2}\right),$$

(80)

$${\lambda }_3^{\left(n\right)}\equiv \frac{i}{4{\lambda }_{\gamma }^{\left(n\right)}}\left(1+\sqrt{1-4{\delta }_n^2}\right)$$

(81)

with

$$E_n\equiv E_0\left[1+(n-1)\Delta z)\right],{\delta }_n\equiv {\xi }_n{\lambda }_{\gamma }^{\left(n\right)}\left(E_0\right)\left(\frac{\mathrm{nG}}{{10}^{11}\mathrm{GeV}}\right),$$

(82)

where ${\xi }_n$ is just $\xi$ as defined by Equation (68) and evaluated in the n-th domain.

5.3. 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_{\gamma \gamma }^{\mathrm{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 ${\left\{{\psi }_n\right\}}_{1\le n\le N_d}$—the whole transfer matrix describing the propagation of the photon/ALP beam is

$$\mathcal{U}\left(E_0,z;{\psi }_1,\dots ,{\psi }_{N_d}\right)=\prod _{n=1}^{N_d}{\mathcal{U}}_n\left(E_n,{\psi }_n\right).$$

(83)

Moreover, the probability that a photon/ALP beam emitted by a blazar at z in the state ${\rho }_1$ will be detected in the state ${\rho }_2$ for the above choice of ${\left\{{\psi }_n\right\}}_{1\le n\le N_d}$ is given by

$$P_{{\rho }_1\to {\rho }_2}\left(E_0,z;{\psi }_1,\dots ,{\psi }_{N_d}\right)=\mathrm{Tr}\left({\rho }_2\mathcal{U}\left(E_0,z;{\psi }_1,\dots ,{\psi }_{N_d}\right){\rho }_1{\mathcal{U}}^{†}\left(E_0,z;{\psi }_1,\dots ,{\psi }_{N_d}\right)\right)$$

(84)

with $\mathrm{Tr}{\rho }_1=\mathrm{Tr}{\rho }_2=1$.

Since the actual values of the angles ${\left\{{\psi }_n\right\}}_{1\le n\le 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

$$P_{{\rho }_1\to {\rho }_2}\left(E_0,z\right)={〈P_{{\rho }_1\to {\rho }_2}\left(E_0,z;{\psi }_1,\dots ,{\psi }_{N_d}\right)〉}_{{\psi }_1,\dots ,{\psi }_{N_d}},$$

(85)

indeed in accordance with the strategy outlined above. In practice, this is accomplished by evaluating the r.h.s. of Equation (84) over a very large number of realizations of the propagation process (we take 5000 realizations) randomly choosing the values of all angles ${\left\{{\psi }_n\right\}}_{1\le n\le N_d}$ for every realization, adding the results and dividing by the number of realizations.

Because the photon polarization cannot be measured at the considered energies, we have to sum the result over the two final polarization states

$${\rho }_x=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),{\rho }_z=\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right),$$

(86)

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

$${\rho }_{\mathrm{unpol}}=\frac{1}{2}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right).$$

(87)

We find in this way the photon survival probability $P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z\right)$

$$P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z\right)=\sum _{i=x,z}{〈P_{{\rho }_{\mathrm{unpol}}\to {\rho }_i}\left(E_0,z;{\psi }_1,\dots ,{\psi }_{N_d}\right)〉}_{{\psi }_1,\dots ,{\psi }_{N_d}}.$$

(88)

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.

6. 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.

• 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.
• 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.
• In order to get rid of evolutionary effects inside blazars we restrict our attention to those with $z\le 0.6$.
• 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\gtrsim 100\mathrm{GeV}$.

We are consequently left with a sample $\mathcal{S}$ of 39 flaring VHE BL Lacs, which are listed in Appendix C.

In first approximation, all observed spectra of the VHE blazars in $\mathcal{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

$${\Phi }_{\mathrm{obs}}(E_0,z)={\widehat{K}}_{\mathrm{obs}}\left(z\right){\left(\frac{E_0}{E_{\mathrm{ref}}}\right)}^{-{\Gamma }_{\mathrm{obs}}\left(z\right)},$$

(89)

where $E_0$ is the observed energy, $E_{\mathrm{ref}}$ is a common reference energy while ${\widehat{K}}_{\mathrm{obs}}\left(z\right)$ and ${\Gamma }_{\mathrm{obs}}\left(z\right)$ denote the normalization constant and the observed slope, respectively, for a source at redshift z. Actually, ${\widehat{K}}_{\mathrm{obs}}\left(z\right)$ 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 rescaling ${\widehat{K}}_{\mathrm{obs}}\left(z\right)\to K_{\mathrm{obs}}\left(z\right)$ in the observed spectrum of the considered blazars in such a way that $K_{\mathrm{obs}}\left(z\right)$ coincides with ${\Phi }_{\mathrm{obs}}(E_0,z)$ at the fiducial energy $E_{0,*}=300\mathrm{GeV}$ for every source in $\mathcal{S}$. Accordingly, Equation (89) becomes

$${\Phi }_{\mathrm{obs}}\left(E_0,z\right)=K_{\mathrm{obs}}\left(z\right){\left(\frac{E_0}{E_{0,*}}\right)}^{-{\Gamma }_{\mathrm{obs}}\left(z\right)}.$$

(90)

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

$${\Phi }_{\mathrm{obs}}(E_0,z)=e^{-{\tau }_{\gamma }(E_0,z)}{\Phi }_{\mathrm{em}}^{\mathrm{CP}}\left(E_0(1+z)\right),$$

(91)

where CP stands throughout this section for conventional physics.

6.1. Conventional Physics

Let us start by deriving the emitted spectrum of every source in $\mathcal{S}$, starting from each observed one. As a preliminary step, thanks to Equation (53) we rewrite Equation (91) as

$${\Phi }_{\mathrm{em}}^{\mathrm{CP}}\left(E_0(1+z)\right)=e^{{\tau }_{\gamma }(E_0,z)}{K}_{\mathrm{obs}}\left(z\right){\left(\frac{E_0}{E_{0,*}}\right)}^{-{\Gamma }_{\mathrm{obs}}\left(z\right)}.$$

(92)

Because of the presence of the exponential in the r.h.s. of Equation (92), ${\Phi }_{\mathrm{em}}^{\mathrm{CP}}\left(E_0(1+z)\right)$ cannot behave as an exact power law. Yet, it turns out to be close to it. Therefore, we best-fit (BF) ${\Phi }_{\mathrm{em}}^{\mathrm{CP}}\left(E_0(1+z)\right)$ to a single power-law expression

$${\Phi }_{\mathrm{em}}^{\mathrm{CP},\mathrm{BF}}\left(E_0(1+z)\right)=K_{\mathrm{em}}^{\mathrm{CP}}\left(z\right){\left(\frac{E_0(1+z)}{E_{0,*}}\right)}^{-{\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)}$$

(93)

over the energy range $\Delta E_0\left(z\right)$ where a source is observed, and so $E_0$ varies inside $\Delta E_0\left(z\right)$ (which changes from source to source). Correspondingly, the resulting values of ${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)$ are plotted in Figure 8.

We proceed by performing a statistical analysis of all values of ${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)$ 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 ${\chi }_{\mathrm{red},\mathrm{CP}}^2$. The values of the ${\chi }_{\mathrm{red},\mathrm{CP}}^2$ 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

$${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)=-5.33z^2-0.66z+2.64.$$

(94)

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 ${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)$ is the exponent of the emitted energy entering ${\Phi }_{\mathrm{em}}^{\mathrm{CP}}\left(E\right)$. Hence, in the two extreme cases $z=0$ and $z=0.6$ we have

$${\Phi }_{\mathrm{em}}^{\mathrm{CP}}(E,0)\propto E^{-2.64},{\Phi }_{\mathrm{em}}^{\mathrm{CP}}(E,0.6)\propto E^{-0.33},$$

(95)

thereby implying that the hardening of the emitted flux progressively increases with the redshift. More generally, we have found a statistical correlation between the $\left\{{\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)\right\}$ and z.

However, this result looks physically absurd. How can the sources get to know their z so as to tune their ${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)$ 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 ${\Gamma }_{\mathrm{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].

6.2. 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_{\mathrm{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) with ${\Phi }_{\mathrm{em}}^{\mathrm{CP}}\left(E_0(1+z)\right)\to {\Phi }_{\mathrm{em}}^{\mathrm{ALP}}\left(E_0(1+z)\right)$, keeping in mind that now ${\tau }^{\mathrm{CP}}\to {\tau }^{\mathrm{ALP}}$. Whence

$$\begin{array}{cc}& {\displaystyle {\Phi }_{\mathrm{em}}^{\mathrm{ALP}}\left(E_0(1+z)\right)={\left(P_{\gamma \to \gamma }^{\mathrm{ALP}}(E_0,z)\right)}^{-1}\times }\\ & {\displaystyle \times K_{\mathrm{obs}}\left(z\right){\left(\frac{E_0}{E_{0,*}}\right)}^{-{\Gamma }_{\mathrm{obs}}\left(z\right)},}\end{array}$$

(96)

Next, we still best-fit ${\Phi }_{\mathrm{em}}^{\mathrm{ALP}}\left(E_0(1+z)\right)$ to a single power law expression

$${\Phi }_{\mathrm{em}}^{\mathrm{ALP},\mathrm{BF}}\left(E_0(1+z)\right)=K_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right){\left(\frac{E_0(1+z)}{E_{0,*}}\right)}^{-{\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)}$$

(97)

over the energy range $\Delta E_0\left(z\right)$ where a source is observed, hence $E_0$ varies within $\Delta E_0\left(z\right)$. Such a best-fitting procedure is performed for every benchmark value of $\xi$ and $L_{\mathrm{dom}}$, namely $L_{\mathrm{dom}}=4\mathrm{Mpc}$, $\xi =0.1,0.5,1,5$, and $L_{\mathrm{dom}}=10\mathrm{Mpc}$, $\xi =0.1,0.5,1,5$.

Moreover, we carry out again the above statistical analysis of the values of ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)$ for all blazars in $\mathcal{S}$, for any benchmark value of $\xi$ and $L_{\mathrm{dom}}$.

Finally, the statistical significance of each fit can be quantified by computing the corresponding ${\chi }_{\mathrm{red},\mathrm{ALP}}^2$, whose values are reported in

Table 1. Values of ${\chi }_{\mathrm{red},\mathrm{CP}}^2$ in the case of conventional physics and ${\chi }_{\mathrm{red},\mathrm{ALP}}^2$ within the ALP scenario, for all blazars belonging to $\mathcal{S}$. The first column indicates the number of fit parameters, the second column concerns conventional physics and the third column refers to the ALP scenario for $L_{\mathrm{dom}}=4\mathrm{Mpc}$ and our benchmark values of $\xi$. The number in boldface corresponds to the minimum of ${\chi }_{\mathrm{red},\mathrm{ALP}}^2$.

# of Fit Parameters	${\mathit{\chi }}_{\mathbf{red},\mathbf{CP}}^2$	${\mathit{\chi }}_{\mathbf{red},\mathbf{ALP}}^2$
		$\xi =0.1$	$\xi =0.5$	$\xi =1$	$\xi =5$
1	2.37	2.29	1.29	1.31	1.43
2	1.49	1.47	1.29	1.31	1.38
3	1.46	1.46	1.32	1.31	1.37

for $L_{\mathrm{dom}}=4\mathrm{Mpc}$, $\xi =0.1,0.5,1,5$, and in

Table 2. Same as Table 1 but for $L_{\mathrm{dom}}=10\mathrm{Mpc}$.

# of Fit Parameters	${\mathit{\chi }}_{\mathbf{red},\mathbf{CP}}^2$	${\mathit{\chi }}_{\mathbf{red},\mathbf{ALP}}^2$
		$\xi =0.1$	$\xi =0.5$	$\xi =1$	$\xi =5$
1	2.37	2.05	1.25	1.39	1.43
2	1.49	1.44	1.26	1.37	1.38
3	1.46	1.46	1.28	1.36	1.37

for $L_{\mathrm{dom}}=10\mathrm{Mpc}$, $\xi =0.1,0.5,1,5$. In both tables the values of ${\chi }_{\mathrm{red},\mathrm{CP}}^2$ are reported for comparison.

The relevance of such a statistical analysis is to single out two preferred situations (corresponding to the minimum of ${\chi }_{\mathrm{red},\mathrm{ALP}}^2$): one for $L_{\mathrm{dom}}=4\mathrm{Mpc}$ and the other for $L_{\mathrm{dom}}=10\mathrm{Mpc}$. In either case, the results are $\xi =0.5$ and a straight best-fit regression line which is exactly horizontal. More in detail, for $L_{\mathrm{dom}}=4\mathrm{Mpc}$ we get ${\chi }_{\mathrm{red},\mathrm{ALP}}^2=1.29$ and ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}=2.54$, while for $L_{\mathrm{dom}}=10\mathrm{Mpc}$ we find ${\chi }_{\mathrm{red},\mathrm{ALP}}^2=1.25$ and ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}=2.60$.

Manifestly, both cases turn out to be very similar. We plot the values of ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)$ in Figure 10 only for the two considered situations.

Because $\xi =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=\mathcal{O}\left(100\right)\mathrm{GeV}$—the ALP mass must be $m=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$, since in Section 5.1 we have seen that ${\omega }_{\mathrm{pl}}<1.92·{10}^{-14}\mathrm{eV}$. Moreover, by recalling Equation (68) with $\xi =0.5$ and the upper bounds on $g_{a\gamma \gamma }$ and B quoted in Section 3.6 we get $2.94·{10}^{-12}{\mathrm{GeV}}^{-1}<g_{a\gamma \gamma }<0.66·{10}^{-10}{\mathrm{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 ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}$ 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!

6.3. 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 ${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)$ by a horizontal straight regression line, at the cost of relaxing the best-fitting requirement. Accordingly, the scatter of the values of ${\Gamma }_{\mathrm{em}}^{\mathrm{CP}}\left(z\right)$ for $95\%$ of blazars belonging to $\mathcal{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 ${\chi }_{\mathrm{red},\mathrm{CP}}^2=2.37$ which is by far too large.

The result obtained in the presence of $\gamma \leftrightarrow a$ oscillations and $\xi =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 ${\chi }_{\mathrm{red},\mathrm{ALP}}^2$ turns out to be considerably smaller. Specifically, the scatter of the values of ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)$ for $95\%$ of the considered blazars is now less than $13\%$ about the mean value set by ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}=2.54$ for $L_{\mathrm{dom}}=4\mathrm{Mpc}$ and less than $13\%$ about the mean value set by ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}=2.60$ for $L_{\mathrm{dom}}=10\mathrm{Mpc}$, namely equal to 0.33 for $L_{\mathrm{dom}}=4\mathrm{Mpc}$ and equal to 0.32 for $L_{\mathrm{dom}}=10\mathrm{Mpc}$. This situation is shown in Figure 12.

We argue that the small scatter in the values of ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)$ implies that the physical emission mechanism is the same for all flaring blazars, with the small fluctuations in ${\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)$ 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_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)$ 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 $\left\{{\Gamma }_{\mathrm{obs}}\left(z\right)\right\}$ distribution (see Appendix C) arises from the small scatter in the $\left\{{\Gamma }_{\mathrm{em}}^{\mathrm{ALP}}\left(z\right)\right\}$ distribution shown in Figure 10? The answer is very simple: most of the scatter in the $\left\{{\Gamma }_{\mathrm{obs}}\left(z\right)\right\}$ distribution arises from the large scatter in the source redshifts.

7. 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.

7.1. Detection of PKS 1222+216 in the VHE Range

The detection of an intense, rapidly varying emission in the energy range $70\mathrm{GeV}–400\mathrm{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 $\gamma$ rays above 20 GeV–30 GeV is expected (see e.g., [187,188,189]). However, photons up to $400\mathrm{GeV}$ have been observed by MAGIC [236]. In addition, the PKS 1222+216 flux doubles in only about 10 min, 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\sim {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\mathrm{GeV}–3\mathrm{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 $\gamma$-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.

7.2. Observations and Setup

The relevant physical parameters for PKS 1222+216 are as follows. We assume a disk luminosity $L_{\mathrm{disk}}\simeq 1.5·{10}^{46}$ erg s${}^{-1}$, a radius of the BLR $R_{\mathrm{BLR}}\simeq 0.23\mathrm{pc}$, and standard values for the cloud number density $n_c\simeq {10}^{10}{\mathrm{cm}}^{-3}$ and the temperature $T_c\simeq {10}^4\mathrm{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\simeq {10}^4{\mathrm{cm}}^{-3}$.

We now estimate the BLR absorption by evaluating the optical depth $\tau \left(E\right)$ of the beam photons inside the BLR according to conventional physics. By following the same procedure developed in [189], the optical depth reads

$$\tau \left(E\right)=\int d\Omega \int dϵ\int dx{n}_{\mathrm{ph}}(ϵ,\Omega ,x){\sigma }_{\gamma \gamma }(E,ϵ,\mu )(1-\mu ),$$

(98)

where E is the energy of a $\gamma$ ray, x is the distance from the center of the BLR, $\mu \equiv cos\theta$ where $\theta$ is the scattering angle between a $\gamma$ ray and a soft photon of energy $ϵ$ of the BLR, $d\Omega =-2\pi d\mu$, while $n_{\mathrm{ph}}(ϵ,\Omega ,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 ${\sigma }_{\gamma \gamma }(E,ϵ,\mu )$ is the pair-production cross-section of Equation (55).

The resulting $\tau \left(E\right)$ 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\mathrm{GeV}–400\mathrm{GeV}$. But MAGIC has detected such photons.

The calculated $\tau \left(E\right)$ is affected by some degree of uncertainty coming from the uncertainty of the input parameters, and in particular the luminosity of the disk $L_{\mathrm{disk}}$. However, since $\tau \propto L_{\mathrm{disk}}^{1/2}$ ($R_{\mathrm{BLR}}\propto L_{\mathrm{disk}}^{1/2}$, see [243] and references therein), the final impact of these uncertainties is moderate. In addition, scattered disk photons [240] show a maximal absorption of $\tau \simeq 0.2$ at about 200 GeV, so that their contribution to the total $\tau \left(E\right)$ can safely be neglected.

7.3. An ALP Model for PKS 1222+216

A natural explanation of the VHE observations arises if $\gamma \leftrightarrow 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 $\gamma \leftrightarrow 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 $\mathbf{B}$ is homogeneous in first approximation, (ii) the large scale jet where $\mathbf{B}$ possesses a smooth y-dependence, (iii) the host galaxy where $\mathbf{B}$ has a domain-like structure (as we shall see) and (iv) the extragalactic space where $\mathbf{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 ${\mathcal{L}}_{\mathrm{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 ${\mathcal{M}}_0$ in Equation (13) two new terms must be added, which read

$${\Delta }_{xx}^{\mathrm{QED}}(E,y)=\frac{2\alpha }{45\pi }{\left(\frac{B_T\left(y\right)}{B_{\mathrm{cr}}}\right)}^{2}E,$$

(99)

and

$${\Delta }_{zz}^{\mathrm{QED}}(E,y)=\frac{7\alpha }{90\pi }{\left(\frac{B_T\left(y\right)}{B_{\mathrm{cr}}}\right)}^{2}E,$$

(100)

respectively, where $B_{\mathrm{cr}}\simeq 4.41·{10}^{13}\mathrm{G}$ is the critical magnetic field. Thus, we can introduce also a high-energy cutoff

$$E_H\equiv \frac{90\pi }{7\alpha }{g}_{a\gamma \gamma }{B}_T{\left(\frac{B_{\mathrm{cr}}}{B_T}\right)}^2,$$

(101)

above which the photon/ALP beam does not propagate within the strong-mixing regime and presents an energy dependent behavior. We will now address $\gamma \leftrightarrow 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=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$ as in Section 5.

7.4. 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_{\mathrm{BLR}}\simeq 0.23\mathrm{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\simeq 0.2\mathrm{G}$ ($B\simeq 2.2\mathrm{G}$ at the base of the jet [244]) and an angle of ${45}^{\circ }$ with the beam direction, since $\gamma \leftrightarrow a$ oscillations vanish if $\mathbf{B}$ is exactly along the beam, while it is maximal for $\mathbf{B}$ transverse to the beam. Whence $B_T=0.14\mathrm{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 $\gamma \leftrightarrow a$ oscillations are relevant well above $E_H$.

We calculate the mean free path inside the BLR as

$${\lambda }_{\gamma }\left(E\right)=\frac{R_{\mathrm{BLR}}}{\tau \left(E\right)},$$

(102)

where $\tau \left(E\right)$ is the optical depth for the $\gamma \gamma \to 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 ${\mathcal{M}}_0$ of Equation (13) are

$${\Delta }_{xx}\left(E\right)=\frac{2\alpha E}{45\pi }{\left(\frac{B_T}{B_{\mathrm{cr}}}\right)}^2+\frac{i\tau \left(E\right)}{2R_{\mathrm{BLR}}}\simeq {10}^{-24}\left[\left(\frac{E}{\mathrm{GeV}}\right)+13.9i\tau \left(E\right)\right]\mathrm{eV},$$

(103)

$${\Delta }_{zz}\left(E\right)=\frac{3.5\alpha E}{45\pi }{\left(\frac{B_T}{B_{\mathrm{cr}}}\right)}^2+\frac{i\tau \left(E\right)}{2R_{\mathrm{BLR}}}\simeq {10}^{-24}\left[1.75\left(\frac{E}{\mathrm{GeV}}\right)+13.9i\tau \left(E\right)\right]\mathrm{eV},$$

(104)

$${\Delta }_{a\gamma }=\frac{1}{2}{g}_{a\gamma \gamma }{B}_T\simeq 1.37·{10}^{-23}\left(g_{a\gamma \gamma }{10}^{11}\mathrm{GeV}\right)\mathrm{eV},$$

(105)

$${\Delta }_{aa}\left(E\right)=0.$$

(106)

It is now possible to evaluate the transfer matrix ${\mathcal{U}}_1(R_{\mathrm{BLR}},0;E)$ in this region by means of the procedure developed in Appendix A.

7.5. Photon-ALP Oscillations in the Large Scale Jet

We now estimate the $\gamma \leftrightarrow a$ oscillations in the jet beyond $R_{\mathrm{BLR}}$, which we call region 2. In this zone $\mathbf{B}$ is believed to possess a toroidal behavior so that $B\left(y\right)\propto y^{-1}$ [245,246] and $B_T\left(y\right)$ reads

$$B_T\left(y\right)\simeq 0.14\left(\frac{R_{\mathrm{BLR}}}{y}\right)\mathrm{G}\simeq 3.22·{10}^{-2}\left(\frac{\mathrm{pc}}{y}\right)\mathrm{G}.$$

(107)

Region 2 extends up to $R_{*}$, which is defined as the distance where $B_T\left(y\right)$ in Equation (107) reaches the value assumed by the strength of turbulent magnetic field in the host elliptical galaxy, whose typical strength is $5\mathsf{\mu }\mathrm{G}$ (more about this, later). Accordingly, we get $R_{*}\simeq 6.7\mathrm{kpc}$.

By employing Equations (37) and (101) with the previous choice of the parameter values and $g_{a\gamma \gamma }=\mathcal{O}\left({10}^{-11}\right){\mathrm{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 ${\mathcal{M}}_0$ in Equation (13) are

$${\Delta }_{xx}\left(E\right)={\Delta }_{zz}\left(E\right)={\Delta }_{aa}\left(E\right)=0,$$

(108)

$${\Delta }_{a\gamma }=\frac{1}{2}{g}_{a\gamma \gamma }{B}_T\simeq 3.1·{10}^{-24}\left(\frac{\mathrm{pc}}{y}\right)\left(g_{a\gamma \gamma }{10}^{11}\mathrm{GeV}\right)\mathrm{eV}.$$

(109)

In the present situation, the transfer matrix can be obtained by analytically solving the beam propagation equation and reads

$$\begin{array}{cc}& {\displaystyle {\mathcal{U}}_2(R_{*},R_{\mathrm{BLR}};E)=}\\ & {\displaystyle =\left(\begin{array}{ccc}1& 0& 0\\ 0& \cos \left(\frac{g_{a\gamma \gamma }{B}_T\left(R_{\mathrm{BLR}}\right)R_{\mathrm{BLR}}}{2}\ln \left(\frac{R_{*}}{R_{\mathrm{BLR}}}\right)\right)& i\sin \left(\frac{g_{a\gamma \gamma }{B}_T\left(R_{\mathrm{BLR}}\right)R_{\mathrm{BLR}}}{2}\ln \left(\frac{R_{*}}{R_{\mathrm{BLR}}}\right)\right)\\ 0& i\sin \left(\frac{g_{a\gamma \gamma }{B}_T\left(R_{\mathrm{BLR}}\right)R_{\mathrm{BLR}}}{2}\ln \left(\frac{R_{*}}{R_{\mathrm{BLR}}}\right)\right)& \cos \left(\frac{g_{a\gamma \gamma }{B}_T\left(R_{\mathrm{BLR}}\right)R_{\mathrm{BLR}}}{2}\ln \left(\frac{R_{*}}{R_{\mathrm{BLR}}}\right)\right)\end{array}\right),}\end{array}$$

(110)

where $B_T\left(R_{\mathrm{BLR}}\right)$ is given by Equation (107).

7.6. Photon-ALP Oscillations in the Host Galaxy

As anticipated, FSRQs are usually located in elliptical galaxies, whose $\mathbf{B}$ is known with great uncertainty. Nevertheless, it is believed that $\mathbf{B}$ possesses a turbulent nature and can be modeled with a domain-like structure, with strength $5\mathsf{\mu }\mathrm{G}$ and domain size equal to $150\mathrm{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_{\mathrm{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 ${\mathcal{M}}_0$ in Equation (13) become simply

$${\Delta }_{xx}\left(E\right)={\Delta }_{zz}\left(E\right)={\Delta }_{aa}\left(E\right)=0,$$

(111)

$${\Delta }_{a\gamma }=\frac{1}{2}{g}_{a\gamma \gamma }{B}_T\simeq 3.4·{10}^{-28}\left(g_{a\gamma \gamma }{10}^{11}\mathrm{GeV}\right)\mathrm{eV},$$

(112)

while the transfer matrix in this region ${\mathcal{U}}_3\left(R_{\mathrm{host}},R_{*};E\right)$ 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 $\gamma \leftrightarrow a$ oscillations in region 3 are totally negligible.

7.7. Photon-ALP Oscillations in Extragalactic Space

Finally, the photon/ALP beam propagates in extragalactic space from $R_{\mathrm{host}}$ to us. We call this region 4. It goes without saying that the treatment of $\gamma \leftrightarrow 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 ${\mathcal{U}}_4\left(D_L,R_{\mathrm{host}};E\right)$, where $D_L$ denotes the luminosity distance of PKS 1222+216 (corresponding to $z=0.432)$.

7.8. Overall Photon-ALP Oscillations

The whole transfer matrix ${\mathcal{U}}_{\mathrm{tot}}\left(D_L,0;E\right)$ 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

$${\mathcal{U}}_{\mathrm{tot}}\left(D_L,0;E\right)={\mathcal{U}}_4\left(D_L,R_{\mathrm{host}};E\right){\mathcal{U}}_3\left(R_{\mathrm{host}},R_{*};E\right){\mathcal{U}}_2(R_{*},R_{\mathrm{BLR}};E){\mathcal{U}}_1(R_{\mathrm{BLR}},0;E).$$

(113)

Now, let us consider the total probability that a photon/ALP beam emitted by the considered FSRQ in the state ${\rho }_1$ will be detected in the state ${\rho }_2$ for a fixed configuration of the $\mathbf{B}$ direction in each domain of the host galaxy and of the extragalactic space. Unless otherwise stated, we suppose that the angles ${\left\{{\phi }_m\right\}}_{1\le m\le N_g}$ and ${\left\{{\psi }_n\right\}}_{1\le n\le N_d}$ of $\mathbf{B}$ in the magnetic domains of the host galaxy and if extragalactic space, respectively are held fixed. Then Equation (49) entails

$$P_{\mathrm{tot},{\rho }_1\to {\rho }_2}\left(D_L,0;E\right)=\mathrm{Tr}\left({\rho }_2{\mathcal{U}}_{\mathrm{tot}}\left(D_L,0;E\right){\rho }_1{\mathcal{U}}_{\mathrm{tot}}^{†}\left(D_L,0;E\right)\right)$$

(114)

with $\mathrm{Tr}{\rho }_1=\mathrm{Tr}{\rho }_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 $\left(N_g{N}_d\right)$-dimensional stochastic process. As already discussed in Section 5, the total photon survival probability is therefore given by

$$P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z\right)=\sum _{i=x,z}{〈P_{{\rho }_{\mathrm{unpol}}\to {\rho }_i}\left(E_0,z;{\phi }_1,\dots ,{\phi }_{N_g};{\psi }_1,\dots ,{\psi }_{N_d}\right)〉}_{{\phi }_1,\dots ,{\phi }_{N_g};{\psi }_1,\dots ,{\psi }_{N_d}},$$

(115)

where for clarity we have explicitly written all angles that are averaged over.

7.9. Results

We now want to clarify the role of the $\gamma \leftrightarrow 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

$$\mathcal{U}\left(R_{*},0;E\right)={\mathcal{U}}_2(R_{*},R_{\mathrm{BLR}};E){\mathcal{U}}_1(R_{\mathrm{BLR}},0;E).$$

(116)

The corresponding photon survival probability in the presence of $\gamma \leftrightarrow a$ oscillations is given by (65), which presently can be written as

$$P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(R_{\mathrm{host}},0;E\right)=e^{-{\tau }_{\gamma }^{\mathrm{ALP}}\left(E\right)}.$$

(117)

We have considered three benchmark cases for B and $g_{a\gamma \gamma }$.

• $B=0.2\mathrm{G},g_{a\gamma \gamma }=1.4·{10}^{-11}{\mathrm{GeV}}^{-1}$.
• $B=0.4\mathrm{G},g_{a\gamma \gamma }=0.7·{10}^{-11}{\mathrm{GeV}}^{-1}$.
• $B=2.0\mathrm{G},g_{a\gamma \gamma }=0.25·{10}^{-11}{\mathrm{GeV}}^{-1}$.

The curves for ${\tau }_{\gamma }^{\mathrm{ALP}}\left(E\right)$ 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 $\gamma \leftrightarrow 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 ${\tau }_{\gamma }^{\mathrm{ALP}}\simeq 4$, which corresponds to a photon survival probability of about $2\%$. Instead, in the optically thin region below ∼30 $\mathrm{GeV}$, the optical depth with $\gamma \leftrightarrow 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

$$\mathcal{P}\equiv \log \left(\frac{P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(R_{*},0;1\mathrm{GeV}\right)}{P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(R_{*},0;300\mathrm{GeV}\right)}\right),$$

(118)

measuring the ratio between $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ at $1\mathrm{GeV}$, which is representative of Fermi/LAT observations, and $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ at $300\mathrm{GeV}$ which accounts for MAGIC ones. In order to get an acceptable shape for the emitted SED at VHE, we need to have $\mathcal{P}$ as low as possible. Indeed, a low value of $\mathcal{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 $\gamma$-ray band. The emitted spectrum ${\Phi }_{\mathrm{em}}\left(E_0(1+z)\right)$ is linked to the observed one ${\Phi }_{\mathrm{obs}}\left(E_0\right)$ by

$${\Phi }_{\mathrm{em}}\left(E_0(1+z)\right)=\frac{{\Phi }_{\mathrm{obs}}(E_0,z)}{P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(D_L,0;E_0\right)},$$

(119)

where z is the redshift of PKS 1222+216 and $P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(D_L,0;E_0\right)$ 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 $\gamma$-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 $\Gamma$, electron normalization K, minimum, break and maximum Lorentz factors ${\gamma }_{\min }$, ${\gamma }_{\mathrm{b}}$, ${\gamma }_{\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 $\gamma$-ray emitting blob we take $r=2.2·{10}^{14}$ cm, $B=0.008\mathrm{G}$, $\Gamma =17.5$, $K=6.7·{10}^9$, ${\gamma }_{\min }=3·{10}^3$, ${\gamma }_{\mathrm{b}}=1.2·{10}^5$, ${\gamma }_{\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_{\gamma }=6·{10}^{48}\mathrm{erg}{\mathrm{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_{\gamma }={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 $\gamma \leftrightarrow 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).

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=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$. Note that our model can be applied to all other VHE FSRQs with analogous results.

8. 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_{\mathrm{osc}}$ is much larger than the size $L_{\mathrm{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_{\mathrm{osc}}$ to decrease as E further increases. Therefore, things completely change whenever $L_{\mathrm{osc}}\lesssim L_{\mathrm{dom}}$, since in this case a full oscillation—or even several oscillations—probe 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 $\mathbf{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_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z;\varphi \left(y\right),\theta \left(y\right)\right)$, where $\varphi \left(y\right)$ and $\theta \left(y\right)$ are the two angles that fix the direction of $\mathbf{B}\left(y\right)$ 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 three-dimensional 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 $\mathbf{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 $B_T={\left(2/3\right)}^{1/2}B$.

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 $\Pi \left(y\right)$ perpendicular to the beam rather than being spread out throughout the full three-dimensional space. As a consequence, $P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z;\varphi \left(y\right),\theta \left(y\right)\right)\to P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E_0,z;\varphi \left(y\right)\right)$, where $\varphi \left(y\right)$ is the angle between ${\mathbf{B}}_T\left(y\right)$ and a fixed fiducial z-direction equal in all domains (namely in all planes $\Pi \left(y\right)$).

8.1. 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

$$\begin{array}{cc}& {\displaystyle \mathcal{M}(E,y)\equiv }\\ & {\displaystyle \equiv \left(\begin{array}{ccc}{\Delta }_{\mathrm{CMB}}\left(E\right)+{\Delta }_{\mathrm{abs}}\left(E\right)+{\Delta }_{\mathrm{pl}}\left(E\right)& 0& {\Delta }_{a\gamma }\sin \varphi \left(y\right)\\ 0& {\Delta }_{\mathrm{CMB}}\left(E\right)+{\Delta }_{\mathrm{abs}}\left(E\right)+{\Delta }_{\mathrm{pl}}\left(E\right)& {\Delta }_{a\gamma }\cos \varphi \left(y\right)\\ {\Delta }_{a\gamma }\sin \varphi \left(y\right)& {\Delta }_{a\gamma }\cos \varphi \left(y\right)& {\Delta }_{aa}\left(E\right).\end{array}\right).}\end{array}$$

(120)

The meaning of the terms of $\mathcal{M}(E,y)$ is as follows. The contribution from photon dispersion on the CMB is ${\Delta }_{\mathrm{CMB}}\left(E\right)=0.522·{10}^{-42}E$ [145], the contribution from the EBL absorption is ${\Delta }_{\mathrm{abs}}\left(E\right)=i/\left(2{\lambda }_{\gamma }\left(E\right)\right)$ where ${\lambda }_{\gamma }\left(E\right)$ 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 ${\Delta }_{\mathrm{pl}}\left(E\right)=-{\omega }_{\mathrm{pl}}^2/\left(2E\right)$ while the remaining terms are ${\Delta }_{a\gamma }=g_{a\gamma \gamma }{B}_T/2$ and ${\Delta }_{aa}\left(E\right)=-m_a^2/\left(2E\right)$, just like in Section 3.2.

The other fact is that we now have three regimes, separated by the low-energy threshold

$$E_L\equiv \frac{|m_a^2-{\omega }_{\mathrm{pl}}^2|}{2g_{a\gamma \gamma }{B}_T}\simeq \frac{2.56}{\xi }\left|{\left(\frac{m_a}{\mathrm{neV}}\right)}^2-{\left(\frac{{\omega }_{\mathrm{pl}}}{\mathrm{neV}}\right)}^2\right|\mathrm{TeV},$$

(121)

which is equal to Equation (37), and the high-energy threshold

$$E_H\equiv 1.92·{10}^{42}{g}_{a\gamma \gamma }{B}_T\simeq 3.74·{10}^2\xi \mathrm{GeV}.$$

(122)

Specifically we have:

• $E<E_L$—This is the low-energy weak-mixing regime, wherein the terms $\propto E^{-1}$ dominate. Correspondingly, we have

  $$L_{\mathrm{osc}}\left(E\right)\simeq \frac{4\pi E}{|m_a^2-{\omega }_{\mathrm{pl}}^2|},$$

  (123)

  and

  $$P_{\gamma \to a}(E,L_{\mathrm{dom}})\simeq {\left(\frac{2g_{a\gamma \gamma }{B}_{T}E}{|m_a^2-{\omega }_{\mathrm{pl}}^2|}\right)}^2{\sin }^2\left(\frac{|m_a^2-{\omega }_{\mathrm{pl}}^2|L_{\mathrm{dom}}}{4E}\right).$$

  (124)
  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=\mathrm{constant}$ term dominates. Accordingly, we obtain

  $$L_{\mathrm{osc}}\simeq \frac{2\pi }{g_{a\gamma \gamma }{B}_T}\simeq 2.05·{10}^2{\xi }^{-1}\mathrm{Mpc},$$

  (125)

  $$P_{\gamma \to a}\left(L_{\mathrm{dom}}\right)\simeq {\sin }^2\left(\frac{g_{a\gamma \gamma }{B}_T{L}_{\mathrm{dom}}}{2}\right)\simeq {\sin }^2\left[1.54·{10}^{-2}\xi \left(\frac{L_{\mathrm{dom}}}{\mathrm{Mpc}}\right)\right].$$

  (126)
  Clearly, $L_{\mathrm{osc}}$ and $P_{\gamma \to a}\left(L_{\mathrm{dom}}\right)$ are independent of all the energy-dependent terms, and $P_{\gamma \to a}\left(L_{\mathrm{dom}}\right)$ 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\gamma \gamma }{B}_T$. Correspondingly, we get

  $$L_{\mathrm{osc}}\left(E\right)\simeq \frac{1.20·{10}^{43}}{E}\simeq 76.15\left(\frac{\mathrm{TeV}}{E}\right)\mathrm{Mpc},$$

  (127)

  $$P_{\gamma \to a}(E,L_{\mathrm{dom}})\simeq 1.39·{10}^{-1}{\xi }^2{\left(\frac{\mathrm{TeV}}{E}\right)}^2{\sin }^2\left[4.12·{10}^{-2}\left(\frac{L_{\mathrm{dom}}}{\mathrm{Mpc}}\right)\left(\frac{E}{\mathrm{TeV}}\right)\right].$$

  (128)
  Manifestly, $L_{\mathrm{osc}}\left(E\right)$ decreases with increasing E and $P_{\gamma \to a}(E,L_{\mathrm{dom}})$ exhibits oscillations in E: this means that the individual realizations of the beam propagation are also oscillating functions of E. Moreover—since $P_{\gamma \to a}(E,L_{\mathrm{dom}})\propto E^{-2}$—as E increases the photon-ALP oscillations become unobservable at some point.

8.2. 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 ${\left\{y_{D,n}\right\}}_{0\le n\le 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\le n\le N_d$) towards the blazar.

Because of our ignorance about the strength of $\mathbf{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 $\mathbf{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 ${\mathbf{B}}_T\left(y\right)$. Therefore, we denote by ${\left\{{\varphi }_n\right\}}_{1\le n\le N_d}$ the set of angles that ${\mathbf{B}}_T\left(y\right)$ forms with the fixed fiducial z-direction in the middle of every domain. Thanks to the previous assumptions also $B_T$—but not ${\mathbf{B}}_T$—can be taken as constant in all domains.

Given the fact that ${\mathbf{B}}_T\left(y\right)$ changes randomly from one domain to the next, in order for ${\mathbf{B}}_T\left(y\right)$ 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 ${\mathbf{B}}_T\left(y\right)$ is homogeneous in the central part, but as the distance from the edge with the $(n+1)$-th domain decreases we assume that ${\mathbf{B}}_T\left(y\right)$ linearly changes thereby becoming equal to ${\mathbf{B}}_T\left(y\right)$ on the same edge but in the $(n+1)$-th domain. Accordingly, the continuity of the components of ${\mathbf{B}}_T\left(y\right)$ 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

$$y_{0,n}\equiv y_{D,n}-\frac{\sigma }{2}\left(y_{D,n}-y_{D,n-1}\right),(1\le n\le N-1),$$

(129)

$$y_{1,n}\equiv y_{D,n}+\frac{\sigma }{2}\left(y_{D,n+1}-y_{D,n}\right),(1\le n\le N-1),$$

(130)

where $\sigma \in [0,1]$ is the smoothing parameter. The interval $[y_{0,n},y_{1,n}]$ is the region where the angle $\varphi \left(y\right)$ changes smoothly from the value ${\varphi }_{0,n}\equiv {\varphi }_n$ in the n-th domain to the value ${\varphi }_{0,n+1}\equiv {\varphi }_{n+1}$ in the ($n+1$)-th domain. Manifestly, for $\sigma =0$ we have $y_{0,n}=y_{1,n}$, and we recover the DLSHE model. On the other hand, for $\sigma =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 $\varphi$ in any domain. The general case is intermediate—represented by a value of $0<\sigma <1$—so that in the central part of a domain the angle is constant (${\varphi }_{0,n}$) and then it linearly joins the value of the constant angle in the next domain (${\varphi }_{1,n}$). Hence, in a generic interval $[y_{1,n-1},y_{1,n}](1\le n\le N_1)$ we have

$$\varphi \left(y\right)=\left\{\begin{array}{cc}{\varphi }_{0,n}=\mathrm{constant},\hfill & y\in [y_{1,n-1},y_{0,n}],\hfill \\ {\varphi }_{0,n}+{\displaystyle \frac{{\varphi }_{0,n+1}-{\varphi }_{0,n}}{y_{1,n}-y_{0,n}}}(y-y_{0,n}),\hfill & y\in [y_{0,n},y_{1,n}].\hfill \end{array}\right.$$

(131)

According to our conventions, the blazar redshift is $z\equiv 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 ${\overline{z}}_n\equiv \left(z_{n-1}+z_n\right)/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 ${\overline{E}}_n\equiv \left(E_{n-1}+E_n\right)/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}$.

8.3. 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_{\mathrm{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 ${\left\{L_{\mathrm{dom},n}\right\}}_{1\le n\le N}$. Owing to the properties of $\mathbf{B}$ at redshift $z=0$, we take for the probability density in question the power law $\propto L_{\mathrm{dom}}^{-1.2}$ inside the range $0.2\mathrm{Mpc}–10\mathrm{Mpc}$, which means that $\langle L_{\mathrm{dom}}\rangle =2\mathrm{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

$${\lambda }_{\gamma ,n}=\frac{L_{\mathrm{dom},\mathrm{n}}}{{\tau }_{\gamma }(E_{N_d},z_{n-1})-{\tau }_{\gamma }(E_{N_d},z_n)},$$

(132)

with ${\tau }_{\gamma }$ again given by the FR model, and

$$B_{T,n}\left(y\right)=\left(B_{T,N_d}\left(y\right)\right){\left(1+{\overline{z}}_n\right)}^2,$$

(133)

where $B_{T,N_d}\left(y\right)$ is the strength of ${\mathbf{B}}_T\left(y\right)$ in the local Universe, namely in the domain closest to the observer ($z=0$).

So, the above mixing matrix $\mathcal{M}(E,y)$ in a single n-th domain is fully determined, and now has $E\to E_n=E_{N_d}\left(1+{\overline{z}}_n\right)$ 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).

8.4. 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

$$\begin{array}{cc}& {\displaystyle {\mathcal{U}}_n\left({\overline{E}}_n;z_n,z_{n-1};{\varphi }_{1,n},{\varphi }_{1,n-1}\right)=}\\ & {\displaystyle ={\mathcal{U}}_{\mathrm{var},n}\left({\overline{E}}_n;z_n,z_{n-1};{\varphi }_{1,n},{\varphi }_{0,n}\right){\mathcal{U}}_{\mathrm{const},n}\left({\overline{E}}_n;z_n,z_{n-1};{\varphi }_{0,n}\right),}\end{array}$$

(134)

for an arbitrary choice of the angle ${\varphi }_{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\to {\overline{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

$$\begin{array}{cc}& {\displaystyle {\mathcal{U}}_T\left(E_{N_d};z;{\left\{{\varphi }_n\right\}}_{1\le n\le N_d}\right)={\mathcal{U}}_{\mathrm{const},N_d}\left(E_{N_d};z_{N_d},z_{N_{d-1}};{\varphi }_{0,N_d}\right)\times }\\ & {\displaystyle \times \prod _{n=1}^{N_{d-1}}{\mathcal{U}}_{\mathrm{var},n}\left({\overline{E}}_n;z_n,z_{n-1};{\varphi }_{1,n},{\varphi }_{0,n}\right){\mathcal{U}}_{\mathrm{const},n}\left({\overline{E}}_n;z_n,z_{n-1};{\varphi }_{0,n}\right).}\end{array}$$

(135)

Note that this product must be ordered in such a way that the transfer matrices with smaller n must be closer to the source.

8.5. 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}\to E_0$, $z_{N_d}\to 0$, $z_0\to z$. Accordingly, we have

$$\begin{array}{cc}& {\displaystyle P_{\gamma \to \gamma ,\mathrm{unp}}^{\mathrm{ALP}}\left(E_0;{\rho }_x,{\rho }_z;z,{\rho }_{\mathrm{unp}};{\left\{{\varphi }_n\right\}}_{1\le n\le N_d}\right)=}\\ & {\displaystyle =\sum _{i=x,z}\mathrm{Tr}\left[{\rho }_i{\mathcal{U}}_T\left(E_0;z;{\left\{{\varphi }_n\right\}}_{1\le n\le N_d}\right){\rho }_{\mathrm{unp}}{\mathcal{U}}_T^{†}\left(E_0;z;{\left\{{\varphi }_n\right\}}_{1\le n\le N}\right)\right]}\end{array}$$

(136)

with

$${\rho }_x\equiv \left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),{\rho }_z\equiv \left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right),{\rho }_{\mathrm{unp}}\equiv \frac{1}{2}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right).$$

(137)

Below, the photon survival probability $P_{\gamma \to \gamma ,\mathrm{unp}}^{\mathrm{ALP}}\left(E_0;{\rho }_x,{\rho }_z;z,{\rho }_{\mathrm{unp}};{\left\{{\varphi }_n\right\}}_{1\le n\le N}\right)$ 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 $\xi =0.5,1,2,5$. In order to simplify notations, we will denote $P_{\gamma \to \gamma ,\mathrm{unp}}^{\mathrm{ALP}}\left(E_0;{\rho }_x,{\rho }_z;z,{\rho }_{\mathrm{unp}};{\left\{{\varphi }_n\right\}}_{1\le n\le N}\right)$ as $P_{\gamma \to \gamma }^{\mathrm{ALP}}(E_0,z)$. Thousand random realizations of the propagation process have been considered for each choice of z, $\xi$, $E_0$. In all figures a random distribution of the domain length $L_{\mathrm{dom}}$ has been taken: a power law distribution function has been chosen with exponent $\alpha =-1.2$ and domain length in the interval between the minimal value $L_{\mathrm{dom}}^{\min }=0.2\mathrm{Mpc}$ and the maximal value $L_{\mathrm{dom}}^{\max }=10\mathrm{Mpc}$. The resulting average domain length is $\langle L_{\mathrm{dom}}\rangle =2\mathrm{Mpc}$. Our results are shown in Figure 17, Figure 18, Figure 19, Figure 20, Figure 21, Figure 22 and Figure 23. We take the smoothing parameter $\sigma =0.2$ for the transition between two adjacent magnetic domains. The dotted-dashed 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 $\xi =0.5$, in the upper-right panel $\xi =1$, in the lower-left panel $\xi =2$ and in the lower-right panel $\xi =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].

9. 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=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$ and $g_{a\gamma \gamma }=\mathcal{O}\left({10}^{-11}\right){\mathrm{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 $\gamma \to a$ inside a blazar and the reconversion $a\to \gamma$ in the Milky Way, neglecting any possible $\gamma \leftrightarrow 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\simeq 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.

9.1. Propagation in the BL Lac Jet

We denote by ${\mathcal{R}}_{\mathrm{VHE}}$ the region where the VHE photons originate inside the BL Lac jet, with $y_{\mathrm{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 ${\mathcal{R}}_{\mathrm{jet}}$ between $y_{\mathrm{VHE}}$ and the end of the jet $y_{\mathrm{jet}}$, which we denote as ${\mathcal{U}}_{{\mathcal{R}}_{\mathrm{jet}}}(E;y_{\mathrm{jet}},y_{\mathrm{VHE}})$.

The region ${\mathcal{R}}_{\mathrm{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_{\mathrm{VHE}}$ from the central SMBH, (2) the transverse magnetic field profile $B_{T,{\mathcal{R}}_{\mathrm{jet}}}\left(y\right)$ from $y_{\mathrm{VHE}}$ to $y_{\mathrm{jet}}$, (3) the electron density profile $n_{e,{\mathcal{R}}_{\mathrm{jet}}}\left(y\right)$ from $y_{\mathrm{VHE}}$ to $y_{\mathrm{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 ${\mathcal{R}}_{\mathrm{VHE}}$ we find $0.1\mathrm{G}\lesssim B_{T,{\mathcal{R}}_{\mathrm{VHE}}}\lesssim 1\mathrm{G}$ and for definiteness we choose $B_{T,{\mathcal{R}}_{\mathrm{VHE}}}=0.5\mathrm{G}$. Moreover, we get $n_{e,{\mathcal{R}}_{\mathrm{VHE}}}\simeq 5·{10}^4{\mathrm{cm}}^{-3}$, leading in turn to a plasma frequency of ${\omega }_{\mathrm{pl}}\simeq 8.25·{10}^{-9}\mathrm{eV}$, thanks to Equation (20). Although there is no direct way to infer a precise value of $y_{\mathrm{VHE}}$, this quantity can be estimated from the size of ${\mathcal{R}}_{\mathrm{VHE}}$—which is assumed as measure of the jet cross-section— thus finding ${10}^{16}\mathrm{cm}\lesssim y_{\mathrm{VHE}}\lesssim {10}^{17}\mathrm{cm}$. For definiteness, we shall take $y_{\mathrm{VHE}}\simeq 3·{10}^{16}\mathrm{cm}$. Once produced, VHE photons propagate unimpeded out to $y_{\mathrm{jet}}\simeq 1\mathrm{kpc}$ where they leave the jet, entering the host galaxy. Within ${\mathcal{R}}_{\mathrm{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

$$B_{T,{\mathcal{R}}_{\mathrm{jet}}}\left(y\right)=B_{T,{\mathcal{R}}_{\mathrm{VHE}}}\left(\frac{y_{\mathrm{VHE}}}{y}\right).$$

(138)

Concerning the electron density profile, due to the conical shape of the jet our expectation is

$$n_{e,{\mathcal{R}}_{\mathrm{jet}}}\left(y\right)=n_{e,{\mathcal{R}}_{\mathrm{VHE}}}{\left(\frac{y_{\mathrm{VHE}}}{y}\right)}^2.$$

(139)

The knowledge of the above quantities allows us to compute the entire propagation process of the photon/ALP beam within the jet, namely ${\mathcal{U}}_{{\mathcal{R}}_{\mathrm{jet}}}(E;y_{\mathrm{jet}},y_{\mathrm{VHE}})$.

It should be kept in mind that in ${\mathcal{R}}_{\mathrm{jet}}$ we consider the photon/ALP beam in a frame co-moving with the jet, so that we must apply the transformation $E\to \gamma E$ to the beam in order to go to a fixed frame—as it will be performed in the next regions—with $\gamma$ being the Lorentz factor. We take $\gamma =15$.

9.2. Propagation in the Host Galaxy

All the three considered blazars are hosted by elliptical galaxies, which we denote by ${\mathcal{R}}_{\mathrm{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 $\gamma \leftrightarrow a$ oscillations is totally negligible. Therefore, denoting by $y_{\mathrm{in},\mathrm{host}}\equiv y_{\mathrm{jet}}$ and by $y_{\mathrm{out},\mathrm{host}}$ the points on the y axis where the beam enters and exits from the host galaxy, respectively, we have ${\mathcal{U}}_{{\mathcal{R}}_{\mathrm{host}}}(E;y_{\mathrm{out},\mathrm{host}},y_{\mathrm{in},\mathrm{host}})=1$.

9.3. Propagation in the Extragalactic Space

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

9.4. Propagation in the Milky Way

We denote by ${\mathcal{R}}_{\mathrm{MW}}$ the region where the photon/ALP beam propagates inside the Milky Way, i.e., from $y_{\mathrm{MW}}$ up to the Earth, whose position is denoted by $y_{\oplus }$.

By closely following the strategy described in [74], we compute ${\mathcal{U}}_{{\mathcal{R}}_{\mathrm{MW}}}(E;y_{\oplus },y_{\mathrm{MW}})$. In order to take into account the structured behaviour of the Galactic magnetic field ${\mathbf{B}}_{\mathrm{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 ${\mathbf{B}}_{\mathrm{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 $\gamma \leftrightarrow 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\simeq 1.1·{10}^{-2}{\mathrm{cm}}^{-3}$, resulting in a plasma frequency ${\omega }_{\mathrm{pl}}\simeq 3.9·{10}^{-12}\mathrm{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 ${\mathcal{U}}_{{\mathcal{R}}_{\mathrm{MW}}}(E;y_{\oplus },y_{\mathrm{MW}})$ for an arbitrary direction of the line of sight to a given blazar.

9.5. Overall Photon Survival Probability

Because all the transfer matrices in each region are now known, the total transfer matrix $\mathcal{U}(E;y_{\oplus },y_{\mathrm{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

$$\begin{array}{cc}& {\displaystyle \mathcal{U}(E;y_{\oplus },y_{\mathrm{VHE}})={\mathcal{U}}_{{\mathcal{R}}_{\mathrm{MW}}}(E;y_{\oplus },y_{\mathrm{MW}})\times }\end{array}$$

(140)

$$\begin{array}{cc}& {\displaystyle \times {\mathcal{U}}_{{\mathcal{R}}_{\mathrm{ext}}}(E;y_{\mathrm{MW}},y_{\mathrm{out},\mathrm{host}}){\mathcal{U}}_{{\mathcal{R}}_{\mathrm{host}}}(E;y_{\mathrm{out},\mathrm{host}},y_{\mathrm{in},\mathrm{host}})\times }\\ \hfill & {\displaystyle \times {\mathcal{U}}_{{\mathcal{R}}_{\mathrm{jet}}}(E;y_{\mathrm{in},\mathrm{host}},y_{\mathrm{VHE}}),}\end{array}$$

(141)

where of course we have $y_{\mathrm{in},\mathrm{host}}\equiv y_{\mathrm{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 $\rho \left(y\right)={(A_x\left(y\right),A_z\left(y\right),a\left(y\right))}^T\otimes (A_x\left(y\right),A_z\left(y\right),a\left(y\right))$. As a consequence, the overall photon survival probability becomes

$$\begin{array}{cc}& {\displaystyle P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E;y_{\oplus },{\rho }_x,{\rho }_z;y_{\mathrm{VHE}},{\rho }_{\mathrm{unp}}\right)=}\\ & {\displaystyle =\sum _{i=x,z}\mathrm{Tr}\left[{\rho }_i\mathcal{U}\left(E;y_{\oplus },y_{\mathrm{VHE}}\right){\rho }_{\mathrm{unp}}{\mathcal{U}}^{†}\left(E;y_{\oplus },y_{\mathrm{VHE}}\right)\right],}\end{array}$$

(142)

where

$${\rho }_x\equiv \left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),{\rho }_z\equiv \left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right),{\rho }_{\mathrm{unp}}\equiv \frac{1}{2}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right).$$

(143)

and

$${\rho }_{\mathrm{unpol}}=\frac{1}{2}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right).$$

(144)

Below—merely for notational convenience—we shall replace $P_{\gamma \to \gamma }^{\mathrm{ALP}}(E;y_{\oplus },{\rho }_x,{\rho }_z;y_{\mathrm{VHE}},$${\rho }_{\mathrm{unp}})$ simply by $P_{\gamma \to \gamma }^{\mathrm{ALP}}\left(E,z\right)$.

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_{\mathrm{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_{\mathrm{osc}}$ versus E is strongly affected by the value of $B_{T,{\mathcal{R}}_{\mathrm{jet}}}\left(y\right)$: as expected (see also [114]), as $B_{T,{\mathcal{R}}_{\mathrm{jet}}}\left(y\right)$ decreases—when the distance from the emission region increases—the maximal value of $L_{\mathrm{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_{\mathrm{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_{\mathrm{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.

9.6. 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 $\nu F_{\nu }$ in the presence of $\gamma \leftrightarrow 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_{\mathrm{cut}}$ as

$${\Phi }_{\mathrm{int}}\left(E\right)={\Phi }_0{\left(\frac{E}{E_0}\right)}^{-\Gamma }{e}^{-E/E_{\mathrm{cut}}},$$

(145)

where ${\Phi }_0$ is a normalization constant accounting for the blazar luminosity, $E_0$ is a reference energy and $\Gamma$ 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 $\gamma \leftrightarrow a$ oscillations are at work. In order to obtain the SED we take $E_{\mathrm{cut}}=10\mathrm{TeV}$, $E_0=1\mathrm{TeV}$ and $\Gamma =1.8$ in Equation (146).
• 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 $\gamma \leftrightarrow 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 $\gamma \leftrightarrow a$ oscillations are present. The SED is obtained by taking in Equation (146) $E_{\mathrm{cut}}=30\mathrm{TeV}$ in the case of conventional physics, and $E_{\mathrm{cut}}=10\mathrm{TeV}$ when also $\gamma \leftrightarrow a$ oscillations are considered, while in either case we choose $E_0=1\mathrm{TeV}$ and $\Gamma =1.4$. Note that $\Gamma$ 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\ge 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_{\mathrm{cut}}=30\mathrm{TeV}$, $E_0=1\mathrm{TeV}$ and $k=1.4$ in Equation (146) for both cases (conventional physics alone, presence of $\gamma \leftrightarrow 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.

9.7. Results

Our results about the SED of the above-considered BL Lacs are exhibited in Figure 25 and Figure 28. Generally speaking, the $\gamma \leftrightarrow 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 $\gamma \leftrightarrow 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_{\gamma \to \gamma }^{\mathrm{ALP}}(E,z)$ to increase with $g_{a\gamma \gamma }$, 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, $\gamma \to a$ conversions in the BL Lac and $a\to \gamma$ back-conversions in the Milky Way help increasing $P_{\gamma \to \gamma }^{\mathrm{ALP}}(E,z)$, but not that much. Suppose instead that both $g_{a\gamma \gamma }$ 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 $\gamma \to a$ and $a\to \gamma$ 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 Figure 27 and Figure 28 around $E\simeq 3\mathrm{TeV}$. As the energy increases, photon dispersion on the CMB becomes dominant, which causes a much smaller number of $\gamma \to a$ and $a\to \gamma$ 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 Figure 27 and Figure 28 around $E=\left(10--30\right)\mathrm{TeV}$. All figures show that—as E progressively increases beyond $70\mathrm{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—${\mathbf{B}}_{\mathrm{ext}}$ orientation angles, domain lengths $L_{\mathrm{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 Figure 27 and Figure 28 around $E\sim 20\mathrm{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 $\gamma \leftrightarrow a$ oscillations naturally matches the data. For 1ES 0229+200 (see Figure 26) the model including $\gamma \leftrightarrow 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 $\gamma \leftrightarrow 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 Figure 27 and Figure 28—as the redshift increases—at high energies the difference between the results from conventional physics alone, and the model including $\gamma \leftrightarrow a$ oscillations becomes more and more dramatic. This is especially true when sizable $\gamma \to 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)\mathrm{TeV}$ as it is evident from Figure 27 and Figure 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 $\mathrm{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_{\mathrm{VHE}}$, $B_{T,{\mathcal{R}}_{\mathrm{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_{\mathrm{ext}}\ll {10}^{-9}\mathrm{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)\mathrm{TeV}$ remains unaffected at high redshift.
• The electromagnetic cascade proposed to mimic $\gamma \leftrightarrow a$ oscillation effects in blazar spectra [266] can work only for $B_{\mathrm{ext}}\lesssim \mathcal{O}\left({10}^{-15}\right)\mathrm{G}$, which is indeed quite close to the $B_{\mathrm{ext}}$ lower limits [208,209,210]. Still, for $B_{\mathrm{ext}}\gtrsim \mathcal{O}\left({10}^{-15}\right)\mathrm{G}$ the charged particles produced in the cascade are deflected by ${\mathbf{B}}_{\mathrm{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 $\gamma +\gamma \to e^{+}+e^{-}$.
• For $E\gtrsim 100\mathrm{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 $\gamma \leftrightarrow a$ oscillations reduce photon absorption, thereby considerably weakening dust absorption.

10. Polarization Effects

As stressed in the Introduction and in Section 2.3, not only the photon-ALP interaction produces $\gamma \leftrightarrow 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.

10.1. 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 $\rho \left(y\right)$ 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

$${\rho }_x=\left(\begin{array}{ccc}1& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right),{\rho }_z=\left(\begin{array}{ccc}0& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right),$$

(146)

respectively, the ALP state reading

$${\rho }_a=\left(\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 1\end{array}\right),$$

(147)

and unpolarized photons represented by

$${\rho }_{\mathrm{unpol}}=\frac{1}{2}\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 0\end{array}\right).$$

(148)

Instead, the polarization density matrix characterizing partially polarized photons shows an intermediate functional expression between Equations (147) and (149).

Once the transfer matrix of the photon-ALP system $\mathcal{U}$ is known, the final polarization density matrix $\rho$ can be computed by means of Equation (48) when the initial polarization density matrix ${\rho }_0$ is specified. It is useful to express the photonic part of $\rho$ in terms of the Stokes parameters as [269]

$${\rho }_{\gamma }=\frac{1}{2}\left(\begin{array}{cc}I+Q& U-iV\\ U+iV& I-Q\end{array}\right).$$

(149)

The photon degree of linear polarization ${\Pi }_L$ is defined as [270]

$${\Pi }_L\equiv \frac{{(Q^2+U^2)}^{1/2}}{I},$$

(150)

while the polarization angle $\chi$ reads [270]

$$\chi \equiv \frac{1}{2}\mathrm{arctg}\left(\frac{U}{Q}\right).$$

(151)

It is simple algebra to express Equations (151) and (152) in terms of the elements of the $\rho$ matrix ${\rho }_{ij}$ with $i,j=1,2$ as

$${\Pi }_L=\frac{{\left[{({\rho }_{11}-{\rho }_{22})}^2+{({\rho }_{12}+{\rho }_{21})}^2\right]}^{1/2}}{{\rho }_{11}+{\rho }_{22}},$$

(152)

and

$$\chi =\frac{1}{2}\mathrm{arctg}\left(\frac{{\rho }_{12}+{\rho }_{21}}{{\rho }_{11}-{\rho }_{22}}\right),$$

(153)

respectively.

The photon degree of linear polarization at emission ${\Pi }_{L,0}$ is linked to the photon conversion and survival probabilities $P_{\gamma \to a}^{\mathrm{ALP}}$ and $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ in the presence of the photon-ALP interaction. In particular, some theorems stated and proven in [127] show what follows.

• When photon absorption is negligible and photons (without initial ALPs) are emitted with initial degree of linear polarization ${\Pi }_{L,0}$, the two conditions $P_{\gamma \to a}^{\mathrm{ALP}}\le (1+{\Pi }_{L,0})/2$ and $P_{\gamma \to \gamma }^{\mathrm{ALP}}\ge (1-{\Pi }_{L,0})/2$ hold. If photons are emitted unpolarized (${\Pi }_{L,0}=0$) we have the two inequalities $P_{\gamma \to a}^{\mathrm{ALP}}\le 1/2$ and $P_{\gamma \to \gamma }^{\mathrm{ALP}}\ge 1/2$.
• Under the previous conditions, ${\Pi }_{L,0}$ can be viewed as the measure of the overlap between the values assumed by $P_{\gamma \to a}^{\mathrm{ALP}}$ and $P_{\gamma \to \gamma }^{\mathrm{ALP}}$. If photons are emitted unpolarized (${\Pi }_{L,0}=0$), then $P_{\gamma \to a}^{\mathrm{ALP}}$ and $P_{\gamma \to \gamma }^{\mathrm{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_{\gamma \to a}^{\mathrm{ALP}}$ and $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ possess an oscillatory behavior in energy, their values oscillate within the bounds ensured by the above theorems and ${\Pi }_{L,0}$ can be inferred. In particular, from an observed astrophysical spectrum ${\Phi }_{\mathrm{obs}}$ we can extract $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ as

$$P_{\gamma \to \gamma }^{\mathrm{ALP}}=\frac{{\Phi }_{\mathrm{obs}}}{{\Phi }_{\mathrm{em}}},$$

(154)

where ${\Phi }_{\mathrm{em}}$ is the emitted spectrum, which is supposed either known or derivable from ${\Phi }_{\mathrm{obs}}$ (for more details see [127]). Moreover, we have $P_{\gamma \to a}^{\mathrm{ALP}}=1-P_{\gamma \to \gamma }^{\mathrm{ALP}}$ since there is no photon absorption. Now, ${\Pi }_{L,0}$ is simply given by the measure of the interval where $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ and $P_{\gamma \to a}^{\mathrm{ALP}}$ overlap, as Figure 29 shows for a typical shape of $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ and $P_{\gamma \to a}^{\mathrm{ALP}}$.

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 ${\Pi }_L$. In addition, only flux observations are needed.

10.2. 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 ${\Pi }_{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 ${\Pi }_L$ and $\chi$. Concerning the parameters of the photon-ALP system, in order to be specific we take $g_{a\gamma \gamma }=0.5·{10}^{-11}{\mathrm{GeV}}^{-1}$ and $m_a={10}^{-10}\mathrm{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 ${\mathbf{B}}^{\mathrm{clu}}$ and the electron number density profile $n_e^{\mathrm{clu}}$ 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_0^{\mathrm{clu}}=15\mu \mathrm{G}$, $k_L=0.2{\mathrm{kpc}}^{-1}$, $k_H=3{\mathrm{kpc}}^{-1}$, $n_{e,0}^{\mathrm{clu}}=0.5·{10}^{-2}{\mathrm{cm}}^{-3}$, ${\eta }_{\mathrm{clu}}=0.75$, ${\beta }_{\mathrm{clu}}=2/3$, $r_{\mathrm{core}}=100\mathrm{kpc}$, and a cluster radius of $1\mathrm{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 ${\mathcal{U}}_{\mathrm{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^{\mathrm{ext}}=1\mathrm{nG}$ with coherence length $L_{\mathrm{dom}}^{\mathrm{ext}}$ in the range $(0.2–10)\mathrm{Mpc}$ and average $\langle L_{\mathrm{dom}}^{\mathrm{ext}}\rangle =2\mathrm{Mpc}$, we compute the transfer matrix ${\mathcal{U}}_{\mathrm{ext}}$ in this region. We recall that $\gamma \leftrightarrow 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 $\gamma \leftrightarrow a$ oscillations. Hence, we are in a position to evaluate the transfer matrix of the photon-ALP system in the Milky Way ${\mathcal{U}}_{\mathrm{MW}}$.

By combining the previous transfer matrices in the correct order, we obtain the whole transfer matrix ${\mathcal{U}}_{\mathrm{tot}}$ associated with the propagation of the photon/ALP beam from the cluster core to us

$${\mathcal{U}}_{\mathrm{tot}}={\mathcal{U}}_{\mathrm{MW}}{\mathcal{U}}_{\mathrm{ext}}{\mathcal{U}}_{\mathrm{clu}},$$

(155)

while the final photon survival probability with $\gamma \leftrightarrow a$ oscillations $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ is given by

$$P_{\gamma \to \gamma }^{\mathrm{ALP}}=\sum _{i=x,z}\mathrm{Tr}\left[{\rho }_i{\mathcal{U}}_{\mathrm{tot}}{\rho }_{\mathrm{in}}{\mathcal{U}}_{\mathrm{tot}}^{†}\right],$$

(156)

with ${\rho }_{\mathrm{in}}\equiv {\rho }_{\mathrm{unpol}}$ reading from Equation (149) and ${\rho }_x$ and ${\rho }_z$ from Equation (147). Then, the corresponding final degree of linear polarization ${\Pi }_L$ emerges from Equation (153).

In the left panel of Figure 30 we report $P_{\gamma \to \gamma }^{\mathrm{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 ${\Pi }_L$, while in the right panel of Figure 30 we present the probability density function $f_{\mathcal{P}}$ of ${\Pi }_L$ associated with several realizations at the benchmark energy $E_0=10\mathrm{MeV}$.

From Figure 30 we observe that the photon/ALP beam propagates in the weak mixing regime, $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ and ${\Pi }_L$ show an oscillatory energy behavior and ${\Pi }_L>0$. Since ${\Pi }_L=0$ is never the most probable result, as the right panel of Figure 30 shows, we can conclude that a signal of ${\Pi }_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_{\gamma \to \gamma }^{\mathrm{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 ${\Pi }_L$ up to so high energies, yet [275].

10.3. 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 ${\Pi }_{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 ${\mathcal{U}}_{\mathrm{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 ${\mathcal{U}}_{\mathrm{host}}$. The transfer matrices of the photon-ALP system in the galaxy cluster ${\mathcal{U}}_{\mathrm{clu}}$, in the extragalactic space ${\mathcal{U}}_{\mathrm{ext}}$ and in the Milky Way ${\mathcal{U}}_{\mathrm{MW}}$ exactly read from the previous section with the same parameters apart from $n_{e,0}^{\mathrm{clu}}=5·{10}^{-2}{\mathrm{cm}}^{-3}$ since we are considering now a CC galaxy cluster. By evaluating the latter three transfer matrices and combining them with ${\mathcal{U}}_{\mathrm{jet}}$ and ${\mathcal{U}}_{\mathrm{host}}$ in the correct order we can calculate the total transfer matrix ${\mathcal{U}}_{\mathrm{tot}}$ associated to the propagation of the photon/ALP beam from the jet base up to us

$${\mathcal{U}}_{\mathrm{tot}}={\mathcal{U}}_{\mathrm{MW}}{\mathcal{U}}_{\mathrm{ext}}{\mathcal{U}}_{\mathrm{clu}}{\mathcal{U}}_{\mathrm{host}}{\mathcal{U}}_{\mathrm{jet}},$$

(157)

while the final photon survival probability in the presence of photon-ALP interaction $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ reads from Equation (157) and the corresponding final degree of linear polarization ${\Pi }_L$ is given by Equation (153).

The left panel of Figure 31 shows $P_{\gamma \to \gamma }^{\mathrm{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 ${\Pi }_L$, while the right panel of Figure 31 reports the probability density function $f_{\Pi }$ of ${\Pi }_L$ associated to several realizations at the benchmark energy $E_0=10\mathrm{MeV}$.

Figure 31 shows that the photon/ALP beam is in the weak mixing regime and $P_{\gamma \to \gamma }^{\mathrm{ALP}}$ presents a pseudo-oscillatory behavior with respect to the energy. Correspondingly, we observe in the central panel of Figure 31 that ${\Pi }_L>0$ with a high variation. From the right panel of Figure 31, we infer that ${\Pi }_L=0$ is never the most probable result. Therefore, we can conclude that a signal of ${\Pi }_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_{\gamma \to \gamma }^{\mathrm{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.

11. 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=\mathcal{O}\left({10}^{-10}\right)\mathrm{eV}$ and $g_{a\gamma \gamma }=\mathcal{O}\left({10}^{-11}\right){\mathrm{GeV}}^{-1}$ have emerged.

An indirect detection can be made by the new generation of gamma-ray observatories, such as CTA [130], HAWC [131], GAMMA-400 [132], LHAASO [133], TAIGA-HiSCORE [134] and HERD [135].

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!