# Axion-photon conversion of GRB221009A

Luohan Wang  Bo-Qiang Ma<EMAIL_ADDRESS>
###### Abstract

The newly observed gamma ray burst GRB221009A exhibits the existence of 10
TeV-scale photons, and the axion-photon conversion has been suggested as a
candidate to explain such energetic features of GRB221009A. In this work we
adopt a model to calculate the conversion probability of the energetic photons
from GRB221009A to the Earth. The result shows that the penetration
probability of photons with energy above $10^{1}$ TeV can be up to $10^{-2}$
to $10^{-4}$ depending on the coupling constant $g_{a\gamma}$ and the axion
mass $m_{a}$, together on the magnetic field parameters of the source galaxy
of GRB221009A. We show that the parameters of the source galaxy, with the
magnetic field handled by a cellular model, contribute a lot of uncertainties
to the penetration probability, so we have more freedom to reconcile a variety
of axionlike particle parameters from other observations with the Large High
Altitude Air Shower Observatory observation. By comparing the results in this
article with the data from Large High Altitude Air Shower Observatory, we can
obtain more precise constraints on the ranges of these parameters.

###### keywords:

axion-like particle, axion-photon conversion, gamma ray burst, very-high-
energy photons, GRB221009A

††journal: arXiv:2304.01819, published in Phys. Rev. D 108 (2023) 023002

[aff1] organization = School of Physics, Peking University, city = Beijing
100871, country = China

[aff2] organization = Center for High Energy Physics, Peking University, city
= Beijing 100871, country = China

[aff3] organization = Collaborative Innovation Center of Quantum Matter, city=
Beijing, country = China

##  1 Introduction

During the propagation of very-high-energy (VHE, $E_{\gamma}\geq 0.1~{}$TeV)
photons in the Universe, the photon flux should be significantly attenuated by
the interaction of these photons with background photons to annihilate into
electron-positron pairs [1, 2, 3]. This mechanism is supposed to strongly
suppress the VHE spectra from distant sources and thus to forbid the detection
of VHE photons. However, with the oscillation of photons into axionlike
particles (ALPs), such suppression can be avoided [4].

ALP is a sort of generalization of the axion, the pseudo-Goldstone boson
associated with Peccei-Quinn symmetry to solve the strong-$CP$ problem [5].
ALPs are very light pseudoscalar bosons represented by $a$. Numerous ALP
couplings to the standard model particles can be considered, including the
case where ALPs couple to two photons. The axions that we consider in this
work are characterized by their coupling with two photons. In the presence of
transverse magnetic field, photons can translate into axions, and vice visa
[6, 7]. For VHE photons from extragalactic sources, the photons can transition
into photon-ALP mixed flux in the source galaxy and the axion part can
propagate unimpeded through the intergalactic medium, while the VHE photons
are absorbed by the extragalactic background photons. When passing through the
Galaxy magnetic field, the back conversion of axions to photons can then bring
back a fraction of VHE photons. Thus, the conclusions from QED are
significantly modified in the presence of ALPs. The above scenario was first
proposed to explain the excess of very-high-energy gamma ray that should
undergo attenuation by the extragalactic background photons [8], later with
source $B$ field or the intergalactic magnetic fields considered [9, 10, 11,
12] and then with the Milky Way magnetic field also included [13].

On October 9, 2022, an extremely powerful gamma-ray burst named GRB221009A,
located at RA=288.282 and Dec=19.495 with z=0.1505, was detected by several
observatories. The Large High Altitude Air Shower Observatory (LHAASO)
detected 64000 of very high energy photons above $0.2~{}$TeV, with the maximum
of energy from above 10 TeV up to 18 TeV [14, 15]. Assuming that the photons
observed by LHAASO indeed came from GRB221009A, the detection of such
energetic photons is forbidden under standard models [16, 17], but is
reasonable under the hypothesis of the photon-ALP oscillation [4, 18, 19, 20,
21]. Alternative explanations, such as the Lorentz invariance violation of
cosmic high-energy photons [16, 17, 22, 23, 24], photonic decay of heavy
neutrinos from the GRB source [25, 26, 27, 28] or large uncertainty in the
energy reconstruction of the observed events [16, 17], are also proposed to
explain the LHAASO observation of above 10 TeV photon events.

In this paper, we analyze the penetration probability of TeV-scale photons
from GRB221009A in the presence of ALPs. The result shows that the penetration
probability of photons with energy above $10^{1}$ TeV can be up to $10^{-2}$
to $10^{-4}$, indicating the viability for using the hypothetical axions to
explain the energetic features of GRB221009A. What is more, this paper gives
the quantitative effects of the source galaxy magnetic field and the choice of
the axion-photon coupling constant $g_{a\gamma}$ and the axion mass $m_{a}$ on
photons with different energy, and calculates the penetration probability for
different parameters. For $g_{a\gamma}\in[0.5,2.1]\times 10^{-11}\
\mathrm{GeV}^{-1}$ and $m_{a}\in[0.01,20]\times 10^{-8}\ \mathrm{eV}$, the
maximum of total penetration probability for 18 TeV photons is
$10^{-2}$-$10^{-4}$, based on the magnetic field parameters chosen at the
source galaxy with a cellular model. The maximum saturation back conversion
probability in the Milky Way is $3.6\times 10^{-2}$, for $g_{a\gamma}$ around
$2\times 10^{-11}\ \mathrm{GeV}^{-1}$. By comparing the result in this article
to the data from LHAASO, we can give more precise constraints on the ranges of
the coupling constant $g_{a\gamma}$ and the mass of axions $m_{a}$, together
on the magnetic field parameters of the source galaxy of GRB221009A.

##  2 The propagation of high-energy photons under the Standard Model

###  2.1 Absorption of VHE photons by background photons

The cosmic microwave background (CMB) lights and extragalactic background
lights (EBL) are two significant background radiations in the Universe.
According to the standard model, there exists an interaction that creates an
electron-positron pair from two photons. Suppose that a high-energy photon
collides with a low-energy background photon and produces an electron-positron
pair $\gamma\gamma\rightarrow e^{+}e^{-}$. Setting $c=\hslash=k_{B}=1$ in the
natural unit, we have

$p^{\mu}p_{\mu}\geq 4m_{e}^{2},$ (1)

where $p^{\mu}=p^{\mu}_{\gamma}+p^{\mu}_{\mathrm{bg}}$, with
$p^{\mu}_{\gamma}$ and $p^{\mu}_{\mathrm{bg}}$ being the four-momentum of the
gamma ray photon and the background photon. Therefore, the threshold energy of
the reaction is

$E_{\mathrm{thr}}=\frac{2m_{e}^{2}}{\epsilon(1-\cos\theta)},$ (2)

where $\epsilon$ is the energy of background photons and $\theta$ is the
colliding angle between the three momentum of two photons.

The cross section for the reaction of $\gamma\gamma\rightarrow e^{+}e^{-}$ is
[29]

$\sigma_{\gamma\gamma}(v)=\frac{\pi}{2}\frac{\alpha^{2}}{m_{e}^{2}}(1-v^{2})\left[(3-v^{4})\ln\frac{1+v}{1-v}-2v(2-v^{2})\right],$
(3)

where
$v(E,\epsilon,\theta)=\sqrt{1-\frac{2m_{e}^{2}}{E\epsilon(1-\cos\theta)}}$ is
the speed of generated electron-positron pair in the center-of-momentum frame.

The optical depth for EBL absorption is

$\begin{split}\tau_{\gamma}(E_{\gamma})=\int_{0}^{z_{0}}\frac{\text{d}x}{\text{d}z}\text{d}z\int_{-1}^{1}\frac{1-\cos\theta}{2}\text{d}\cos\theta\\\
\int_{\epsilon_{\mathrm{thr}}(E_{\gamma}(1+z),\theta)}^{\infty}\text{d}\epsilon\
\sigma_{\gamma\gamma}(v(E_{\gamma}(1+z),\epsilon,\theta))\frac{\text{d}n(\epsilon)}{\text{d}\epsilon},\end{split}$
(4)

where $E_{\gamma}$ is the observed energy of gamma-ray photons, $n(\epsilon)$
is the number density of EBL photons and
$\epsilon_{\mathrm{thr}}={2m_{e}^{2}}/{E_{\gamma}(1+z)}(1-\cos\theta)$ denotes
the threshold energy of background photons reacting with gamma photons of
observed energy $E_{\gamma}$ at redshift $z$, and

$\frac{\text{d}x}{\text{d}z}=\frac{c}{H_{0}\sqrt{\Omega_{m}(1+t)^{3}+\Omega_{\Lambda}}},$
(5)

where $H_{0}=73\ \text{km}\cdot\text{Mpc}^{-1}\cdot\text{s}^{-1}$ is the
Hubble constant, $\Omega_{m}=0.24$ is the energy density and
$\Omega_{\Lambda}=1-\Omega_{m}$ is the dark energy density (in flat
cosmology).

Hence, without the assumption of ALPs, the observed energy intensity is
$I_{\mathrm{ob}}(E_{\gamma})=I_{0}(E_{0})e^{-\tau_{\gamma}(E_{\gamma})}$,
where $I_{0}$ is the emitted spectrum with initial photon energy
$E_{0}=E_{\gamma}(1+z)$.

###  2.2 Interaction with EBL photons

For CMB photons, the mean energy is $\epsilon\approx 6.35\times 10^{-4}$ eV.
Therefore, the corresponding threshold energy is
$E_{\mathrm{thr}}^{\text{CMB}}\approx 411$ TeV, which is not in our
consideration since the VHE photons detected by LHAASO in GRB221009A have the
energy of the order around 10 TeV. EBL photons have the energy between
$10^{-3}$ to 1 eV. Consequently, the corresponding threshold energies for
collision with EBL photons are between 261 GeV to 261 TeV. Thus, EBL is the
dominant contribution to the suppression of 10 TeV-scale photons. The solid
curve in Fig. 1 is the EBL penetration probability $(e^{-\tau_{\gamma}})$ for
TeV-scale photons from $z=0.1505$, which is the redshift of GRB221009A.
Apparently, VHE photons are strongly suppressed by EBL, and photons with high
energy can hardly survive, as already shown in Refs. [16, 17].

Figure 1: The EBL penetration probability $(e^{-\tau_{\gamma}(E_{\gamma})})$
as a function of observed gamma-ray energy of TeV photons from $z=0.1505$
under the standard model. The EBL model is taken from Ref. [30] and
http://side.iaa.es/EBL/. One sees that 10 TeV scale photons are strongly
suppressed by EBL, as already shown in Refs. [16, 17].

The EBL model picked for the above calculation has been compared with other
EBL models in the original EBL model paper in Ref. [30] and discussed in Refs.
[16, 17], from which we know that the background EBL photons with larger
energy are more suppressed than other EBL models, so the EBL model we adopted
predicts a larger penetration ability for high energy photons to reach the
Earth than other EBL models. Choice of other model calculations would offer
more strong conflict between the observation of 18 TeV photons with
predictions from the standard model [16, 17].

##  3 Axions coupling with photons

In this section we introduce the principle for the $\gamma\gamma a$ coupling
process and give the method of calculating conversion probability.

###  3.1 General theory

The axions are characterized by their coupling with two photons, which plays a
key role in most searches. The Lagrangian of the interaction is

$\mathcal{L}=-\frac{g_{a\gamma}}{4}F_{\mu\nu}\tilde{F}^{\mu\nu}a=g_{a\gamma}\vec{E}\cdot\vec{B}a,$
(6)

where $F$ is the electromagnetic field-strength tensor and
$\tilde{F}_{\mu\nu}=\frac{1}{2}\epsilon_{\mu\nu\rho\sigma}F^{\rho\sigma}$ is
its dual, and $g_{a\gamma}$ is the coupling constant of the axion-photon
oscillation.

By considering very relativistic axions ($m_{a}<<\omega$), the equation of the
photon and axion amplitudes reduces to a first-order propagation equation. If
a beam of polarized monochromatic light with energy $E$ travels along the
$z$-direction, then we have

$\displaystyle(E-i\partial_{z}+\mathcal{M})\begin{pmatrix}A_{x}\\\ A_{y}\\\
a\end{pmatrix}=0,$ (7)

where $E$ is the energy of gamma ray photons, $A_{x}$ and $A_{y}$ are the
photon polarization amplitudes along the $x$-axis and the $y$-axis, $a$ is the
ALP amplitude, and $\mathcal{M}$ denotes the photon-axion mixing matrix. The
mixing matrix is

$\displaystyle\mathcal{M}=\begin{pmatrix}\Delta_{xx}&\Delta_{xy}&\Delta_{a\gamma}\cos\theta\\\
\Delta_{yx}&\Delta_{yy}&\Delta_{a\gamma}\sin\theta\\\
\Delta_{a\gamma}\cos\theta&\Delta_{a\gamma}\sin\theta&\Delta_{a}\end{pmatrix},$
(8)

where $\Delta_{a}=-{m_{a}^{2}}/{2E}$,
$\Delta_{a\gamma}={g_{a\gamma}B_{T}}/{2}$, $B_{T}$ is the transverse magnetic
field perpendicular to the $z$ direction, and $\theta$ is the angle between
$B_{T}$ and $x$ axis, i.e. $B_{x}=B_{T}\cos\theta,B_{y}=B_{T}\cos\theta$. The
term $\Delta_{ij}$ mixing $A_{x}$ and $A_{y}$ is determined both by the
properties of the medium and the QED vacuum polarization effect.

If we set the $y$ axis along the direction of $B_{T}$, then we can decouple
the $x$ component, and Eq. (8) reduces to

$\displaystyle\mathcal{M}=\begin{pmatrix}\Delta_{\perp}&\Delta_{R}&0\\\
\Delta_{R}&\Delta_{\parallel}&\Delta_{a\gamma}\\\
0&\Delta_{a\gamma}&\Delta_{a}\end{pmatrix},$ (9)

where
$\Delta_{\perp,\parallel}=\Delta_{\mathrm{pl}}+\Delta^{\mathrm{CM}}_{\perp,\parallel}$.
In a plasma, the photons acquire an effective mass given by
$\omega_{\mathrm{pl}}^{2}={4\pi\alpha n_{e}}/{m_{e}}$, leading to
$\Delta_{\mathrm{pl}}=-{\omega_{\mathrm{pl}}^{2}}/{2E}$. $m_{e}$ is the mass
of electron, $n_{e}$ is the electron density in the medium, and $\alpha$ is
the fine structure constant. The $\Delta^{\mathrm{CM}}_{\perp,\parallel}$ term
describes the Cotton-Mouton effect, i.e., the birefringence of fluids in the
presence of a transverse magnetic field where
$|\Delta^{\mathrm{CM}}_{\perp}-\Delta^{\mathrm{CM}}_{\parallel}|\propto
B_{T}^{2}$. These terms are of little importance to the following topics and
can be neglected [31]. $\Delta_{R}$ is the Faraday rotation term. It depends
on the energy and the longitudinal component $B_{z}$ and causes the coupling
between $A_{x}$ and $A_{y}$. Since Faraday rotation angle is proportion to
$E^{-2}$ and is negligible for VHE photons, $\Delta_{R}$ can be omitted.

For the relevant parameters, we can numerically calculate

$\displaystyle\Delta_{a\gamma}$ $\displaystyle\simeq 3.1\times
10^{-2}\left(\frac{g_{a\gamma}}{2\times 10^{-11}\
\text{GeV}^{-1}}\right)\left(\frac{B_{T}}{\mathrm{\mu
G}}\right)\mathrm{kpc}^{-1},$ (10) $\displaystyle\Delta_{a}$
$\displaystyle\simeq-7.8\times 10^{-3}\left(\frac{m_{a}}{10^{-8}\
\mathrm{eV}}\right)^{2}\left(\frac{E}{\mathrm{TeV}}\right)^{-1}\mathrm{kpc}^{-1},$
(11) $\displaystyle\Delta_{\mathrm{pl}}$ $\displaystyle\simeq-1.1\times
10^{-10}\left(\frac{n_{e}}{10^{-3}\
\mathrm{cm^{3}}}\right)\left(\frac{E}{\mathrm{TeV}}\right)^{-1}\mathrm{kpc}^{-1}.$
(12)

The coupling constant $g_{a\gamma}\leq 2.1\times 10^{-11}\ \text{GeV}^{-1}$
for $m_{a}\leq 2\times 10^{-7}\ \text{eV}$ is based on the result in Ref.
[32]. Reference [20] also gives the condition for $g_{a\gamma}$ and $m_{a}$ to
detect 18 TeV photons. In this paper, we use the combination of conditions
from the two articles. In the following section, we take $g_{a\gamma}=2\times
10^{-11}\ \mathrm{GeV}^{-1}$ and $m_{a}=10^{-8}\ \mathrm{eV}$ as an
approximation for the parameters.

###  3.2 Conversion of $\gamma\gamma\to a$ in the source galaxy

The gamma ray burst is possibly generated by the collapsing of an extremely
massive stellar or a newly born, very rapidly rotating and highly magnetized
neutron stars [33, 34]. However, in both models, collimated relativistic jets
pummel through the surface of its progenitor and radiate as gamma rays. The
jet of gamma rays is strongly centered and propagates through the Universe as
a narrow focused beam. The process of $\gamma\rightarrow a\rightarrow\gamma$
has three steps. First, after going out of the compact region of the source,
the photon flux translates into photon-axion mixture due to the interaction
with the magnetic filed in the source galaxy. Then the flux propagates through
intergalactic area, while VHE photons are absorbed by background photons and
the corresponding axions are left undisturbed. Finally, the beam arrived at
the Milky Way, and the high energy axions can converse back into VHE photons.
The analyses for the second and third steps are in subsection 3.3 and
subsection 3.4. In this section, we focus on the conversion of
$\gamma\rightarrow a$ in the source galaxy.

Due to the poor knowledge we have about the magnetic field near the GRB
source, for simplicity, we can assume that the magnetic field in the source
galaxy is composed of many small regions with the same scale $s$. The
direction of $\vec{B}$ is the same in each region and changes randomly from
one region to another. Such a model can be called the cellular model (see Ref.
[35]). In fact, we will show that the structure of the cellular model is
analogy to our Galaxy magnetic field, so it is reasonable to assume that the
source galaxy magnetic field has a cellular structure.

The magnetic field of our Galaxy (i.e., the Milky Way) is relatively well
known and can be handled by the Jansson model (see Ref. [36]). The Galactic
magnetic field consists of a regular component and a turbulent component. The
regular component maintains its strength and direction within several kpc, and
the direction of the magnetic field can change abruptly at the interface
between different regions, as described in the Jansson model. From the
viewpoint of the cellular model, we can treat the Galactic magnetic field as a
composition of many small regions with the direction and strength of each
region magnetic field adjusted by the Jansson model. If we only concern about
an averaged effect such as the propagation of a particle through a number of
small regions, we can also adopt the cellular model to treat the magnetic
field of each region with random arranged direction and fixed strength. We
will show that the regular component of the Milky Way magnetic field can be
effectively treated by the cellular model with a scale of $s=4$ kpc to handle
the back conversion of axion to photon in subsection 3.4. The turbulent
component, with the scale of $10^{-2}$ pc, which dominates on small scale,
however, is proved to be irrelevant to the conversion based on the discussion
at the end of subsection 3.2. Therefore, if we assume that the magnetic field
of the source galaxy holds similar properties to the Galactic magnetic field,
the assumption of a cellular model is reasonable.

Therefore, in each small region, we can decouple the $x$ component by setting
the $y$ axis parallel to the transverse magnetic component $B_{T}$, and Eq.
(7) reduces to a two-dimension equation consequently,

$\displaystyle(E-i\partial_{z}+\mathcal{M}_{2})\begin{pmatrix}A_{2}\\\
a\end{pmatrix},\hfill$ (13)

with

$\displaystyle\mathcal{M}_{2}=\begin{pmatrix}\Delta_{\mathrm{pl}}&\Delta_{a\gamma}\\\
\Delta_{a\gamma}&\Delta_{a}\end{pmatrix}.\hfill$ (14)

We can diagonalize $\mathcal{M}_{2}$ by a rotation angle

$\theta=\frac{1}{2}\arctan(\frac{2\Delta_{a\gamma}}{\Delta_{\mathrm{pl}}-\Delta_{a}}).$
(15)

Then we can arrive at the conversion probability ($\gamma$ to $a$) in each
small region

$\displaystyle P_{0}(\gamma\rightarrow a)$
$\displaystyle=|\bra{A_{2}(0)}{a(s)}\rangle|^{2},$ (16)
$\displaystyle=\sin^{2}(2\theta)\sin^{2}(\frac{\Delta_{\mathrm{osc}}s}{2}),$
(17)
$\displaystyle=(\Delta_{a\gamma}s)^{2}\frac{\sin^{2}(\frac{\Delta_{\mathrm{osc}}s}{2})}{(\frac{\Delta_{\mathrm{osc}}s}{2})^{2}},$
(18)

where the oscillation wave number is

$\Delta_{\mathrm{osc}}^{2}=(\Delta_{\mathrm{pl}}-\Delta_{a})^{2}+4\Delta^{2}_{a\gamma}.$
(19)

If the flux passes through a magnetic field with length $L$ composed of $n$
randomly arranged small regions, then the total conversion probability is

$\displaystyle P_{\gamma\rightarrow a}$
$\displaystyle=\frac{1}{3}\left[1-(1-\frac{3}{2}P_{0})^{n}\right],$ (20)
$\displaystyle=\frac{1}{3}\left(1-e^{-\frac{3P_{0}n}{2}}\right),\ (P_{0}\
\mathrm{is\ small}),$ (21)

the detailed derivation can be seen in Refs. [35, 31]. Apparently,
$P_{\gamma\rightarrow a}$ reaches its maximum of 1/3 when $P_{0}n$ is big
enough.

If $2\Delta_{a\gamma}<<|\Delta_{\mathrm{pl}}-\Delta_{a}|$, i.e., when $E$ is
of GeV scale, $P_{0}$ is proportional to
$(\frac{g_{a\gamma}B_{T}}{m_{a}^{2}})^{2}\sin^{2}(\frac{m_{a}^{2}}{4E})E^{2}$,
but is relatively small. We can thereby define a critical energy $E_{c}$
satisfying

$|\Delta_{a}(E=E_{c})-\Delta_{\mathrm{pl}}(E=E_{c})|\equiv 2\Delta_{a\gamma},$
(22)

below which $P_{0}$ is negligible. Then $P_{0}$ can be written as

$P_{0}=(\Delta_{a\gamma}s)^{2}\frac{\sin^{2}(\Delta_{a\gamma}s\sqrt{1+(E_{c}/E)^{2}})}{(\Delta_{a\gamma}s\sqrt{1+(E_{c}/E)^{2}})^{2}}.$
(23)

Note that $E_{c}\propto m_{a}^{2}/g_{a\gamma}B_{T}$ since
$\Delta_{\mathrm{pl}}$ is negligible compared with $\Delta_{a\gamma}$ and
$\Delta_{a}$. For the parameters from Eq. (12), and $g_{a\gamma}\sim 10^{-11}\
\mathrm{GeV}^{-1}$, $m_{a}\sim 10^{-8}\ \mathrm{eV},E_{c}$ is around 0.1 TeV.
Since we expect $\gamma\rightarrow a$ conversion in the source galaxy and back
conversion in the Milky Way for VHE photons, $E_{c}$ is expected to be less
than 1 TeV, therefore we have the relation

$(\frac{m_{a}}{10^{-8}\ \mathrm{eV}})^{2}(\frac{g_{a\gamma}}{10^{-11}\
\mathrm{GeV}})^{-1}(\frac{B}{\mathrm{\mu}G})^{-1}\leq 2\times 10^{2}.$ (24)

When $2\Delta_{a\gamma}>>|\Delta_{\mathrm{pl}}-\Delta_{a}|$ (i.e., $E>>E_{c}$)
or in any case $\Delta_{\mathrm{osc}}s\rightarrow 0$, the conversion
probability $P_{0}$ is independent of energy, which is just the case in this
article. What is more, based on the parameter chosen in this article, the
oscillation length $2\pi/\Delta_{\mathrm{osc}}$ is much larger than the domain
size $s$. This means the discontinuity at the interface of two adjacent
domains is not felt by the oscillation [37]. In this case, the application of
the cellular model is successful despite the abrupt change of $B_{T}$.

The cosmic accelerators producing photons may host strong magnetic fields
$B\sim O(1\to 10)\ \mathrm{\mu G}$ (see Refs. [38, 39]). Suppose that the
regular component of magnetic field has $s\simeq 3\ \mathrm{kpc}$, $L\simeq
30\ \mathrm{kpc}$, and that $n_{e}\simeq 1.0\times 10^{-3}\ \mathrm{cm^{3}}$.
We can plot the conversion probability of $\gamma$ to $a$ in the source galaxy
by the regular component as functions of photon energy $E$ as shown in Fig. 2.

Figure 2: The conversion probability of $\gamma$ to $a$ in the source galaxy
caused by the regular component with different strength. $B_{T}\simeq
O(\mathrm{\mu G})$. For photons with energy above 0.1 TeV, $P$ is a constant.
For photons with energy below 0.1 TeV, $P$ is an oscillation term dominated by
a small coefficient.

The conversion probability of $\gamma$ to $a$ for VHE photons is between 0.05
and 0.33 for $B_{T}$ from $1\ \mathrm{\mu G}$ to $10\ \mathrm{\mu G}$. We will
take $B_{T}=5\ \mathrm{\mu G}$ as an average field strength [40] in the
analysis below. We also take $s\simeq 1\ \mathrm{kpc}$, $L\simeq 10\
\mathrm{kpc}$ as the minimum typical lengths of the source galaxy, which are
actually small parameters for most galaxies.

As to the turbulent component, $s\simeq 0.01\ \mathrm{pc}$, causing
$P_{0}\sim\sin^{2}(\Delta_{a\gamma}s)\sim(\Delta_{a\gamma}s)^{2}$ to vanish
for the chosen parameters. For example, in the case where $B_{T}=10\
\mathrm{\mu G}$, $s=0.01\ \mathrm{pc}$, $L=30\ \mathrm{kpc}$, the conversion
probability is $10^{-5}$. Therefore, the regular component of magnetic field
is the dominant factor in the photon-axion oscillation.

###  3.3 Propagation in the intergalactic area

As has been discussed in Ref. [4], the suggestion of adopting the axion-photon
conversion as a mechanism to explain the high energy features of GRB 221009A
requires the following conditions: (i) photon-to-axion conversion upon leaving
the source galaxy or cluster; (ii) non-significant re-conversion along the
intergalactic medium (IG) in the path to observer; (iii) axion-to-photon
conversion upon entry into the Milky Way.

The intergalactic magnetic field is very weak and has an order of magnitude
varying from $10^{-11}$ to several nG. Assume that the strength of magnetic
field in intergalactic area is $B_{T}\sim 1~{}\mathrm{nG}$. Therefore, the
critical energy $E_{c}\sim 128$ TeV based on the chosen parameters, is far
above the detection feasibility. In fact, we expect $E_{c}$ to be larger than
18 TeV, which is the highest energy for photons detected by LHAASO, and there
is

$(\frac{m_{a}}{10^{-8}\ \mathrm{eV}})^{2}(\frac{g_{a\gamma}}{10^{-11}\
\mathrm{GeV}})^{-1}(\frac{B}{\mathrm{n}G})^{-1}\geq 3.5\times 10^{-2}.$ (25)

As a consequence, the photon-axion conversion process for photons with energy
$E<E_{c}$ is negligible in the intergalactic area.

Using the combination of Eq. (24) and Eq. (25), we can arrive at the same
range of axion parameters required by the detection of a 18 TeV photon as
shown in Fig. 1 in Ref. [20].

During the propagation in the intergalactic area, the intensity of flux is
suppressed as shown in Fig. 3. The intensity of the flux arriving at the
Galaxy is

$\displaystyle I_{\gamma}=(1-P_{\gamma\rightarrow a})e^{-\tau_{\gamma}}I_{0},$
(26) $\displaystyle I_{a}=P_{\gamma\rightarrow a}I_{0},$ (27)

where $I_{0}$ is the initial flux intensity at the source.

By ignoring the polarization of photons, the probability amplitude of the
mixed flux at the edge of the Milky Way is

$\begin{pmatrix}A_{x}\\\ A_{y}\\\
a\end{pmatrix}=\begin{pmatrix}\sqrt{\frac{I_{\gamma}}{2}}\\\
\sqrt{\frac{I_{\gamma}}{2}}\\\ \sqrt{I_{a}}\end{pmatrix}.$ (28)

###  3.4 Back conversion in the Milky Way

As shown in subsection 3.2, we can neglect the turbulent component in the
Galaxy and discuss the effect of the regular component specifically.

Measurements from Faraday rotation show that the regular component is
basically parallel to the disk, with a small vertical component. By choosing
the Jansson model (see Ref. [36]) as the Galaxy magnetic field model, the
regular component can be split into a disk part, a toroidal halo part and a
X-shaped part.

The disk component is composed of eight logarithmic spirals located on the
Galaxy disk, and the field strength falls off as $r^{-1}$ in each spiral. In
the vertical direction, the field strength falls off as a function of $z$,
with typical length $h\simeq 0.4$ kpc. The electron density inside the Milky
Way disk is $n_{e}\simeq 1.1\times 10^{-3}\ \mathrm{cm^{3}}$, resulting in
$\omega_{\mathrm{pl}}\simeq 4.1\times 10^{-12}\ \mathrm{eV}$. The toroidal
halo component is completely azimuthal. It functions mainly in the area
slightly above the disk, and has a typical size of 9.22 kpc ($z>0$) or 16.7
kpc ($z<0$). The X-shaped component has a “X” shape on the $r$-$z$ plane, and
the detailed properties can be seen in Ref. [36].

In the Galaxy center, there is a possible dipole component with stronger $B$.
However, since photons from GRB221009A do not pass through that area, we can
focus on the outer regions. What is more, we suppose that the magnetic field
strength is zero for $r>20$ kpc.

The propagation of photons in the Milky Way is actually a three dimensional
problem, which is more complicated in calculation than simply estimated by the
cellular model as we have done in subsection 3.2. However, we can still reduce
$\mathcal{M}$ to the following expression, which is similar to Eq. (9)

$\displaystyle\mathcal{M}=\begin{pmatrix}\Delta_{\mathrm{pl}}&0&\frac{g_{a\gamma}}{2}B_{x}\\\
0&\Delta_{\mathrm{pl}}&\frac{g_{a\gamma}}{2}B_{y}\\\
\frac{g_{a\gamma}}{2}B_{x}&\frac{g_{a\gamma}}{2}B_{y}&\Delta_{a}\end{pmatrix}.$
(29)

Since $E\sim 10^{12}\ \text{eV}>>g_{a\gamma}B\sim 10^{-28}\ \text{eV}$, we can
take $\mathcal{M}$ as a perturbation of $E$, $A=A^{(0)}+A^{(1)}$, where
$A^{(0)}$ is the zero order term and $A^{(1)}$ is the first order term. The
zero order term of Eq. (7) is

$i\partial_{z}A^{(0)}=EA^{(0)}.$ (30)

Therefore,
${\color[rgb]{0,0,0}{A^{(0)}(z)}}=\mathcal{U}_{0}(z)A^{(0)}{\color[rgb]{0,0,0}{(0)}}$,
where

$\mathcal{U}_{0}(z)=\exp(-i\int_{0}^{z}\text{d}z^{\prime}E).$

In the interaction representation, we have
$A_{\mathrm{int}}=\mathcal{U}_{0}^{\dagger}A$,
$\mathcal{M}_{\mathrm{int}}=\mathcal{U}^{\dagger}_{0}\mathcal{M}\
\mathcal{U}_{0}=\mathcal{M}$, where $\mathcal{M}$ from Eq. (29) and Eq. (7)
can be written as

$i\frac{\partial A_{\mathrm{int}}}{\partial
z}=\mathcal{M}_{\mathrm{int}}A_{\mathrm{int}},$ (31)

and

$A_{\mathrm{int}}(z)=A_{\mathrm{int}}(0)-i\int_{0}^{z}\text{d}z^{\prime}\mathcal{M}_{\mathrm{\mathrm{int}}}(z^{\prime})A_{\mathrm{int}}(z^{\prime}),$
(32)

where
$A^{(0)}_{\mathrm{int}}(z)=\mathcal{U}_{0}^{\dagger}A^{(0)}(z)=\mathcal{U}_{0}^{\dagger}\
\mathcal{U}_{0}A^{(0)}(0)=A^{(0)}(0)$.

Preserving the zero order term of $A_{\mathrm{int}}$, we have
$A_{\mathrm{int}}(0)=A^{(0)}_{{\mathrm{int}}}(0)=A^{(0)}(0)$. We can
numerically calculate Eq. (32) based on the Jansson and Farrar model, with
mixing matrix given in Eq. (29) along the galactic line of sight of
GRB221009A. The initial probability amplitude is given in Eq. (28) and the
location of GRB221009A is given in Ref. [41].

Since there is only a phase shift between $A(z)$ and $A_{\mathrm{int}}(z)$,
the total penetration probability for photons from GRB221009A to the Earth is
$|A_{\mathrm{int},x}(l)|^{2}+|A_{\mathrm{int},y}(l)|^{2}$, where $l\simeq
23.95$ kpc is the distance from the edge of the Milky Way magnetic field to
the Earth. We can plot the total penetration probability as a function of
$E_{\gamma}$. For example, assuming that $g_{a\gamma}=2\times 10^{-11}\
\text{GeV}^{-1}$, $m_{a}=10^{-8}\ \text{eV}$, with varying $B$, $L$, $s$ at
the source, we have Fig. 3.

Figure 3: The total penetration probability with $g_{a\gamma}=2\times
10^{-11}\ \text{GeV}^{-1},m_{a}=10^{-8}\ \text{eV}$ for different parameters
at the source galaxy. For $s\leq 1\ \mathrm{kpc}$, $n={L}/{s}=10$,
$B_{s}\in[1,10]\ \mathrm{\mu G}$, the error is shown in the figure. The dot-
dashed line is the penetration probability for $B_{T}=5\ \mathrm{\mu G}$,
$s=3\ \mathrm{kpc}$, $L=30\ \mathrm{kpc}$, while the dotted line is the
penetration probability for $B_{T}=1\ \mathrm{\mu G}$, $s=1\ \mathrm{kpc}$,
$L=10\ \mathrm{kpc}$.

The upper limit of the probability does not grow with $s>3$ kpc since Eq. (21)
has converged to the value of 1/3. For different parameters given in the
source, the penetration probability ranges from $10^{-4}$ to $10^{-2}$,
depending on the properties of magnetic field at the source galaxy.
Apparently, the existence of axion-photon oscillations can strongly enlarge
the penetration probability of VHE photons.

The back conversion probability in the Milky Way also converges to a constant
as $E_{\gamma}\rightarrow\infty$, which can be seen intuitively in Fig. 4,
where $p$ denotes the back conversion probability at
$E_{\gamma}=18~{}\mathrm{TeV}$. Defining
$g_{11}=g_{a\gamma}/10^{-11}\mathrm{GeV}^{-1}$, we take $g_{11}\geq 0.5$ so
that the back conversion in the Milky Way is not too weak. However, note that
the real penetration probability for photons with energy above $E_{c}$ in the
intergalactic magnetic field is not a constant, since these axions may convert
into photons and be absorbed by EBL. The correction to the standard model will
only be significant for photons with energy between
$\mathrm{max}(E_{c,\mathrm{source}},E_{c,\mathrm{MW}})$ and
$E_{c,\mathrm{IG}}$, where $E_{c,\mathrm{source}}$, $E_{c,\mathrm{MW}}$, and
$E_{c,\mathrm{IG}}$ denote the critical energy in the source galaxy, the Milky
way, and the intergalactic area respectively.

Figure 4: The back conversion probability in the Milky Way with different
axion parameters $g_{a\gamma}$ and $m_{a}$. In the figure $g_{11}$ denotes
$g_{a\gamma}/(10^{-11}~{}\text{GeV}^{-1})$ and $p$ denotes the corresponding
probability at $E_{\gamma}=18~{}\mathrm{TeV}$. The maximum of back conversion
probability is $3.6\times 10^{-2}$.

What is more, we can calculate the penetration probability with varying
$g_{a\gamma}$ and $m_{a}$ for photons with energy of 18 TeV as shown in Fig.
5. We take the values of $g_{a\gamma}$ and $m_{a}$ from the combined limit of
Refs. [18, 19, 20, 21] to explain the LHAASO observation and Refs. [43, 44,
45] to meet constraints from other astrophysical observations:
$g_{a\gamma}\in[0.5,2.1]\times 10^{-11}\ \mathrm{GeV}^{-1}$ and
$m_{a}\in[0.01,20]\times 10^{-8}\ \mathrm{eV}$. From Fig. 5, we notice that
the calculated penetration probability is sensible to the axion-photon
coupling $g_{a\gamma}$ but not so sensitive to the axion mass $m_{a}$.

Figure 5: The penetration probability for a 18 TeV photon with varying axion
parameters $g_{a\gamma}\in[0.5,2.1]\times 10^{-11}\ \mathrm{GeV}^{-1}$ and
$m_{a}\in[0.01,20]\times 10^{-8}\ \mathrm{eV}$. The picture above has
$B_{T}=10\ \mathrm{\mu G}$, $L=30\ \text{kpc}$, $s=3\ \text{kpc}$
corresponding to the maximum of penetration, while the lower picture has
$B_{T}=1\ \mathrm{\mu G}$, $L=10\ \text{kpc}$, $s=1\ \text{kpc}$ corresponding
to the minimum of penetration.

If we restrict $m_{a}<2\times 10^{-8}\ \mathrm{eV}$, in the range given in
Ref. [32], then the penetration probability for 18 TeV photons from GRB221009A
is $10^{-2}$ for $g_{11}=2.1$.

We still need the detected data for LHAASO to determine the specific values of
$g_{a\gamma}$ and $m_{a}$, however, the analysis in this article shows that
the existence of axions gives the penetration probability of VHE photons
within $10^{-2}$ to $10^{-4}$, which does explain the detection of VHE
photons, and gives the criterion to judge the ranges of $g_{a\gamma}$ and
$m_{a}$.

According to Fig. 5, if we know the order of magnitude of the penetration
probability, we can basically constrain the parameters at the source, and
focus on the choice of $g_{a\gamma}$ and $m_{a}$. Then we can limit the values
of $g_{a\gamma}$ and $m_{a}$ to more narrow ranges.

The ALP parameters can explain the LHAASO observation variances with
significant differences in literature [18, 19, 20, 21]:
$g_{a\gamma}\in[0.5,1]\times 10^{-11}\ \mathrm{GeV}^{-1}$ and
$m_{a}\in[0.01,20]\times 10^{-8}\ \mathrm{eV}$. There are also constraints
from other astrophysical observations, such as these from the analysis of
cosmic sub-PeV gamma rays from the Crab Nebula with $g_{a\gamma}\leq 1.8\times
10^{-10}\ \mathrm{GeV}^{-1}$ and $m_{a}\leq 2\times 10^{-7}$– $6\times
10^{-7}\ \mathrm{eV}$ [42], from galactic sub-PeV gamma rays [43] with
$g_{a\gamma}\leq 2.1\times 10^{-11}\ \mathrm{GeV}^{-1}$ for $m_{a}\leq 2\times
10^{-7}\ \mathrm{eV}$, from optical polarisation measurements of magnetic
white dwarfs [44] with $g_{a\gamma}\leq 5.4\times 10^{-12}\ \mathrm{GeV}^{-1}$
and $m_{a}\leq 3\times 10^{-7}\ \mathrm{eV}$, and from gamma-ray spectra of
flat-spectrum radio quasars [45] with $g_{a\gamma}\leq 5\times 10^{-12}\
\mathrm{GeV}^{-1}$. We show that by adjusting the parameters of the source
galaxy, we can reproduce the required magnitude to explain the LHAASO
observation of above 10 TeV photons with a variety of
$g_{a\gamma}\in[0.5,2.1]\times 10^{-11}\ \mathrm{GeV}^{-1}$ and
$m_{a}\in[0.01,20]\times 10^{-8}\ \mathrm{eV}$, as shown in Fig. 4 and Fig. 5.
The novelty of the present work, in comparison with previous literature [18,
19, 20, 21], reveals the large uncertainties introduced by the parameters of
the source galaxy to reconcile a variety of axion parameters with the LHAASO
observation. The parameters of the cellular model about the source galaxy can
be further constrained to narrower ranges when there will be more information
about the source and the magnetic field of the source galaxy.

In our above calculations, the results are obtained by adopting the realistic
Jansson model to treat the magnetic field of the Milky Way. If we adopt the
cellular model to calculate the back conversion of axion to photon conversion,
we find that we need effective parameters $s=4$ kpc, $B_{T}=1\ \mathrm{\mu
G}$, $L=20$ kpc for getting the same results obtained with the Jansson model.
This proves the efficiency of adopting the cellular model to handle the
magnetic field of the Milky Way. Therefore we can use the simple cellular
model of galactic magnetic field to handle both the source galaxy and our
local Galaxy flexibly.

##  4 Conclusions

In this work, we use axionlike particles, which appear in many extensions of
the standard model, to explain the energetic features of the newly detected
GRB221009A. In the presence of the axion-photon conversion process
$\gamma\rightarrow a\rightarrow\gamma$, we calculate the conversion
probability, back conversion probability and the total penetration probability
of VHE photons with varying parameters at the source galaxy, the coupling
constant $g_{g\gamma}$, and the mass of axions $g_{a\gamma}$. The presence of
axions can significantly increase the penetration probability of VHE photons,
with a maximum of $10^{-2}$. Even with the smallest conversion probability
assumed in the source galaxy, the penetration probability can reach $10^{-4}$.

Though we have given the penetration probability for different $g_{a\gamma}$
and $m_{a}$, however, as shown in Fig. 5, the parameters in the source galaxy,
with the magnetic field handled by a cellular model, contribute a lot of
uncertainties to the penetration probability. This provides us more freedom to
reconcile a variety of axionlike paricle parameters from other observations
with the LHAASO observation of above 10 TeV photons. More precise ranges of
the parameters can be obtained by confronting our theoretical prediction with
the data from LHAASO. Further studies should be based on more experimental
results, especially VHE photons from galaxy clusters or other extragalactic
sources, in order to reach more precise constraints on $g_{a\gamma}$ and
$m_{a}$, together on the parameter space of the source galaxy. Alternative
possibilities, such as the Lorentz invariance violation or heavy neutrino
decay, should be also checked by confronting various features of experimental
observations with theoretical predictions comprehensively.

## Acknowledgments

This work is supported by National Natural Science Foundation of China (Grant
No. 12075003).

## References

  * [1] R. Gould and G. Schréder, Opacity of the Universe to High-Energy Photons, Phys. Rev. Lett. 16, no.6, 252-254 (1966) doi:10.1103/PhysRevLett.16.252.
  * [2] G. G. Fazio and F. W. Stecker, Predicted high energy break in the isotropic gamma-ray spectrum: A test of cosmological origin, Nature 226, 135-136 (1970) doi:10.1038/226135a0.
  * [3] R. J. Protheroe and H. Meyer, An infrared background TeV gamma-ray crisis?, Phys. Lett. B 493, 1-6 (2000) doi:10.1016/S0370-2693(00)01113-8 [arXiv:astro-ph/0005349 [astro-ph]].
  * [4] See, e.g., G. Zhang and B.-Q. Ma, Axion-Photon Conversion of LHAASO Multi-TeV and PeV Photons, Chinese Phys. Lett. 40, 011401 (2023) doi:10.1088/0256-307X/40/1/011401 [arXiv:2210.13120 [hep-ph]], and references therein.
  * [5] R. D. Peccei and H. R. Quinn, CP Conservation in the Presence of Instantons, Phys. Rev. Lett. 38, 1440-1443 (1977) doi:10.1103/PhysRevLett.38.1440.
  * [6] P. Sikivie, Phys. Rev. Lett. 51, 1415-1417 (1983) [erratum: Phys. Rev. Lett. 52, 695(E) (1984)], doi:10.1103/PhysRevLett.51.1415.
  * [7] E. Masso and R. Toldra, On a light spinless particle coupled to photons, Phys. Rev. D 52, 1755-1763 (1995) doi:10.1103/PhysRevD.52.1755 [arXiv:hep-ph/9503293 [hep-ph]].
  * [8] 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) doi:10.1103/PhysRevD.76.121301 [arXiv:0707.4312 [astro-ph]].
  * [9] A. De Angelis, O. Mansutti, M. Persic and M. Roncadelli, Photon propagation and the VHE gamma-ray spectra of blazars: how transparent is really the Universe?, Mon. Not. Roy. Astron. Soc. 394, L21-L25 (2009) doi:10.1111/j.1745-3933.2008.00602.x [arXiv:0807.4246 [astro-ph]].
  * [10] 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) doi:10.1103/PhysRevD.79.123511 [arXiv:0905.3270 [astro-ph.CO]].
  * [11] A. Mirizzi and D. Montanino, Stochastic conversions of TeV photons into axion-like particles in extragalactic magnetic fields, JCAP 12 (2009) 004, doi:10.1088/1475-7516/2009/12/004 [arXiv:0911.0015 [astro-ph.HE]].
  * [12] A. Dominguez, M. A. Sanchez-Conde and F. Prada, Axion-like particle imprint in cosmological very-high-energy sources, JCAP 11 (2011) 020, doi:10.1088/1475-7516/2011/11/020 [arXiv:1106.1860 [astro-ph.CO]].
  * [13] M. Simet, D. Hooper, and P. D. Serpico, Milky Way as a kiloparsec-scale axionscope, Phys. Rev. D 77, 063001 (2008) doi:10.1103/PhysRevD.77.063001 [arXiv:0712.2825 [astro-ph]].
  * [14] Yong Huang, Shicong Hu, Songzhan Chen, et al., LHAASO observed GRB 221009A with more than 5000 VHE photons up to around 18 TeV, GCN Circ. 32677 (2022).
  * [15] Z. Cao et al. [LHAASO Collaboration], A tera–electron volt afterglow from a narrow jet in an extremely bright gamma-ray burst, Science doi:10.1126/science.adg9328 (2023) doi:10.1126/science.adg9328 [arXiv:2306.06372].
  * [16] H. Li and B.-Q. Ma, Lorentz invariance violation induced threshold anomaly versus very-high energy cosmic photon emission from GRB 221009A, Astropart. Phys. 148, 102831 (2023) doi:10.1016/j.astropartphys.2023.102831 [arXiv:2210.06338 [astro-ph.HE]].
  * [17] H. Li, B.-Q. Ma, Searching Lorentz invariance violation from cosmic photon attenuation, Eur. Phys. J. C 83, 192 (2023) doi:10.1140/epjc/s10052-023-11334-z [arXiv:2210.05563].
  * [18] G. Galanti, M. Roncadelli and F. Tavecchio, [arXiv:2210.05659 [astro-ph.HE]].
  * [19] A. Baktash, D. Horns and M. Meyer, [arXiv:2210.07172 [astro-ph.HE]].
  * [20] S.V. Troitsky, Parameters of axion-like particles required to explain high-energy photons from GRB 221009A, JETP Letters (Pis’ma ZhETF) 116, 745 (2022) doi:10.1134/S0021364022602408 [arXiv:2210.09250 [astro-ph.HE]].
  * [21] S. Nakagawa, F. Takahashi, M. Yamada and W. Yin, Axion dark matter from first-order phase transition, and very high energy photons from GRB 221009A, Phys. Lett. B 839, 137824 (2023) doi:10.1016/j.physletb.2023.137824 [arXiv:2210.10022 [hep-ph]].
  * [22] J. D. Finke and S. Razzaque, Possible Evidence for Lorentz Invariance Violation in Gamma-Ray Burst 221009A, Astrophys. J. Lett. 942, no.1, L21 (2023) doi:10.3847/2041-8213/acade1 [arXiv:2210.11261 [astro-ph.HE]].
  * [23] J. Zhu and B.-Q. Ma, Light speed variation from GRB 221009A, J. Phys. G 50, no.6, 06LT01 (2023) doi:10.1088/1361-6471/accebb [arXiv:2210.11376 [astro-ph.HE]].
  * [24] H. Li and B.-Q. Ma, Revisiting Lorentz invariance violation from GRB 221009A, [arXiv:2306.02962 [astro-ph.HE]].
  * [25] K. Cheung, [arXiv:2210.14178 [hep-ph]].
  * [26] V. Brdar and Y. Y. Li, Neutrino origin of LHAASO’s 18 TeV GRB221009A photon, Phys. Lett. B 839, 137763 (2023) doi:10.1016/j.physletb.2023.137763 [arXiv:2211.02028 [hep-ph]].
  * [27] J. Huang, Y. Wang, B. Yu and S. Zhou, Invisible neutrino decays as origin of TeV gamma rays from GRB221009A, JCAP 04, 056 (2023) doi:10.1088/1475-7516/2023/04/056 [arXiv:2212.03477 [hep-ph]].
  * [28] S. Y. Guo, M. Khlopov, L. Wu and B. Zhu, [arXiv:2301.03523 [hep-ph]].
  * [29] G. Breit and J. A. Wheeler, Collision of two light quanta, Phys. Rev. 46 (12), 1087-1091 (1934) doi:10.1103/PhysRev.46.1087.
  * [30] A. Domínguez, J. R. Primack, D. J. Rosario, F. Prada, R. C. Gilmore, S. M. Faber, D. C. Koo, R. S. Somerville, M. A. Pérez-Torres, P. Pérez-González, J.-S. Huang, M. Davis, P. Guhathakurta, P. Barmby, C. J. Conselice, M. Lozano, J. A. Newman, M. C. Cooper, Extragalactic background light inferred from AEGIS galaxy-SED-type fractions, Monthly Notices of the Royal Astronomical Society 410 (4), 2556--2578 (2011)  doi:10.1111/j.1365-2966.2010.17631.x.
  * [31] A. Mirizzi, G. Raffelt, P. Serpico, Photon-Axion Conversion in Intergalactic Magnetic Fields and Cosmological Consequences, Lecture Notes in Physics 741, 115–134 (2008) doi:10.1007/978-3-540-73518-2$\\_$7, Berlin Heidelberg: Springer; 2008.
  * [32] L. Mastrototaro, P. Carenza, M. Chianese, D. F. G. Fiorillo, G. Miele, A. Mirizzi and D. Montanino, Constraining axion-like particles with the diffuse gamma-ray flux measured by the Large High Altitude Air Shower Observatory, Eur. Phys. J. C 82 (11), 1012 (2022)  doi:10.1140/epjc/s10052-022-10979-6.
  * [33] J. Granot, T. Piran, O. Bromberg, J. L. Racusin, F. Daigne, Gamma-ray bursts as sources of strong magnetic fields, Space Science Reviews 191 (1-4), 471--518 (2015)  doi:10.1007/s11214-015-0191-6.
  * [34] A. I. MacFadyen, S. E. Woosley, Collapsars: Gamma-ray bursts and explosions in “failed supernovae”, The Astrophysical Journal 524 (1), 262--289 (1999)  doi:10.1086/307790.
  * [35] Y. Grossman, S. Roy and J. Zupan, Effects of initial axion production and photon axion oscillation on type Ia supernova dimming, Phys. Lett. B 543, 23-28 (2002) doi:10.1016/S0370-2693(02)02448-6 [arXiv:hep-ph/0204216 [hep-ph]].
  * [36] R. Jansson, G. R. Farrar, A new model of the galactic magnetic field, The Astrophysical Journal 757 (1), 14 (2012)  doi:10.1088/0004-637X/757/1/14.
  * [37] G. Galanti, M. Roncadelli,Behavior of axionlike particles in smoothed out domainlike magnetic fields, Phys. Rev. D 98, 043018 (2018)  doi:10.1103/PhysRevD.98.043018.
  * [38] G. Galanti, F. Tavecchio, M. Roncadelli, C. Evoli, Blazar VHE spectral alterations induced by photon–ALP oscillations, Monthly Notices of the Royal Astronomical Society 487 (1), 123--132 (2019)  doi:10.1093/mnras/stz1144.
  * [39] F. Tavecchio, M. Roncadelli, G. Galanti, G. Bonnoli, Evidence for an axion-like particle from PKS 1222+216?, Physical Review D 86 (8), 085036 (2012)  doi:10.1103/physrevd.86.085036.
  * [40] A. Fletcher, Magnetic fields in nearby galaxies (2011), published in The Dynamic Interstellar Medium: A Celebration of the Canadian Galactic Plane Survey, Astronomical Society of the Pacific Conference Series Vol. 438, 197-210 (2011),  arXiv:1104.2427.
  * [41] R. Pillera, E. Bissaldi, N. Omodei, G. La Mura, F. Longo, Fermi-LAT. team, GRB 221009a: Fermi-LAT refined analysis, GRB Coordinates Network 32658 (2022).
  * [42] X. J. Bi, Y. Gao, J. Guo, N. Houston, T. Li, F. Xu and X. Zhang, Axion and dark photon limits from Crab Nebula high energy gamma-rays, Phys. Rev. D 103, no.4, 043018 (2021) doi:10.1103/PhysRevD.103.043018 [arXiv:2002.01796 [astro-ph.HE]].
  * [43] C. Eckner, F. Calore, First constraints on axionlike particles from galactic sub-PeV gamma rays, Physical Review D 106 (8), 083020 (2022)  doi:10.1103/physrevd.106.083020.
  * [44] C. Dessert, D. Dunsky and B. R. Safdi, Upper limit on the axion-photon coupling from magnetic white dwarf polarization, Phys. Rev. D 105, 103034 (2022)  doi:10.1103/PhysRevD.105.103034 [arXiv:2203.04319 [hep-ph]].
  * [45] L. Dirson and D. Horns, Phenomenological modelling of the Crab Nebula’s broadband energy spectrum and its apparent extension, Astron. Astrophys. 671, A67 (2023)  doi:10.1051/0004-6361/202243578 [arXiv:2203.11502 [astro-ph.HE]].