# Strain effects in phosphorous bound exciton transitions in silicon

Teemu Loippo  Department of Physics and Nanoscience Center, University of
Jyväskylä, P.O. Box 35, FI-40014 University of Jyväskylä, Finland  Antti
Kanniainen  Department of Physics and Nanoscience Center, University of
Jyväskylä, P.O. Box 35, FI-40014 University of Jyväskylä, Finland  Juha T.
Muhonen<EMAIL_ADDRESS>Department of Physics and Nanoscience Center,
University of Jyväskylä, P.O. Box 35, FI-40014 University of Jyväskylä,
Finland

###### Abstract

Donor spin states in silicon are a promising candidate for quantum information
processing. One possible donor spin readout mechanism is the bound exciton
transition that can be excited optically and creates an electrical signal when
it decays. This transition has been extensively studied in bulk, but in order
to scale towards localized spin readout, microfabricated structures are needed
for detection. As these electrodes will inevitably cause strain in the silicon
lattice, it will be crucial to understand how strain affects the exciton
transitions. Here we study the phosphorous donor bound exciton transitions in
silicon using hybrid electro-optical readout with microfabricated electrodes.
We observe a significant zero-field splitting as well mixing of the hole
states due to strain. We can model these effects assuming the known asymmetry
of the hole g-factors and the Pikus-Bir Hamiltonian describing the strain. In
addition, we describe the laser power, electric-field and light polarization
dependence of the transitions. Importantly, the hole-mixing should not prevent
donor electron spin readout and using our measured parameters and numerical
simulations we anticipate that hybrid spin readout in a silicon-on-insulator
platform should be possible, allowing integration to silicon photonics
platforms.

##  I Introduction

Donor spin qubits in silicon are a promising candidate for processing quantum
information [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]. Their application potential is
however somewhat constrained by the lack of an optical interface, and
considerable attention has been recently focused on spin-photon interfaces in
silicon [11, 12, 13, 14, 15]. One possible readout avenue is the donor bound
exciton transition [16, 17, 18, 19, 20], which however decays mostly via Auger
recombination hence not enabling an efficient optical readout protocol. Hybrid
electro-optical readout, where the spin-selective transition is excited
optically but readout electrically, is possible and has been demonstrated with
both direct electrical and capacitive readout[21] in bulk. Scaling the
electro-optical readout towards single-spin level will, however, require
moving to microfabricated structures as well as most likely taking advantage
of silicon photonics platforms. So far, only one demonstration using
microfabricated electrodes exists [22] and the integration to photonics
platforms remains an open challenge. The latter will require moving to
silicon-on-insulator (SOI) substrates where photoluminescence measurements
have shown the existence of the bound exciton transition [23] but no
electrical readout has been demonstrated.

One open question when moving towards microfabricated structures on SOI
substrates is the role of strain in the exciton readout. At low temperatures
the different thermal expansion coefficients of silicon and metal cause sharp
strain profiles around the electrodes. This might be especially important when
all the spins are located in the surface layer in an SOI structure (although
see Section VIII for further discussion on this point). Several recent studies
[24, 25, 26] have highlighted how strain can substantially change the donor
spin states and their decoherence properties. Substantially less attention in
this context has so far been paid to strain effects of the hole and exciton
states.

Here we report experiments demonstrating the electro-optical donor exciton
response using microfabricated electrodes and focused light on doped natural
silicon. We have performed experiments on both silicon grown with the
Czochralski method (CZ) and float-zone method (FZ). All data shown is from the
FZ silicon unless otherwise mentioned. The strain from the on-chip electrodes
creates significant effects but importantly does not create any fundamental
obstacles for localized spin readout. We found strain-induced heavy-hole
light-hole splitting to be significant in all samples and we see evidence of
avoided crossings between the heavy-hole and light-hole transitions indicating
a coupling term in the Hamiltonian. We can model these effects with remarkable
agreement assuming the known asymmetry of the hole g-factors and the Pikus-Bir
Hamiltonian describing the strain, but we also have unresolved questions
regarding electric field dependence and light polarization response which is
complicated by the mixing of the different hole states. Nevertheless, using
the extracted parameters we predict with numerical simulations how the exciton
signal should behave in an SOI structure where the strain effects must be
considered with care.

##  II Experimental setup

At low temperatures phosphorous donors in silicon lattice can bind an electron
around them in a Coulomb potential forming a system that resembles a hydrogen
atom but with a Bohr radius of around 1.8 nm [27]. The spin of this bound
electron, in combination with the phosphorous nuclear spin, has been shown to
be a promising quantum computing unit [1, 5]. Unfortunately for optical
applications, the loose binding potential means that the ”atomic” transitions
of this system are in the technologically inconvenient terahertz regime. There
is nevertheless a bound exciton transition at a more technologically
convenient 1078 nm wavelength. The electron of the exciton will form a singlet
pair with the donor electron, leaving the hole 3/2-spin as the spin degree of
freedom in the bound exciton state. Hence, there are six allowed transitions
between the original donor spin state and the exciton state in a finite
magnetic field, these are depicted in Fig. 1(b-c). We neglect the hyperfine
coupling between the donor electron spin and the nuclear spin throughout the
paper as our linewidths are not narrow enough to resolve it. This could be
ameliorated by moving to isotopically pure silicon [17, 6]. As any single
transition originates from a particular donor electron spin state, the
existence or non-existence of these transitions can be used for electron spin
readout. The decay of the exciton happens mostly via Auger recombination,
ending up with an ionized donor and an extra ”hot” electron now in the
conduction band. This causes a change in the conductivity of the silicon,
which we detect.

The measurement scheme is presented in Fig. 1(a). The samples studied are
uniformly doped silicon chips onto which two gold electrodes with length of
100 $\mu$m and separated by 60 $\mu$m (CZ sample) or 100 $\mu$m (FZ sample)
are fabricated using electron beam lithography and liftoff techniques. For the
measurements the samples are placed inside a dry dilution fridge (Bluefors
LD400), which however is mostly operated at 3.3 K temperatures for the data
presented here. The laser (tunable diode laser 1045-1085 nm, Newport TLB-6722)
is guided via free space from outside the cryostat to a fixed lens (NA=0.83)
above the sample chip. The laserspot diameter is roughly 60 $\mu$m in the data
presented, and hence far from the diffraction limit. The sample is glued on a
piezo stage (Attocube ANPxyz101/2) allowing us to align and focus the laser
spot with the readout electrodes using an infrared camera outside the
cryostat. The electrodes connect to a voltage (CZ measurements) or a current
source (FZ measurements). Current and voltage are measured over the sample
with a low noise current amplifier (Femto DLPCA-200), and a low noise voltage
preamplifier (DL instruments model 1201), respectively. The amplified signals
are read out with a digital oscilloscope (Rohde & Schwarz RTB2004). We also
have LEDs (peak wavelength 1020 nm) inside the cryostat to provide light with
energy above the silicon bandgap. The magnetic field is created by a homemade
superconducting magnet in the Helmholtz configuration. The magnetic field
direction is perpendicular to the light field direction, meaning we are using
the Voigt geometry here, allowing us to address all the transitions with
linearly polarized light. The magnet is calibrated with a Hall sensor. We have
a small uncertainty in the location of the sample with regards to the
calibration location, from which we estimate a conservative error range of
$\pm 5$ % for the magnetic field calibration.

Figure 1: a) A schematic of our sample and measurement geometry and the
optical setup. Light is guided to the sample in free space via a fixed lens.
The magnetic field is in-plane and parallel to the direction of measured
current and the [011] crystal direction. b) Energy level diagram showing the
Zeeman splitting of donor electron ($E_{g}$) and hole ($E_{D^{0}X}$) states.
c) Arrows showing the six allowed exciton transitions in a non-zero magnetic
field. d) Simulated current signal of the exciton transitions as a function of
magnetic field in strain-free silicon. Transition brightness corresponds to
the Clebsch-Gordan coefficient of the transition.

Figure 1(b-c) shows schematically the six allowed bound exciton transitions
and Fig.1(d) depicts the signal we would expect to get in different magnetic
fields. The actual measured signal is shown in Fig. 2(c). There are several
main differences which we assign to the effects of strain.

##  III Magnetic field dependency and strain effects

As explained above, the transitions we observe happen essentially between a
single-electron spin-state (D0) and a hole spin-state (D0X). Hence, we need to
write down the Hamiltonians for both the originating and the final state of
the system, while taking strain effects into consideration. We note that we do
not try to model the center frequency of the transition but rather only the
magnetic field dependency. Our center transition frequency at zero field
($\sim$ 1078.02 nm) is somewhat less than what has been reported before [28,
18, 29], we assign this to the substrate material that is not especially free
of impurities.

The magnetic field dependency of the neutral donor state is assumed to simply
follow

$H_{D}^{Z}=\frac{1}{2}\mu_{B}g_{d}\mathbf{B}\cdot\bm{\sigma},$ (1)

where $\mu_{B}$ is the Bohr magneton, $g_{d}=1.9985$ is the donor electron
$g$-factor, $\mathbf{B}$ is the magnetic field vector and $\bm{\sigma}$ is a
vector containing the Pauli spin matrices. We assume isotropic $g$-factor for
the donor electron, although the isotropy might be broken by strain [30, 25].
Nevertheless, the expected anisotropy is of the order of 10-3[31], much
smaller than the strain effects for the holes as we will show below, and can
be neglected here. We also neglect the hyperfine interaction since we cannot
resolve it. There are also strain variation in the electron energy levels
coming from the variation of the hyperfine interaction as a function of strain
but also these are expected to be small [30, 24, 25, 26] compared to the
effects we study here.

For the magnetic field dependence of the hole we use the anisotropic
$g$-factor model [32]

$H^{Z}_{B}=\mu_{B}\left(g_{1}\mathbf{J}\cdot\mathbf{B}+g_{2}\mathbf{J}^{3}\cdot\mathbf{B}\right),$
(2)

where $\mathbf{J}$ is a vector containing 3/2-spin matrices in $x,\>y\>\text{
and }z$ directions (which we take to correspond to the crystal directions of
[100] and equivalent in order to match the strain directions in Eq. (3)
below), $\mathbf{B}$ is the magnetic field vector, and $g_{1}$ and $g_{2}$ are
the isotropic and anisotropic $g$-factors, respectively.

Neglecting the strain effects, that is, only using Hamiltonians described by
Eq. (1) for the donor and Eq. (2) for the hole, we will get exciton transition
lines shown in Fig. 1(d). Comparing this to our actual measurement data shown
in Fig. 2(c), reveals several features not described by the model. The most
prominent discrepancy is the existence of two peaks at zero magnetic field. We
assign this zero field splitting (ZFS) to the splitting of the heavy-hole and
light-hole states caused by strain. The strain effects affecting hole states
in silicon are conventionally modeled with the Pikus-Bir strain Hamiltonian
[33]

$\displaystyle H_{PB}(\bm{\epsilon})=$ $\displaystyle
a\mathrm{Tr}(\bm{\epsilon})+b\sum_{i=x,y,z}\left(J^{2}_{i}-\frac{\mathbf{J}^{2}}{3}\right)\epsilon_{ii}+$
$\displaystyle\frac{d}{\sqrt{3}}\sum_{i\neq
j}(J_{i}J_{j}+J_{j}J_{i})\epsilon_{ij},$ (3)

where $\bm{\epsilon}$ is the strain matrix and $a,b$ and $d$ are deformation
potentials. The parameter $a$ only changes the center frequency of the
transition and is not important for our model here. The other factors $b$ and
$d$ we use as fit parameters. The extracted values are listed in Table 1. For
strain, we use values from a Comsol simulation of the strain caused by the
different thermal coefficients of silicon and the metallic electrodes and the
stage, see section IX.2 for details. Note that as the strain and parameters
$b$ and $d$ are always multiplied, any discrepancy in the strain parameters
will then directly affect these values.

Table 1: Extracted values for deformation potentials $b$ and $d$ along with hole g-factors $g_{1}$ and $g_{2}$ from fitting Eqs. (3) and (2) to the data shown in Fig.2(c). Reference values for $b$ and $d$ are taken from [34]. Reference values for $g_{1}$ and $g_{2}$ from [16]. $g_{LH}$ and $g_{HH}$ were calculated using Eqs. (6–9). |  $b$ (eV) |  $d$ (eV) | $g_{1}$ | $g_{2}$ | $g_{LH}$ | $g_{HH}$
---|---|---|---|---|---|---
This work [011] | -8 | -4 | 0.83 | 0.22 | 1.40 | 1.28
Reference | -2.2 | -5.1 | - | - | - | -
Reference [100] | - | - | 0.8 | 0.24 | 0.86 | 1.34
Reference [111] | - | - | 0.86 | 0.21 | 1.57 | 1.27

By itself the Pikus-Bir strain Hamiltonian simply splits the states at zero-
field. However, the anisotropic component of the Zeeman Hamiltonian ($g_{2}$
in Eq. (2)) causes coupling of the hole states and we get the avoided
crossings between the transitions (Fig. 2(a-b)) which we clearly also see in
the measured data. Using just these equations we can get a remarkable
agreement with the data, as shown in Fig. 2(c). We use the QuTiP Python
package [35, 36] to calculate the transition energies numerically from the
Hamiltonians above.

Notably, we find that the data shows eight transition lines instead of six.
The two extra lines are the outermost and less bright lines. These transitions
are seemingly violating the transition selection rules $\Delta
m=m_{h}-m_{e}=-1,0,1$, where $m_{h}$ and $m_{e}$ are the hole and electron
spins, respectively. However, due to the mixing of the hole states at low
fields (close to the avoided crossings) these transitions are not strictly
forbidden. As the magnetic field increases the forbidden transition lines are
getting fainter, which is the expected result, since the hole states become
less mixed and the spin-values of the hole states are getting better defined
at higher magnetic fields.

Figure 2: a) Strain split hole eigenenergies and b) exciton transition
energies calculated using Eqs. (1), (2) and (3). The magnetic field lies in
the [011] crystal direction. The labels show the dominating state at high
magnetic fields. c) Measured transition energies as a function of the magnetic
field, in the [011] direction, overlaid with the calculated transition
energies presented in (b). The center transition frequency is adjusted to
match the data. Blue lines indicate transitions involving the donor electron
spin up state and red lines transitions coming from donor electron spin down
state. Transitions that are forbidden without the hole mixing are shown with
transparent lines.

From the data in Fig. 2(c) we can also extract the g-factors
($g_{d},g_{1},g_{2}$) defined in Eqs. (1) and (2). The donor g-factor $g_{d}$
anisotropy and strain dependence are well studied and we use the literature
value 1.9985 for all data. The hole $g$-factors are known to have strong
anisotropy between [100] and [111] crystal orientations [16, 37]. Here we
determine the hole $g$-factors in [011] orientation, and fits to the data give
$g_{1}$ and $g_{2}$ parameters shown in Table 1. Our $g_{1}$ and $g_{2}$ are
well in line with the measurements in other crystal directions. Using these
values we can also extract the g-factors for heavy-holes and light-holes (the
slopes of the transitions if there would be no strain induced crossings)
$g_{LH}$ and $g_{HH}$. These values are also shown in Table 1 with known
literature values for other crystal orientations. Details of analytical
calculation of $g_{LH}$ and $g_{HH}$ are shown in Section IX.3. They can also
be extracted numerically.

##  IV Polarization dependence

In all optical transitions the polarization of the absorbed or emitted light
is intimately tied to the spin angular momentum change of the transition. Thus
we would also expect a polarization dependence of the transitions, based on
the change in the spin quantum number $\Delta m$. As can be seen from the
energy level diagram in Fig. 1(c), there are two transitions with $\Delta
m=-1,\,0,\,1$ each. To convert these to the expected light polarization, one
must also consider the relative direction of the light with regards to the
magnetic dipole moment. In our geometry where the magnetic field is
perpendicular to the light travel direction (Voigt geometry), we would expect
to be able to excite the $\Delta_{m}=0$ transitions ($\pi$-transitions) with
one linear polarization component and the $\Delta m=\pm 1$ components
($\sigma$-transitions) with the orthogonal linear polarization.

The hole-mixing changes this picture as the hole spin states are not anymore
well defined and hence also the polarization dependence of the transitions is
relaxed. And indeed we see only a very small polarization dependence in our
signal as can be seen in Fig. 3, where the polarization dependence of the
signal is depicted at 427 mT. We see all seven transitions at all
polarizations even at this field, with only a very slight variation in the
amplitude. We also can still distinguish both forbidden transitions (the
higher energy one being significantly brighter), the appearance of which shows
that there is still a significant amount of hole spin state mixing. (This is
also apparent in our numerical modelling based on the parameters above.)

Figure 3:  (a) Half-waveplate (HWP) angle dependence on the exciton signal at
427 mT. (b) Peak height profiles taken along the dashed lines in (a). The
absence of proper polarization control is due to hole mixing. There is
possibly a small trend visible where the middle transitions ($\pi$-transitions
in absence of hole mixing) get weaker and other transitions
($\sigma$-transitions in absence of hole mixing) get stronger from left to
right as would be expected, but the trend is within the noise of the
measurement.

##  V Above band gap light, laser power dependency and linewidths

We now turn to the more technical details of the hybrid electro-optic readout.
First, for good signal-to-noise ratio it seems to be crucial to provide above
bandgap light. In Fig. 4(a) we show the zero magnetic field data measured at
different optical powers, with and without the LEDs providing the above
bandgap light. It is obvious that the effect of the LEDs is considerable in
increasing the signal. This is due to the extra electrons in the conduction
band allowing for faster charge neutralization rate of the donors, increasing
the ”recycling” rate of electrons and hence the measured current, as has also
been reported before [18, 19, 22].

Figure 4: a) Laser power dependence of the signal at zero magnetic field.
Current values are offset to show the difference in the peak heights. The base
current difference between LED on and LED off measurement is $\approx$3.5 nA.
b) Exciton transition at zero applied magnetic field for CZ and FZ silicon.
Current values are normalized and CZ data is offset to illustrate the
difference in the linewidths. CZ data was captured without LED illumination
and FZ data with LEDs. The laser power was the same in both measurements.

Interestingly, the shape of the zero-field data changes when the above bandgap
light is either on or off. When the LEDs are on the zero-field data does not
have a dependency on the laser power, as shown in Fig. 4(a) showing a constant
ZFS of roughly 23 $\mu$eV. However, when the LEDs are off, we see a
significant change in the shape of the data. Naively fitting still a two-peak
function to the LEDs off data would lead to a ZFS to roughly 15 $\mu$eV. If
this would imply a ZFS change due to heating from the LEDs changing strain
sample would be heating up to 100 K according to Comsol simulations, which is
not possible as at these temperatures neither the donor electrons nor excitons
would be bound.

Figure 4(b) shows comparison data at zero applied magnetic field in CZ and FZ
samples. As is clearly visible the linewidth improves considerably when moving
to FZ-silicon, as is expected from the lower concentration of oxygen and
carbon impurities [28]. In the CZ sample the full-width-half-maximums (FWHM)
are roughly 25 $\mu$eV whereas in FZ we see linewidth of 8 $\mu$eV for the
heavy-hole (lower energy) peak and 5 $\mu$eV for the light-hole (higher
energy) peak. When fitting to data acquired at a high magnetic field we find
FWHMs ranging 3–5 $\mu$eV, which is well in line with values reported before
for natural silicon [29].

##  VI Millikelvin measurements

In addition to 3.3 K measurements, we also performed measurements at
millikelvin temperatures. However, we did not see a significant difference to
the 3.3 K data when high laser power of 30 $\mu$W was used, as shown in Fig.
5(a) and (c) (sample stage thermometer temperature 115 mK). We assume that the
heating from the laser power was enough to basically keep the sample at the
same temperature as before (it also raised the fridge temperature). When the
laser power was reduced to one tenth, the signal became very weak and we
observe a very different pattern, see Fig. 5(b) and (d) (sample stage
thermometer temperature 45 mK). LEDs providing above bandgap light were on in
both measurements. We cannot explain the low power data with strain effects,
since the thermal contraction at low temperatures should be negligible.
Somewhat similar ’additional’ zero field splitting has been observed before at
1.4 K temperatures [18] and was attributed to other charge centers interacting
with the spin state.

Figure 5: Measured magnetic field dependence on the spectrum at millikelvin
temperatures with a) 30 $\mu$W (sample stage temperature 115 mK) and b) 3
$\mu$W of laser power (sample stage temperature 45 mK). Figure (c) and (d)
show crosscuts at zero magnetic field.

##  VII Stark shift

We also performed measurements sweeping the biasing voltage at zero magnetic
field in order to characterize the Stark shift in our devices in both 3.3 K
and in 50 mK temperatures. Peculiarly, in the CZ device we found a quadratic
Stark shift of the transition signal at millikelvin temperatures and a linear
Stark shift at 3.3 K. In addition, at millikelvin temperatures we saw the
higher energy peak only become visible after certain applied electric field.
Unfortunately, in the FZ sample we only managed to see a reliable signal at
3.3 K temperature and confirmed a similar linear shift there as was measured
in CZ samples.

In order to quantify the shift $\Delta E$, we fit the data using

$\frac{\Delta E(\bf{F})}{E(0)}=\eta|\bf{F}|^{2},$ (4)

for the quadratic shift and

$\frac{\Delta E(\bf{F})}{E(0)}=\xi|\bf{F}|,$ (5)

for the linear case. The $\bf{F}$ is the applied external electric field and
$\eta$ and $\xi$ are the fitting parameters.

In Fig. 6(a) and (b) we show the data from the CZ grown sample measured in 3.3
K and in mK conditions. The blue curves are the best fits. The 3.3 K data
shows a linear Stark shift and we find a good fit with $\xi=-1\cdot 10^{-4}$
$\mu$m/V. For the mK data the fit yields a quadratic Stark shift parameter
$\eta=7\cdot 10^{-3}$ $\mu\text{m}^{2}/\text{V}^{2}$. For context this can be
compared with the value $-2.65\cdot 10^{-3}$ $\mu\text{m}^{2}/\text{V}^{2}$
calculated for donor electron hyperfine shift [38]. Using a simple hydrogen
atom model, our measured value would imply a Bohr radius of 8 nm.

Figure 6: a) An electric field sweep of the CZ grown sample at 3.3 K and b) at
100 mK. Linear Stark shift was fitted to a) where blue line follows Eq. (5)
and $\xi=-1\cdot 10^{-4}$ $\mu$m/V. For b) we fitted quadratic Stark shift
presented in Eq. (4) where $\eta=7\cdot 10^{-3}$
$\mu\text{m}^{2}/\text{V}^{2}$. Panel c) is the same measurement performed on
FZ grown silicon at 3.3 K. For c) we find $\xi=1.5\cdot 10^{-4}$ $\mu$m/V. We
were not able to capture comparable millikelvin data from the FZ sample.

Figure 6(c) shows similar data measured from the FZ sample, yielding a similar
result to the CZ sample at 3.3 K, namely a linear Stark shift with
$\xi=1.5\cdot 10^{-4}$ $\mu$m/V. At high biases we see artefacts that require
further investigations. It is however interesting to note that the data
exhibits a feature where the ZFS gets considerably smaller. Hence, we
speculate that it could be possible that one could use electric fields to
partially overcome the effects of strain induced ZFS.

##  VIII Discussion and future outlook

We have characterized the phosphorous bound exciton transitions in silicon
using hybrid electro-optical readout. The inevitable strain coming from the
electrodes needed for this readout method causes hole mixing which both
complicates the analysis of the transitions and prevents polarization
dependent addressing. Importantly, it does not however fundamentally prevent
electron spin readout.

The strain can be seen as a non-wanted property, when dealing with sensitive
exciton transitions, but it could also be useful for tuning them. However the
tuning needs to take place already in the design of the system. Strain
engineering in micro-electro-mechanical-systems is already widely studied
topic [39], so the methods for the strain tuning are already available. This
could be useful for e.g. purposefully matching the transition frequency for
photonic components.

We have also directly measured the heavy-hole and light-hole $g$-factors in
[011] orientation, and our results confirm the expectation from the $g_{1}$
and $g_{2}$ measured before. Indeed, with only slight adjusting of the
simulated strain parameters we get a remarkable agreement between theory and
experiment for the magnetic field dependency using the Pikus-Bir strain
Hamiltonian. It is notable that the $g$-factor anisotropy of holes also allows
one to tune the exciton transition position just by rotating the sample, which
could be done in-situ with a rotating stage.

For widespread applications, it would be beneficial to integrate the hybrid
readout with silicon photonics components. This would require integrating the
readout electronics with SOI devices, where strain effects can be increased
due to the thin device layer. On the other hand, in a SOI substrate it is easy
to suspend the device layer. If one then places the readout electronics on the
non-suspended part, the strain caused by the electrodes around the donors is
actually decreased compared to the bulk case (see strain simulations in
Section IX.2). More importantly, with our simulation parameters all the strain
components have completely flat profiles in the suspended parts, since we
leave 10 $\mu$m space between the electrode and the suspended Si film. This
uniformity in the strain field will be crucial for any larger scale
architecture in order to avoid inhomogenous broadning of the exciton
transition.

If in addition isotopically pure silicon is used, the ensemble linewidths
could be extremely narrow. Previously it has been reported that the exciton
transitions in 99.991 % isotopically enriched 28Si have a FWHM of 150 neV
[18]. We show expected (simulated) signals from the SOI samples made with
natural silicon and enriched 28Si in Fig. 7 using the measured parameters and
strain values from the SOI simulation (assuming strain free substrate at room
temperature). From these it is expected that hybrid electro-optical spin
readout should be possible in SOI platforms at low magnetic fields.

Figure 7: Simulated exciton transitions on SOI platform in a) natural silicon
and in b) 99.991% enriched 28Si. Simulation parameters are presented in Table
1 and strain is taken from simulations shown in Sec. IX.2. Transition
linewidth in a) is the bulk value we have extracted from natural silicon 4
$\mu$eV, and in b) 150 neV as reported in [18]. The white rings in a) show
possible spin readout conditions in natural silicon at just below 200 mT,
where the transition is isolated enough to be addressed optically, and
similarly for the white line in (b) at 20 mT.

##  IX Methods

###  IX.1 Samples

The CZ wafer used was a 380 $\mu$m thick phosphorous doped wafer acquired from
Okmetic Oy with crystal direction (100). The wafer resistivity range at room
temperature was specified to be 0.2-0.25 $\Omega\,\text{cm}$, meaning doping
levels of 2.4-3.2e16 $\text{cm}^{-3}$. The FZ wafer used was a 200 $\mu$m
thick phosphorous doped wafer acquired from Sil’Tronix ST with crystal
direction (100).The wafer resistivity range at room temperature was specified
to be 1-5 $\Omega\,\text{cm}$, meaning doping levels of 0.9-4.9e15
$\text{cm}^{-3}$.

The electrodes were fabricated using electron beam lithography with PMMA
resist and UHV electron beam evaporation. The electrodes consist of a 5 nm
titanium layer for adhesion to silicon surface and 50 nm of gold.

###  IX.2 Strain simulations

The strain in our samples was estimated by simulating the thermal expansion
mismatches in Comsol Multiphysics software. In the simulation a silicon chip
rests on a copper stage and has gold electrodes on the top. We assume zero-
strain at room temperature and then calculate the induced strain (by both the
electrodes and the stage) at low temperatures. We made custom interpolation
tables for the linear thermal expansion coefficients [40, 41, 42], since at
low temperatures the coefficients are considerably different compared to the
room temperature case. The temperature is swept down to 1 K and the strain is
calculated for each temperature step. An anisotropic model and parameters were
used for silicon. Then all strain components were plotted on the surface (top
part in Fig. 8) and at 10 $\mu$m depth (bottom part in Fig. 8). In calculating
the Pikus-Bir Hamiltonian Eq. (3) we use an average of each strain component
along a line connecting the electrodes at 10 $\mu$m depth.

Figure 8: Simulated strain tensor elements on the sample surface (top) and at
10 $\mu$m depth (bottom) in the bulk sample. The dark blue colored area in the
panel showing $\epsilon_{11}$ is the gold electrode. The black arrows are
showing the coordinate directions. Figure 9: Simulated strain tensor elements
at 100 nm depth in the suspended 220 nm silicon film. The dark blue colored
area is under the gold electrode in the panel showing $\epsilon_{11}$. The
black arrows are showing the coordinate directions. The stripes between the
electrodes are narrow holes that allow the release of the Si film by HF
etching.

We repeated the simulation also for SOI material that has a 220 nm natural Si
film, 3 $\mu$m SiO2 layer and 750 $\mu$m Si substrate. In the simulations we
suspend the area between the electrodes, leaving a $10$ $\mu$m gap between the
suspension and the electrode edge. This allows the strain, created by the
electrodes, to relax before the suspended silicon film starts and considerably
both lessens the strain and relaxes any gradients in the suspended area as can
be seen in the simulated strain components at 100 nm depth that are shown in
Fig. 9.

###  IX.3 Extracting hole g-factors

The Zeeman energy splitting of a 3/2-hole can be calculated from [33]

$\displaystyle E_{1,2}=$
$\displaystyle\pm\mu_{0}\left\\{\frac{\mathbf{B}^{2}}{8}\left[9\left(g_{1}+\frac{9}{4}g_{2}\right)^{2}+\left(g_{1}+\frac{g_{2}}{4}\right)^{2}\right]+\right.$
$\displaystyle\left.\left(g_{1}+\frac{7}{4}g_{2}\right)\left[\left(g_{1}+\frac{13}{4}g_{2}\right)^{2}\mathbf{B}^{4}-\right.\right.$
$\displaystyle\left.\left.9g_{2}\left(g_{1}+\frac{5}{2}g_{2}\right)\left(B^{2}_{x}B^{2}_{y}+B^{2}_{x}B^{2}_{z}+B^{2}_{y}B^{2}_{z}\right)\right]^{\frac{1}{2}}\right\\}^{\frac{1}{2}}$
(6)

and

$\displaystyle E_{3,4}=$
$\displaystyle\pm\mu_{0}\left\\{\frac{\mathbf{B}^{2}}{8}\left[9\left(g_{1}+\frac{9}{4}g_{2}\right)^{2}+\left(g_{1}+\frac{g_{2}}{4}\right)^{2}\right]-\right.$
$\displaystyle\left.\left(g_{1}+\frac{7}{4}g_{2}\right)\left[\left(g_{1}+\frac{13}{4}g_{2}\right)^{2}\mathbf{B}^{4}-\right.\right.$
$\displaystyle\left.\left.9g_{2}\left(g_{1}+\frac{5}{2}g_{2}\right)\left(B^{2}_{x}B^{2}_{y}+B^{2}_{x}B^{2}_{z}+B^{2}_{y}B^{2}_{z}\right)\right]^{\frac{1}{2}}\right\\}^{\frac{1}{2}}$
(7)

where $E_{1,2}$ and $E_{3,4}$ give the heavy-hole and light-hole splitting
energies, respectively. By setting the magnetic field direction using notation
$[xyz]$, referring to the crystal orientation, one can calculate the effective
$g$-factors along any direction. The heavy hole and light hole $g$-factors,
that are observed at the high field asymptote, can be solved from

$\displaystyle\frac{1}{2}\mu_{0}g_{LH}\mathbf{B}=E_{3,4}$ (8)
$\displaystyle\frac{3}{2}\mu_{0}g_{HH}\mathbf{B}=E_{1,2}.$ (9)

###### Acknowledgements.

We acknowledge useful discussions with Jarryd Pla and Joe Salfi, and Henri
Lyyra for proofreading the manuscript. This project has received funding from
the European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No 852428), from Academy of
Finland Grant No 321416 and from Jenny and Antti Wihuri Foundation.

## References

  * Kane [1998] B. E. Kane, A silicon-based nuclear spin quantum computer, Nature 393, 133 (1998).
  * Morton  _et al._ [2008] J. J. L. Morton, A. M. Tyryshkin, R. M. Brown, S. Shankar, B. W. Lovett, A. Ardavan, T. Schenkel, E. E. Haller, J. W. Ager, and S. A. Lyon, Solid-state quantum memory using the 31p nuclear spin, Nature 455, 1085 (2008).
  * Morello  _et al._ [2010] A. Morello, J. J. Pla, F. A. Zwanenburg, K. W. Chan, K. Y. Tan, H. Huebl, M. Möttönen, C. D. Nugroho, C. Yang, J. A. van Donkelaar, A. D. C. Alves, D. N. Jamieson, C. C. Escott, L. C. L. Hollenberg, R. G. Clark, and A. S. Dzurak, Single-shot readout of an electron spin in silicon, Nature 467, 687 (2010).
  * Tyryshkin  _et al._ [2012] A. M. Tyryshkin, S. Tojo, J. J. L. Morton, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, T. Schenkel, M. L. W. Thewalt, K. M. Itoh, and S. A. Lyon, Electron spin coherence exceeding seconds in high-purity silicon, Nature Materials 11, 143 (2012).
  * Pla  _et al._ [2012] J. J. Pla, K. Y. Tan, J. P. Dehollain, W. H. Lim, J. J. L. Morton, D. N. Jamieson, A. S. Dzurak, and A. Morello, A single-atom electron spin qubit in silicon, Nature 489, 541 (2012).
  * Steger  _et al._ [2012] M. Steger, K. Saeedi, W. Thewalt M. L., L. Morton J. J., H. Riemann, V. Abrosimov N., P. Becker, and H.-J. Pohl, Quantum information storage for over 180 s using donor spins in a 28si “semiconductor vacuum”, Science 336, 1280 (2012).
  * Saeedi  _et al._ [2013] K. Saeedi, S. Simmons, J. Salvail, P. Dluhy, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, J. J. Morton, and M. L. Thewalt, Room-temperature quantum bit storage exceeding 39 minutes using ionized donors in silicon-28, Science 342, 830 (2013).
  * Muhonen  _et al._ [2014] J. T. Muhonen, J. P. Dehollain, A. Laucht, F. E. Hudson, R. Kalra, T. Sekiguchi, K. M. Itoh, D. N. Jamieson, J. C. McCallum, A. S. Dzurak, and A. Morello, Storing quantum information for 30 seconds in a nanoelectronic device, Nature Nanotechnology 9, 986 (2014).
  * He  _et al._ [2019] Y. He, S. K. Gorman, D. Keith, L. Kranz, J. G. Keizer, and M. Y. Simmons, A two-qubit gate between phosphorus donor electrons in silicon, Nature 571, 371 (2019).
  * Madzik  _et al._ [2022] M. T. Madzik, S. Asaad, A. Youssry, B. Joecker, K. M. Rudinger, E. Nielsen, K. C. Young, T. J. Proctor, A. D. Baczewski, A. Laucht, V. Schmitt, F. E. Hudson, K. M. Itoh, A. M. Jakob, B. C. Johnson, D. N. Jamieson, A. S. Dzurak, C. Ferrie, R. Blume-Kohout, and A. Morello, Precision tomography of a three-qubit donor quantum processor in silicon, Nature 601, 348 (2022).
  * Songtao  _et al._ [2020] C. Songtao, R. Mouktik, C. M. Phenicie, S. Ourari, and J. D. Thompson, Parallel single-shot measurement and coherent control of solid-state spins below the diffraction limit, Science 370, 592 (2020).
  * Hollenbach  _et al._ [2020] M. Hollenbach, Y. Berencen, U. Kentsch, M. Helm, and G. V. Astakhov, Engineering telecom single-photon emitters in silicon for scalable quantum photonics, Opt. Express 28, 26111 (2020).
  * Bergeron  _et al._ [2020] L. Bergeron, C. Chartrand, A. T. K. Kurkjian, K. J. Morse, H. Riemann, N. V. Abrosimov, P. Becker, H.-J. Pohl, M. L. W. Thewalt, and S. Simmons, Silicon-integrated telecommunications photon-spin interface, PRX Quantum 1, 020301 (2020).
  * Redjem  _et al._ [2020] W. Redjem, A. Durand, T. Herzig, A. Benali, S. Pezzagna, J. Meijer, A. Y. Kuznetsov, H. S. Nguyen, S. Cueff, J.-M. Gérard, I. Robert-Philip, B. Gil, D. Caliste, P. Pochet, M. Abbarchi, V. Jacques, A. Dréau, and G. Cassabois, Single artificial atoms in silicon emitting at telecom wavelengths, Nature Electronics 3, 738 (2020).
  * Durand  _et al._ [2021] A. Durand, Y. Baron, W. Redjem, T. Herzig, A. Benali, S. Pezzagna, J. Meijer, A. Y. Kuznetsov, J.-M. Gérard, I. Robert-Philip, M. Abbarchi, V. Jacques, G. Cassabois, and A. Dréau, Broad diversity of near-infrared single-photon emitters in silicon, Physical Review Letters 126, 083602 (2021).
  * Kaminskil  _et al._ [1980] A. S. Kaminskil, V. A. Karasyuk, and Y. E. Pokrovskil, Luminescence of excitons bound to phosphorus atoms insilicon subjected to a magnetic field, Sov. Phys. JETP 52, 211 (1980).
  * Karaiskaj  _et al._ [2001] D. Karaiskaj, M. L. W. Thewalt, T. Ruf, M. Cardona, H.-J. Pohl, G. G. Deviatych, P. G. Sennikov, and H. Riemann, Photoluminescence of isotopically purified silicon: How sharp are bound exciton transitions?, Physical Review Letters 86, 6010 (2001).
  * Yang  _et al._ [2006] A. Yang, M. Steger, D. Karaiskaj, M. L. W. Thewalt, M. Cardona, K. M. Itoh, H. Riemann, N. V. Abrosimov, M. F. Churbanov, A. V. Gusev, A. D. Bulanov, A. K. Kaliteevskii, O. N. Godisov, P. Becker, H.-J. Pohl, J. W. Ager, and E. E. Haller, Optical detection and ionization of donors in specific electronic and nuclear spin states, Physical Review Letters 97, 227401 (2006).
  * Steger  _et al._ [2011] M. Steger, T. Sekiguchi, A. Yang, K. Saeedi, M. E. Hayden, M. L. W. Thewalt, K. M. Itoh, H. Riemann, N. V. Abrosimov, P. Becker, and H.-J. Pohl, Optically-detected nmr of optically-hyperpolarized 31p neutral donors in 28si, Journal of Applied Physics 109, 102411 (2011).
  * Gullans and Taylor [2015] M. J. Gullans and J. M. Taylor, Optical control of donor spin qubits in silicon, Physical Review B 92, 195411 (2015).
  * Ross  _et al._ [2019] P. Ross, B. C. Rose, C. C. Lo, M. L. Thewalt, A. M. Tyryshkin, S. A. Lyon, and J. J. Morton, Electron spin resonance of p donors in isotopically purified $\mathrm{Si}$ detected by contactless photoconductivity, Phys. Rev. Applied 11, 054014 (2019).
  * Lo  _et al._ [2015] C. C. Lo, M. Urdampilleta, P. Ross, M. F. Gonzalez-Zalba, J. Mansir, S. A. Lyon, M. Thewalt, and J. J. L. Morton, Hybrid optical-electrical detection of donor electron spins with bound excitons in silicon, Nature Materials 14, 490 (2015).
  * Sumikura  _et al._ [2011] H. Sumikura, K. Nishiguchi, Y. Ono, A. Fujiwara, and M. Notomi, Bound exciton photoluminescence from ion‑implanted phosphorus in thin silicon layers, Opt. Express 19, 25255 (2011).
  * Pla  _et al._ [2018] J. J. Pla, A. Bienfait, G. Pica, J. Mansir, F. A. Mohiyaddin, Z. Zeng, Y. M. Niquet, A. Morello, T. Schenkel, J. J. L. Morton, and P. Bertet, Strain-induced spin-resonance shifts in silicon devices, Phys. Rev. Applied 9, 044014 (2018).
  * Mansir  _et al._ [2018] J. Mansir, P. Conti, Z. Zeng, J. J. Pla, P. Bertet, M. W. Swift, C. G. Van de Walle, M. L. W. Thewalt, B. Sklenard, Y. M. Niquet, and J. J. L. Morton, Linear hyperfine tuning of donor spins in silicon using hydrostatic strain, Phys. Rev. Lett. 120, 167701 (2018).
  * Ranjan  _et al._ [2021] V. Ranjan, B. Albanese, E. Albertinale, E. Billaud, D. Flanigan, J. J. Pla, T. Schenkel, D. Vion, D. Esteve, E. Flurin, J. J. L. Morton, Y. M. Niquet, and P. Bertet, Spatially resolved decoherence of donor spins in silicon strained by a metallic electrode, Phys. Rev. X 11, 031036 (2021).
  * Smith  _et al._ [2017] J. S. Smith, A. Budi, M. C. Per, N. Vogt, D. W. Drumm, L. C. L. Hollenberg, J. H. Cole, and S. P. Russo, Ab initio calculation of energy levels for phosphorus donors in silicon, Scientific Reports 7, 10.1038/s41598-017-06296-8 (2017).
  * Safonov  _et al._ [1996] A. N. Safonov, G. Davies, and E. C. Lightowlers, Line shape of the no-phonon luminescence of excitons bound to phosphorus in carbon-doped silicon, Phys. Rev. B 54, 4409 (1996).
  * Yang  _et al._ [2009] A. Yang, M. Steger, T. Sekiguchi, M. L. W. Thewalt, J. W. Ager, and E. E. Haller, Homogeneous linewidth of the p31 bound exciton transition in silicon, Applied Physics Letters 95, 122113 (2009).
  * Wilson and Feher [1961] D. K. Wilson and G. Feher, Electron spin resonance experiments on donors in silicon. iii. investigation of excited states by the application of uniaxial stress and their importance in relaxation processes, Phys. Rev. 124, 1068 (1961).
  * Usman  _et al._ [2018] M. Usman, H. Huebl, A. R. Stegner, C. D. Hill, M. S. Brandt, and L. C. L. Hollenberg, Measurements and atomistic theory of electron $g$-factor anisotropy for phosphorus donors in strained silicon, Phys. Rev. B 98, 035432 (2018).
  * Bir  _et al._ [1963] G. Bir, E. Butikov, and G. Pikus, Spin and combined resonance on acceptor centres in ge and si type crystals—i, Journal of Physics and Chemistry of Solids 24, 1467 (1963).
  * Bir and Pikus [1974] G. Bir and G. Pikus,  _Symmetry and Strain-induced Effects in Semiconductors_, A Halsted Press book (Wiley, 1974).
  * Blacha  _et al._ [1984] A. Blacha, H. Presting, and M. Cardona, Deformation potentials ofk = 0 states of tetrahedral semiconductors, physica status solidi (b) 126, 11 (1984).
  * Johansson  _et al._ [2012] J. Johansson, P. Nation, and F. Nori, Qutip: An open-source python framework for the dynamics of open quantum systems, Computer Physics Communications 183, 1760 (2012).
  * Johansson  _et al._ [2013] J. Johansson, P. Nation, and F. Nori, Qutip 2: A python framework for the dynamics of open quantum systems, Computer Physics Communications 184, 1234 (2013).
  * Beck [2019] M. Beck, Optical spectroscopy of donor bound excitons and spin relaxation of donor electrons in isotopically enriched silicon (2019).
  * Ma  _et al._ [2022] H. Ma, Y.-L. Hsueh, S. Monir, Y. Jiang, and R. Rahman, Ab-initio calculations of shallow dopant qubits in silicon from pseudopotential and all-electron mixed approach, Communications Physics 5, 165 (2022).
  * Shchepetov  _et al._ [2013] A. Shchepetov, M. Prunnila, F. Alzina, L. Schneider, J. Cuffe, H. Jiang, E. I. Kauppinen, C. M. S. Torres, and J. Ahopelto, Ultra-thin free-standing single crystalline silicon membranes with strain control, Applied Physics Letters 102, 192108 (2013).
  * Rodriguez  _et al._ [2018] J. Rodriguez, S. A. Chandorkar, G. M. Glaze, D. D. Gerrard, Y. Chen, D. B. Heinz, I. B. Flader, and T. W. Kenny, Direct detection of anchor damping in MEMS tuning fork resonators, Journal of Microelectromechanical Systems 27, 800 (2018).
  * Pamato  _et al._ [2018] M. G. Pamato, I. G. Wood, D. P. Dobson, S. A. Hunt, and L. Vočadlo, The thermal expansion of gold: point defect concentrations and pre-melting in a face-centred cubic metal, Journal of Applied Crystallography 51, 470 (2018).
  * White [1964] G. White, Thermal expansion of silica at low temperatures, Cryogenics 4, 2 (1964).


# External Inverse-Compton and Proton Synchrotron Emission from the Reverse
Shock as the Origin of VHE Gamma-Rays from the Hyper-Bright GRB 221009A

B. Theodore Zhang (张兵)11affiliation: Center for Gravitational Physics and
Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto, Kyoto 606-8502, Japan Kohta Murase22affiliation: Department
of Physics, The Pennsylvania State University, University Park, PA 16802, USA
33affiliation: Department of Astronomy & Astrophysics, The Pennsylvania State
University, University Park, PA 16802, USA 44affiliation: Center for
Multimessenger Astrophysics, Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA 55affiliation:
School of Natural Sciences, Institute for Advanced Study, Princeton, NJ 08540,
USA 11affiliation: Center for Gravitational Physics and Quantum Information,
Yukawa Institute for Theoretical Physics, Kyoto University, Kyoto, Kyoto
606-8502, Japan Kunihito Ioka11affiliation: Center for Gravitational Physics
and Quantum Information, Yukawa Institute for Theoretical Physics, Kyoto
University, Kyoto, Kyoto 606-8502, Japan Deheng Song11affiliation: Center for
Gravitational Physics and Quantum Information, Yukawa Institute for
Theoretical Physics, Kyoto University, Kyoto, Kyoto 606-8502, Japan Chengchao
Yuan (袁成超)66affiliation: Deutsches Elektronen-Synchrotron DESY, Platanenallee
6, 15738 Zeuthen, Germany Péter Mészáros22affiliation: Department of Physics,
The Pennsylvania State University, University Park, PA 16802, USA
33affiliation: Department of Astronomy & Astrophysics, The Pennsylvania State
University, University Park, PA 16802, USA 44affiliation: Center for
Multimessenger Astrophysics, Institute for Gravitation and the Cosmos, The
Pennsylvania State University, University Park, PA 16802, USA

###### Abstract

The detection of the hyper-bright gamma-ray burst (GRB) 221009A enables us to
explore the nature of the GRB emission and the origin of very-high-energy
(VHE) gamma-rays. We analyze the Fermi-LAT data of this burst and investigate
the GeV-TeV emission in the framework of the external reverse shock model. We
show that the early $\sim 1-10$ GeV emission can be explained by the external
inverse-Compton mechanism via upscattering MeV gamma-rays by electrons
accelerated at the reverse shock, in addition to the synchrotron self-Compton
component. The predicted early optical flux could have been brighter than that
of the naked-eye GRB 080319B. We also show that proton synchrotron emission
from accelerated ultra-high-energy cosmic rays (UHECRs) is detectable, and
could potentially explain $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\rm TeV$ photons detected by LHAASO or
constrain the UHECR acceleration mechanism. Our model suggests that the
detection of $\mathcal{O}(10\rm~{}TeV)$ photons with energies up to $\sim 18$
TeV is possible for reasonable models of the extragalactic background light
without invoking new physics, and predicts anti-correlations between MeV
photons and TeV photons, which can be tested with the LHAASO data.

##  1 Introduction

Gamma-ray bursts (GRBs) are among the most luminous explosions in the Universe
(Mészáros, 2006; Kumar & Zhang, 2014). In 2019, the detection of two TeV
bursts, GRB 190114C (MAGIC Collaboration, 2019; MAGIC Collaboration et al.,
2019) and GRB 180720B (Abdalla et al., 2019), has opened a new window in the
VHE ($\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}0.1\rm~{}TeV$) band for studying GRBs,
providing us with new opportunities to investigate the nature of GRBs (see
Miceli & Nava, 2022; Gill & Granot, 2022, for reviews).

On October 9 2022, GRB 221009A was triggered by the Fermi Gamma-Ray Burst
Monitor (GBM) at $T_{0}$ = 13:16:59.99 UT (Veres et al., 2022). The Swift
Burst Alert Telescope (BAT) also triggered GRB 221009A around one-hour later
(Dichiara et al., 2022). GRB 221009A is an extraordinarily bright and
energetic GRB with isotropic-equivalent energy $\mathcal{E}_{k}\sim 3\times
10^{55}\rm~{}erg$ for a radiative efficiency of 10% (Frederiks et al., 2022)
at a redshift $z=0.15$ (de Ugarte Postigo et al., 2022). The Fermi-LAT
reported the detection of $>100\rm~{}MeV$ gamma-rays with the maximum photon
energy reaching 99.3 GeV (Pillera et al., 2022). Remarkably, it was reported
by the Large High Altitude Air Shower Observatory (LHAASO) that there are more
than 5000 gamma-rays with energy beyond 500 GeV from GRB 221009A detected, and
the highest-energy gamma-ray energy reaches 18 TeV (Huang et al., 2022).

The production of VHE gamma-rays from GRB up to $\sim\rm TeV$ has been widely
discussed in the standard afterglow model via synchrotron self-Compton (SSC)
process (e.g., Mészáros & Rees, 1994; Dermer et al., 2000; Sari & Esin, 2001;
Zhang & Mészáros, 2001; Ren et al., 2022) or external inverse-Compton (EIC)
process (e.g., Wang et al., 2006; Murase et al., 2010; Toma et al., 2011;
Murase et al., 2011; He et al., 2012; Veres & Mészáros, 2014; Kimura et al.,
2019; Zhang et al., 2021a, b). Proton synchrotron emission has been proposed
as one alternative mechanism to generate VHE gamma-rays from GRBs (e.g.,
Totani, 1998; Zhang & Mészáros, 2001; Murase et al., 2008; Asano et al., 2009;
Isravel et al., 2022), which usually requires protons to be accelerated to the
ultra-high-energy (UHE) range. The proton synchrotron emission has an
advantage in the generation of $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}10\rm~{}TeV$ gamma-rays, which is usually
difficult for SSC and EIC processes due to the limitation of the Klein-Nishina
effect. Motivated by the detection of $\mathcal{O}(10\rm~{}TeV)$ gamma-rays
from GRB 221009A by the LHAASO observatory, it is reasonable to study the
proton synchrotron radiation process of GRB 221009A in detail.

It has been proposed that the acceleration of UHECRs is possible in the
internal shock model (e.g., Waxman, 1995; Vietri, 1996) or reverse shock model
of GRBs (e.g., Waxman & Bahcall, 2000; Murase, 2007; Murase et al., 2008;
Zhang et al., 2018), while the acceleration of UHECRs in the forward shock
region via the diffusive shock acceleration mechanism is difficult due to the
low magnetic field strength of the external medium (Gallant & Achterberg,
1999; Murase et al., 2008; Sironi et al., 2015). GRB internal shocks may occur
at much smaller radii, where the escape of $\mathcal{O}(10\rm~{}TeV)$ gamma-
rays is difficult unless the Lorentz factor is extremely large (e.g., Murase
et al., 2022). One of the possibilities is the UHECR proton synchrotron
emission process in the reverse shock model. The observed prompt emission in
the MeV band indicates that GRB 221009A is a long-lasting GRB with $T_{\rm
90}$ of at least $\sim 600$ seconds, implying a thick ejecta shell case where
the long-lasting reverse shock can be expected. Also, during the early phase
of the reverse shock, the extraordinarily bright prompt photons can be
upscattered by the high-energy electrons accelerated in the reverse shock
region to produce high-energy gamma-rays which may be related to Fermi-LAT
observations (e.g., Beloborodov, 2005).

Throughout this work, we use $Q_{x}/Q=10^{x}$ in CGS units.

##  2 GeV-TeV gamma-rays from a reverse shock model

Figure 1: The flux light curve observed by the Fermi-LAT between $T_{0}$ \+
203 s and $T_{0}$ \+ 1000 s. The first time interval overlaps with the LAT Bad
Time Interval due to photon pile-up and therefore is only included for
reference purposes. The black-dashed lines indicate the three peak times of
the prompt emission observed by the Konus-Wind.

We analyze the Fermi-LAT data collected between 203 s and 1000 s after the
Fermi GBM trigger time of $T_{0}$ = 13:16:59.99 UT (Veres et al., 2022;
Pillera et al., 2022). Our region of interest (ROI) is a 10∘ by 10∘ region
centered around the location of the GRB as reported by Pillera et al. 2022 (RA
= 288.28∘, DEC = 19.49∘). We select the events from the P8R3_TRANSIENT020
class with energies between 0.1 GeV and 10 GeV. We consider 3 time intervals
after $T_{0}$: 203 s - 294 s, 294 s - 410 s, and 410 s - 1000 s. In each time
interval, we perform unbinned likelihood analysis by varying a point source
placed at the center of the ROI. We also include the Galactic interstellar
emission model gll_iem_v07.fits and the isotropic diffuse emission model
iso_P8R3_TRANSIENT020_V3_v1.txt, and their normalizations are allowed to vary
freely. The first time interval overlaps with the LAT Bad Time Interval due to
photon pile-up (Omodei et al., 2022a, b). We only include it for reference
purposes. The LAT photons were detected in the time interval $T_{0}$ \+ 294s
to $T_{0}$ \+ 410 s, and the integrated photon flux between 0.1 GeV and 10 GeV
is $(4.0\pm 0.4)\times 10^{-3}$ ph cm-2 s-1 with a photon index of $1.76\pm
0.09$. We do not find photons from the direction of the GRB during the time
interval $T_{0}$ \+ 410 s to $T_{0}$ \+ 1000 s. The 95% upper limit of the
flux is $3.5\times 10^{-3}$ ph cm-2 s-1. The flux light curve is shown in Fig.
1.

In Fig. 1, we also show the peak time of the observed pulse of GRB prompt
emission detected by Konus-Wind at 13:17:01.648 UT (Frederiks et al., 2022).
We can see there was active prompt emission during the Fermi-LAT observations
within the time interval $T-T_{0}\sim 200-400\rm~{}s$, while the prompt
emission becomes weaker at later times up to $\sim 600\rm~{}s$, which may also
be regarded as the flare phase (Frederiks et al., 2023; An et al., 2023).
Considering the long-lasting prompt emission of GRB 221009A, in this work, we
assume that the GRB ejecta is thick with duration time $600\rm~{}s$ (e.g.,
Ioka et al., 2005).

We consider two episodes, where the first episode (Episode I) is strongly
affected by the active prompt emission phase, while the second episode
(Episode II) is much less affected by the prompt emission until the reverse
shock finishes crossing the ejecta. However, the detailed modeling of the
dynamical evolution of the reverse shock can be complicated (Sari & Piran,
1995; Sari, 1997). For the purpose of this work, we analytically derive the
characteristic radius and time of the reverse shock at the shock crossing time
$t_{\times}$ (Panaitescu & Kumar, 2004; Murase, 2007), where the reverse shock
finishes crossing the ejecta.

When the ultrarelativistic thick ejecta propagates into the external medium,
two types of shocks are formed: a reverse shock (RS), which propagates back
into the ejecta shell increasing the internal energy, and a forward shock
(FS), which propagates into the external medium and energizes the swept-up
matter. The typical crossing time $t_{\times}$ for the RS to complete its
crossing through the ejecta depends on the width of the ultrarelativistic
ejecta shell $\Delta$ that represents the geometrical thickness measured in
the stellar frame. In the thick ejecta shell regime, the ejecta width can be
estimated as $\Delta\approx c\delta T\simeq 1.9\times 10^{13}~{}\delta
T_{2.8}\rm~{}cm$, where ${\color[rgb]{0,0,0}\delta}T$ is the duration of the
GRB ejecta released by their source measured in the GRB frame. Then the
crossing radius of the reverse shock is $r_{\times}\simeq 3.8\times
10^{17}~{}\mathcal{E}_{k,55}^{1/4}n^{-1/4}\delta T_{2.8}^{1/4}\rm~{}cm$ where
$\mathcal{E}_{k}$ is the isotropic-equivalent kinetic energy and $n$ is the
constant density of the external medium. The crossing time is
$t_{\times}(\delta T)\approx 0.71\delta T(1+z)\simeq
540~{}{\color[rgb]{0,0,0}\delta}T_{2.8}\rm~{}s$, and hereafter the dependence
on $z$ is neglected for simplicity. The Lorentz factor of the shocked ejecta
at the crossing radii $r_{\times}$ can be written as $\Gamma_{\times}\simeq
83~{}\mathcal{E}_{k,55}^{1/8}n^{-1/8}\delta T_{2.8}^{-3/8}$ and $\Gamma_{\rm
rel}\approx 0.5(\Gamma_{\times}/\Gamma_{0}+\Gamma_{0}/\Gamma_{\times})\simeq
1.7\Gamma_{0,2.4}$ when measured in the stellar frame and in the frame of
unshocked ejecta, respectively. Here $\Gamma_{0}$ is the initial Lorentz
factor of GRB ejecta. The reverse shock is relativistic before it crosses the
shell. The contribution to the shell width due to the velocity spread can be
estimated as $r_{\times}/2\Gamma_{0}^{2}\simeq 3.1\times
10^{12}\mathcal{E}_{k,55}^{1/4}n^{-1/4}\delta
T_{2.8}^{1/4}\Gamma_{0,2.4}^{-2}\rm~{}cm<\Delta$ as expected in the thick
ejecta shell regime (Sari & Piran, 1995).

For the earlier evolution phase of the reverse shock, especially for Episode
I, we assume that only a fraction of the total ejecta energy
$\mathcal{E}_{k,\rm I}\approx({\color[rgb]{0,0,0}\delta T_{\rm
I}/\delta}T)\mathcal{E}_{k}$ carried by the outer edge of the ejecta shell is
transferred to the external medium during ${\color[rgb]{0,0,0}t_{\times,{\rm
I}}\equiv t_{\times}(\delta T_{\rm I})<t_{\times}}$. Thus, we can use the same
method as above to estimate the characteristic radius and Lorentz factor of
the shocked ejecta after the reverse shock crosses. For Episode II, we adopt
the same assumptions and use the same method, except that $\delta T_{\rm
II}\sim\delta T$ and $\mathcal{E}_{k,\rm II}\sim\mathcal{E}_{k}$.

Now we proceed to derive the emission properties of the shocked ejecta at the
shock crossing radius for both episodes. Note that for Episode I, the physical
quantities in the following equations are normalized to the numerical values
closer to those expected for Episode II. For Episode II, the magnetic field
strength in the comoving frame is estimated as
$B_{\times}=[32\pi\epsilon_{B}n_{\rm ej}m_{p}c^{2}(\Gamma_{\rm
rel}-1)(\Gamma_{\rm rel}+3/4)]^{1/2}\simeq
8.6{\rm~{}G}\epsilon_{B,-1}^{1/2}\mathcal{E}_{k,55}^{1/4}\Gamma_{0,2.4}^{-1}\delta
T_{2.8}^{-3/4}n^{1/4}(g(\Gamma_{\rm rel})/1.7)^{1/2}$, where $\epsilon_{B}$ is
the energy fraction of internal energy that is converted into the magnetic
energy, the proton number density of the unshocked ejecta is $n_{\rm ej}\simeq
3.2\times 10^{3}~{}\mathcal{E}_{k,55}^{1/2}\Gamma_{0,2.4}^{-2}\delta
T_{2.8}^{-3/2}n^{1/2}\rm~{}cm^{-3}$ and $g(\Gamma_{\rm rel})\equiv(\Gamma_{\rm
rel}-1)(\Gamma_{\rm rel}+3/4)$. The total number of electrons energized by the
reverse shock is $N_{e}^{r}=\mathcal{E}_{k}/\Gamma_{0}m_{p}c^{2}\simeq
2.8\times 10^{55}~{}\mathcal{E}_{k,55}\Gamma_{0,2.4}^{-1}$. The minimum
electron Lorentz factor after shock acceleration can be estimated as
$\gamma_{m}\approx(\epsilon_{e}/f_{e})[(s_{e}-2)/(s_{e}-1)](m_{p}/m_{e})(\Gamma_{\rm
rel}-1)\simeq 4.2\times 10^{2}\epsilon_{e,-1}f_{e,-1}^{-1}$ for $s_{e}=2.6$
and $\Gamma_{\rm rel}=1.7$, where the energy fraction $\epsilon_{e}$ of the
post-shock internal energy is converted into electron non-thermal energy,
$f_{e}$ is the number fraction of accelerated electrons and $s_{e}$ is the
electron spectral index. The main difference between the treatment of Episode
I and Episode II are the values adopted for $\epsilon_{B}$ and $\epsilon_{e}$
(see Table I). We then derive the steady-state electron energy distribution
considering various cooling processes, including adiabatic cooling,
synchrotron cooling, SSC cooling, and EIC cooling, using the iteration method
described in Murase et al. (2011); Zhang et al. (2021a). The comoving frame
non-thermal proton energy density is determined by
$U_{p}\approx\epsilon_{p}\mathcal{E}_{k}/(4\pi
r_{\times}^{2}\Gamma_{\times}^{2}\Delta$), where $\epsilon_{p}$ is the
fraction of downstream energy transferred to the non-thermal protons. The
comoving frame minimum proton energy is $\varepsilon_{p,\rm
min}\approx\Gamma_{\rm rel}m_{p}c^{2}$. The maximum proton energy achieved
under the confinement condition $t_{\rm acc}<t_{\rm dyn}$ is $E_{\rm
max,dyn}\approx\eta^{-1}eB_{\times}r_{\times}\simeq 1.0\times
10^{21}~{}{\rm~{}eV}~{}\eta^{-1}\mathcal{E}_{k,55}^{1/2}\Gamma_{0,2.4}^{-1}\delta
T_{2.8}^{-1/2}\epsilon_{B,-1}^{1/2}(g(\Gamma_{\rm rel})/1.7)^{1/2}$. Here, we
define the acceleration timescale as $t_{\rm acc}=\eta t_{L}$, where $t_{L}$
is the Larmor time, and $\eta$ is a coefficient which is $\sim$a few in the
Bohm limit (e.g., Sironi et al., 2015). The production of UHECRs at the GRB
reverse shock is possible (e.g., Murase et al., 2008; Zhang et al., 2018),
where the maximum proton energy can be limited by various cooling processes,
e.g., for synchrotron cooling we have $E_{\rm max,syn}\simeq 5.5\times
10^{20}~{}{\rm~{}eV}~{}Z^{-3/2}\eta^{-1/2}\epsilon_{B,-1}^{-1/4}\Gamma_{0,2.4}^{1/2}n^{-1/4}(g(\Gamma_{\rm
rel})/1.7)^{-1/4}$. We also include the effect of photomeson cooling on the
maximum energy of protons which depends on the energy density of the target
photon fields. The effect of overlapping of the prompt emission in the reverse
shocked region is considered, where the comoving prompt photon energy density
measured in the reverse shocked region is $U_{\rm GRB\gamma}\approx L_{\rm
GRB\gamma}^{\rm iso}/4\pi r_{\times}^{2}\Gamma_{\times}^{2}c$. The Band
function is used for modeling the energy spectra of the prompt emission, where
the low-energy power law index is $\alpha=1.1$, the high-energy power law
index is $\beta=2.6$ and the observed peak energy is $E_{\rm pk}=1\rm~{}MeV$
(Frederiks et al., 2022).

##  3 Results

The predicted multi-wavelength energy spectra for the two episodes are shown
in Fig. 2, where the corresponding parameters are summarized in Table 1.

Table 1: Physical parameters used in the reverse shock model.

Parameter | Episode I | Episode II
---|---|---
$\Gamma_{0}$ | 250 | 250
$\mathcal{E}_{k}$a [$\rm erg$] | $2\times 10^{55}$ | $2\times 10^{55}$
$n_{\rm ex}$ [$\rm cm^{-3}$] | 1 | 1
$\delta T[$s$]$ |  ${\color[rgb]{0,0,0}300}$b | ${\color[rgb]{0,0,0}600}$
$\epsilon_{B}$ | 0.05 | 0.5
$\epsilon_{e}$ | 0.35 | 0.02
$f_{e}$ | 0.8 | 0.01
$s_{e}$ | 2.6 | 2.6
$\epsilon_{p}$ | 0.1 | 0.08
$s_{p}$ | 2.0 | 2.0
$L_{\rm GRB\gamma}^{\rm iso}$ [$\rm erg~{}s^{-1}$] | $2\times 10^{52}$ | $2\times 10^{50}$

  * a

The dissipated energy during Episode I is only a fraction of the total kinetic
energy $\mathcal{E}_{k,\rm I}\approx(\delta T_{\rm I}/\delta
T)\mathcal{E}_{k}$ as explained in the main text. For Episode II, we have
$\mathcal{E}_{k,\rm II}\approx\mathcal{E}_{k}$.

  * b

This is the duration of the GRB ejecta during Episode I $\delta T_{\rm I}$.

Figure 2: Multi-wavelength energy spectra up to the VHE energy range in our
reverse shock model in Episode I (upper panel) and Episode II (lower panel),
respectively. The black dashed line in the upper panel is the prompt spectrum
with an exponential cutoff at $E_{\gamma,\rm max}=20\rm~{}MeV$. The orange
line in the upper panel indicates the averaged energy spectrum observed by
Fermi-LAT within $T-T_{0}=294-410\rm~{}s$, extrapolated to the energy of 300
GeV. The black line in the low panel is the Fermi-LAT 95% UL within
$T-T_{0}=410-1000\rm~{}s$.

###  3.1 Episode I - Upscattered prompt emission?

In the upper panel of Fig. 2, we show the multi-wavelength energy spectrum
from various processes for Episode I. The orange dashed region represents the
average energy spectrum observed by Fermi-LAT within
$T-T_{0}=294-410\rm~{}s$111We do not consider gamma-rays detected in the first
time interval $T-T_{0}=203-294\rm~{}s$ which overlaps with the LAT Bad Time
Interval., extrapolated to 300 GeV. The dot-dashed curve is the upscattered
prompt emission by the non-thermal electrons accelerated in the reverse shock
region in Episode I which dominates energy flux beyond $\sim 1$ GeV, while the
dashed curve is the SSC component, which dominates the energy flux from 0.1
GeV to 1 GeV. The synchrotron component is indicated as the blue solid curve,
while the proton synchrotron component is marked as the brown solid curve. We
also overlay the prompt emission assuming Band function with an exponential
cutoff at $E_{\gamma,\rm max}=20\rm~{}MeV$. The corresponding parameters are
summarized in Table. 1.

At Episode I, the microphysical parameters satisfy $\epsilon_{e}>\epsilon_{B}$
and the inverse-Compton component dominates the gamma-ray emission. We find
that the EIC spectra can explain the hard energy spectrum observed by Fermi-
LAT above $\sim 1$ GeV. Note that the detected $\sim 99.3\rm~{}GeV$ photon by
Fermi-LAT at $T-T_{0}=240\rm~{}s$ may be explained by the EIC component
considering the enhancement of the prompt photons during the first time
interval mentioned in Fig. 1. However, the proton synchrotron component is not
important in Episode I because the dominance of prompt MeV gamma-rays in the
energy density limits the emission power of the synchrotron component. We note
that the predicted early optical flux could have been brighter than that of
the naked-eye GRB 080319B (Racusin et al., 2008).

###  3.2 Episode II - Proton synchrotron emission?

In the lower panel of Fig. 2, we show the energy spectrum from various
processes for Episode II. Due to the decrease of the prompt luminosity by at
least two orders of magnitude, the EIC component is no longer important, while
the proton synchrotron emission is more prominent at Episode II given that the
magnetic fields are strong. As shown in Table 1, we assume that the
microphysical parameter $\epsilon_{B}=0.5$ is larger than the value adopted in
Episode I. Such a change of $\epsilon_{B}$ is possible for long-lasting
magnetic energy-dominated GRB ejecta. Initially, the value of $\epsilon_{B}$
is very high for magnetic energy-dominated GRB ejecta. During Episode I, most
of the GRB ejecta energy is released in the form of the extraordinarily strong
prompt emission pulse with abundant pair production, which may effectively
suppress $\epsilon_{B}$ especially if the prompt gamma-rays are produced by
magnetic dissipation, and the reverse shock is stronger due to the weak
magnetic energy in the GRB ejecta. Thus, we could expect strong inverse-
Compton emission at Episode I. Later, the reverse shock continues to cross the
inner region of the GRB ejecta, which is still magnetic energy dominated, and
it becomes weak, consistent with the observations. We can expect higher values
of $\epsilon_{B}$ in Episode II, and the proton synchrotron emission is very
efficient. The above physical process could also explain the decrease of
$\epsilon_{e}$ and $\epsilon_{p}$ in Episode II.

The peak energy from proton synchrotron emission can reach $\sim 10\rm~{}TeV$
without extragalactic background light (EBL) absorption. Remarkably, the
corresponding spectral index of the energy flux from proton synchrotron
emission is $(3-s_{p})/2=1/2$ for $s_{p}=2$, which is larger than the spectral
index $\sim 0.2$ inferred from Fermi-LAT observations in Episode I. The harder
spectral index enhances the fraction of $\mathcal{O}(10\rm~{}TeV)$ photons in
the total observed gamma-rays, even though it is undergoing EBL-induced
attenuation during the propagation from the source to Earth (e.g., Baktash et
al., 2022; Zhao et al., 2022).

###  3.3 High-energy neutrino production

Figure 3: Neutrino fluences emitted during the time window $100\rm~{}s$
(Episode I) and $300\rm~{}s$ (Episode II).

In Fig. 3, we show the predicted neutrino fluences from Episode I and Episode
II, respectively. The physical parameters used for the corresponding
calculations are identical to those in Table 1. The neutrino energy spectrum
predicted in Episode I has two bumps, where the low-energy bump with peak
energy at $\sim\rm~{}PeV$ energy range is due to the photomeson production
interaction between high-energy protons with prompt target photon fields,
while the high-energy bump with peak energy at $\sim 10\rm~{}EeV$ is the
results between the interaction of UHE protons with the lower-energy
synchrotron photons from the reverse shock. The neutrino energy spectrum
predicted in Episode II only shows one bump, which is consistent with the
assumption that the effect of prompt emission is no longer important in
Episode II. The maximum proton energy is also higher in Episode II due to the
strong magnetic fields, which is consistent with the neutrino spectrum where
the neutrino has higher energy in Episode II. Note that the cooling of
secondary muons and pions is neglected when calculating the neutrino flux in
our model, which will reduce the flux of the higher energy part of the
observed neutrino spectrum, especially for Episode II (Murase, 2007). The
predicted neutrino fluence is consistent with the upper limits reported by
IceCube (The IceCube Collaboration, 2022; Murase et al., 2022; Ai & Gao,
2022).

##  4 Discussion and implications

###  4.1 Anti-correlation

The $\gamma\gamma$ optical depth $\tau_{\gamma\gamma}$ due to the interaction
with prompt photons at the shocked ejecta region can be estimated with the
following formula (e.g., Murase et al., 2022),

$\tau_{\gamma\gamma}^{\rm
prompt}\simeq{\color[rgb]{0,0,0}11}\frac{\eta_{\gamma\gamma,-1}L_{\rm
GRB\gamma,{\color[rgb]{0,0,0}53.5}}^{\rm
iso}}{r_{17.3}\Gamma_{\times,2}^{2}E_{\rm pk,\rm
MeV}}\begin{cases}(E_{\gamma}/\tilde{E}_{\gamma,b})^{\beta-1},E_{\gamma}<\tilde{E}_{\gamma,b}\\\
(E_{\gamma}/\tilde{E}_{\gamma,b})^{\alpha-1},E_{\gamma}>\tilde{E}_{\gamma,b}\end{cases},$
(1)

where $\eta_{\gamma\gamma}\sim 0.1$ (Svensson, 1987),
$\tilde{E}_{\gamma,b}\approx\Gamma_{\times}^{2}m_{e}^{2}c^{4}/E_{\rm pk,\rm
MeV}\simeq 2.6~{}\Gamma_{\times,2}^{2}\rm~{}GeV$ represents the typical energy
of high-energy $\gamma$-rays that interact with target photons of peak energy
$E_{\rm pk,MeV}$ in the observer frame, and $L_{\rm GRB\gamma}^{\rm iso}$
represents the isotropic-equivalent luminosity in the Konus-Wind band where
the bolometric correction has been accounted for. The optical depth at 1 TeV
is $\tau_{\gamma\gamma}(1\rm~{}TeV)\sim{\color[rgb]{0,0,0}20}$ for
$\alpha=1.1$. In addition, the synchrotron photons become dominant in the
energy range below $\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}\rm keV$ in Episode I compared to the
prompt target photons as shown in Fig. 2. Then the optical depth can be
estimated as

$\tau_{\gamma\gamma}^{\rm syn}\simeq
2300\frac{\eta_{\gamma\gamma,-1}L_{\gamma,50}^{\rm
syn}}{r_{17.3}\Gamma_{\times,2}^{2}E_{\gamma,m}^{\rm
syn}}\begin{cases}(E_{\gamma}/\tilde{E}_{\gamma,m})^{\frac{s_{e}}{2}},E_{\gamma}<\tilde{E}_{\gamma,m}\\\
(E_{\gamma}/\tilde{E}_{\gamma,m})^{\frac{1}{2}},\tilde{E}_{\gamma,c}>E_{\gamma}>\tilde{E}_{\gamma,m}\end{cases}$
(2)

where $L_{\gamma}^{\rm syn}$ is the luminosity of synchrotron emission at the
synchrotron peak energy of the reverse shock component, $E_{\gamma,m}^{\rm
syn}\simeq 1\rm~{}eV$ ($E_{\gamma,c}^{\rm syn}\simeq 0.1\rm~{}eV$) is the
observed synchrotron peak (cooling) energy in the fast cooling regime,
$\tilde{E}_{\gamma,m}\approx\Gamma_{\times}^{2}m_{e}^{2}c^{4}/E_{\gamma,m}^{\rm
syn}\simeq 2600~{}\Gamma_{\times,2}^{2}\rm~{}TeV$ and
$\tilde{E}_{\gamma,c}\approx\Gamma_{\times}^{2}m_{e}^{2}c^{4}/E_{\gamma,c}^{\rm
syn}\simeq 26000~{}\Gamma_{\times,2}^{2}\rm~{}TeV$. Then the optical depth at
10 TeV is $\tau_{\gamma\gamma}(10\rm~{}TeV)\sim 2$ for $s_{e}=2.6$. Note that
the analytical calculation of the optical depth in Eq. 2 is only valid for
$E_{\gamma}<\tilde{E}_{\gamma,c}$, as the differential spectral index drops
below unity ($2/3<1$) when the photon energy is below $E_{\gamma,c}^{\rm
syn}$. Thus, in our model, we predict the anti-correlation of the
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\rm TeV$ photons during the strongest
prompt GRB emission phase, with the former escaping at later times. For
Episode II, due to the lower luminosity of prompt emission and the synchrotron
photons, we can expect the escape of $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\rm TeV$ photons
simultaneously with low-energy photons.

###  4.2 Implications for LHAASO detection

Figure 4: The predicted detection (cumulative) number for the LHAASO detector
with the energy spectrum predicted from our model using three different EBL
models with the total integration time $\Delta T=100$ s for Episode I and
$\Delta T=300$ s for Episode II. The vertical lines indicate the observed
photon energy, and the shaded region corresponds to the energy range
$18\times(100\pm 40\%)\rm~{}TeV$ with a $\sim 40\%$ energy resolution at
$18\rm~{}TeV$ energy band of LHAASO-KM2A (Ma et al., 2022).

In Fig. 4, we show the number of photons predicted in our model that can be
detected by LHAASO. Considering the uncertainties on the measured EBL models,
we adopt three representative EBL models, namely Franceschini08 (Franceschini
et al., 2008), Gilmore12 (Gilmore et al., 2012) and Finke10 (Finke et al.,
2010). We adopt the effective area of the LHAASO-WCDA for photons provided in
Wang et al. (2022) for zenith window $15^{\circ}<\theta<30^{\circ}$ and the
effective area from LHAASO-KM2A (Ma et al., 2022). The integrated number of
events can be estimated as

$N(>E)=\int_{E}^{E_{\rm max}}dEF(E)A_{\rm eff}^{\gamma}(E,\theta)\Delta T,$
(3)

where $F(E)$ is the predicted energy spectrum, $A_{\rm
eff}^{\gamma}(E,\theta)$ is the effective area including both LHAASO-WCDA and
LHAASO-KM2A, and $\Delta T$ is the integration time. Even though the LHAASO-
WCDA dominates the effective area at low energy, the effective area of LHAASO-
KM2A is comparable to LHAASO-WCDA above $\sim 10\rm~{}TeV$ (Ma et al., 2022).

If we adopt a typical integration time $300\rm~{}s$ for Episode II, then the
number of photons detected is $N(>0.5\rm~{}TeV)\sim{\color[rgb]{0,0,0}4000}$
for photons with energy larger than 500 GeV,
$N(>1\rm~{}TeV)\sim{\color[rgb]{0,0,0}1500}$ for photons with energy larger
than 1 TeV, and $N(>10\rm~{}TeV)\sim{\color[rgb]{0,0,0}7}$ for photons with
energy larger than 10 TeV. The expected cosmic ray background for LHAASO-WCDA
can be estimated with the same method used in Wang et al. (2022). We find that
the photons with energies larger than 500 GeV (1 TeV) could be detected with
LHAASO-WCDA with a significance level $\sim 25\rm~{}s.d.$ ($\sim
15\rm~{}s.d.$) estimated through Eq. 17 of Li & Ma (1983) using fitted cosmic
ray spectrum in Particle Data Group et al. (2020).

Even though we expect $N(>0.5\rm~{}TeV)\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$<$}}}{\color[rgb]{0,0,0}300}$ photons detected
in Episode I, the detection of photons with energy above 1 TeV may be
challenging for LHAASO-WCDA in Episode I, where the significance level is
expected to be less than $\sim 2\rm~{}s.d.$ in our fiducial case. TeV photons
during Episode I mainly come from SSC emission in our model. However, there is
significant uncertainty in the SSC component, and if more
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\rm TeV$ gamma-rays have been detected by
LHAASO in Episode I, it would indicate that the SSC component becomes more
important than the fiducial case shown in this work, or the proton synchrotron
emission is enhanced because of the larger value of $\epsilon_{B}$ as in the
case of Episode II. The VHE gamma-rays detected by LHAASO in Episode I may be
primarily contributed by the SSC component resulting from non-thermal
electrons that are accelerated by the forward shock (e.g., Ren et al., 2022;
Sato et al., 2022). The fitting of the low-frequency radio to optical spectra
at the earlier stages indicates that additional components, possibly resulting
from synchrotron emission from ejecta swept by the reverse shock, are required
beyond the standard forward shock model (O’Connor et al., 2023). However,
further investigation using multi-wavelength data is needed to determine the
relative contribution of the observed VHE gamma-rays from forward and reverse
shock.

In the thick shell case, it takes time for the shocked ejecta to get adjusted
to the Blandford-McKee (BM) profile after the reverse shock crosses the
ejecta. As in synchrotron emission from primary electrons accelerated at the
reverse shock (Kobayashi, 2000), we expect that the proton synchrotron light
curve will decay with a temporal index of $-(73s_{p}+21)/96\sim-1.7$ for
$s_{p}=2$, which is typically steeper than the forward shock case. However,
the details of the reverse shock light curve depend on the profile of the
ejecta, for long-lasting reverse shocks especially from $t_{\times,\rm I}$ to
$t_{\times}$, which can be shallower than the standard case.

The LHAASO-KM2A is extremely effective for suppressing the cosmic ray
background (Ma et al., 2022), so the detection of dozens of photons with
energy $\mathcal{O}(10\rm~{}TeV)$ is reasonable considering the energy
resolution. The energy resolution at $\mathcal{O}(10\rm~{}TeV)$ is $\Delta
E/E\sim 45\%$ for LHAASO-KM2A and $\Delta E/E\sim 60\%$ for LHAASO-WCDA. Our
results indicate that even without considering the effect of new physics on
the propagation of VHE gamma-rays, e.g., photon-ALP mixing (Galanti et al.,
2022; Baktash et al., 2022; Nakagawa et al., 2022; Troitsky, 2022) and Lorentz
invariance violation (LIV) (Li & Ma, 2022; Zhu & Ma, 2022; Finke & Razzaque,
2022), the detection of $\sim 18\rm~{}TeV$ photons by LHAASO-KM2A can be
explained for reasonable EBL models. Another possible explanation for the
detection of $\sim 18\rm~{}TeV$ is the intergalactic electromagnetic cascade
due to the propagation of UHECRs (Batista, 2022; Das & Razzaque, 2022;
Mirabal, 2022). Both models also require the efficient production of UHECRs
from GRBs. But in the intergalactic cascade scenario, the time delay of the
photon arrival during the propagation of UHECRs and the development of the
electromagnetic cascade depends strongly on the magnetic field structure of
the host galaxy, the host galaxy cluster, and the intergalactic medium, and
the time delay is longer given that the magnetic field strength is larger than
$10^{-17}~{}\rm G$ (Takahashi et al., 2008; Murase et al., 2009, 2012;
Mirabal, 2022).

The default value of the coefficient in the acceleration timescale formula
$t_{\rm acc}=\eta t_{L}$ is set to $\eta=1$, which may be optimistic. If we
adopt a more considerable value of $\eta=10$, then the proton synchrotron
component would be reduced, which makes the detection of $\mathrel{\hbox{\hbox
to0.0pt{\hbox{\lower 4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\rm TeV$ gamma-rays
more challenging in Episode II as shown in Fig. 4. Thus, our results imply
that the reverse shock of the long-lasting GRB 221009A should be a very
efficient accelerator of non-thermal particles. If most of the
$\mathrel{\hbox{\hbox to0.0pt{\hbox{\lower
4.0pt\hbox{$\sim$}}\hss}\hbox{$>$}}}\rm TeV$ gamma-rays are detected after
Episode II, which is $T\sim 600\rm~{}s$ in this work, then this would mean
that the GRB ejecta duration is a bit longer and the reverse shock finishes
crossing the GRB ejecta at later times.

##  5 Summary

The detection of GRB 221009A in the VHE gamma-ray band up to
$\mathcal{O}(10\rm~{}TeV)$ by LHAASO provides us with an opportunity to study
the radiative processes of GRBs in the highest energy range. In this work, we
studied the origin of VHE gamma-rays from GRB 221009A in the framework of the
reverse shock model where non-thermal electrons and protons are expected to be
accelerated. The reverse shock could last up to thousands of seconds after the
start of GRB, which is consistent with the long-lasting prompt emission
observed from GRB 221009A.

We considered two episodes, where the emission from the reverse shock in the
first episode (Episode I) is strongly affected by the strong prompt emission,
while the effect of prompt emission becomes weak in the second episode
(Episode II). In addition, the microphysical parameters, e.g., $\epsilon_{B}$
and $\epsilon_{e}$, could be different from Episode I to Episode II causing
different behavior in the emission processes.

Our results show that the upscattered prompt MeV photons by the non-thermal
electrons accelerated in the reverse shock region in Episode I mainly
contribute to the energy flux observed by Fermi-LAT above $\sim\rm GeV$, in
addition to the SSC component. The emission from the proton synchrotron
process is not important in Episode I because of the presence of MeV gamma-
rays and weak magnetic field strength.

We found that the proton synchrotron process can dominate the output in the
VHE band in Episode II, where the magnetic field strength is strong enough to
increase the proton synchrotron emission significantly. Our rough estimates
show that $\sim{\color[rgb]{0,0,0}4000}$ photons with energy larger than 0.5
TeV can be detected by LHAASO in Episodes I and II, which may be consistent
with the number of photons detected by LHAASO within 2000 seconds. Due to the
hard spectral index of the proton synchrotron emission compared to the
inverse-Compton process as inferred from Fermi-LAT data, it is plausible to
detect dozens of $\mathcal{O}(10\rm~{}TeV)$ photons using reasonable EBL
models without invoking new physics.

Note that the forward shock has a similar energy to the reverse shock, and it
may also contribute to the VHE emission. The forward shock could enhance the
number of $\sim\rm~{}TeV$ gamma-rays detected by LHAASO, but difficult for
$\mathcal{O}(10\rm~{}TeV)$ gamma-rays, which emphasizes the role of proton
synchrotron emission from the reverse shock as proposed in this work.

In the future, we can expect more GRBs to be detected in the VHE gamma-ray
band, especially with the Cherenkov Telescope Array (Inoue et al., 2013;
Kakuwa et al., 2012). Our work suggests that the observation of GRB in the VHE
gamma-ray band can be used for constraining the particle acceleration and
radiative processes of non-thermal electrons and protons in the reverse shock
model.

The work was partly supported by the NSF Grants No. AST-1908689, No.
AST-2108466 and No. AST-2108467, and KAKENHI No. 20H01901 and No. 20H05852
(K.M.) and No. 22H00130, 20H01901, 20H01904, 20H00158, 18H01215, 17H06357,
17H06362 (K.I.).

## References

  * Abdalla et al. (2019) Abdalla, H., Adam, R., Aharonian, F., et al. 2019, Nature, 575, 464.  http://www.nature.com/articles/s41586-019-1743-9
  * Ai & Gao (2022) Ai, S., & Gao, H. 2022, Model constraints based on the IceCube neutrino non-detection of GRB 201009A, arXiv, arXiv:2210.14116 [astro-ph], doi:10.48550/arXiv.2210.14116.  http://arxiv.org/abs/2210.14116
  * An et al. (2023) An, Z.-H., Antier, S., Bi, X.-Z., et al. 2023, Insight-HXMT and GECAM-C Observations of the Brightest-of-All-Time GRB 221009A, arXiv, arxiv:arXiv:2303.01203
  * Asano et al. (2009) Asano, K., Inoue, S., & Mészáros, P. 2009, The Astrophysical Journal, 699, 953.  https://iopscience.iop.org/article/10.1088/0004-637X/699/2/953
  * Baktash et al. (2022) Baktash, A., Horns, D., & Meyer, M. 2022, Interpretation of multi-TeV photons from GRB221009A, arXiv, arXiv:2210.07172 [astro-ph], doi:10.48550/arXiv.2210.07172.  http://arxiv.org/abs/2210.07172
  * Batista (2022) Batista, R. A. 2022, GRB 221009A: a potential source of ultra-high-energy cosmic rays, arXiv, arXiv:2210.12855 [astro-ph], doi:10.48550/arXiv.2210.12855.  http://arxiv.org/abs/2210.12855
  * Beloborodov (2005) Beloborodov, A. M. 2005, ApJL , 618, L13
  * Das & Razzaque (2022) Das, S., & Razzaque, S. 2022, Ultrahigh-energy cosmic-ray signature in GRB 221009A, arXiv, arXiv:2210.13349 [astro-ph, physics:hep-ph], doi:10.48550/arXiv.2210.13349.  http://arxiv.org/abs/2210.13349
  * de Ugarte Postigo et al. (2022) de Ugarte Postigo, A., Izzo, L., & Pugliese, G. 2022, GCN CIRCULAR, 32648, 1
  * Dermer et al. (2000) Dermer, C. D., Chiang, J., & Mitman, K. E. 2000, Astrophys. J., 537, 785
  * Dichiara et al. (2022) Dichiara, S., Gropp, J. D., Kennea, J. A., et al. 2022, The Astronomer’s Telegram, 15650, 1, aDS Bibcode: 2022ATel15650….1D.  https://ui.adsabs.harvard.edu/abs/2022ATel15650....1D
  * Finke & Razzaque (2022) Finke, J. D., & Razzaque, S. 2022, Possible Evidence for Lorentz Invariance Violation in Gamma-ray Burst 221009A, arXiv, arXiv:2210.11261 [astro-ph, physics:hep-ph], doi:10.48550/arXiv.2210.11261.  http://arxiv.org/abs/2210.11261
  * Finke et al. (2010) Finke, J. D., Razzaque, S., & Dermer, C. D. 2010, The Astrophysical Journal, 712, 238, arXiv:0905.1115 [astro-ph].  http://arxiv.org/abs/0905.1115
  * Franceschini et al. (2008) Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, Astronomy & Astrophysics, 487, 837.  http://www.aanda.org/10.1051/0004-6361:200809691
  * Frederiks et al. (2022) Frederiks, D., Lysenko, A., Ridnaia, A., & Svinkin, D. 2022, GCN CIRCULAR, 32668, 1
  * Frederiks et al. (2023) Frederiks, D., Svinkin, D., Lysenko, A. L., et al. 2023, Properties of the Extremely Energetic GRB~221009A from Konus-\textit{}WIND{} and \textit{}SRG{}/ART-XC Observations, arXiv, arxiv:arXiv:2302.13383
  * Galanti et al. (2022) Galanti, G., Roncadelli, M., & Tavecchio, F. 2022, Axion-like particles explain the very-high energy emission from GRB221009A, , , arXiv:2210.05659 [astro-ph, physics:hep-ph], doi:10.48550/arXiv.2210.05659.  http://arxiv.org/abs/2210.05659
  * Gallant & Achterberg (1999) Gallant, Y. A., & Achterberg, A. 1999, Monthly Notices of the Royal Astronomical Society, 305, L6.  https://academic.oup.com/mnras/article-lookup/doi/10.1046/j.1365-8711.1999.02566.x
  * Gill & Granot (2022) Gill, R., & Granot, J. 2022, Galaxies, 10, 74
  * Gilmore et al. (2012) Gilmore, R. C., Somerville, R. S., Primack, J. R., & Domínguez, A. 2012, Monthly Notices of the Royal Astronomical Society, 422, 3189.  https://academic.oup.com/mnras/article-lookup/doi/10.1111/j.1365-2966.2012.20841.x
  * He et al. (2012) He, H.-N., Zhang, B.-B., Wang, X.-Y., Li, Z., & Mészáros, P. 2012, Astrophys. J., 753, 178
  * Huang et al. (2022) Huang, Y., Hu, S., Chen, S., et al. 2022, GCN CIRCULAR, 32677.  https://gcn.gsfc.nasa.gov/gcn3/32677.gcn3
  * Inoue et al. (2013) Inoue, S., Granot, J., O’Brien, P. T., et al. 2013, Astroparticle Physics, 43, 252.  https://linkinghub.elsevier.com/retrieve/pii/S0927650513000054
  * Ioka et al. (2005) Ioka, K., Kobayashi, S., & Zhang, B. 2005, The Astrophysical Journal, 631, 429\.  https://iopscience.iop.org/article/10.1086/432567
  * Isravel et al. (2022) Isravel, H., Pe’er, A., & Begue, D. 2022, Proton Synchrotron Origin of the Very High Energy Emission of GRB 190114C, arXiv, arXiv:2210.02363 [astro-ph], doi:10.48550/arXiv.2210.02363.  http://arxiv.org/abs/2210.02363
  * Kakuwa et al. (2012) Kakuwa, J., Murase, K., Toma, K., et al. 2012, Monthly Notices of the Royal Astronomical Society, 425, 514.  https://academic.oup.com/mnras/article-lookup/doi/10.1111/j.1365-2966.2012.21490.x
  * Kimura et al. (2019) Kimura, S. S., Murase, K., Ioka, K., et al. 2019, Astrophys. J. Lett., 887, L16
  * Kobayashi (2000) Kobayashi, S. 2000, The Astrophysical Journal, 545, 807
  * Kumar & Zhang (2014) Kumar, P., & Zhang, B. 2014, Phys. Rept., 561, 1
  * Li & Ma (2022) Li, H., & Ma, B.-Q. 2022, Lorentz invariance violation induced threshold anomaly versus very-high energy cosmic photon emission from GRB20221009A, , , arXiv:2210.06338 [astro-ph, physics:gr-qc, physics:hep-ph], doi:10.48550/arXiv.2210.06338.  http://arxiv.org/abs/2210.06338
  * Li & Ma (1983) Li, T. P., & Ma, Y. Q. 1983, The Astrophysical Journal, 272, 317, aDS Bibcode: 1983ApJ…272..317L.  https://ui.adsabs.harvard.edu/abs/1983ApJ...272..317L
  * Ma et al. (2022) Ma, X.-H., Bi, Y.-J., Cao, Z., et al. 2022, Chinese Physics C, 46, 030001.  https://iopscience.iop.org/article/10.1088/1674-1137/ac3fa6
  * MAGIC Collaboration (2019) MAGIC Collaboration. 2019, Nature, 575, 455.  http://www.nature.com/articles/s41586-019-1750-x
  * MAGIC Collaboration et al. (2019) MAGIC Collaboration, Veres, P., Bhat, P. N., et al. 2019, Nature, 575, 459.  http://www.nature.com/articles/s41586-019-1754-6
  * Mészáros (2006) Mészáros, P. 2006, Reports on Progress in Physics, 69, 2259.  https://iopscience.iop.org/article/10.1088/0034-4885/69/8/R01
  * Mészáros & Rees (1994) Mészáros, P., & Rees, M. J. 1994, Mon. Not. Roy. Astron. Soc., 269, L41
  * Miceli & Nava (2022) Miceli, D., & Nava, L. 2022, Galaxies, 10, 66, arXiv:2205.12146 [astro-ph].  http://arxiv.org/abs/2205.12146
  * Mirabal (2022) Mirabal, N. 2022, Secondary GeV-TeV emission from ultra-high-energy cosmic rays accelerated by GRB 221009A, arXiv, arXiv:2210.14243 [astro-ph, physics:hep-ex], doi:10.48550/arXiv.2210.14243.  http://arxiv.org/abs/2210.14243
  * Murase (2007) Murase, K. 2007, Physical Review D, 76, 123001.  https://link.aps.org/doi/10.1103/PhysRevD.76.123001
  * Murase et al. (2012) Murase, K., Dermer, C. D., Takami, H., & Migliori, G. 2012, Astrophys. J., 749, 63
  * Murase et al. (2008) Murase, K., Ioka, K., Nagataki, S., & Nakamura, T. 2008, Physical Review D, 78, 023005.  https://link.aps.org/doi/10.1103/PhysRevD.78.023005
  * Murase et al. (2022) Murase, K., Mukhopadhyay, M., Kheirandish, A., Kimura, S. S., & Fang, K. 2022, Neutrinos from the Brightest Gamma-Ray Burst?, arXiv, arXiv:2210.15625 [astro-ph, physics:hep-ph], doi:10.48550/arXiv.2210.15625.  http://arxiv.org/abs/2210.15625
  * Murase et al. (2011) Murase, K., Toma, K., Yamazaki, R., & Mészáros, P. 2011, The Astrophysical Journal, 732, 77.  https://iopscience.iop.org/article/10.1088/0004-637X/732/2/77
  * Murase et al. (2010) Murase, K., Toma, K., Yamazaki, R., Nagataki, S., & Ioka, K. 2010, MNRAS , 402, L54
  * Murase et al. (2009) Murase, K., Zhang, B., Takahashi, K., & Nagataki, S. 2009, Mon. Not. Roy. Astron. Soc., 396, 1825
  * Nakagawa et al. (2022) Nakagawa, S., Takahashi, F., Yamada, M., & Yin, W. 2022, Axion dark matter from first-order phase transition, and very high energy photons from GRB 221009A, arXiv, arXiv:2210.10022 [astro-ph, physics:hep-ph], doi:10.48550/arXiv.2210.10022.  http://arxiv.org/abs/2210.10022
  * O’Connor et al. (2023) O’Connor, B., et al. 2023, arXiv:2302.07906
  * Omodei et al. (2022a) Omodei, N., Bruel, P., & Bregeon, J. 2022a, GCN CIRCULAR, 32760, 1
  * Omodei et al. (2022b) —. 2022b, GCN CIRCULAR, 32916, 1
  * Panaitescu & Kumar (2004) Panaitescu, A., & Kumar, P. 2004, Monthly Notices of the Royal Astronomical Society, 353, 511.  https://academic.oup.com/mnras/article-lookup/doi/10.1111/j.1365-2966.2004.08083.x
  * Particle Data Group et al. (2020) Particle Data Group, Zyla, P. A., Barnett, R. M., et al. 2020, Progress of Theoretical and Experimental Physics, 2020, 083C01.  https://academic.oup.com/ptep/article/doi/10.1093/ptep/ptaa104/5891211
  * Pillera et al. (2022) Pillera, R., Bissaldi, E., Omodei, N., La Mura, G., & Longo, F. 2022, The Astronomer’s Telegram, 15656, 1, aDS Bibcode: 2022ATel15656….1P.  https://ui.adsabs.harvard.edu/abs/2022ATel15656....1P
  * Racusin et al. (2008) Racusin, J. L., Karpov, S. V., Sokolowski, M., et al. 2008, Nature, 455, 183.  http://www.nature.com/articles/nature07270
  * Ren et al. (2022) Ren, J., Wang, Y., & Zhang, L.-L. 2022, Very High Energy Afterglow Emission of GRB~221009A: Lessons Learned from the Brightest Long Gamma-ray Burst in a Wind Environment, arXiv, arXiv:2210.10673 [astro-ph], doi:10.48550/arXiv.2210.10673.  http://arxiv.org/abs/2210.10673
  * Sari (1997) Sari, R. 1997, The Astrophysical Journal, 489, L37.  https://iopscience.iop.org/article/10.1086/310957
  * Sari & Esin (2001) Sari, R., & Esin, A. A. 2001, Astrophys. J., 548, 787
  * Sari & Piran (1995) Sari, R., & Piran, T. 1995, The Astrophysical Journal, 455, doi:10.1086/309835.  https://iopscience.iop.org/article/10.1086/309835
  * Sato et al. (2022) Sato, Y., Murase, K., Ohira, Y., & Yamazaki, R. 2022, arXiv:2212.09266
  * Sironi et al. (2015) Sironi, L., Keshet, U., & Lemoine, M. 2015, Space Sci. Rev. , 191, 519
  * Svensson (1987) Svensson, R. 1987, Monthly Notices of the Royal Astronomical Society, 227, 403.  https://academic.oup.com/mnras/article-lookup/doi/10.1093/mnras/227.2.403
  * Takahashi et al. (2008) Takahashi, K., Murase, K., Ichiki, K., Inoue, S., & Nagataki, S. 2008, Astrophys. J. Lett., 687, L5
  * The IceCube Collaboration (2022) The IceCube Collaboration, I. 2022, GCN CIRCULAR, 32665, 1
  * Toma et al. (2011) Toma, K., Wu, X.-F., & Mészáros, P. 2011, Monthly Notices of the Royal Astronomical Society, 415, 1663.  https://academic.oup.com/mnras/article-lookup/doi/10.1111/j.1365-2966.2011.18807.x
  * Totani (1998) Totani, T. 1998, The Astrophysical Journal, 502, L13.  https://iopscience.iop.org/article/10.1086/311489
  * Troitsky (2022) Troitsky, S. V. 2022, Parameters of axion-like particles required to explain high-energy photons from GRB 221009A, arXiv, arXiv:2210.09250 [astro-ph, physics:hep-ph], doi:10.48550/arXiv.2210.09250.  http://arxiv.org/abs/2210.09250
  * Veres et al. (2022) Veres, P., Bissaldi, E., Meegan, C., & Fermi GBM Team. 2022, GRB Coordinates Network, 31816, 1, aDS Bibcode: 2022GCN.31816….1V.  https://ui.adsabs.harvard.edu/abs/2022GCN.31816....1V
  * Veres & Mészáros (2014) Veres, P., & Mészáros, P. 2014, Astrophys. J., 787, 168
  * Vietri (1996) Vietri, M. 1996, Monthly Notices of the Royal Astronomical Society, 278, L1.  https://academic.oup.com/mnras/article-lookup/doi/10.1093/mnras/278.1.L1
  * Wang et al. (2022) Wang, X.-Y., Bi, X.-J., Cao, Z., et al. 2022, Chinese Physics C, 46, 030003.  https://iopscience.iop.org/article/10.1088/1674-1137/ac3fa9
  * Wang et al. (2006) Wang, X.-Y., Li, Z., & Mészáros, P. 2006, Astrophys. J. Lett., 641, L89
  * Waxman (1995) Waxman, E. 1995, Physical Review Letters, 75, 386.  https://link.aps.org/doi/10.1103/PhysRevLett.75.386
  * Waxman & Bahcall (2000) Waxman, E., & Bahcall, J. N. 2000, The Astrophysical Journal, 541, 707.  https://iopscience.iop.org/article/10.1086/309462
  * Zhang & Mészáros (2001) Zhang, B., & Mészáros, P. 2001, Astrophys. J., 559, 110
  * Zhang et al. (2018) Zhang, B. T., Murase, K., Kimura, S. S., Horiuchi, S., & Mészáros, P. 2018, Physical Review D, 97, 083010.  https://link.aps.org/doi/10.1103/PhysRevD.97.083010
  * Zhang et al. (2021a) Zhang, B. T., Murase, K., Veres, P., & Mészáros, P. 2021a, The Astrophysical Journal, 920, 55.  https://iopscience.iop.org/article/10.3847/1538-4357/ac0cfc
  * Zhang et al. (2021b) Zhang, B. T., Murase, K., Yuan, C., Kimura, S. S., & Mészáros, P. 2021b, The Astrophysical Journal Letters, 908, L36.  https://iopscience.iop.org/article/10.3847/2041-8213/abe0b0
  * Zhao et al. (2022) Zhao, Z.-C., Zhou, Y., & Wang, S. 2022, Standard physics is still capable to interpret $\sim$18 TeV photons from GRB~221009A, arXiv, arXiv:2210.10778 [astro-ph], doi:10.48550/arXiv.2210.10778.  http://arxiv.org/abs/2210.10778
  * Zhu & Ma (2022) Zhu, J., & Ma, B.-Q. 2022, Light speed variation from GRB 221009A, arXiv, arXiv:2210.11376 [astro-ph, physics:gr-qc, physics:hep-ph], doi:10.48550/arXiv.2210.11376.  http://arxiv.org/abs/2210.11376