# Probing $\mu$eV ALPs with future LHAASO observation of AGN $\gamma$-ray
spectra

Guangbo Long  Siyu Chen  Shuo Xu  Hong-Hao Zhang<EMAIL_ADDRESS>School of Physics, Sun Yat-sen University, Guangzhou, GuangDong, People’s
Republic of China

###### Abstract

Axion-like particles (ALPs) are predicted in some well-motivated theories
beyond the Standard Model. The TeV gamma-rays from active galactic nuclei
(AGN) suffers attenuation by the pair production interactions with the cosmic
infrared background light (EBL/CMB) during its travel to the earth. The
attenuation can be circumvented through photon-ALP conversions in the AGN and
Galaxy magnetic-field, and a flux enhancement is expected to arise in the
observed spectrum. In this work, we study the potential of the AGN gamma-ray
spectrum for energy up to above 100 TeV to probe ALP-parameter space at around
$\mu$eV, where the coupling $g_{a\gamma}$ is so far relatively weak
constrained. We find the nearby and bright sources, Mrk 501, IC 310 and M 87,
are suitable for our objective. Assuming an intrinsic spectrum exponential
cutoff energy at $E_{\rm c}$=100 TeV, we extrapolate the observed spectra of
these sources up to above 100 TeV by the models with/without ALPs. For
$g_{a\gamma}\gtrsim 2\times$$10^{-11}\rm GeV^{-1}$ with $m_{a}\lesssim
1\,\mu$eV, the flux at around 100 TeV predicted by the ALP model can be
enhanced more than an order of magnitude than that from the standard
absorption, and could be detected by LHAASO. Our result is subject to the
uncertainty from the intrinsic spectrum above tens of TeV, which require
further observations on these sources by the forthcoming CTA, LHAASO, SWGO and
so on.

##  I Introduction

Several extensions of the Standard Model suggest the existence of very light
pseudoscalar bosons called axion-like particles (ALPs) Svrcek2006 ;
Jaeckel2010 . These spin-0 neutral particles are a sort of generalization of
the axion, which was originally proposed to solve the strong CP problem
naturally PQ1977 ; Weinberg1978 ; Wiczek1978 , and they are also a promising
dark-matter candidate Preskill1983 ; Abbott1983 ; Dine1983 ; Marsh2011 . One
of their characteristics is their coupling to photons by
$g_{a\gamma}\mathbf{E}\cdot\mathbf{B}a$, with $g_{a\gamma}$ being the coupling
strength, $\mathbf{E}$ the electric field of the photons, $\mathbf{B}$ an
external magnetic field, and $a$ the ALP field strength Raffelt1988 ;
Sikivie1983 . As a consequence, the phenomenon of photon-ALP mixing take
place, and lead to photon-ALP oscillations (or conversions) Raffelt1988 ;
Sikivie1983 ; Hochmuth2007 ; Angelis2008 . To reach efficient conversions,
they should take place above a critical energy given by Hooper2007 ;
Angelis2007 ; Mirizzi2009

$E_{\rm crit}\sim 20(\frac{m_{a}}{10^{-6}\rm eV})^{2}(\frac{10^{-5}\rm
G}{B})(\frac{2.5\cdot 10^{-11}{\rm GeV}^{-1}}{g_{a\gamma}})\,\rm TeV,$ (1)

where $B$ is the homogeneous magnetic-field component transverse to the
propagation direction. In contrast to axion, ALP mass $m_{a}$ is independent
of $g_{a\gamma}$. Around the critical energy, oscillatory features that depend
on the configuration of magnetic field are expected to occur Wouters2012 .

Many laboratory experiments and astronomical observations are being carried
out to search for ALPs via the effect mentioned above. The representative
experiments are photon-regenerated experiments, such as“Light shining through
a wall” experiments ALPS ALPS2013 , solar ALPs experiments CAST CAST2011  and
dark-matter haloscopes ADMX ADMX2006 .

Owing to the universal presence of magnetic fields along the line of sight to
active galactic nuclei (AGN), photon-ALP oscillations can lead to distinctive
signatures in AGN spectra Angelis2007 ; Simet2008 ; Hooper2007 ; Angelis2008 .
Thus, ALP-photon coupling can be detected through the observations of AGN
(see, e.g. Refs Angelis2007 ; Conde2009 ; Angelis2011 ; Dominguez2011 ;
Tavecchio2012 ).

On the one hand, it is possible to detect the ALP-induced observational
effects on the $\gamma$-rays transparency of the Universe Angelis2007 ;
Simet2008 ; Conde2009 ; Angelis2011 ; Dominguez2011 ; Meyer2013 ; Meyer2014 ;
Troitsky2016 ; Montanino2017 ; Buehler2020 . The very-high-energy (VHE, above
100 GeV) $\gamma$-rays from the extragalactic sources suffer attenuation by
pair-production interactions with the background (extragalactic background
light, EBL; or cosmic microwave background, CMB) photons during the
propagation Nikishov1962 ; Hauser2001 ; HESS2006 ; Dwek2013 ; Costamante2013 .
The attenuation increases with the distance to the source and the energy of
the VHE photons Dwek2013 . If the photon-ALP conversions exist for a
sufficiently large coupling, the emitted photons convert into ALPs and then
these ALPs reconvert back into photons before arriving in the Earth, i.e. ALPs
circumvent pair production. Thus, the opacity of the Universe for VHE gamma-
rays is reduced and the observed flux is enhanced significantly (i.e. causing
a hardening of the spectra above $E_{\rm crit}$, see e.g. Refs. Angelis2007 ;
Mirizzi2009 ; Dominguez2011 ; Angelis2013 ; Meyer2013 ; Troitsky2016 ;
hardening2 ; Galanti2018 ; Galanti2019 ). The range of the parameters where
ALPs would increase the $\gamma$-ray transparency of the Universe (for 1.3
times the optical depth of Franceschini _et al._ EBL model Franceschini2008 )
is constrained from VHE $\gamma-$rays observations of blazar (AGN with jet
closely aligned to the line of sight) Meyer2013 . Data from the Fermi-LAT
observations of distant (redshift $z>$0.1) blazar limit
$g_{a\gamma}<10^{-11}\rm GeV^{-1}$ for $m_{a}<$3 neV assuming a value of the
intergalactic magnetic field strength of 1 nG Buehler2020 .

On the other hand, taking seriously the irregularities of AGN gamma-ray
spectra produced by the oscillations at energies around $E_{\rm crit}$, strong
bounds on $g_{a\gamma}$ have been derived CTA2020alp . In particular, for 0.4
neV$<m_{a}<$100 neV, the strongest bounds on $g_{a\gamma}$ are derived from
the absence of irregularities in H.E.S.S. and Fermi-LAT observations as well
as Fermi-LAT observations of AGN Hess2013alp ; Fermi2016alp ; Zhang2018alp ;
Li2020alp . It is worth emphasizing that this method highly depend on the
configuration of magnetic fields adopted Libanov2020 .

So far, the coupling $g_{a\gamma}<6.6\times 10^{11}\,\rm GeV^{-1}$ for 0.2
$\mu$eV$\lesssim m_{a}\lesssim 2\,\mu$eV, containing viable ALP dark matter
parameter space (i.e. $g_{a\gamma}\lesssim 2\times 10^{11}\,\rm GeV^{-1}$ for
$m_{a}\sim\mu$eV) Arias2012 , almost have not been limited (see e.g. Figure 5
of Ref Buehler2020 ), although they are expected to be probed by future
experiments (e.g. ALPS II ALPS2013 , IAXO IAXO2019 ) or radio-astronomical
observations (e.g. Refs Sigl2017 ; Edwards2020 ; Caputo2019 ; Ghosh2020 ).
According to Eq. (1), effective oscillations for ALPs at this parameter range
takes place when the energy of photons is larger than $\sim$20 TeV for the
value of AGN magnetic field $B\sim 10^{-5}$ G Kohri2017 ; Meyer2014magnetic .
Therefore, photons with energy larger than 20 TeV are required to be detected
if probing these ALPs through the VHE $\gamma-$rays observations of AGN.
However, only very few photons of such high energy from extragalactic sources
have been observed by past and present telescopes TEVCAT .

Thanks to the upcoming Large High Altitude Air Shower Observatory (LHAASO
LHAASO2019 ) with the ability to survey the TeV sky continuously, it is
expected to reach sensitivities above 30 TeV about 100 times higher than that
of the current VHE instruments (e.g. H.E.S.S. HESSp , MAGIC MAGICp , VERITAS
VERITASp ) LHAASO2019 . Furthermore, the conversions in the intergalactic
magnetic field (IGMF) can be neglected. With current upper limits on the IGMF
strength of $\lesssim 10^{-9}$ G and on $g_{a\gamma}<6.6\times 10^{-11}\,\rm
GeV^{-1}$ IGMF ; CAST2011 , Eq.(1) give that $E_{\rm crit}\lesssim$100 TeV
only for $m_{a}\lesssim 35$ neV. Obviously, the unprecedented sensitivities of
LHAASO to detect TeV AGN provide a good chance to probe $\mu$ eV ALPs.

In this paper, we assume ALPs converted from the gamma-rays photons in the
AGN’s magnetic field travel unhindered through extragalactic space, and then
these ALPs partially back-convert into photons in galactic magnetic field
(GMF), see Fig. 1. We investigate the LHAASO sensitivity to detect the ALP-
induced flux enhancement at the highest energies by using the extrapolated
observations of suitable AGNs for energy up to above 100 TeV, and estimate the
corresponding ALP-parameter space.

The paper is structured as follows. In section II we give the formula for
evaluating the photon survival probability along the line of sight. The sample
selection is described and the data analysis is introduced in section III
before presenting our results in section IV. We discuss our model assumptions
in section V and conclude in section VI.

##  II Photon survival probability

Figure 1: Cartoon of the formalism adopted in this article, where the TeV
photon ($\gamma$, purple line) /ALP (a, black line) beam propagates from the
AGN $\gamma$-ray source to the Earth. The interaction $\gamma+\gamma_{\rm
EBL/CMB}\rightarrow\rm e^{\pm}$ takes place during the photon propagation. The
photon-ALP conversions (green line) take place in the magnetic field around
the gamma-ray emission region and GMF respectively, leading to an improvement
of photon survival probability. There are two $\gamma\rightarrow\gamma$
channels: $\gamma\rightarrow\gamma(\rm e^{\pm})\rightarrow\gamma$;
$\gamma\rightarrow a\rightarrow\gamma$, and the latter is dominant in this
situation.

During the propagation from $\gamma$-ray emission regions to the Earth, we
assume the emitted photons mix with ALPs in the AGN $B$-field and in the GMF
respectively, and undergo the pair production with EBL/CMB in extragalactic
space, shown in Fig. 1. The photon survival probability
$P_{\gamma\rightarrow\gamma}$ will be calculated under these effects.

###  II.1 Photon-ALP conversion

The probability of the conversion from an unpolarized photon to an ALP after
passing through a homogeneous magnetic field $\mathbf{B}$ over a distance of
length $r$, is expressed as Tavecchio2012 ; Angelis2013 ; Masaki2017

$P_{\gamma\rightarrow
a}=\frac{1}{2}\Big{[}\frac{g_{a\gamma}B}{\bigtriangleup_{\rm
osc}(E)}\Big{]}^{2}{\rm sin}^{2}\Big{[}\frac{\bigtriangleup_{\rm
osc}(E)r}{2}\Big{]}$ (2)

with $\bigtriangleup_{\rm osc}(E)$=$g_{a\gamma}B\sqrt{1+(\frac{E_{\rm
crit}}{E}+\frac{E}{E_{\rm H}})^{2}}$, where $E_{\rm crit}$ is the critical
energy defined in Eq. 1 and $E_{\rm H}$ is derived from the Cotton-Mouton (CM)
effect accounting for the photon one-loop vacuum polarization given by
Tavecchio2012

$E_{\rm H}=2.1\times 10^{5}\big{(}\frac{10^{-5}\rm
G}{B}\big{)}\big{(}\frac{g_{a\gamma}}{10^{-11}\rm GeV^{-1}}\big{)}\,\,\rm
GeV.$ (3)

The factor “1/2” in Eq. 2 results from the average over the photon helicities
Kohri2017 .

$P_{a\rightarrow\gamma}$ (=$2P_{\gamma\rightarrow a}$) tends to be sizable and
constant, when $E_{\rm crit}<E<E_{\rm H}$ and $1\lesssim g_{a\gamma}Br/2$. For
the case of the Galaxy, the latter can be expressed as

$1\lesssim(\frac{r}{10\,\rm kpc})(\frac{B}{1.23\,\mu\rm
G})(\frac{g_{a\gamma}}{5\cdot 10^{-11}{\rm GeV}^{-1}}).$ (4)

Here, we refer to Ref Han2017  for the values corresponding to $r$ and $B$.

###  II.2 Magnetic field assumption

The magnetic fields around the AGN gamma-ray source commonly include those in
the jet, the radio lobes, and the host galaxy. The jet $B$ is believed to
decrease as the distance to the central engine along the jet axis
Tavecchio2015 ; Zheng2017 ; Meyer2014magnetic . At the VHE emission region,
the typical value of $B$ is in the interval of 0.1-5 G Tavecchio2015 ;
Zheng2017 ; Kang2014 ; Meyer2014magnetic ; Xue2016 ; Sullivan2009 . The
typical value of the magnetic field in the radio lobes is $10\,\mu$G for a
coherence length of 10 kpc IAXO2019 ; Pudritz2012 . The magnetic field in the
host galactic is poorly known and its strength is approximately equal to
$\mu$G with coherence lengths of the order of 0.1 to 0.2 kpc Widrow2002 ;
Meyer2014magnetic ; IAXO2019 . A part of VHE $\gamma-$ray AGNs are located in
galaxy clusters hardening2 , where the typical $B$ value is 5 $\mu$G with
coherence lengths of 10-100 kpc IAXO2019 ; Widrow2002 ; Fermi2016alp .

Our method takes advantage of the ALP-induced flux enhancement. It is more
sensitive to the average size of magnetic fields than the complicated magnetic
field configuration and detail with considerable uncertainty Meyer2013 ;
Meyer2014magnetic ; Fermi2016alp . Therefore, for simplicity, we assume the
“source” magnetic field is homogeneous within a region of 10 kpc, with
strength $B_{\rm s}$=10 $\mu$G, following Refs. Kohri2017 ; Hooper2007 .
Similarly, an average value of 1.23 $\mu$G for the Galactic magnetic field
(GMF) $B_{\rm GMF}$ with 10 kpc is adopted in our model Han2017 .

According to Eq. 4, the minimum coupling $g_{a\gamma}$ reaching the
significant transformation is about 2.5 $g_{11}$ (=$10^{-11}$ $\rm GeV^{-1}$)
for the GMF. As the larger $B_{\rm s}$, the minimum $g_{a\gamma}$
corresponding to the conversion in the source can be low to $g_{11}$.
Analogously, the critical energy for the conversion in the GMF, $\sim$160 TeV
for $g_{a\gamma}\simeq 2.5g_{11}$, $m_{a}\simeq 1\,\mu$eV, is higher than that
(the critical energy) for the source. However, the CM effect can suppress the
source conversion above 210 TeV for $g_{a\gamma}=g_{11}$ due to the larger
$B_{\rm s}$, as seen in Eq. 3, and it can be neglected in the case of GMF for
the energy considered in this paper.

We neglect the CM effect contributed by the CMB photons Dobrynina2015 ;
Montanino2017 , which is important when $E\gtrsim 460$ TeV for the B-fields
and the coupling ($g_{a\gamma}\gtrsim g_{11}$) considered in this work.

###  II.3 Photon survival probability

We employ the EBL model of Ref. Gilmore2012  to account the VHE photon
absorption onto the EBL. This recent EBL model has been tested repeatedly and
is generally consistent with the VHE $\gamma$-ray observations (e.g., Ref.
Fermi2012 ; HESS2013 ; VERITAS2015 ; Biteau2015 ; MAGIC2016 ; Armstrong2017 ;
Yuan2012 ; long2020 ). Furthermore, the infrared EBL intensity from this model
is in the mid-level among several recent EBL models Dominguez2011EBL ;
Finke2010 ; Franceschini2008 , which are more inconsistent but basically match
the direct measurements at infrared band HESS2017 . Hence, choosing this model
to account the EBL optical depth is helpful to reduce EBL uncertainty. Above
140 TeV, the CMB optical depth of TeV photons becomes dominant as its
intensity is much stronger than the EBL’s at the wavelength longer than 400
$\mu$m.

After obtaining the EBL/CMB spectrum, we can further estimate the optical
depth $\tau_{\gamma\gamma}(E,z)$ for the photon with energy $E$ from the
source of redshift $z$ Gilmore2012 ; Gong2013 . Then, the photon survival
probability on the whole path from the source to the earth can be derived

$P_{\gamma\rightarrow\gamma}=P_{\gamma\rightarrow\gamma}^{\rm S}{\rm
exp}(-\tau_{\gamma\gamma})P_{\gamma\rightarrow\gamma}^{\rm
G}+P_{\gamma\rightarrow a}^{\rm S}P_{a\rightarrow\gamma}^{\rm G},$ (5)

where $P_{\gamma\rightarrow a}^{\rm S}$ and $P_{a\rightarrow\gamma}^{\rm S}$
are the conversion probabilities from photons/ALPs to ALPs/photons in the
source respectively. There is a relation $P_{\gamma\rightarrow\gamma}^{\rm
S}=1-P_{\gamma\rightarrow a}^{\rm S}$; Similarly, the variables with a
superscript “G” represent those for the GMF. The derivation of Eq. 5 can be
illustrated vividly in Fig. 1: the first term corresponds to
$\gamma\rightarrow\gamma(\rm e^{\pm})\rightarrow\gamma$ channel suffered from
EBL/CMB absorption. the second is related to $\gamma\rightarrow
a\rightarrow\gamma$ channel unaffected by the absorption.

Figure 2: The photon survival probability $P_{\gamma\rightarrow\gamma}(E,z)$
on the whole path from the $\gamma-$ray source to the earth, where the ALP
mass $m_{a}=1\,\mu$eV and coupling $g_{a\gamma}=3g_{11}$. The meaning of each
colored curve is annotated in the diagram. $P_{1}=P_{\gamma\rightarrow a}^{\rm
S}P_{a\rightarrow\gamma}^{\rm G}$ and $P_{2}=P_{\gamma\rightarrow\gamma}^{\rm
S}{\rm exp}(-\tau_{\gamma\gamma}(E,0.01))P_{\gamma\rightarrow\gamma}^{\rm G}$,
i.e., the second and first (redshift-independent) term in Eq. 5. They
correspond to the channels of $\gamma\rightarrow a\rightarrow\gamma$ and
$\gamma\rightarrow\gamma(\rm e^{\pm})\rightarrow\gamma$ shown in Fig. 1,
respectively

Fig. 2 shows the change of $P_{\gamma\rightarrow\gamma}(E,z)$ with energy $E$
for different $z$, where the ALP mass $m_{a}=1\,\mu$eV and coupling
$g_{a\gamma}=3g_{11}$. $P_{2}=P_{\gamma\rightarrow\gamma}^{\rm S}{\rm
exp}(-\tau_{\gamma\gamma}(E,0.01))P_{\gamma\rightarrow\gamma}^{\rm G}$ and
$P_{1}=P_{\gamma\rightarrow a}^{\rm S}P_{a\rightarrow\gamma}^{\rm G}$, i.e.,
the first and second term in Eq. 5. In the low energy region, the conversion
is noneffective and $P_{\gamma\rightarrow\gamma}(E,z)$ is dominated by the
absorption term $P_{2}\sim{\rm exp}(-\tau_{\gamma\gamma})$. As the energy
turns to the higher region, ${\rm exp}(-\tau_{\gamma\gamma})\rightarrow 0$,
while the channel $\gamma\rightarrow a\rightarrow\gamma$ is getting “wider”,
and hence $P_{\gamma\rightarrow\gamma}\simeq P_{1}$, which is independent of
$z$. As a consequence, the curves of $P_{\gamma\rightarrow\gamma}(E,z)$ for
different $z$ at high energy region show v-shaped lines, and converge to
$P_{1}$. When $E>E_{\rm crit}\thickapprox$150 TeV, $P_{1}$ is getting closer
and closer to its maximum, and when $E>E_{\rm H}\thickapprox$630 TeV, the CM
effect suppresses the source conversion. Thus the peak appears at the highest
energy band.

Note that since we do not concern the spectral irregularities, we take the
average value for the square of the sine function in Eq. 2 when the phase is
larger than 1 rad, e.g., according to Refs. Kohri2017 ; Hooper2007 ;
Mirizzi2007 , smearing out the rapid-oscillatory features of the probability
function. The average value is approximatively taken 2/3 rather than 1/2 to
match the saturation-conversion probability ($P_{\gamma\rightarrow a}$) of
1/3, which corresponds to a more realistic scenario the beam propagates
through many domains of randomly oriented magnetic fields with constant size
$B$, for instance, in Refs. Mirizzi2007 ; Meyer2014 .

In the limit of saturated conversion $E_{\rm crit}\ll E\ll E_{\rm H}$ and
$1<g_{a\gamma}Br/2$, about

$P_{\gamma\rightarrow a}^{\rm S}P_{a\rightarrow\gamma}^{\rm
G}=\frac{1}{3}\times\frac{2}{3}$ (6)

of the original photons survive through $\gamma\rightarrow a\rightarrow\gamma$
channel. Obviously, for the ALP-parameter value applied to Fig. 2, the
condition of the saturated conversion does not match, for example
$g_{a\gamma}B_{\rm GMF}r_{\rm GMF}/2<1$.

Figure 3: Top panel and left of the bottom panels: fitting and extrapolating
the observations of M 87, IC 310 and Mrk 501. The blue and red lines
represents respectively the PLC or LPC fit with the EBL (Gilmore)/CMB-
absorption correction and that considering further the photon-ALP conversion
and the CM effect. $E_{c}$ is the cutoff energy assumed for the intrinsic
spectrum. The meaning of other symbols are indicated in the legend. Right of
the bottom panel: expected ALP limits based on our model as well as the future
LHAASO (or+CTA) observations of M 87 (blue dash line), IC 310 (green dash
line) and Mrk 501 (red dash line). For comparison, limits (black line) and 5
$\sigma$ sensitivities of future experiments (black dashed line) are also
shown.

##  III method

###  III.1 Sample selection

In our method, the VHE $\gamma-$ray observations are utilized to model the
intrinsic spectrum of the emitted source. So far about 75 VHE AGN have been
detected by the VHE instruments TEVCAT . In principle, most of these sources
may be used to search for the ALP-induced flux boost since it is independent
of the redshift above $E_{\rm crit}$, as shown in Fig. 2. But we should
acquire as many data below $E_{\rm crit}$ as possible, which are expected to
be slightly affected by the ALPs, so that a more realistic spectrum at higher
energy could be extrapolated by the observations together with the assumed
model (see Eq. 7). Hence, we preferentially consider the nearby sources whose
$P_{\gamma\rightarrow\gamma}$ curve can show a “shallow valley” due to the
relatively slight $\gamma\gamma$-absorption (see Fig. 2).

Based on the study of Franceschini et al. Franceschini2019 , the adjacent
sources of M 87, IC 310 and Mkn 501 would likely be detected by LHAASO up to
75 TeV, 50 TeV, 25 TeV respectively, when taking into account the standard
EBL-absorption. Thus, we predict that they could provide more detectable data
below $E_{\rm crit}$. In this paper, we will fit and extrapolate the spectral
data of these sources to the highest VHE energies.

_M 87_($z$=0.004)—a giant radio galaxy of Fanaroff-Riley type I with kpc radio
jet, located in the Virgo Cluster. It has been detected by almost all the
Imaging Air Cherenkov telescopes (IACTs) M87MAGIC2019 . Strong and rapid flux
variability in gamma-ray band is shown, but no significant spectral changes
with a typical photon index of 2.2 M87MAGIC2019 ; M87VERITAS2011 ; M87HESS2006
. We adopt H.E.S.S. data taken during 21 h of effective observation, in the
2005 12. Feb.-15. May high state (see Fig. 3).

_IC 310_($z$=0.019)—seems to be a transitional AGN between a low-luminosity
HBL (high-frequency peaked BL Lac, namely blazar with weak optical emission
lines Urry1995 ) and a radio galaxy Franceschini2019 , located on the
outskirts of the Perseus galaxy cluster. An extraordinary TeV flare in 2012
Nov. 12-13 and then a high state during several of the following months was
detected by MAGIC IC310MAGIC2017 . We use the observed spectrum with photon
spectral index 1.9 during 3.7 h of observation in the flare state (see Fig.
3).

_Mrk 501_($z$=0.034)—the next-closest known HBL. It is known for showing the
spectral variability at VHE band. During a famous outburst in 1997, the source
shown a very hard intrinsic spectrum with no softening up to the highest-
energy detected photons of 20 TeV Franceschini2019 . We choose the spectrum
detected by Fermi-LAT and MAGIC during the 4.5 month long multifrequency
campaign (2009 March 15 - August 1 during its relatively low activity)
Mrk501fermi2011  (see Fig. 3). Note that this BL Lac is also supposed to be
located in galaxy clusters hardening2 .

###  III.2 Theoretical and intrinsic spectra

We model VHE gamma-ray spectra with

$\psi_{0}=\rm e^{-\tau_{\gamma\gamma}}\phi\,\,\rm
or\,\,\psi_{1}=P_{\gamma\rightarrow\gamma}\phi,$ (7)

where $P_{\gamma\rightarrow\gamma}$ and $\tau_{\gamma\gamma}$ are defined in
Eq. 5. $\phi$ represents the intrinsic spectrum assumed for the sources. The
model with ALP has two additional free parameters, $g_{a\gamma}$ and $m_{a}$,
relative to the traditional model.

$\phi$ is assumed as one of the two common models HESS2013 : power-law with
exponential cut-off (PLC), and log-parabola with exponential cut-off (LPC).
The PLC spectrum is described by three parameters: $\phi_{\rm
PLC}=\phi_{0}(E/E_{0})^{-\alpha}{\rm exp}(-E/E_{\rm c})$, where $E_{\rm c}$ is
the cut-off energy, $\alpha$ is the photon spectral index constrained by the
particle acceleration theory as $\alpha\geq 1.5$, $\phi_{0}$ is the flux
normalization, and $E_{0}$ is the fixed reference energy. While the LPC
spectrum has additional curvature parameter $t>0$: $\phi_{\rm
LPC}=\phi_{0}(E/E_{0})^{-s-t\,{\rm log}(E/E_{0})}{\rm exp}(-E/E_{\rm c})$ and
also $\langle s+t\,{\rm log}(E/E_{0})\rangle\geq 1.5$.

Since the highest energy of detected photons in our samples is $\lesssim$ 20
TeV and the spectra are hard with no sign of convergence, the constrain on
$E_{\rm c}$ by the observations is very weak. If the parent particles
responsible for the VHE emission are electrons, the cutoff can be derived from
the Klein-Nishina suppression, energy loss of the electrons and pair
attenuation in the VHE emission region, see e.g., Lefa2012 ; Stawarz2008 ;
Lewis2019 ; Warren2020 ; Lemoine2020 ; Mrk501fermi2011 . Here, we uniformly
take $E_{\rm c}$=100 TeV for our samples Franceschini2019 , though $E_{\rm c}$
could be higher if the VHE $\gamma$-ray emission is of a hadronic origin
LHAASO2019 ; Xue2019 .

###  III.3 Fitting and extrapolating the observations

To simulate the observations at the highest energies, we firstly fit the three
observed spectra with $\psi_{0}$ and extrapolate it to hundreds of TeV
energies, respectively. Meanwhile, the form (PLC or LPC) of $\phi$ to be
chosen for each source is determined in terms of its average chi-square value
per degree of freedom. Then, we use $\psi_{1}$ containing the determined
$\phi$ and $P_{\gamma\rightarrow\gamma}$ with given $g_{a\gamma}$ and $m_{a}$
to fit and extrapolate the observations of each source.

We assume if the ALP-induced flux enhancement $\frac{\psi_{1}}{\psi_{0}}$ is
more than one order of magnitude and the predicted spectra $\psi_{1}$ is over
the equipment sensitivity, then the given ALP could be constrained. As the
continuity and (approximative) monotonicity of $P_{\gamma\rightarrow\gamma}$,
we only need to test the ALP parameters with small $g_{a\gamma}$ and that with
large $m_{a}$ to obtain the constrained region.

##  IV Results

Fig. 3 and Fig. 5 in the appendix report the predicted sensitivity limits for
future LHAASO/CTA observations of the three most promising nearby AGNs. The
blue line represents the standard absorption fit including an extrapolation
above 100 TeV. The results corresponding to the minimum allowable coupling
$g_{a\gamma}$ for different mass $m_{a}$ are shown by the red line. The 50-h 5
$\sigma$ sensitivity limits for CTA, 5-year and 1-year 5 $\sigma$ limits for
LHAASO are shown with the black dotted line, the black dashed line and the
black solid line, respectively. Combined observations from the instruments
will reach sensitivities of a few times $10^{-11}$ GeV $\rm cm^{-1}$ $\rm
s^{-1}$.

_M 87_. The befitting intrinsic spectrum for the VHE observations is PLC with
a photon index of 2.1. The predicted spectra extrapolated by the best-fit
model ($\psi_{0}$ and $\psi_{1}$) are above the 5-year sensitivity of LHAASO
up to 100 TeV (this will allow measurements of the M 87 spectrum up to about
70 TeV), which is beneficial to constrain the intrinsic spectrum. Above 100
TeV, the photon-ALP conversion gradually become important so that the photons
survive mainly thought $\gamma\rightarrow a\rightarrow\gamma$ channel.
Consequently, the flux enhancement is over an order of magnitude and the flux
is over the LHAASO sensitivity around 100 TeV for three given ALP-parameter
values. The line of $\psi_{1}$ shows a very or no “shallow valley”, as the
transition of survival probability $P_{\gamma\rightarrow\gamma}\approx P_{1}$
to $P_{\gamma\rightarrow\gamma}\approx P_{2}$ is smooth due to the low
redshift of M 87. At the highest energy, the intrinsic-spectrum cutoff makes
the curve go down.

_IC 310_. A PLC-intrinsic spectrum is used to model the intrinsic spectrum.
The predicted spectra extrapolated by the best-fit model ($\psi_{0}$ and
$\psi_{1}$) are above the sensitivity of LHAASO up to about 30 TeV, which
theoretically will allow measurements of the IC 310 spectrum up to 60 TeV.
Above $\sim$70 TeV, the photon-ALP conversion gradually become important, so
that the flux enhancement is over an order of magnitude and could be detected
by LHAASO around 100 TeV. But the intrinsic-spectrum cutoff or it together
with the CM effect makes the curve turn down at the highest energy band.

_Mrk 501_. A LPC-intrinsic spectrum is chosen. The predicted absorption-
corrected spectrum is over the sensitivity of LHAASO up to about 30 TeV, which
theoretically will give a relatively weak constraint on the intrinsic
spectrum. The very prominent enhanced flux at around 100 TeV is over the
LHAASO sensitivity.

We estimate the ALP parameter space that will be possible probed by the future
LHAASO (or+CTA) observations in the last picture of Fig. 3, where other limits
and sensitivity projections are also given for comparing. For M 87, LHAASO
would be able to explore $g_{a\gamma}$ down to about 2$g_{11}$ for
$m_{a}<0.3\,\mu$eV. For IC 310 and Mrk 501, a lower value of
$g_{a\gamma}\simeq g_{11}$ for $m_{a}<0.7\,\mu$eV and $m_{a}<0.4\,\mu$eV would
be explored respectively, some of which is invoked to explain the cold dark
matter Fermi2016alp . The results corresponding to IC 310 give stronger
exploitable bound on the coupling for $m_{a}\,\lesssim 1\mu$eV:
$g_{a\gamma}\simeq 2g_{11}$. In the case of M 87, a relatively weak bound is
given, as its observed flux is lower and its lower redshift leads to that
higher energy is required to achieve the same enhancement.

##  V Discussion

In this section, we will discuss our model assumptions.

As one of the effective conversion conditions requires the photon energy to
satisfy $E_{\rm crict}<E<E_{\rm H}$ and depends on the $B-$field, particularly
the highest energy conversions in the source $B-$field is prone to be
suppressed by the CM effect. Therefore, the uncertainty of $B_{\rm s}$ can
translate into an uncertainty of the photon survival probability and our
result. Fig. 4 shows how $B_{\rm s}$ affects $P_{\gamma\rightarrow\gamma}$ for
a fixed redshift $z=0.005$ and ALP parameters $m_{a}=1\,\mu$eV,
$g_{a\gamma}=2g_{11}$. For $B_{\rm s}$ in the range from several $\mu$G to 20
$\mu$G, the photon survival probabilities are close at around 200 TeV. It
means that provided $B_{\rm s}$ is between several $\mu$G and 20 $\mu$G the
ALP-induced flux enhancement could be achieved and comparable as done with
$B_{\rm s}=10\,\mu$G above. From this perspective, our result appears to be
robust.

We assume the exponential cutoff energy $E_{\rm c}=100$ TeV for the intrinsic
spectrum, and the ALP-induced flux would be detected for $g_{a\gamma}\gtrsim
2\times$$g_{11}$ with $m_{a}\lesssim 1\,\mu$eV. This result is sensitive to
$E_{\rm c}$, as the predicted photons around 100 TeV are mainly from the
conversion channel $\gamma\rightarrow a\rightarrow\gamma$. We therefore
investigate how the limits change if we alter the fiducial assumptions for the
$E_{\rm c}$ calculation. We find if $E_{\rm c}=200$ TeV, the limited ALP-
parameter region expands to $g_{a\gamma}\gtrsim 1\times$$g_{11}$ with
$m_{a}\lesssim 1\,\mu$eV, and $E_{\rm c}=50$ TeV, it reduces to
$g_{a\gamma}\gtrsim 3\times$$g_{11}$ with $m_{a}\lesssim 1\,\mu$eV.

The reasons we take $E_{\rm c}=100$ TeV is as follows. First, the spectrum
from our sample show very hard and have no sign of cutoff up to the highest
energy of about $10-20$ TeV. Second, if the spectrum is dominated by emission
of leptonic origin (with evidence that most of the rapid variable emission has
a leptonic origin), the cutoff above 100 TeV is possible. The recent
observation from the Crab Nebula with energy beyond 100 TeV show no
exponential cutoff below 100 TeV, which is usually interpreted in the
framework of leptonic models Amenomori2019 ; MAGIC2019 ; HAWC2019 . As
powerful cosmic particle accelerators Kotera2011 , that may happen on some
extreme TeV AGNs, too. Thirdly, AGNs are excellent candidates as Ultra-High-
Energy Cosmic Rays sources Mbarek2019 , and the hadronic cosmic rays are
capable of producing spectrum without cutoff below 100 TeV if the VHE emission
is dominated by hadronic origin CTA2017 . To determine the magnitude of
$E_{\rm c}$ without ambiguity, we need to further research on the intrinsic
physics (including parent particle species and its spectral energy
distribution, the radiation mechanism, and pair attenuation in the emission
region) of the $\gamma-$ray sources, as well as the forthcoming observations
above tens of TeV by CTA, LHAASO, SWGO and so on.

The spectrum of IC 310 we adopted was observed during a huge flare state that
lasted for only a few days by MAGIC, and nobody knows how many times such
strong flare or more occurs in one year. Even so, we simulate the observations
with 1 yr or 5 yr sensitivities of LHAASO, which will be able to continuously
survey every day the TeV sky to assess their suitability to constrain ALPs. In
this sense, the extrapolated flux from this spectrum should correspond to a
very high predicted one and thus its ALP limit should be treated as an
optimistic estimate.

Figure 4: The photon survival probabilities $P_{\gamma\rightarrow\gamma}$ for
different values of the source field magnetic $B_{\rm s}$. The redshift
$z=0.005$ and field magnetic region $r_{s}$=10 kpc. For $B_{\rm s}$ in the
range from several $\mu$G to 20 $\mu$G, the photon survival probabilities are
close at around 200 TeV. The CM effect almost completely suppresses the
photon-ALP conversion for $B_{\rm s}=50\,\mu$G.

##  VI Conclusion

In this article, we have discussed the potential of the gamma-ray spectrum of
AGN for energy up to above 100 TeV to probe ALP parameter space at around
$\mu$eV, where the coupling $g_{a\gamma}$ is so far relatively weak
constraint.

In case of conventional physics, most of the photons above tens of TeV emitted
from distant (distance$>$10 Mpc) AGN would be absorbed by the EBL/CMB during
its travel to the earth (see Fig 1 and 2). But more such photons, no matter
how far away, could survive, if we assume that the photon-ALP conversions
($\gamma\rightarrow a\rightarrow\gamma$) take place in the homogeneous source
and Galaxy magnetic-field for typical values of $l$=10 kpc, $B_{\rm s}$=10
$\mu$G, and $B_{\rm GMF}$=1.23 $\mu$G. Consequently, a very significant ALP-
induced flux enhancement, shaped as a peak, is expected to arise in the
observed spectrum above tens of TeV (see Fig 2). This provides the upcoming
LHAASO a good chance to detect the enhancement as its unprecedented
sensitivity above 30 TeV.

In order to acquire as many observations at tens of TeV as possible and thus
reduce the uncertainty from the intrinsic spectrum, the nearby and bright
sources, such as Mrk 501, IC 310 and M 87, are recommended to constrain the
ALPs around $\mu$eV. Assuming an intrinsic spectrum with exponentially
truncated at a fixed cutoff energy $E_{\rm c}$=100 TeV, we have extrapolated
the observed spectra of our sample up to above 100 TeV by the models
with/without ALPs. For $g_{a\gamma}\gtrsim 2\times$$10^{-11}\rm GeV^{-1}$ with
$m_{a}\lesssim 1\,\mu$eV, the flux at around 100 TeV predicted by the ALP
model can be more than an order of magnitude larger than that from the
standard absorption, and the enhanced flux could be detected by LHAASO (see
Fig 3 and 5).

Our result is subject to the uncertainty from the extrapolation of intrinsic
spectrum above tens of TeV. This will require further research on these
sources (Mrk 501, IC 310 and M 87) based on the forthcoming observations by
CTA, LHAASO, SWGO and so on.

###### Acknowledgements.

We would like to thank Weipeng Lin, P. H. T. Tam, Chengfeng Cai, Yu-Zhao Yu,
Seishi Enomoto, Yi-Lei Tang and Yu-Pan Zeng for useful discussions and
comments. This work is supported by the National Natural Science Foundation of
China (NSFC) under Grant No. 11875327, the Fundamental Research Funds for the
Central Universities, China, and the Sun Yat-Sen University Science
Foundation.

## References

  * (1) P. Svrcek, E. Witten, Axions In String Theory. J. High Energy Phys. 0606. 051 (2006).
  * (2) J. Jaeckel and A. Ringwald, The Low-Energy Frontier of Particle Physics. Ann. Rev. Nucl. Part. Sci. 60, 405 (2010).
  * (3) R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Pseudoparticles. Phys. Rev. Lett. 38, 1440 (1977);
  * (4) S. Weinberg, A New Light Boson?. Phys. Rev. Lett. 40, 223 (1978).
  * (5) F. Wiczek, Problem of Strong P and T Invariance in the Presence of Instantons. Phys. Rev. Lett. 40, 279 (1978).
  * (6) J. Preskill, M. B. Wise, and F. Wilczek, Cosmology of the invisible axion. Phys. Lett. B. 120, 127 (1983).
  * (7) L. F. Abbott and P. Sikivie, A cosmological bound on the invisible axion. Phys. Lett. B. 120, 133 (1983).
  * (8) M. Dine and W. Fischler, The Not So Harmless Axion. Phys. Lett. B. 120, 137 (1983).
  * (9) J. E. Marsh, Axiverse extended: Vacuum destabilization, early dark energy, and cosmological collapse. Phys. Rev. D. 83, 123526 (2011).
  * (10) G. Raffelt and L. Stodolsky, Mixing of the photon with low-mass particles. Phys. Rev. D. 37, 1237 (1988)
  * (11) P. Sikivie, Phys. Rev. Lett. 51, 1415 (1983).
  * (12) K. A. Hochmuth and G. Sigl, Effects of axion-photon mixing on gamma-ray spectra from magnetized astrophysical sources. Phys. Rev. D. 76, 123011 (2007).
  * (13) A. de Angelis, O. Mansutti, and M. Roncadelli, Axion-like particles, cosmic magnetic fields and gamma-ray astrophysics. Phys. Lett. B. 659, 847 (2008)
  * (14) D. Hooper and P. D. Serpico, Detecting Axionlike Particles with Gamma Ray Telescopes. Phys. Rev. Lett. 99, 231102 (2007).
  * (15) A. De Angelis, M. Roncadelli and O. Mansutti, Evidence for a new light spin-zero boson from cosmological gamma-ray propagation?. Phys. Rev. D. 76, 121301(R) (2007).
  * (16) A. Mirizzi, D. Montanino, and J. Cosmol. Stochastic conversions of TeV photons into axion-like particles in extragalactic magnetic fields. J. Cosmol. Astropart. Phys. 12 (2009) 004.
  * (17) D. Wouters and P. Brun, Irregularity in gamma ray source spectra as a signature of axionlike particles. Phys. Rev. D. 86, 043005 (2012).
  * (18) R. Bahre _et al_. (ALPS), Any light particle search II¡ªTechnical Design Report. Journal of Instrumentation. 8, T09001 (2013).
  * (19) S. Aune _et al_. (CAST), CAST search for sub-eV mass solar axions with 3He buffer gas. Phys. Rev. Lett. 107, 261302 (2011).
  * (20) L. D. Duffy, P. Sikivie, D. B. Tanner, S. J. Asztalos, C. Hagmann, D. Kinion, L. J. Rosenberg, K. van Bibber, D. B. Yu, and R. F. Bradley (ADMX), High resolution search for dark-matter axions. Phys. Rev. D. 74, 012006 (2006).
  * (21) M. Simet, D. Hooper, and P. D. Serpico, Milky Way as a kiloparsec-scale axionscope. Phys. Rev. D. 77, 063001 (2008).
  * (22) M. A. Sanchez-Conde, D. Paneque, E. Bloom, F. Prada, and A. Dominguez, Hints of the existence of Axion-Like-Particles from the gamma-ray spectra of cosmological sources. Phys. Rev. D. 79, 123511 (2009).
  * (23) A. De Angelis, G. Galanti, and M. Roncadelli, Relevance of axionlike particles for very-high-energy astrophysics. Phys. Rev. D. 84, 105030 (2011).
  * (24) A. Dominguez, M. A. Sanchez-Conde, and F. Prada, Axion-like particle imprint in cosmological very-high-energy sources. J. Cosmol. Astropart. Phys. 11 (2011) 020.
  * (25) A. De. Angelis, G. Galanti, and M. Roncadelli, Relevance of axion-like particles for very-high-energy astrophysics. Phys. Rev. D. 87, 109903(E) (2013).
  * (26) F. Tavecchio, M. Roncadelli, G. Galanti et al., Evidence for an axion-like particle from PKS 1222+216? Phys. Rev. D. 86, 085036(E) (2012).
  * (27) M. Meyer, D. Horns, and M. Raue, First lower limits on the photon-axion-like particle coupling from very high energy gamma-ray observations. Phys. Rev. D. 87, 035027 (2013).
  * (28) M. Meyer, and J. Conrad, Sensitivity of the Cherenkov Telescope Array to the detection of axion-like particles at high gamma-ray opacities. J. Cosmol. Astropart. Phys. 12 (2014) 016.
  * (29) S. Troitsky, Towards discrimination between galactic and intergalactic axion-photon mixing. Phys. Rev. D. 93, 045014 (2016).
  * (30) D. Montanino, F. Vazza, A. Mirizzi et al., Enhancing the Spectral Hardening of Cosmic TeV Photons by Mixing with Axionlike Particles in the Magnetized Cosmic Web. Phys. Rev. Lett. 119, 101101 (2017).
  * (31) R. Buehler, G. Gallardo, G. Maier, A. Dominguez, M. Lopez and M. Meyer. Search for the imprint of axion-like particles in the highest-energy photons of hard $\gamma$-ray blazars. J. Cosmol. Astropart. Phys. 09 (2020) 027.
  * (32) A. I. Nikishov, Sov. Phys. Pair production by a constant external field. JETP 393, 14(1962).
  * (33) M, G, Hauser and E, Dwek, THE COSMIC INFRARED BACKGROUND: Measurements and Implications. Annu. Rev. Astron. Astrophys. 39, 249 (2001).
  * (34) F. Aharonian et al. (H.E.S.S. Collaboration), A low level of extragalactic background light as revealed by $\gamma$-rays from blazars. Nature(London) 440,1018(2006).
  * (35) E. Dwek and A, Kusenko, The Extragalactic Background Light and the Gamma-ray Opacity of the Universe. Astropart. Phys. 43, 112 (2013).
  * (36) L. Costamante, Gamma-rays from Blazars and the Extragalactic Background Light. Int. J. Mod. Phys. D. 22, 1330025(2013).
  * (37) D. Horns, L. Maccione, M. Meyer, A. Mirizzi, D. Montanino, and M. Roncadelli, “Hardening of TeV gamma spectrum of active galactic nuclei in galaxy clusters by conversions of photons into axionlike particles,” Phys Rev D, 86, 075024(2012).
  * (38) G. Galanti, M. Roncadelli, Behavior of axionlike particles in smoothed out domainlike magnetic fields. Phys. Rev. D. 98, 043018 (2018).
  * (39) G. Galanti, F. Tavecchio, M. Roncadelli and C. Evoli, Blazar VHE spectral alterations induced by photon-ALP oscillations. Mon. Not. R. Astron. Soc. 487, 123 (2019).
  * (40) A. Franceschini, G. Rodighiero, and M. Vaccari, Extragalactic optical-infrared background radiation, its time evolution and the cosmic photon-photon opacity. Astron. Astrophys. 487, 837(2008).
  * (41) H. Abdalla et al. (CTA), Sensitivity of the Cherenkov Telescope Array for probing cosmology and fundamental physics with gamma-ray propagation. arXiv:2010.01349.
  * (42) A. Abramowski, F. Acero, F. Aharonian, F. Ait Benkhali, A. G. Akhperjanian, E. Anguner et al., Constraints on axionlike particles with H.E.S.S. from the irregularity of the PKS 2155-304 energy spectrum, Phys. Rev. D. 88, 102003 (2013).
  * (43) M. Ajello, A. Albert, B. Anderson, L. Baldini, G. Barbiellini, D. Bastieri et al., Search for Spectral Irregularities due to Photon-Axionlike-Particle Oscillations with the Fermi Large Area Telescope, Phys. Rev. Lett. 116,161101 (2016).
  * (44) C. Zhang, Y. F. Liang, S. Li et al., New bounds on axionlike particles from the Fermi Large Area Telescope observation of PKS 2155-304. Phys. Rev. D. 97, 063009 (2018).
  * (45) H. J. Li, J. G. Guo, X. J. Bi, S. J. Lin and P. F. Yin, Limits on axion-like particles from Mrk 421 with 4.5-years period observations by ARGO-YBJ and Fermi-LAT. arXiv:2008.09464.
  * (46) M. Libanov and S. Troitsky, On the impact of magnetic-field models in galaxy clusters on constraints on axion-like particles from the lack of irregularities in high-energy spectra of astrophysical sources. Phys. Lett. B. 802, 135252 (2020).
  * (47) P. Arias, D. Cadamuro, M. Goodsell, J. Jaeckel, J. Redondo and A. Ringwald, WISPy cold dark matter. J. Cosmol. Astropart. Phys. 6 (2012) 013.
  * (48) E. Armengaud et al. (IAXO), Physics potential of the International Axion Observatory (IAXO). J. Cosmol. Astropart. Phys. 06 (2019) 047.
  * (49) G. Sigl, Phys. Rev. D. Astrophysical haloscopes. 96, 103014 (2017).
  * (50) T. D. P. Edwards, M. Chianese, B. J. Kavanagh, S. M. Nissanke and C. Weniger, Unique Multimessenger Signal of QCD Axion Dark Matter. Phys. Rev. Lett. 124, 161101 (2020).
  * (51) A. Caputo, C. PenaGaray and S. J. Witte, Detecting the Stimulated Decay of Axions at Radio Frequencies. Looking for axion dark matter in dwarf spheroidal galaxies. Phys. Rev. D. 98, 083024 (2018)
  * (52) O. Ghosh, J. Salvado and J. Miralda-Escude, Axion Gegenschein: Probing Back-scattering of Astrophysical Radio Sources Induced by Dark Matter. arXiv:2008.02729.
  * (53) K. Kohri and H. Kodama, Axion-like particles and recent observations of the cosmic infrared background radiation. Phys. Rev. D. 96, 051701(R)(2017).
  * (54) M. Meyer, D. Montanino and J. Conrad, J. Cosmol. Astropart. On detecting oscillations of gamma rays into axion-like particles in turbulent and coherent magnetic fields. Phys. 09 (2014) 003.
  * (55) see http://tevcat.uchicago.edu/.
  * (56) X. Bai et al. (LHAASO collaboration), The large high altitude air shower observatory (LHAASO) science white paper. arXiv:1905.02773.
  * (57) F. Aharonian et al. (HESS collaboration), Observations of the Crab Nebula with H.E.S.S, Astron. Astrophys. 899, 457(2006).
  * (58) MAGIC collaboration, J. Aleksic et al., Performance of the MAGIC stereo system obtained with Crab Nebula data, Astropart. Phys. 435, 35 (2012).
  * (59) J. Holder, V.A. Acciari, E. Aliu, T. Arlen, M. Beilicke et al., Status of the VERITAS Observatory, AIP Conf. Proc. 657, 1085 (2009) [arXiv:0810.0474].
  * (60) M. S. Pshirkov, P. G. Tinyakov, and F. R. Urban, New limits on extragalactic magnetic fields from rotation measures. Phys. Rev. Lett. 116, 191302(2016).
  * (61) E. Masaki, A. Aoki, and Jiro Soda, Photon-Axion Conversion, Magnetic Field Configuration, and Polarization of Photons. Phys. Rev. D. 96, 043519 (2017).
  * (62) J. L. Han, Annu. Rev. Astron. Astrophys. 55, 1(2017).
  * (63) F. Tavecchio, M. Roncadelli, and G. Galanti, Photons to axion-like particles conversion in Active Galactic Nuclei. Phys. Lett. B. 744, 375 (2015).
  * (64) Y. G. Zheng, C. Y. Yang, L. Zhang, and J. C. Wang, Discerning the Gamma-Ray-emitting Region in the Flat Spectrum Radio Quasars. Astrophys. J. Suppl. Ser. 228, 1(2017).
  * (65) S. J. Kang, L. Chen, and Q. w. Wu, Constraints on the minimum electron Lorentz factor and matter content of jets for a sample of bright Fermi blazars. Astrophys. J. Suppl. Ser. 215, 5(2014).
  * (66) R. Xue, D. Luo, Z. R. Wang et al., Curvature of the spectral energy distribution, the inverse Compton component and the jet in Fermi 2LAC blazars. Mon. Not. R. Astron. Soc. 463, 3038 (2016).
  * (67) S. P. O. Sullivan and D.C. Gabuzda, Magnetic field strength and spectral distribution of six parsec-scale active galactic nuclei jets. Mon. Not. Roy. Astron. Soc. 400, 26 (2009).
  * (68) R. E. Pudritz, M. J. Hardcastle and D. C. Gabuzda, Magnetic fields in astrophysical jets: from launch to termination. Space Sci. Rev. 169, 27 (2012).
  * (69) L. M. Widrow, Reviews of Modern Physics. Origin of Galactic and Extragalactic Magnetic Fields. 74 775 (2002).
  * (70) A. Dobrynina, A. Kartavtsev, and G. Raffelt, Photon-photon dispersion of TeV gamma rays and its role for photon-ALP conversion. Phys. Rev. D. 91, 083003 (2015).
  * (71) R. C. Gilmore, R. S. Sommerville, J. R. Primack et al., Semi-analytic modelling of the extragalactic background light and consequences for extragalactic gamma-ray spectra. Mon. Not. R. Astron. Soc. 422, 3189 (2012).
  * (72) M. Ackermann et al. (Fermi-LAT Collaboration), The Imprint of The Extragalactic Background Light in the Gamma-Ray Spectra of Blazars. Science. 338, 1190(2012).
  * (73) A. Abramowski et al. (H.E.S.S. Collaboration). Measurement of the extragalactic background light imprint on the spectra of the brightest blazars observed with H.E.S.S. Astron. Astrophys. 550, A4 (2013).
  * (74) A. U. Abeysekara et al. (VERITAS Collaboration), Gamma-rays from the quasar PKS 1441+ 25: Story of an escape. Astrophys. J. 815, L22(2015).
  * (75) J. Biteau and D. A. Williams, The extragalactic background light, the Hubble constant, and anomalies: conclusions from 20 years of TeV gamma-ray observations. Astrophys. J. 812, 60(2015).
  * (76) M. L. Ahnen et al. (MAGIC Collaboration), MAGIC observations of the February 2014 flare of 1ES 1011+496 and ensuing constraint of the EBL density. Astron. Astrophys. 590, A24(2016).
  * (77) T. Armstrong, A. M. Brown, and P. M. Chadwick, Fermi-LAT high z AGN and the Extragalactic Background Light. Mon. Not. R. Astron. Soc. 470, 4089(2017).
  * (78) Q. Yuan, H. L. Huang, X. J Bi, and H. H. Zhang, Measuring the extragalactic background light from very high energy gamma-ray observations of blazars. arXiv:1212.5866.
  * (79) G. B. Long, W. P. Lin, P. H. T. Tam, and W. S. Zhu, Testing the CIBER cosmic infrared background measurements and axionlike particles with observations of TeV blazars. Phys. Rev. D. 101, 063004(2020).
  * (80) A. Dominguez, J. R. Primack, D. J. Rosario et al., Extragalactic background light inferred from AEGIS galaxy-SED-type fractions. Mon. Not. R. Astron. Soc. 410, 2556(2011).
  * (81) J. D. Finke, S. Razzaque, and C. D. Dermer, Modeling the extragalactic background light from stars and dust. Astrophys. J. 712, 238 (2010).
  * (82) H. Abdalla et al. (H.E.S.S. Collaboration), Measurement of the EBL spectral energy distribution using the VHE $\gamma$-ray spectra of H.E.S.S. blazars. Astron. Astrophys. 606, A59(2017).
  * (83) Y. Gong and A. Cooray, Astrophys. The extragalactic background light from the measurements of the attenuation of high-energy gamma-ray spectrum. J. Lett. 772, L12(2013).
  * (84) A. Mirizzi, G.G. Raffelt and P.D. Serpico, Signatures of axion-like particles in the spectra of TeV gamma-ray sources, Phys. Rev. D 76, 023001(2007).
  * (85) A. Franceschini, L. Foffano, E. Prandini, and F. Tavecchio, Very high-energy constraints on the infrared extragalactic background light. 629, A2 (2019).
  * (86) MAGIC collaboration, V. A. Acciari et al., Mon. Not. R. Astron. Soc. 492, 5354 (2020).
  * (87) VERITAS collaboration, E. Aliu et al., VERITAS Observations of Day-scale Flaring of M 87 in 2010 April. Astrophys. J. 746, 141 (2011).
  * (88) H.E.S.S. collaboration, F. Aharonian. Fast variability of TeV $\gamma$-rays from the radio galaxy M87. Sience 314,1424(2006).
  * (89) C. M. Urry and P. Padovani, Unified Schemes for Radio-Loud Active Galactic Nuclei. Publ. Astron. Soc. Jan. 107, 803(1995).
  * (90) MAGIC collaboration, M. L. Ahnen et al., First multi-wavelength campaign on the gamma-ray-loud active galaxy IC 310, Astron. Astrophys. 603, A25 (2017).
  * (91) Fermi-LAT, MAGIC and VERITAS collaboration. Insights into the high-energy gamma-ray emission of Markarian 501 from extensive multifrequency observations in the Fermi era, Astrophys. J. 727, 129 (2011).
  * (92) E. Lefa, S. R. Kelner, and F. A. Aharonian. On the spectral shape of radiation due to Inverse Compton Scattering close to the maximum cut-off. Astrophys. J. 753, 1176 (2012).
  * (93) L, Stawarz and V Petrosian. On the momentum diffusion of radiating ultrarelativistic electrons in a turbulent magnetic field. Astrophys. J. 681, 1725 (2008).
  * (94) T. R. Lewis, P. A. Becker and J. D. Finke. Electron Acceleration in Blazars: Application to the 3C 279 Flare on 2013 December 20. Astrophys. J. 884, 116 (2019).
  * (95) D. C. Warren, C. A. A. Beauchemin, M. V. Barkov, and S. Nagataki. The maximum energy of shock-accelerated electrons in a microturbulent magnetic field. arXiv:2010.06234.
  * (96) M. Lemoine and M. A. Malkov. Powerlaw spectra from stochastic acceleration. Mon. Not. R. Astron. Soc. 499, 4972 (2020).
  * (97) R. Xue, R.Y. Liu, M. Petropoulou, F. Oikonomou, Z. R.Wang, K. Wang, and X.Y. Wang. Powerlaw spectra from stochastic acceleration. Astrophys. J. 886, 23 (2019).
  * (98) M. Amenomori et al. (Tibet AS$\gamma$ Collaboration), First Detection of Photons with Energy beyond 100 TeV from an Astrophysical Source. Phys. Rev. Lett. 123, 051101 (2019).
  * (99) MAGIC collaboration, V. A. Acciari et al., MAGIC very large zenith angle observations of the Crab Nebula up to 100 TeV. Astron. Astrophys. 635, A158 (2019).
  * (100) HAWC collaboration, A. U. Abeysekara et al., Measurement of the Crab Nebula Spectrum Past 100 TeV with HAWC. Astrophys. J. 881, 134 (2019).
  * (101) K. Kotera and A. V. Olinto. The Astrophysics of Ultrahigh-Energy Cosmic Rays. Annu. Rev. Astron. Astrophys.49, 119 (2011).
  * (102) R. Mbarek, D. Caprioli. Bottom-up Acceleration of Ultra-High-Energy Cosmic Rays in the Jets of Active Galactic Nuclei. Astrophys. J. 886, 8 (2019).
  * (103) CTA collaboration, B.S. Acharya et al., Science with the Cherenkov Telescope Array. arXiv:1709.07997.

##  Appendix A fitting Results for the other 12 spectra

Figure 5: Fitting and extrapolating the observations of M 87, IC 310 and Mrk
501. The meanings of various symbols in the diagram are the same as Fig. 3.