Search for the imprint of axion-like particles in the highest-energy photons of hard $\gamma$-ray blazars

Axion-like particles (ALPs), predicted in theories beyond the Standard Model, can have observational effects on the transparency of the Universe to $\gamma$ rays in the presence of magnetic fields. In this work, we search for effects compatible with the existence of ALPs with 80 months of data from the ${\it Fermi}$ Large Area Telescope, by comparing the distributions of observed highest energy photons from sources beyond redshifts of z $\geq$ 0.1 with theoretical predictions in the presence of ALPs. We find no evidence for an increased $\gamma$-ray transparency due to ALPs and therefore we set limits on the ALPs parameters assuming a value of the intergalactic magnetic field strength of 1 nG. Photon-ALP couplings above $10^{-11}$ GeV$^{-1}$ are excluded for ALP masses $m_{a}$ $\lesssim3.0$ neV. These limits exclude a region of the parameter space not covered by other $\gamma$-ray telescopes and are compatible with constraints imposed by other experiments.


Introduction
Axion-like particles (ALPs), very light pseudo-scalar bosons predicted by multiple extensions of the Standard Model [1][2][3][4][5][6], could be detected through their coupling to photons in the presence of external magnetic fields, L aγ = − 1 4 g aγ F µνF µν a = g aγ E · Ba, (1.1) latter plays a key role in the observation of the γ-ray sky, being the fundamental source of opacity for the Universe to γ rays. It is also relevant for re-ionization models in cosmology and galaxy formation and evolution [14,[34][35][36]. Direct measurements of the EBL are challenging due to the presence of other backgrounds in the solar system and our galaxy, such as the bright zodiacal light and the galactic emission [37]. The main constraints come from the integrated light of discrete extragalactic sources and from γ-ray observations [38][39][40][41], for which blazars are good probes [42]. Many different approaches have been used to model the intensity and spectral shape of the EBL. In this work, we use the observationally-based Domínguez et al. model [43] to derive all the results and the Finke et al. [44] model for comparison.
The survival probability, or attenuation factor, is described by a decreasing exponential law, P γγ = exp [−τ (E, z)]. It depends on the optical depth parameter τ (E, z), which is an increasing function of the photon energy and the distance to the source. The γ-ray spectra of sources are then described by, where φ obs and φ int are the observed and intrinsc spectra, respectively. The latter would be the observed spectrum of the source if there was no EBL absorption. The cosmic γ-ray horizon (CGRH) [45,46] is defined as the energy E 0 at which the optical depth becomes unity τ (E 0 , z) = 1. The solid line in Fig. 13 of Ref. [47] is the CGRH derived with the Domínguez et al. model. Above the CGRH curve, the Universe is more opaque to γ rays due to larger values of τ that translate into smaller photon survival probabilities, whereas in the region below the curve, the survival probabilities are larger and the Universe is more transparent to γ rays. For a given redshift, the CGRH quantifies the maximum energy of photons that survive the EBL. If there were modifications of the canonical γ-ray propagation, the observed HEP event for each source should change correspondingly. This is what we use in order to search for ALPs effects.

Photon-ALP mixing
The photon-ALP system is described by the lagrangian of Eq. 1.1. Given a homogeneous magnetic field B over a distance of length s and a polarized photon beam, the photon-axion conversion probability P γ→a can be calculated analytically [48], where B T is the transverse component to the direction of propagation, g γa is the coupling constant, E is the photon energy and E c is the critical energy around which we expect the ALPs effects to occur, defined in Eq. 1.2. Equation 2.2 is only valid if the photon-ALP beam is in a pure polarization state and the conversion takes place in a homogeneous magnetic field. Since photon polarization is not measured in the γ-ray band, we have to treat the photon-ALP beam as unpolarized and take into account general magnetic field morphologies.
The IGMF is often modeled as a domain-like structure in redshift, with a homogeneous magnetic field strength in each cell and a random orientation in different domains [26]. The size of each cell is the so-called coherence length, a distance upon which the magnetic field is homogeneous. Unfortunately, there are only upper limits available for the strength and coherence length of large-scale magnetic fields, with B ≤ 1 nG and s ∼ Mpc [17,49]. These upper limits are used as model parameters for the analysis performed in this work, B = 1 nG and s = 1 Mpc. We take n e ∼ 10 −7 cm −3 as the electron density in the intergalactic medium [50]. With this value of n e , the plasma frequency is ω pl ∼ 1.17 · 10 −14 eV.
The propagation of the beam in a domain-like structure is a stochastic process. For this reason, the effects of a single trajectory of the beam cannot be measured. An average over a large number of realizations of the process is required in order to evaluate the photon survival probability. We compute the photon-ALP oscillation probability in the transfermatrix formalism as described in, e.g., Refs. [51][52][53], using the gammaALPs code 2 . A summary is provided in Appendix A. The systematic uncertainties associated with the magnetic field parameters B T and s have the biggest impact on the results, changing the exclusion region size up to 35%, as discussed in Section 5 and in Appendix B. Additionally, instead of using a single magnetic field realization for all the sources, we use different realizations for the different lines of sight.

Data and simulations
We use the energy of the HEP events from the Second Catalog of Hard 3 Fermi-LAT sources (2FHL). The catalog reports the properties of 360 sources significantly detected by the LAT above 50 GeV from August 2008 to April 2015 [47]. The characterization of such sources at these energies was made possible with the Pass 8 event-level analysis [54] and long telescope exposure. More than 80% of the sources are extragalactic, of which 75% are AGN.
We chose this set of sources because the critical energy of the ALPs mixing in the IGMF lies within the energy range of the 2FHL catalog, for the IGMF field parameters, couplings g 11 and masses m a,neV tested in this work. Out of these available sources we only use source with redshifts z ≥ 0.1. For smaller values of redshift the effects of ALPs are too small to be detected. The catalog also offers the probabilities of each photon event belonging to each one of the sources in the region of interest. We only take sources for which the probability of the HEP assigned to them is P ≥ 0.99. This cut improves the background rejection but it also entails a reduction in statistics, from a sample of 96 sources to 79 sources. With this cut, results are more conservative and the excluded region decreases by 15%, see Appendix B for more details.
An overview of the method goes as follows. We require two main elements, namely, the HEP events from the 2FHL catalog sources and their corresponding probability distribution functions (p.d.f.s). The former are taken directly from the catalog results, whereas the latter are derived from simulations. Then, we compute the probability or likelihood that each HEP has with its corresponding distribution, and we combine them because the sources are independent.
For each source, we simulate a HEP probability distribution from which we expect the measured HEP to come. We use 40 logarithmically spaced energy bins from 50 GeV to 2 TeV, the same range as in the catalog, and we compute the expected number of events between energies E 1 and E 2 through (3.1) The first term in the equation, P γγ (E, z, θ), is the photon survival probability from Eq. 2.1. The parameter θ = (m a , g aγ ) represents the mass of the axion and the coupling constant for the ALPs model of Section 2.2. The second term, φ(E), is the intrinsic spectrum of the source, which has to be derived in a region with negligible EBL attenuation, below the range of the catalog. Therefore, we make use of a recent re-analysis of the 2FHL sources, carried out in Ref. [55], which extends the analysis spectral range down to 300 MeV. Following a χ 2 minimization procedure, we fit these spectral data to power laws and logarithmic parabolas. This is possible, as γ-ray sources typically have smooth spectra over a limited energy range and blazars are well described by these functions [35,41]. The fits take spectral points until τ (z, E) ∼ 0.1, energies in which EBL effects cease to be negligible. A list of the sources and the spectral fits can be found in Appendix C. An example of the spectral fit and absorption models on the spectrum of a source is displayed in Fig. 1. The last term in Eq. 3.1, (E), is the exposure map of the Fermi-LAT, an integral of the total instrument response function over the entire region of interest of each source, provided in the analysis files of Ref. [47]. The probability of detecting c counts in the i-th bin is given by Poisson statistics, i.e. p ij = N c ij exp(−N ij )/c!, where the index j stands for the source and N ij is given by Eq. 3.1. Using these probability functions, we generate events for each source and energy bin. The last non-empty energy bin, c ≥ 1, is taken as the bin with the HEP. For each source j and attenuation model θ, we build a histogram of HEPs by repeating these pseudo-experiments 10 4 times. These normalized histograms are the HEP p.d.f.s. An example is displayed in Fig. 2.

Likelihood analysis
For a random variable x distributed according to a p.d.f. f (x, θ), where θ is any parameter of the function, the probability for a measurement x i to be in [ Assuming N independent observations of x, the joint likelihood function is [56], In our work, the random variable is the HEP of each source, E j , which is the maximum energy event with a probability larger than 99% to be associated with the j-th source. Since all the sources are independent, the joint likelihood is the product of likelihoods for each individual source. Each likelihood is computed using the p.d.f.s. simulated following the procedure of Section 3. For the null hypothesis, with only EBL, the parameter θ is set to θ 0 = (m a , g aγ ) = (0, 0). For the alternative hypothesis θ 1 = (m a , g aγ ) takes values in the ALPs parameter space for which the critical energy remains within the Fermi-LAT energy range. The joint likelihood function is given by, With the joint likelihood functions we can define different test-statistics (TS) and perform a statistical hypothesis testing between models. We use the following TS, Λ, defined as the log-likelihood ratio test where θ 1 now stands for a composite ALPs domain consisting of a set of points within the main region, taken between 0.1 ≤ m a ≤ 10 neV and 0.5 ≤ g 11 ≤ 7.0. In order to draw any conclusions, we need to compare the observed value of TS, Λ obs , to different acceptance or rejection thresholds integrated from the null and alternative TS distributions. The null Λ distribution, f Λ (H 0 ), is derived by generating 10 4 Monte-Carlo events under the HEP probability distributions simulated with the null hypothesis. We compute the detection threshold Λ thr by integrating 95% of this distribution. If Λ obs < Λ thr , observations are compatible with the null hypothesis. On the contrary, if Λ obs > Λ thr , the TS is too large to be compatible with the null hypothesis and we could claim a 2σ significance signal discovery. The alternative TS distributions, f Λ (H 1 (m a , g aγ )), are derived in the same way, but the events are now simulated under the alternative hypotheses. The exclusion thresholds for the ALPs hypotheses, Λ exc (m a , g aγ ), are computed by integrating 5% of these distributions. If Λ obs < Λ exc (m a , g aγ ), the set of ALPs parameters would be rejected. This sub-domain of parameters is chosen a posteriori, testing different sets of points until the medians of the alternative distributions are larger than the rejection thresholds for the null hypothesis.

Results and conclusions
The null and alternative TS distributions, f Λ (H 0 ) and f Λ (H 1 (m a , g aγ )) for m a = 1.3 neV and g 11 = 5.2 respectively, are shown in Fig. 3. The observed TS found with the HEP data is shown as well, Λ obs = −4.7. The 2σ detection threshold, derived by integrating the null distribution, is Λ thr = 3.2. As can be seen, Λ obs < Λ thr , therefore the results are compatible with the null hypothesis and no evidence for ALPs was found in these data. An upper limit is set by computing the 95% exclusion thresholds Λ exc (m a , g aγ ), also displayed in Fig. 3, and testing whether Λ obs < Λ exc (m a , g aγ ). This is repeated for each point tested in the sub-domain. The resulting upper limits can be seen in Figs. 4 and 5. For B T = 1 nG and s = 1 Mpc, photon-ALP couplings between 1.0 g 11 7.0 are excluded for masses below m a 3.0 neV. On the right side, the contour follows the constant critical energy diagonal from Eq. 1.2, whereas the horizontal line around g 11 ∼ 1 depends upon the product B T · s. This horizontal line extends to arbitrarily small masses since B T · s does not depend on m a . When m 2 a < ω 2 pl,neV , the effective mass takes the value of the plasma frequency of the medium. The fluctuations in the contours are due to the limited number of magnetic field realizations and pseudo experiments in the simulation.
The dominant systematic uncertainties are related to the choice of magnetic field parameters. Only upper limits exist for the strength and coherence length of the IGMF [49]. We repeat the simulation and likelihood analysis by decreasing the field strength to B = 0.5 nG and B = 0.1 nG. The former value yields upper limits that are smaller in area of the excluded region by ∼ 30% compared to the initial case, while for the latter we cannot set any upper limits. We also increase the coherence length to s = 5 Mpc. In this scenario, the resulting upper limits increase by roughly ∼ 30%. These results, seen in Fig. 4, are consistent with the ALPs mixing equations. Other sources of systematic uncertainties, such as the EBL model, energy dispersion effects, a different set of magnetic field realizations and HEP probability cuts are smaller than ∼ 15%and are discussed in Appendix B. The limits derived in this work are compatible with other limits and sensitivities of future ALPs experiments, as shown in Fig. 5. Together with the SN 1987A γ burst experiment [57] and with previous Fermi-LAT [21] results, our limits strongly constrain part of the parameter space in which ALPs can modify the opacity of the Universe to γ rays. Our limits also constrain a part of the unexplored parameter space that the Fermi-LAT work using observations of the central AGN in the Perseus galaxy cluster could not cover before, i.e., the hole around g 11 ∼ 3 and m a ∼ 3 neV. They are also within the planned sensitivities of ALPS II [58] and IAXO [59]. None of these limits constrain the region where ALPs could compose the entirety of dark matter content of the Universe [60] which is below the dashed black line in Fig. 5.
Magnetic field morphologies in the interstellar and intergalactic space are not fully understood yet. Better observations of cosmic magnetic fields are needed in order to reduce the systematic uncertanties associated to these fields, which are crucial for the photon-ALP beam propagation. Future experiments like JVLA [61], ALMA [62], and SKA [63] will be able to improve these current limitations [20]. Recent EBL results, such as Ref. [64], may also be used in future ALPs analyses.
A part of the parameter space where ALPs could affect the transparency of the Universe to γ rays remains to be explored. For larger values of m a , the critical energy increases and the maximal conversion probability takes place in the TeV range. The conversion is further enhanced by the inclusion of the galactic magnetic field in the propagation of the beam. Current Cherenkov telescopes, with energy ranges within the TeV range can reach higher masses in the parameter space and improve the limits derived in this work. The future Cherenkov Telescope Array [65] can probe even higher masses. Hence, the combined likelihood analysis of many sources presented in this work can be extended to these instruments. The number of events within a given energy range can be measured and compared to simulations that include ALPs models. Other magnetic field scenarios may also be included to improve the results.

Acknowledgments
The Fermi LAT Collaboration acknowledges generous ongoing support from a number of agencies and institutes that have supported both the development and the operation of the LAT as well as scientific data analysis. These include the National Aeronautics and Space Administration and the Department of Energy in the United States, the Commissariat à l'Energie Atomique and the Centre National de la Recherche Scientifique / Institut National

A Survival probability
In this Appendix we provide a brief summary of the survival probability derivation, which is treated in detail in Ref. [51]. From Eq. 1.1, the photon-ALP propagation for a monochromatic system along the y-axis, in case that the photon energy is E m a , can be described by [66] (i∂ y + E + M 0 ) where E is the photon energy, M 0 is the mixing matrix and the second factor is ψ(y), the beam state vector. The photon polarization amplitudes along the x-and z-axis are denoted by A x (y) and A z (y), respectively, while the ALP field amplitude is a(y). The mixing matrix is real and symmetric, and it involves different terms. The general form is given by: The ∆ aγ -terms and the ∆ aa -term represent the mixing of photons with ALPs from Eq. 1.1 and ALPs self-interactions, respectively. The remaining terms depend on the properties of the medium, namely QED vacuum effects and absorption mechanisms. The former come from the Heisenberg-Euler-Weisskopf (HEW) effective Lagrangian for the photon one-loop vacuum polarization under an external magnetic field [67]. Vacuum QED terms can be ignored at high energies. The other contribution is due to background particles in the medium that may annihilate the primary photon, such as EBL.
The photon-ALP beam has to be treated as unpolarized, therefore it is described by a polarization density matrix, ρ(y) = ψ(y) ⊗ ψ(y) † , which obeys the Von Neumann equation of non-relativistic quantum mechanics [51], The solution of the propagation of the beam is given by the transfer matrix, The transition probability from one state to another state is given by the trace of the projection of both states, For the intergalactic medium, the magnetic region is split into N domains with a homogeneous field in each one. For a homogeneous field, the problem is simplified because we can always align B T with the z-axis. In this case, this can no longer be done across all the domains, but the problem can still be simplified through similarity transformations, where V (φ) is the rotation matrix in the x−z plane perpendicular to the propagation direction and φ is the angle B T forms with the z-axis in each domain. The full transfer matrix of the system across N domains is, The transition probability between two states is computed with Eqs. A.6 and A.8. The propagation is a stochastic process due to the random orientations of the field in each domain, therefore only the mean properties of the beam can be evaluated. Moreover, we have to sum over the two final polarization states because of the lack of polarization measurements. The photon survival probability is, The modified photon survival probabilities for each grid point are computed with this equation.

B Other uncertainties
In order to evaluate uncertainties associated to model parameters and analysis choices, we derive different sets of limits repeating the simulation and analysis procedures of Sections 3 and 4. For each test, we extract the percentage changed in area from one case to another for the particular grid explored in this study, between 0.1 ≤ m a ≤ 10 neV and 0.5 ≤ g 11 ≤ 7.0. We refer as area the visual region in logarithmic scale that is excluded for each case.  In Section 3, we discussed the AGN data sample and took sources based on the HEP probability to belong to a source. All of the sources in the 2FHL catalog have a HEP with P ≥ 0.85, with most of them above P ≥ 0.99, due to the low background of the LAT at high energies. Selecting sources with higher values of P allows us to reduce the events that come from background. However, this also entails a reduction of statistics in our sample. We tested the effects of different HEP probability cuts within one realization, resulting in contours with area changes smaller than ∼ 10%, as displayed in Fig. 7. Our results were derived with the Domínguez et al. EBL model. We test the effects of choosing a different model by repeating the analysis with the Finke et al. model [44]. The upper limits increase by ∼ 15%, as seen in Fig. 8. Finally, we did not consider energy dispersion effects. The reason for this is that, above 1 GeV, these effects are below 10% at 68% confidence and thus we do not expect a large change in the results. All the results were computed with the mean values of the spectral fit parameters, which have their own statistical uncertainties that we did not take into account. A way to do this would be to bias the fluxes of the sources to produce more photons at high energies, then the probability distributions would shift to the right, changing the likelihood values for both models and all sources. We did not perform this for two reasons. First, this is similar to modifying the opacity of the Universe with a different propagation model, which we already carried out with the Finke et al. model, resulting in a limits difference that is still smaller than the dominant uncertainty imposed by the magnetic field. Second, the parameters of the spectra would take random values for each source that might cancel out in the combined likelihood analysis of all the sources. Additionally, some of these spectra accept different model representations and are not bound to only power laws or only logarithmic parabolas, hence we tested the effects of choosing different sets of parameters for the sources. The results were still smaller than 10% compared to the main limits we obtained in this work. Finally, the effective areas have also associated systematic uncertainties at high energies. It is usually evaluated through a bracketing procedure during the LAT data analysis. Its effects at the energies considered in our analysis are estimated to yield differences in photon counts between 5% and 15%. This is still below the dominant systematic uncertainty from the magnetic field. Table 1 summarizes the spectral fits for all the sources in the 2FHL catalog with z ≥ 0.1. The first column is the spectral shape, which can be a power law (pow) or a logarithmic parabola (log). The second column is the χ 2 value from the minimization procedure. The third column, ndf, is the number of degrees of freedom. The last column denotes the p-value. Table 2 summarizes the fit parameters for all the sources with a pivot energy of 100 GeV. The column indices follow the logarithmic parabola equation parameters and its corresponding errors: