# Multi-messenger model for the prompt emission from GRB 221009A

Annika Rudolph Niels Bohr International Academy and DARK, Niels Bohr
Institute, University of Copenhagen, Blegdamsvej 17,
2100, Copenhagen, Denmark  Deutsches Elektronen-Synchrotron DESY,
Platanenallee 6, 15738 Zeuthen, Germany  Maria Petropoulou Department of
Physics, National and Kapodistrian University of Athens, University Campus
Zografos, GR 15783, Athens, Greece  Institute of Accelerating Systems &
Applications, University Campus Zografos, GR 15783, Athens, Greece  Walter
Winter Deutsches Elektronen-Synchrotron DESY, Platanenallee 6, 15738 Zeuthen,
Germany  Željka Bošnjak Faculty of Electrical Engineering and Computing,
University of Zagreb, Unska ul. 3, 10000 Zagreb, Croatia

###### Abstract

We present a multi-messenger model for the prompt emission from GRB 221009A
within the internal shock scenario. We consider the time-dependent evolution
of the outflow with its impact on the observed light curve from multiple
collisions, and the self-consistent generation of the electromagnetic spectrum
in synchrotron and inverse Compton-dominated scenarios. Our leptohadronic
model includes UHE protons potentially accelerated in the outflow, and their
feedback on spectral energy distribution and on the neutrino emission. We find
that we can roughly reproduce the observed light curves with an engine with
varying ejection velocity of ultra-relativistic material, which has an
intermediate quiescent period of about 200 seconds and a variability timescale
of $\sim 1$ s. We consider baryonic loadings of 3 and 30 that are compatible
with the hypothesis that the highest-energetic LHAASO photons might come from
UHECR interactions with the extragalactic background light, and the paradigm
that energetic GRBs may power the UHECR flux. For these values and the high
dissipation radii considered we find consistency with the non-observation of
neutrinos and no significant signatures on the electromagnetic spectrum.
Inverse Compton-dominated scenarios from the prompt emission are demonstrated
to lead to about an order of magnitude higher fluxes in the HE-range; this
enhancement is testable by its spectral impact in the Fermi-GBM and LAT
ranges.

gamma-ray burst: GRB 221009A – cosmic rays – neutrinos – radiation mechanisms:
non-thermal

##  1 Introduction

Gamma-ray bursts (GRBs) are extremely energetic explosions involving the
collapse of a massive star or the merger of two compact objects. As their name
suggests, GRBs release most of their electromagnetic output in $\gamma$ rays
within a short period of time (ranging from tens of milliseconds to hundreds
of seconds). This brief and variable emission, known as the prompt GRB phase,
is followed by the afterglow, a long-lasting multi-wavelength emission (for a
review see, e.g. Kumar & Zhang, 2015). The prompt emission is thought to be
produced in a relativistic collimated plasma outflow (jet) launched by the
central engine via some dissipative mechanism. For kinetically dominated jets,
a leading scenario involves energy dissipation and particle acceleration at
internal shocks that are produced when portions of the jet are moving outwards
with varying Lorentz factors (Rees & Meszaros, 1994; Kobayashi et al., 1997;
Daigne & Mochkovitch, 1998). Part of the remaining jet energy can be later
transformed to non-thermal radiation at a relativistic blast wave sweeping up
material from the circumburst medium and powering the afterglow (Rees &
Mészáros, 1992; Chiang & Dermer, 1999).

GRBs are one of the prime targets of multi-messenger astronomy, as they have
been detected in gravitational waves (GW 170817/GRB 170817, Abbott et al.,
2017) and have been proposed as candidate sources of ultra-high-energy cosmic
rays (UHECRs) and astrophysical neutrinos (e.g. Waxman, 1995; Vietri, 1995;
Waxman & Bahcall, 1997). While the contribution of typical luminosity GRBs
during their prompt phase to the diffuse neutrino flux measured by IceCube has
been constrained to $\lesssim 1\%$ of the diffuse astrophysical neutrino flux
(Abbasi et al., 2022), the hypothesis that GRBs are UHECR sources cannot be
ruled out yet (for a recent review see Kimura, 2022). In fact, GRBs populating
the high-end of the isotropic $\gamma$-ray energy distribution
($E_{\gamma,{\rm iso}}>10^{54}$ erg) may provide the necessary energy output
per event for powering UHECRs, requiring only a moderate baryonic loading
(defined as the energy injected into non-thermal protons versus electrons) and
without violating existing neutrino limits (Rudolph et al., 2022). In
addition, energetic bursts may be detected as a single source (e.g. Gao et
al., 2013).

On October 9th 2022, a very bright GRB was observed at redshift $z=0.151$ (de
Ugarte Postigo et al., 2022). The burst triggered the Gamma-Ray Burst Monitor
(GBM) on board Fermi at 2022-10-09 13:16:59.000 UT (Veres et al., 2022), about
an hour before the detection of a hard X-ray transient by the Burst Alert
Telescope (BAT) of the Neil Gehrels Swift satellite (S. Lesage et al., 2022;
Dichiara et al., 2022). The prompt phase of GRB 221009A consisted of a
precursor (at about 10 s), followed by an extremely bright emission period
about 200 s post GBM trigger. Overall, the prompt emission period lasted
roughly 327 s and was composed of several peaks (S. Lesage et al., 2022). The
preliminary Konus-Wind light curve showed several peaks of roughly 40 s
duration (D. Frederiks et al., 2022); the initial peaks were separated from
the late-time peak by a quiescent period of about 220 s. Some short-timescale
variability on the order of seconds might be also visible in the 0.4-100 MeV
AGILE MCAL light curve (A. Ursi et al., 2022) and INTEGRAL SPI-ACS light curve
(J. Rika et al., 2022). The burst was also observed in high energies (HE,
$\gtrsim 0.1$ GeV) by the Fermi Large Area Telescope (LAT), starting about 200
s after the GBM trigger (i.e. during the prompt phase) and extending up to
$\sim 25$ ks into the afterglow phase (R. Pillera et al., 2022). In addition,
very high energy (VHE, $>100$ GeV) photons were detected by LHAASO, but their
association to the prompt phase is not clear given that these photons were
observed within a period of up to 2000 s after trigger (Yong Huang et al.,
2022). It has been speculated that the highest energy photons (up to 18 TeV)
might come from UHECR interactions with the extragalactic background light
(EBL), since such energetic photons escaping the source would be otherwise
attenuated by the EBL (Das & Razzaque, 2022; Alves Batista, 2022; Mirabal,
2022). These scenarios require a significant amount of energy carried by UHE
protons, which might also leave signatures in the electromagnetic spectrum.

The extreme brightness of this burst caused pile-up in almost all GRB
detectors, namely Fermi-GBM and LAT, KONUS-Wind, and AGILE. For this reason,
detailed spectral analysis was unavailable at the time of writing. Preliminary
analysis of LAT data (100 MeV - 1 GeV) for $200-800$ s after the GBM trigger
provided an estimate of the photon index ($-1.87\pm 0.04$) and the photon flux
$(6.2\pm 0.4)\times 10^{-3}$ ph cm-2 s-1 (Pillera et al., 2022), but thet time
interval excluded from the analysis was recently extended to 300 s (Omodei et
al., 2022). Hence, the relevance of these results to the main GRB episode
detected by GBM is not yet clear. Moreover, preliminary analysis of Konus data
during the brightest phase of the event (i.e. $\sim 180-200$ s after the Konus
trigger) produced a rest-frame peak energy of $E_{\mathrm{peak}}=1060$ keV and
$E_{\gamma,\mathrm{iso}}\simeq 3\times 10^{54}$ erg. No associated muon-
neutrino track was detected by IceCube in a time range of [-1 hour, +2 hours]
from the initial GBM trigger, which resulted in an upper limit on the muon
neutrino fluence of $3.9\times 10^{-2}$ GeV cm-2 assuming an $E^{-2}$ neutrino
spectrum (The IceCube Collaboration, 2022). The inferred $\gamma$-ray
isotropic energy, the proximity of this event to Earth, and the lack of prompt
neutrino detection make GRB 221009A a unique case for the study of multi-
messenger signatures from GRBs.

In this Letter, we present a multi-messenger model for the prompt emission of
GRB 221009A. Under the assumption that protons and electrons are accelerated
at the fastest possible rate in internal shocks occurring at different radii
within the jet, we compute the multi-messenger emission from each collision
while taking into account the varying physical conditions in the outflow and
the UHECR feedback on both the photon and neutrino emissions. Our goal is to
test the hypothesis of UHECR acceleration in the GRB jet by comparing the
self-consistently computed broadband photon spectrum and the accompanying
neutrino flux with available observational information. Since details on the
prompt spectrum are not yet available, and the mentioned pile-up effects
likely introduce some degeneracy, our model also has some predictive power. In
this work we indicatively use the following observables: (1) the peak energy
$E_{\mathrm{peak}}=1060$ keV (as fitted for the onset of the bright emission
period by Konus), (2) the estimated fluence in the GBM band $\mathcal{F}_{\rm
GBM}=2.91\cdot 10^{-2}$ erg cm-2 ($1-1000$ keV) and (3) the approximate light
curve structure observed by Konusand INTEGRAL SPI-ACS. Peak energy and fluence
are reproduced within $\pm 5$ %.

Table 1: Fireball characteristics and microphysics parameters for the
different scenarios.

|  | “$R_{\mathrm{16}}$”-scenario | “$R_{\mathrm{17}}$”-scenario
---|---|---|---
Quantity | Symbol | SYN-dom. | IC-dom. | SYN-dom. | IC-dom.
Engine time for main emission period |  $t_{\mathrm{main}}$ [s] | 74
Engine quiescent time |  $t_{\mathrm{quiet}}$ [s] | 213
Engine time for late-time activity |  $t_{\mathrm{late}}$ [s] | 8
Variability timescale of engine activity |  $\delta t_{\mathrm{var}}$ [s] | 1.4
Averaged $\Gamma$ at beginning | $\langle\Gamma_{\mathrm{ini}}\rangle$ | 265 | 663
Averaged $\Gamma$ at end | $\langle\Gamma_{\mathrm{fin}}\rangle$ | 228 | 570
Averaged $\Gamma$ of emitting plasma | $\langle\Gamma_{\mathrm{em}}\rangle$ | 293 | 731
Averaged radius of emitting plasma |  $\langle R_{\mathrm{Coll}}\rangle$ [cm] | $1.2\cdot 10^{16}$ | $2.0\cdot 10^{17}$
Total energy transferred to non-thermal electrons |  $E_{\mathrm{e,NT}}$ [erg] | $4.6\cdot 10^{54}$ | $9.3\cdot 10^{54}$ | $4.8\cdot 10^{54}$ | $9.3\cdot 10^{54}$
Initial fireball kinetic energy |  $E_{\mathrm{kin,ini}}$ [erg] | $5.7\cdot 10^{56}$ | $9.2\cdot 10^{56}$ | $3.8\cdot 10^{57}$ | $7.1\cdot 10^{57}$
Averaged maximal proton energy |  $\langle E_{\mathrm{p,max}}\rangle$ [$10^{11}$ GeV] | 21.9 | 2.7 | 23.7 | 1.1
Emitted gamma-ray energy (1 - $10^{4}$ keV) |  $E_{\mathrm{\gamma,iso}}$ [erg] | $2.9\cdot 10^{54}$ | $3.1\cdot 10^{54}$ | $2.9\cdot 10^{54}$ | $2.9\cdot 10^{54}$
Relative fraction of energy transferred to magnetic field | $f_{\mathrm{B/e}}=\epsilon_{\mathrm{B}}/\epsilon_{\mathrm{e}}$ | $1$ | $10^{-3}$ | $1$ | $10^{-3}$
Relative fraction of energy transferred to acc. protons | $f_{\mathrm{p/e}}=\epsilon_{\mathrm{p}}/\epsilon_{\mathrm{e}}$ | 3 | 3 | 30 | 30
Power-law index of accelerated electrons | $p_{\mathrm{e}}$ | 2.2
Minimum Lorentz factor of accelerated electrons |  $\gamma_{\mathrm{e,min}}^{\prime}$ [$10^{4}$] | 3 | 9 | 6 | 11
Power-law index of accelerated protons | $p_{\mathrm{p}}$ | 2.0
Minimum Lorentz factor of accelerated electrons | $\gamma_{\mathrm{p,min}}^{\prime}$ | 10

  * •

Notes. – All listed energies refer to the isotropic equivalent values.
Energies and times are reported in the source (engine) frame. Averaged
quantities are computed weighing with the shell mass
($\langle\Gamma_{\mathrm{ini}}\rangle$ and
$\langle\Gamma_{\mathrm{fin}}\rangle$) or with the dissipated energy of a
collision ($\langle\Gamma_{\mathrm{em}}\rangle$, $\langle
R_{\mathrm{Coll}}\rangle$, $\langle E_{\mathrm{p,max}}\rangle$). The maximal
proton energies are computed considering synchrotron and adiabatic losses and
equally refer to the engine frame.

##  2 Model and implementation

Our model is fully described in an accompanying paper (Rudolph et al., 2022)
to which we refer the interested readers for details. Here we summarize the
key points of the model and present its main parameters in Tab. 1. In short, a
relativistic outflow is discretized as shells ejected with different Lorentz
factors from a central engine operating over a time $t_{\mathrm{eng}}$ (as
measured in the engine rest frame). Shells catch up with each other and merge
at a radius $R_{\mathrm{Coll}}$ from the central emitter, where energy is
dissipated and particles are accelerated. Shells continue to propagate, merge
and dissipate energy until their velocity distribution is such that no more
collision occur. The dissipated energy in each inelastic collision, which can
be obtained from energy and momentum conservation, is distributed into non-
thermal electrons, protons and magnetic fields (for simplicity we assume that
no energy remains in thermal plasma). We introduce the partition parameters
$f_{\rm p/e}$ (baryonic loading) and $f_{\rm B/e}$ (magnetic loading)
describing the ratio between proton/electron and magnetic/electron energy
density, respectively (here only referring to accelerated, non-thermal
particles). Assuming that the sub-MeV prompt spectrum is predominately
produced by synchrotron radiation of accelerated electrons, we normalize the
fireball kinetic energy to the total energy transferred to non-thermal
electrons $E_{\mathrm{e,NT}}^{\mathrm{tot}}$ that is needed to produce a given
$E_{\mathrm{\gamma,iso}}$. Based on our simulations the energy dissipation
efficiency111The fireball radiative efficiency is then found by multiplying
$\varepsilon_{\rm diss}$ with the fraction of energy carried by primary
electrons and the synchrotron radiative efficiency. of the fireball is
$\varepsilon_{\mathrm{diss}}\simeq 0.04$, for the specific Lorentz factor
distribution assumed here. The required fireball kinetic energy (isotropic
equivalent) is then indirectly related to the energy partition parameters and,
in particular, increases with increasing baryonic loading as
$E_{\mathrm{kin,ini}}=\varepsilon_{\rm
diss}^{-1}E_{\mathrm{e,NT}}^{\mathrm{tot}}(1+f_{\mathrm{p/e}}+f_{\mathrm{B/e}})$.
In contrast to Rudolph et al. (2022), we normalize all models to a similar
$E_{\mathrm{\gamma,iso}}$, leading to higher $E_{\mathrm{kin,ini}}$ for cases
with a low electron synchrotron efficiency.

Figure 1: Left: Distribution of Lorentz factors of plasma shells launched by
the central engine, here $R_{\mathrm{ini}}$ is the initial radius of a plasma
shell. Two cases are shown leading to high (left axis) and low (right axis)
collision radii. The parameter $t_{\mathrm{eng}}$ represents the total
activity time of the engine, i.e.,
$t_{\mathrm{eng}}=t_{\mathrm{main}}+t_{\mathrm{quiet}}+t_{\mathrm{late}}$.
Middle and right: Two-dimensional phase space of the collision radius and the
total energy carried by non-thermal electrons (middle) or the Lorentz factor
of the emitting plasma (right). Results are shown for the case of low
dissipation radii and SYN-dominated electron cooling (Table 1). The collisions
(each represented by a circle) are separated in color by the time interval
during which they will be observed: purple (0–36.5 s), blue (36.5–71.5 s),
aquamarine (71.5–129.5 s) and green (from 308.5 s). Grouped collisions
correspond to four distinctive pulses of the light curve shown in Fig. 2.

To reproduce the structure of the observed light curve we use a varying
Lorentz factor profile, see Fig. 1 (left panel), here without addressing the
question which engine properties would lead to such a profile. The peaks of
the light curve will be reproduced by the overall variation of the initial
Lorentz factor $\Gamma_{\rm ini}$, and the short-term flux variability by
superimposed fluctuations on $\Gamma_{\rm ini}$ on a scale $\delta
t_{\mathrm{var}}\simeq 1$s. The relative peak heights will be controlled by
the initial Lorentz factor and kinetic energy distribution of the shells.
Assuming that the power of the engine is constant (while it is active), the
initial shells will have equal kinetic energies in the source frame. Note also
that a long quiescent engine period (about 200 s in the engine frame) is
needed. The engine activity is thus three time intervals: The activity time of
the main emission period $t_{\mathrm{main}}$, the engine quiescent time
$t_{\mathrm{qiet}}$ and activity time for the late-time emission
$t_{\mathrm{late}}$.

We present two different scenarios for the initial Lorentz factor distribution
– denoted as “$R_{\mathrm{16}}$” and “$R_{\mathrm{17}}$” – that produce
collisions at an average radius $\langle R_{\mathrm{Coll}}\rangle\sim 10^{16}$
cm and $\langle R_{\mathrm{Coll}}\rangle\sim 2\cdot 10^{17}$ cm respectively,
and probe higher and lower typical plasma densities. We verified that in both
scenarios the bulk of energy dissipation occurs below the estimated
deceleration radius for typical parameters of the circumburst medium. These
scenarios are motivated by estimates for the optical thickness of the emitting
plasma to $\gamma\gamma$ pair production (Murase et al., 2022) and the
requirement that the bulk of energy is dissipated below the approximate
deceleration radius. In principle, high $\Gamma$ factors are also expected
from the empirical $E_{\gamma,\mathrm{iso}}-\Gamma$ relationships (Ghirlanda
et al., 2018). As an example, we show for “$R_{\mathrm{16}}$” the two-
dimensional distributions of the non-thermal electron energy and the Lorentz
factor of the emitting shell with collision radius in the middle and right
panels of Fig. 1.

We self-consistently evaluate the emission from non-thermal electrons and
protons, accelerated in shocks from the collision of shells, and injected with
power-law distributions into the radiation zone, which is the hot plasma of
the merged shell. The full-burst emission is then integrated over all
collisions, taking into account the curvature of the emitting surface. The
maximal electron and proton energies are limited by the dominating energy loss
processes assuming efficient acceleration (operating at the Bohm limit). For
each scenario, the minimum electron Lorentz factor
$\gamma_{\mathrm{e,min}}^{\prime}$ is set such the peak energy of 1060 keV is
reproduced; The minimal proton energies are generally fixed to
$\gamma^{\prime}_{\mathrm{p,min}}=10$. We assume a power-law slope of primary
protons of $p_{\mathrm{p}}=2.0$ as in Rudolph et al. (2022). In our model the
primary electron power-law slope would usually be determined by the high-
energy photon index. As there was no reliable measurement for the high-energy
slope at the time of writing, we adapt $p_{\mathrm{e}}=2.2$ as suggested for
mildly relativistic shocks (Crumley et al., 2019). We then numerically
solve222The simulations are performed with the proprietary code AM3 (Gao et
al., 2017). the coupled system of integrodifferential equations describing the
temporal evolution of particle distributions (electrons/positrons, photons,
protons, neutrons, pions, muons, neutrinos) including the relevant processes,
such as synchrotron emission and absorption, inverse Compton scattering (ICS),
photo-pair and photo-pion production, $\gamma\gamma$ pair production,
adiabatic cooling, and escape. Our approach fully captures the electromagnetic
cascade induced by photo-hadronic and $\gamma\gamma$ pair production in each
merged shell333We do not consider interactions of particle populations between
shells..

As discussed in detail in Rudolph et al. (2022), the electromagnetic spectrum
will be, even in the leptohadronic case, dominated by the primary leptonic
emission: a synchrotron component peaking in the MeV band and an inverse
Compton component emerging at higher energies (GeV band) for low enough
$f_{\rm B/e}$ values. We therefore discuss a scenarios synchrotron
(“SYN”)-dominated and inverse Compton (“IC”)-dominated one with different
values of $f_{\rm B/e}$ (see Tab. 1), characterized by different photon
spectra and light curves. We also include the effects of EBL attenuation using
the model of Dominguez et al. (2011), calculated with the open-source gammapy-
package (Deil et al., 2018; Nigro et al., 2019).

A crucial parameter in every multi-messenger model is the baryonic loading.
Rudolph et al. (2022) showed that $f_{\rm p/e}\gtrsim 3$ is required for
$E_{\gamma,\mathrm{iso}}\simeq 3\cdot 10^{54}\,\mathrm{erg}$, if such
energetic GRBs ought to power the UHECRs. Much higher baryonic loadings in
combination with low $R_{\mathrm{Coll}}$ lead to spectral distortions of the
photon spectrum (even in the GBM and LAT bands), and efficient neutrino
production, which might be in tension with multi-messenger limits (for
details, see Sec. 6 in Rudolph et al., 2022). A different argument comes from
the observed VHE photons in LHAASO, if these were produced by interactions of
UHECRs with the EBL (Das & Razzaque, 2022; Alves Batista, 2022; Mirabal,
2022). Das & Razzaque (2022) derive an isotropic-equivalent
$E_{p,\mathrm{iso}}\simeq 3.9\cdot 10^{54}\,\mathrm{erg}$ in the range 0.1-100
EeV, which yields $E_{p,\mathrm{iso}}\gtrsim 10^{55}\,\mathrm{erg}$
bolometrically corrected for our considered energy range, if all UHECRs can
escape. This corresponds to our $f_{\rm p/e}\gtrsim 3$. For Alves Batista
(2022), the given fraction into UHECRs is lower, but the bolometric correction
factor (only UHECRs at the highest energies are considered) is much higher. In
what follows, we choose $f_{\rm p/e}=3$ for “$R_{\mathrm{16}}$” and $f_{\rm
p/e}=30$ for “$R_{\mathrm{17}}$”. The former is compatible with the above
estimates, while the latter is a more aggressive version (still compatible
with neutrino bounds) that will challenge the GRB energetics in terms of the
kinetic energy required. Purely leptonic results will not be shown, as we
found no substantial modifications due to hadronic contributions in the
observable photon spectra for the baryonic loadings considered.

##  3 Results

Figure 2: Synthetic light curves for three energy ranges obtained in the
“$R_{\mathrm{16}}$” scenario for the IC- and SYN-dominated cases. In the left
panel we show only the IC-dominated case as it resembles the SYN-dominated
one. We also add observed light curves of Konus and INTEGRAL SPI-ACS that
operate in a comparable energy range, re-normalised to match the same scale
and shifted in time to match the same peak times (-175 s for Konus, and -225 s
for INTEGRAL SPI-ACS). We set $T_{\mathrm{obs}}=0$ for the synthetic light
curve as the minimum collision time in the observer’s frame.

Light curves. By construction of our model, the rough structure of the light
curves shown in Fig. 2 reflects the initial Lorentz factor distribution in
Fig. 1. Each pulse of the synthetic light curves is formed by shells colliding
over a wide range of radii (as can be inferred from the colour coding in Fig.
1). If all these collisions produced the same amount of radiated power then
time-lags between different energy ranges even within the same pulse would be
possible.

In the left panel, we indicate the observed light curves of INTEGRAL SPI-ACS
and Konus, shifted to match the beginning of our synthetic light curve that
does not include the precursor emission. The general structure of the observed
light curve is reproduced by our model, but our results should also be
representative for light curves with a similar general structure. The relative
intensity of the two bright pulses in the observed light curve may strongly be
impacted by pile-up effects in the detector. Here, we assumed that the first
bright peak is intrinsically much more energetic than the second one. In our
model the relative height of the peak is controlled by the ratio of Lorentz
factors of the colliding shells (e.g. a bright peak is produced when this
ratio is large). Available $\gamma$-ray light curves at the time of writing
indicate short-time variability on the second timescale, which was modeled by
stochastic variations of $\Gamma_{\rm ini}$ on a timescale $\delta
t_{\mathrm{var}}=1.4$ s (see Fig. 1). However, if refined analysis of the
light curves revealed variability on a shorter timescale, this would shift the
typical collision radius inwards where the densities are higher. As shown in
Rudolph et al. (2022), this would lead to stronger signatures of secondary
particles, such as secondary leptons and neutrinos, for the same baryonic
loading.

The synthetic HE light curve (shown in the middle panel of Figure 2) is very
similar to the keV-MeV light curve. Preliminary analysis of LAT data using
gtburst444https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gtburst.html
showed no evidence for a fourth peak in the $20~{}\rm MeV-300~{}\rm GeV$ light
curve. If this is later verified by detailed LAT analysis, then late-time
collisions with low Lorentz factors will be needed in our model in order to
suppress the late-time HE emission due to high internal opacity to
$\gamma\gamma$ pair production. While the shape of the HE light curve is
similar in the SYN- and IC-dominated scenarios, the HE flux is higher in the
latter case, suggesting that primary electrons are cooling more efficiently
via ICS (see also next paragraph). Non-negligible VHE emission between 1-10
TeV is also expected in the IC-dominated scenario (right-hand panel), with
similar light curve as the in the lower energy bands.

Figure 3: Modelled spectra $E_{\mathrm{obs}}\mathcal{F}_{\mathrm{E_{obs}}}$
for the (left) “$R_{\mathrm{16}}$” and “$R_{\mathrm{17}}$” (right) scenarios,
in both cases including a SYN- and a IC-dominated scenario (for parameters see
Table 1). Transparent curves correspond to the fluence without taking EBL
attenuation into account, dashed curves to the per-flavour neutrino fluences.
The latter are compared with the IceCube upper limit reported in The IceCube
Collaboration (2022). The inset shows a zoom-in to the spectrum around the
peak in the Fermi GBM/LAT range (indicated as a shaded region in all plots),
with a dotted vertical line at 1060 keV to indicate the position of the
observed peak. For this energy range we further indicate the photon indices,
that are defined of the slope of d$N_{\mathrm{ph}}$/d$E_{\mathrm{obs}}$.

Photon spectra. Fig. 3 shows the predicted time-integrated photon fluences as
a function of observed energy; the inset is a zoom-in on the spectra around
the peak including the photon index (that is defined as the spectral slope of
d$N_{\mathrm{ph}}$/d$E_{\mathrm{obs}}$).

We first discuss the results for the SYN-dominated cases (blue lines in both
panels). The spectrum around the peak is by construction independent of the
typical collision radius and given by the standard synchrotron fast-cooling
predictions: approaching a photon index $-3/2$ below the peak and $-(p_{\rm
e}+2)/2$ above the peak. Differences between “$R_{\mathrm{17}}$” and
“$R_{\mathrm{16}}$” are visible at low and high energies with respect to the
MeV peak. For “$R_{\mathrm{16}}$”, where the densities are higher, emission of
secondaries increases the fluence in the eV range, leading to a low-energy
spectral break. However, the spectrum in the LAT band is still dominated by
the synchrotron emission of fast-cooling primaries. For “$R_{\mathrm{17}}$”,
the power-law spectrum with photon index $-1.5$ extends to the lowest energies
(till synchrotron self-absorption becomes important), while higher maximal
synchrotron energies yield a harder high-energy spectrum.

In the “$R_{\mathrm{16}}$” IC-dominated case (red curves in left panel), the
spectral slope below the peak is still $-1.5$, but the HE emission (in the LAT
range) is modified by the ICS emission of primary and secondary leptons that
creates an almost flat spectrum (i.e. photon index $\sim-2$). In this case,
the peak of the broadband spectrum is shifted to $\sim 1$ GeV. In the
“$R_{\mathrm{17}}$” IC-dominated case (red curves in right panel) the photon
index below the MeV peak (but not within the GBM band) becomes asymptotically
$\sim-1.25$, which is the highest value that can be produced in our model.
This is indicative of electron cooling via ICS in the Klein-Nishina regime
(Daigne et al., 2011). In the LAT band the spectrum is a power-law that
extends to about 10 TeV (where the Klein-Nishina ICS emission of primary
electrons with Lorentz factor $\gamma^{\prime}_{\rm e,min}$ dominates).
However, due to EBL attenuation this spectral feature is washed out. We note
that in the “$R_{\mathrm{17}}$” IC-dominated case the contribution of
secondary leptons in the GBM and LAT bands is negligible. For a more detailed
discussion on spectra and their decomposition, we point to Sec. 4 and 5 in
Rudolph et al. (2022).

In addition to the full time-integrated spectra, we computed the spectra for
the three main emission periods, namely the first two bright pulses and the
late-time pulse (not explicitly shown). Finding little difference in the
spectra for these three pulses we predict no significant spectral evolution
during these three emission periods.

Neutrinos. We include in Fig. 3 the predicted neutrino spectra (per
flavor)555The synchrotron cooling of pions and muons is taken into account,
see more detailed discussion in Rudolph et al. (2022). and the corresponding
IceCube limits (The IceCube Collaboration, 2022) for the different scenarios
as dashed curves. We also compute the number of expected neutrino events in
IceCube with the appropriate point-source effective area for the declination
range 18-21 degrees (Abbasi et al., 2021), and find $n_{\nu_{\mu}}=0.012$
(0.006) for the “$R_{\mathrm{17}}$” SYN- (IC-) dominated case, and
$n_{\nu_{\mu}}=0.17$ (0.29) for the “$R_{\mathrm{16}}$” SYN- (IC-) dominated
case. The predicted neutrino fluences are thus below the IceCube limits in all
cases and consistent with non-detection, as a result of the relatively large
$R_{\mathrm{Coll}}$ paired with the chosen baryonic loadings in consistency
with the findings in one zone models (Ai & Gao, 2022; Murase et al., 2022).
For “$R_{\mathrm{16}}$”, however, a baryonic loading of $f_{\mathrm{p/e}}\gg
3$ is expected to be in tension with neutrino limits, and if observations
eventually favor a variability timescale $\delta t_{\mathrm{var}}\lesssim 1$ s
the baryonic loading would be limited to an even smaller value.

Contrary to Rudolph et al. (2022), we normalise the initial engine kinetic
energy to achieve the same $\gamma$-ray isotropic energy, which for the IC-
dominated scenarios increases the required energy (see Tab. 1). This increases
the energy transferred to non-thermal protons, and subsequently the neutrino
fluences. The effect can be noticed clearly for the
“$R_{\mathrm{16}}$”-scenario. For the “$R_{\mathrm{17}}$”-scenario the
neutrino production efficiency is limited by the low(er) maximal proton
energies, due to the lower magnetic fields obtained in the IC-dominated case.
It is interesting that the peak neutrino energies of $10^{17}$ to
$10^{19}\,\mathrm{eV}$ exceed the expectation of the standard neutrino model
for GRBs ($10^{15}\,\mathrm{eV}$, see e.g. Hummer et al. (2012)) by at least
two orders of magnitude in energy. This is a result of the synchrotron-cooling
dominated spectral indices below the peak (the photon number density peaks at
lower energies) paired with weak magnetic field effects on the secondaries as
a consequence of large $R_{\mathrm{Coll}}$. Therefore, energetic GRBs may be a
target for future radio detection experiments.

##  4 Discussion

There are a number of effects which can be included in order to enrich model.
For example, the different pulses may be produced by collisions of shells with
very different Lorentz factor ranges or even microphysics parameters. This can
cause a combination of SYN- and IC-dominated scenarios in different peaks (see
also Zhang et al., 2022, for the reverse shock), or suppression of VHE
emission in others (by low Lorentz factors enhancing the $\gamma\gamma$
optical thickness, and also the neutrino production). One may speculate that a
non-observation of the late-term peak by Fermi-LAT or a contribution to the
LHAASO signal could be produced by such effects. Within each peak, low Lorentz
factors and a strong correlation between collision radius and observation time
can cause an early suppression VHE photons, which may be interpreted as a
delay (Bustamante et al., 2017). While our spectral index above the peak can
be adjusted by the electron injection index and the efficiency of ICS, our
model may not accommodate very hard low-energy photon indices of $\sim-1$,
unless additional components, such as photospheric thermal emission, are
considered. In lepto-hadronic models these components would increase the
number of target photons available for photo-hadronic interactions. We note
that the highest photon index of $-1.25$ was obtained in the
“$R_{\mathrm{17}}$” IC-dominated scenario, but well below the GBM band.
Whether such a high dissipation radius is indeed realistic should be confirmed
by observations of the variability timescale and estimates of the Lorentz
factor.

The GRB central engine may be a newly formed accreting black hole or magnetar.
For energetic bursts, such as GRB 221009A, the latter scenario can be excluded
because the available energy is limited by the magnetar’s rotational energy to
$\sim 2\cdot 10^{52}$ erg (Usov, 1992; Thompson et al., 2004). The rotational
energy, $E_{\rm rot}$, of an accreting black hole can be extracted via
electromagnetic fields, provided there is a strong large-scale magnetic field
threading the black hole horizon (Blandford-Znajek mechanism Blandford &
Znajek, 1977). The isotropic equivalent energy of the jet can be written as
$E_{\rm jet,iso}=\eta_{j}f_{\rm b}^{-1}E_{\rm rot}$, where $\eta_{j}<1$ is the
fraction of rotational energy ending up in the jet, $f_{\rm
b}=1-\cos(\theta_{j})\approx\theta_{j}^{2}/2$ is the beaming factor,
$\theta_{j}$ is the jet half-opening angle, $E_{\rm rot}=f(a)M_{\rm BH}c^{2}$,
$f(a_{*})=1-\sqrt{(1+\sqrt{1-a_{*}^{2}})/2}$, and $a_{*}$ is the dimensionless
black hole spin. Adopting $\theta_{j}=3.5$ deg (D’Avanzo et al., 2022)666This
was obtained using a prompt-phase efficiency of 0.2. In our model this
efficiency is lower, which would require either a larger density of the
surrounding medium or yield a smaller opening angle. and $\eta_{j}=0.5$ we
find that a maximally spinning ($a_{*}=1$) black hole with $M_{\rm
BH}=10~{}M_{\odot}$ can produce $E_{\rm jet,iso}\simeq 1.4\cdot 10^{57}$ erg,
which is comfortably larger than $E_{\rm kin,ini}$ required for the
“$R_{\mathrm{16}}$” scenario. Moderate spins ($a_{*}\sim 0.5$) would require
$M_{\rm BH}\sim 40~{}M_{\odot}$ to meet the model’s energetic requirements.
The formation of such massive black holes in the collapsar scenario for long-
duration GRBs (with stellar masses $<40~{}M_{\odot}$ in the zero-age main
sequence) is not expected (Woosley, 1993) – see, however, Siegel et al. (2021)
for very massive collapsars. The “$R_{\mathrm{17}}$”-scenario with baryonic
loading 30 corresponds to $E_{\rm kin,ini}\geq(3-5)\cdot 10^{57}$ erg which
would need $a_{*}=1$ and $M_{\rm BH}\geq 20~{}M_{\odot}$. Hence, this scenario
is unlikely on energetic grounds.

The energy budget is further typically constrained from the afterglow
brightness. Here it is noteworthy to mention that a part of the energy in the
afterglow being dissipated in VHE may relax the requirements on the prompt-
phase efficiency, especially for proton-synchrotron models for the VHE
emission (Isravel et al., 2022). A part of the afterglow energy going into
thermal particles may equally increase the allowed kinetic energy of the
blastwave after the prompt phase (but should be taken into account for the
modeling, see Warren et al. (2018); Margalit & Quataert (2021)). Afterglow
observations are also used to infer typical parameters of the GRB, for example
Ren et al. (2022) find a Lorentz factor $\Gamma_{0}=190$ (with room for
slightly higher $\Gamma_{0}$, see their Fig. 3) for the afterglow. This would
favour our “$R_{\mathrm{16}}$”-scenario, that has a typical Lorentz factor
$\langle\Gamma_{\mathrm{fin}}\rangle\sim 230$ after the prompt phase.

Our obtained (averaged) maximal proton energies are in the range between
$10^{20}$ and $2\cdot 10^{21}\,\mathrm{eV}$ under the assumption of efficient
particle acceleration, see Table 1. Higher proton energies are expected in the
SYN-dominated scenarios than the IC-dominated ones because the acceleration
rate is higher due to higher magnetic fields. Our maximal energies are
compatible with the values used e.g. in Das & Razzaque (2022) (100 EeV) to
describe the LHAASO VHE photons from EBL interactions. We note that the time
delay induced by the extragalactic magnetic fields (EGMFs) requires extremely
low field values paired with large proton energies, which means that this
challenge can be somewhat mitigated in our SYN-dominated scenario by the high
proton energies. We also note that the EGMF induced delay is very large
(assumed to be limited to the observed LHAASO window of 2000s) even under
aggressive assumptions for the magnetic field, which means that the protons
must be accelerated in the prompt phase of the GRB and further delays induced
by the afterglow cannot be accommodated. We find that most of the UHECR
protons are emitted within the first 100s in our model, which is compatible
with this picture. Our obtained maximal proton energies are also higher than
the maximum corresponding rigidity $R_{\mathrm{max}}\simeq 1-3\,\mathrm{EV}$
required to describe UHECR data (Heinze et al., 2019) – where details depend
on the assumed cutoff shape.

##  5 Summary and conclusions

In this letter we have presented a state-of-the-art multi-messenger emission
model for the prompt emission of GRB 221009A. In this model, plasma shells are
ejected from a central engine with varying Lorentz factors; these eventually
catch up and energy is dissipated in internal shocks. Our radiation model
includes the effect of UHECR protons self-consistently, such as the
electromagnetic cascade in each shell, generated from secondary electrons,
positrons and photons produced in photo-hadronic interactions. Our assumptions
for the baryonic loading (3 and 30) have been motivated by the paradigm that
energetic GRBs, such as GRB 221009A, could be sources of the UHECRs; they are
also consistent with the hypothesis that the highest-energy LHAASO photons
come from interactions of UHECRs with the EBL.

We have demonstrated that an intermittent engine can reproduce the observed
prompt light curves if a quiescent period of about 200 seconds is included and
assuming a variability timescale of $\sim 1$s. We have implemented relatively
large Lorentz factors and therefore collision radii supported by various
arguments (such as $\gamma\gamma$ optical thickness, neutrino non-
observation). Our predicted electromagnetic spectra exhibit synchrotron fast
cooling dominated spectral indices below the MeV peak, except for an ICS
dominated scenario with energy dissipation at high radii that yielded a softer
spectral index. Above the peak, ICS can affect the spectral index very
strongly, an effect which can extend up to the highest energies and enhance
the VHE fluxes. On the other hand, hadronic contributions did not
significantly alter the photon spectra. Since pile-up effects affect the
analysis for nearly all instruments due to the high brightness of the GRB, the
predictive power of our model may be useful. The predicted neutrino emission
was consistent with the non-observation of neutrinos by IceCube due to high
predicted peak energies in the range interesting for radio neutrino
telescopes. For lower emission radii (due to variability on shorter timescales
and/or lower Lorentz factors) the baryonic loading would likely be constrained
to lower values than the ones assumed here. Our findings are consistent with
the available rotational energy which can be extracted from a maximally
spinning black hole with a mass of the order of $10\,M_{\odot}$ if the
baryonic loading is not too far away from energy equipartition; our standard
assumption of a baryonic loading of 3 is consistent with this picture.

We conclude that GRB 221009A is an interesting object to test the internal
shock model and the paradigm that energetic GRBs could be the sources of
UHECRs. While direct signatures of cosmic rays (such as neutrinos) have not
been seen, the LHAASO observation of TeV photons could point towards UHE
proton acceleration. Our model connects the different messengers for the
prompt phase of this GRB.

Note. – During completion of this work Liu et al. (2022) appeared. In contrast
to their paper, we self-consistently model the photon spectra and account for
several emission regions along the jet. Their low(er) emission radii are the
consequence of a variability timescale of $82\,\mathrm{ms}$ that was derived
from an analysis of the GBM light curve up to $219\,\mathrm{s}$, whereas we
inferred the variability timescale from the INTEGRAL SPI-ACS light curve also
including the bright emission period most relevant for the spectra. Overall,
their limits on the baryonic loading are compatible with our findings.

Acknowledgements

We would like to thank Irene Tamborra, Iftach Sadeh and Marc Klinger for
reading the manuscript and useful comments and Sylvia Zhu and Andrew Taylor
for discussions around GRB 221009A. A.R. received funding from the Carlsberg
Foundation (CF18-0183). M.P. acknowledges support from the MERAC Fondation
through the project THRILL and from the Hellenic Foundation for Research and
Innovation (H.F.R.I.) under the “2nd call for H.F.R.I. Research Projects to
support Faculty members and Researchers” through the project UNTRAPHOB
(Project ID 3013).

Software: astropy(Robitaille et al., 2013; Price-Whelan et al., 2018), gammapy
(Deil et al., 2018; Nigro et al., 2019) and Python v3.9.

## References

  * A. Ursi et al. (2022) A. Ursi et al. 2022, GCN Circ. 32650
  * Abbasi et al. (2021) Abbasi, R., et al. 2021, arXiv e-prints, doi: 10.21234/CPKQ-K003
  * Abbasi et al. (2022) —. 2022, arXiv e-prints.  https://arxiv.org/abs/2205.11410
  * Abbott et al. (2017) Abbott, B. P., Abbott, R., Abbott, T. D., et al. 2017, Phys. Rev. Lett., 119, 161101, doi: 10.1103/PhysRevLett.119.161101
  * Ai & Gao (2022) Ai, S., & Gao, H. 2022, arXiv e-prints.  https://arxiv.org/abs/2210.14116
  * Alves Batista (2022) Alves Batista, R. 2022, arXiv e-prints.  https://arxiv.org/abs/2210.12855
  * Blandford & Znajek (1977) Blandford, R. D., & Znajek, R. L. 1977, Mon. Not. Roy. Astron. Soc., 179, 433
  * Bustamante et al. (2017) Bustamante, M., Murase, K., Winter, W., & Heinze, J. 2017, Astrophys. J., 837, 33, doi: 10.3847/1538-4357/837/1/33
  * Chiang & Dermer (1999) Chiang, J., & Dermer, C. D. 1999, The Astrophysical Journal, 512, 699, doi: 10.1086/306789
  * Crumley et al. (2019) Crumley, P., Caprioli, D., Markoff, S., & Spitkovsky, A. 2019, Mon. Not. Roy. Astron. Soc., 485, 5105, doi: 10.1093/mnras/stz232
  * D. Frederiks et al. (2022) D. Frederiks et al. 2022, GCN Circ. 32668
  * Daigne et al. (2011) Daigne, F., Bosnjak, Z., & Dubus, G. 2011, Astron. Astrophys., 526, A110, doi: 10.1051/0004-6361/201015457
  * Daigne & Mochkovitch (1998) Daigne, F., & Mochkovitch, R. 1998, Mon. Not. Roy. Astron. Soc., 296, 275, doi: 10.1046/j.1365-8711.1998.01305.x
  * Das & Razzaque (2022) Das, S., & Razzaque, S. 2022, arXiv e-prints.  https://arxiv.org/abs/2210.13349
  * D’Avanzo et al. (2022) D’Avanzo, P., Ferro, M., Brivio, R., et al. 2022, GRB Coordinates Network, 32755, 1
  * de Ugarte Postigo et al. (2022) de Ugarte Postigo, A., Izzo, L., Pugliese, G., et al. 2022, GRB Coordinates Network, 32648, 1
  * Deil et al. (2018) Deil, C., et al. 2018, PoS, ICRC2017, 766, doi: 10.22323/1.301.0766
  * Dichiara et al. (2022) Dichiara, S., Gropp, J. D., Kennea, J. A., et al. 2022, The Astronomer’s Telegram, 15650, 1
  * Dominguez et al. (2011) Dominguez, A., et al. 2011, Mon. Not. Roy. Astron. Soc., 410, 2556, doi: 10.1111/j.1365-2966.2010.17631.x
  * Gao et al. (2013) Gao, S., Kashiyama, K., & Mészáros, P. 2013, Astrophys. J. Lett., 772, L4, doi: 10.1088/2041-8205/772/1/L4
  * Gao et al. (2017) Gao, S., Pohl, M., & Winter, W. 2017, Astrophys. J., 843, 109, doi: 10.3847/1538-4357/aa7754
  * Ghirlanda et al. (2018) Ghirlanda, G., Nappo, F., Ghisellini, G., et al. 2018, Astron. Astrophys., 609, A112, doi: 10.1051/0004-6361/201731598
  * Heinze et al. (2019) Heinze, J., Fedynitch, A., Boncioli, D., & Winter, W. 2019, Astrophys. J., 873, 88, doi: 10.3847/1538-4357/ab05ce
  * Hummer et al. (2012) Hummer, S., Baerwald, P., & Winter, W. 2012, Phys. Rev. Lett., 108, 231101, doi: 10.1103/PhysRevLett.108.231101
  * Isravel et al. (2022) Isravel, H., Pe’er, A., & Begue, D. 2022.  https://arxiv.org/abs/2210.02363
  * J. Rika et al. (2022) J. Rika et al. 2022, GCN Circ. 32685
  * Kimura (2022) Kimura, S. S. 2022, arXiv e-prints, arXiv:2202.06480.  https://arxiv.org/abs/2202.06480
  * Kobayashi et al. (1997) Kobayashi, S., Piran, T., & Sari, R. 1997, Astrophys. J., 490, 92, doi: 10.1086/512791
  * Kumar & Zhang (2015) Kumar, P., & Zhang, B. 2015, Phys. Rep., 561, 1, doi: 10.1016/j.physrep.2014.09.008
  * Liu et al. (2022) Liu, R.-Y., Zhang, H.-M., & Wang, X.-Y. 2022, arXiv e-prints.  https://arxiv.org/abs/2211.14200
  * Margalit & Quataert (2021) Margalit, B., & Quataert, E. 2021, Astrophys. J. Lett., 923, L14, doi: 10.3847/2041-8213/ac3d97
  * Mirabal (2022) Mirabal, N. 2022, arXiv e-prints.  https://arxiv.org/abs/2210.14243
  * Murase et al. (2022) Murase, K., Mukhopadhyay, M., Kheirandish, A., Kimura, S. S., & Fang, K. 2022, arXiv e-prints.  https://arxiv.org/abs/2210.15625
  * Nigro et al. (2019) Nigro, C., et al. 2019, Astron. Astrophys., 625, A10, doi: 10.1051/0004-6361/201834938
  * Omodei et al. (2022) Omodei, N., Bruel, P., Bregeon, J., et al. 2022, GRB Coordinates Network, 32916, 1
  * Pillera et al. (2022) Pillera, R., Bissaldi, E., Omodei, N., et al. 2022, GRB Coordinates Network, 32658, 1
  * Price-Whelan et al. (2018) Price-Whelan, A. M., et al. 2018, Astron. J., 156, 123, doi: 10.3847/1538-3881/aabc4f
  * R. Pillera et al. (2022) R. Pillera et al. 2022, GCN Circ. 32658
  * Rees & Meszaros (1994) Rees, M. J., & Meszaros, P. 1994, Astrophys. J., 430, L93, doi: 10.1086/187446
  * Rees & Mészáros (1992) Rees, M. J., & Mészáros, P. 1992, Monthly Notices of the Royal Astronomical Society, 258, 41P, doi: 10.1093/mnras/258.1.41P
  * Ren et al. (2022) Ren, J., Wang, Y., & Zhang, L.-L. 2022.  https://arxiv.org/abs/2210.10673
  * Robitaille et al. (2013) Robitaille, T. P., et al. 2013, Astron. Astrophys., 558, A33, doi: 10.1051/0004-6361/201322068
  * Rudolph et al. (2022) Rudolph, A., Petropoulou, M., Bosnjak, Z., & Winter, W. 2022, in preparation.
  * S. Lesage et al. (2022) S. Lesage et al. 2022, GCN Circ. 32642
  * Siegel et al. (2021) Siegel, D. M., Agarwal, A., Barnes, J., et al. 2021, arXiv e-prints, arXiv:2111.03094.  https://arxiv.org/abs/2111.03094
  * The IceCube Collaboration (2022) The IceCube Collaboration. 2022, GCN Circ. 32665
  * Thompson et al. (2004) Thompson, T. A., Chang, P., & Quataert, E. 2004, ApJ, 611, 380, doi: 10.1086/421969
  * Usov (1992) Usov, V. V. 1992, Nature, 357, 472, doi: 10.1038/357472a0
  * Veres et al. (2022) Veres, P., Burns, E., Bissaldi, E., et al. 2022, GRB Coordinates Network, 32636, 1
  * Vietri (1995) Vietri, M. 1995, Astrophys. J., 453, 883, doi: 10.1086/176448
  * Warren et al. (2018) Warren, D. C., Barkov, M. V., Ito, H., Nagataki, S., & Laskar, T. 2018, Mon. Not. Roy. Astron. Soc., 480, 4060, doi: 10.1093/mnras/sty2138
  * Waxman (1995) Waxman, E. 1995, Phys. Rev. Lett., 75, 386, doi: 10.1103/PhysRevLett.75.386
  * Waxman & Bahcall (1997) Waxman, E., & Bahcall, J. N. 1997, Phys. Rev. Lett., 78, 2292, doi: 10.1103/PhysRevLett.78.2292
  * Woosley (1993) Woosley, S. E. 1993, ApJ, 405, 273, doi: 10.1086/172359
  * Yong Huang et al. (2022) Yong Huang et al. 2022, GCN Circ. 32677
  * Zhang et al. (2022) Zhang, B. T., Murase, K., Ioka, K., et al. 2022.  https://arxiv.org/abs/2211.05754


††thanks: These authors contributed equally to this work††thanks: These
authors contributed equally to this work

# Multiple magnetic transitions and complex magnetic structures in Fe2SiSe4
with the sawtooth lattice

Feihao Pan  Laboratory for Neutron Scattering and Beijing Key Laboratory of
Optoelectronic Functional Materials and MicroNano Devices, Department of
Physics, Renmin University of China, Beijing 100872, China  Xunwu Hu  State
Key Laboratory of Optoelectronic Materials and Technologies, Guangdong
Provincial Key Laboratory of Magnetoelectric Physics and Devices, Center for
Neutron Science and Technology, School of Physics, Sun Yat-sen University,
Guangzhou, 510275, China  Jiale Huang  Bingxian Shi  Jinchen Wang  Juanjuan
Liu  Hongxia Zhang  Daye Xu  Laboratory for Neutron Scattering and Beijing Key
Laboratory of Optoelectronic Functional Materials and MicroNano Devices,
Department of Physics, Renmin University of China, Beijing 100872, China
Hongliang Wang  Lijie Hao  China Institute of Atomic Energy, PO Box-275-30,
Beijing 102413, China  Peng Cheng<EMAIL_ADDRESS>Laboratory for Neutron
Scattering and Beijing Key Laboratory of Optoelectronic Functional Materials
and MicroNano Devices, Department of Physics, Renmin University of China,
Beijing 100872, China  Dao-Xin Yao<EMAIL_ADDRESS>State Key
Laboratory of Optoelectronic Materials and Technologies, Guangdong Provincial
Key Laboratory of Magnetoelectric Physics and Devices, Center for Neutron
Science and Technology, School of Physics, Sun Yat-sen University, Guangzhou,
510275, China

###### Abstract

The sawtooth lattice shares some structural similarities with the kagome
lattice and may attract renewed research interest. Here, we report a
comprehensive study on the physical properties of Fe2SiSe4, an unexplored
member in the olivine chalcogenides with the sawtooth lattice of Fe. Our
results show that Fe2SiSe4 is a magnetic semiconductor with band gap of 0.66
eV. It first undergoes an antiferromagnetic transition at Tm1=110 K, then an
ferrimagnetic-like one at Tm2=50 K and finally a magnetic transition at Tm3=25
K which is likely driven by the thermal populations of spin-orbit manifold on
the Fe site. Neutron diffraction analysis reveals a non-collinear
antiferromagnetic structure with propagation vector $\mathbf{q_{1}}$=(0,0,0)
at Tm2$\textless$T$\textless$Tm1. Interestingly, below Tm2, an additional
antiferromagnetic structure with $\mathbf{q_{2}}$=(0,0.5,0) appears and
Fe2SiSe4 exhibits a complex double-$\mathbf{q}$ magnetic structure which has
never been observed in sawtooth olivines. Density functional theory
calculations suggest this complex noncollinear magnetic structure may
originate from the competing antiferromagnetic interactions for both intra-
and inter-chain in the sawtooth lattice. Furthermore, band structural
calculations show that Fe2SiSe4 has quasi-flat band features near the valence
and conduction bands. Based on the above results, we propose Fe2SiSe4 as a new
material platform to condensed matter researches.

##  I Introduction

Materials with the sawtooth lattice have attracted attentions from the
condensed matter physics community for several reasons. Firstly, the sawtooth
antiferromagnetic (AFM) chain with corner-sharing triangles of spins
represents one of the fundamental models of geometrically frustrated quantum
magnetism as that in triangular and kagome latticesHagemann  _et al._ (2000);
Lau  _et al._ (2006). Secondly, the sawtooth lattice also exhibits flat-band
featureHuber and Altman (2010) which may give rise to high thermoelectric
performanceGudelli  _et al._ (2015), quantum topological phaseBlundell and
Núñez-Regueiro (2003); Yang  _et al._ (2020, 2021), flat-band spin dynamics
and phonon anomaliesGnezdilov  _et al._ (2019). Besides, a two-dimensional
kagome lattice can be viewed as the combination of one-dimensional sawtooth
chains. Therefore the recent discoveries of novel physical properties in
kagome materialLin  _et al._ (2018); Wang  _et al._ (2018); Ortiz  _et al._
(2020); Nie  _et al._ (2022); Yin  _et al._ (2021) may generate renewed
interest in sawtooth material since these two kinds share structural
similarities and connections.

The A2BX4 (A=Mn, Fe, Ni; B=Si, Ge; X=O,S,Se,Te) olivines represent a large
material family where the transitional-metal atoms form a sawtooth lattice.
Most members in this series are ferrimagnets or antiferromagnets whose spin
structure could be described by magnetic propagation vector
$\mathbf{q}$=(0,0,0)Ehrenberg and Fuess (1993); Baron  _et al._ (1999). Among
them, magnetic frustration is observed for Mn2SiSe4Nhalil  _et al._ (2019).
Mn2SiS4 and Mn2GeS4 are reported to exhibit anomalous magnetic properties
result from the quantum fluctuations near a spin-flop bicritical pointOhgushi
and Ueda (2005). Mn2GeO4 is identified as a functional material that exhibit
coupled ferromagnetism and ferroelectricityWhite  _et al._ (2012). Fe2GeS4 and
Fe2GeSe4 are theoretically proposed as promising candidates with good
thermoelectric performance due to the quasi-flat-band featureGudelli  _et al._
(2015). These findings highlight the exotic physical properties which are
closely related to the peculiar sawtooth lattice.

In A2BX4 series, we noticed that although the magnetic properties of Fe2SiO4
and Fe2SiS4 have been reportedAronson  _et al._ (2007); Baron  _et al._
(1999), Fe2SiSe4 remains an unexplored member whose chemical phase is even
missing in the current inorganic crystal structure database (ICSD). So it
would be interesting to check whether the Fe2SiSe4 phase with sawtooth lattice
really exists. In this paper, we report the synthesize of Fe2SiSe4 single
crystals which is characterized by sawtooth lattice of Fe. Using
magnetization, heat capacity, neutron scattering techniques and band structure
calculations, Fe2SiSe4 is identified as a magnetic semiconductor with multiple
magnetic transitions, non-collinear double-$\mathbf{q}$ magnetic structures
and quasi-flat-band. The underlying physical mechanism for the magnetic
properties and potential applications of Fe2SiSe4 are discussed combined with
the results of density functional theory (DFT) calculations.

##  II Experimental methods

Single crystals of Fe2SiSe4 were grown by chemical vapor transport method. The
pure powder of Fe, Si and Se were mixed in molar ratio 1:2:4 (total mass 1 g),
put into a quartz tube (inner diameter of 12 mm, length 10 cm) with 100 mg
iodine. Then the quartz tube was evacuated to 3.0$\times$10-3Pa and sealed
before put into a two-zone tube furnace. The quartz tube was heated to
590$\,{}^{\circ}\mathrm{C}$ in the raw material end and
660$\,{}^{\circ}\mathrm{C}$ in the other end in 750 minutes, then maintained
at the temperatures for 5760 minutes. The next step is a so-called
”temperature reversing process” in which the two ends switch temperatures in
70 minutes and held in the new temperatures for 12.5 days before cooled with
the furnace. The shining black single crystals of Fe2SiSe4 appeared in the
final cold end with typical dimensions of 1 mm$\times$2 mm$\times$0.5 mm.

X-ray diffraction (XRD) patterns of powder samples were collected from a
Bruker D8 Advance X-ray diffractometer using Cu Kα radiation. Magnetization
measurements were carried out in Quantum Design MPMS3. Resistivity and heat
capacity of the samples were measured on Quantum Design PPMS-14T. The powder
neutron diffraction experiments were carried out on Xingzhi cold neutron
triple-axis spectrometer at the China Advanced Research Reactor (CARR). For
neutron experiments on Xingzhi, the incident neutron energy was fixed at $16$
meV with a neutron velocity selector used upstream to remove higher order
neutronsCheng  _et al._ (2016). About 3 g Fe2SiSe4 powders (crushed from
single crystals) were used in neutron experiments. The program FullProf Suite
package was used in the Rietveld refinement of neutron powder diffraction
dataRodríguez-Carvajal (1993); Rietveld (1969).

The electronic structure and magnetic properties calculations were performed
using the DFT as implemented in the Vienna ab initio simulation package (VASP)
codeKresse and Hafner (1993); Kresse and Furthmüller (1996). The generalized
gradient approximation (GGA) in the form of Perdew-Burke-Ernzerhof (PBE)Perdew
_et al._ (1996) is used for exchange-correlation functional. The projector
augmented-wave (PAW) methodBlöchl (1994) with a 300 eV plane-wave cutoff
energy is employed. The valence electrons configurations for each atom are
$3d^{7}4s^{1}$ for Fe, $3s^{2}3p^{2}$ for Si, and $3d^{10}4s^{2}4p^{4}$ for
Se. A $\Gamma$-centered 3 $\times$ 2 $\times$ 5 k-points mesh within the
Monkhorst-Pack scheme is used for Brillouin zone sampling. The Hubbard U
Dudarev  _et al._ (1998) for the 3d electrons of Fe is chosen as 2 eV to
reproduce the experimental magnetic moments and bandgap of Fe2SiSe4. A 1
$\times$ 2 $\times$ 1 supercell is adopted due to the non-collinear AFM
structure. The direction of the magnetic moment is constrained, that is, the
superimposed double-$\mathbf{q}$ magnetic structure is adopted. All
calculations are performed using the experimental structural parameters.
Convergence criteria employed for both the electronic self-consistent
iteration and the ionic relaxation are set to $10^{-6}$ eV and 0.01
eV/$\mathring{A}$, respectively. Heisenberg exchange interactions are
calculated by the four-states methodXiang  _et al._ (2011). A 1 $\times$ 2
$\times$ 2 supercell is adopted to avoid spurious interaction by the periodic
boundary conditions.

Figure 1: (a) Crystal structure of Fe2SiSe4. (b) The x-ray reflections from
the $ab$-plane of Fe2SiSe4 single crystal. The inset shows the photo of one
crystal. (c) Room temperature XRD patterns of powders (crushed from single
crystals) and the Rietveld refinement results.

##  III Results and discussions

###  III.1 Crystal structure and magnetization

In Fig. 1(c), the XRD patterns on powder samples (crushed from single
crystals) and Rietveld refinement confirm Fe2SiSe4 adopts the orthorhombic
symmetry with space group $Pnma$ (No.62), which belongs to the olivine-type
structure, same as Fe2SiS4 and Fe2SiO4Ehrenberg and Fuess (1993); Baron  _et
al._ (1999). The obtained lattice parameters are $a=13.032$Å, $b=7.549$Å and
$c=6.123$Å. The x-ray reflection pattern from the $ab$-plane of Fe2SiSe4
single crystal in Fig. 1(b) also confirms the above results. As shown in Fig.
1(a), The Fe atoms form infinite sawtooth chains along the $b$-axis.

The Laue diffraction patterns allow us to distinguish the major crystal axes
with respect to the three-dimensional shape of Fe2SiSe4 single crystal. So the
temperature dependent magnetization along $a$-, $b$\- and $c$-axis were
performed respectively, the results are shown in Fig. 2. Three successive
magnetic transitions could be identified. At Tm1=110 K, a cusp appears which
is most prominent under $H\parallel b$ indicating an AFM transition. At lower
temperatures, a second magnetic transition occurs at Tm2=50 K. Its temperature
dependent magnetization behavior seems to be AFM-like under $H\parallel b$
while ferromagnetic-like along other field directions. Then at Tm3=25 K, the
magnetization anomalies along all three directions suggest the existence of a
third magnetic transition.

Figure 2: Temperature dependent magnetizations with applied field along the
$a$, $b$ and $c$-axis are shown respectively in (a), (b) and (c).

The temperature dependent magnetization shows strong anisotropic behavior
along three different directions. Even in the paramagnetic region, the
magnetization under $H\parallel a$ or $H\parallel c$ exhibits typical Curie-
Weiss (CW) behavior while that under $H\parallel b$ shows an anomalous linear
temperature dependence up to 400 K. The CW fit of the high-temperature
magnetization data yields $\mu_{eff}/Fe=4.02~{}\mu_{B}$ and $\theta_{CW}$=-66
K for $H\parallel a$, $\mu_{eff}/Fe=4.08~{}\mu_{B}$ and $\theta_{CW}$=-51 K
for $H\parallel c$. This result indicates dominant antiferromagnetic
correlations in Fe2SiSe4. On the other hand, the T-linear magnetization under
$H\parallel b$ is abnormal. Similar behavior has been observed in Fe-based
superconductors and parent compoundWang  _et al._ (2009); Yan  _et al._
(2008); Li  _et al._ (2018), which is explained as a result of AFM spin
fluctuations in the normal stateZhang  _et al._ (2009). Further theoretical
works are needed to check if this explanation applies to Fe2SiSe4.

###  III.2 Magnetic structure

Since the T-dependent magnetization indicates the existence of three possible
magnetic transitions, powder neutron diffraction on Fe2SiSe4 was performed at
150 K, 80 K, 35 K and 3.5 K to determine the magnetic structure. The emergence
of new peaks and large enhancement of the intensities for some nuclear peaks
could be observed in the neutron diffraction spectra at lower temperature
comparing with that at 150 K (Fig. S1). This allows us to identify plenty of
magnetic Bragg peaks and determine the magnetic structures at different
temperatures through Rietveld refinement. The detailed neutron diffraction
data, representational analysis and refinement process are presented in
Supplemental MaterialSup . The main results are shown in Fig. 3.

Figure 3: (a) Schematics of refined magnetic structure of Fe2SiSe4 at 80 K.
(b) The temperature-dependent intensities of magnetic peaks (0,0,1) and
(1,0.5,1). (c) The refined magnetic structures of Fe2SiSe4 at 3.5 K with
different propagation vectors are illustrated respectively.

First of all, as shown in Fig. 4(b), the intensities of Bragg peak (0,0,1)
could be served as the magnetic order parameter for the first AFM transition
at Tm1. Similarly, the temperature dependent intensity of (1,0.5,1) signifies
new magnetic order develops below Tm2. It should be noted there is almost no
difference for the diffraction patterns at 35 K and 3.5 K within our
instrumental resolution. Therefore the magnetic transition Tm3 may not change
the zero-field magnetic structure.

Secondly, at 80 K, the non-collinear AFM structure with magnetic wave-vector
$\mathbf{q_{1}}$=(0,0,0) is identified and illustrated in Fig. 3(a). The Fe
ions occupy two inequivalent Wyckoff positions. We label the ones at $4a$ site
as Fe1 (black spheres) and that at $4c$ sites as Fe2 (gray spheres). The
magnetic moments all lie in the $ab$-plane and the data refinement yields a
total moment of 1.45$\mu_{B}$/Fe for Fe1 and 4.04$\mu_{B}$/Fe for Fe2. The
latter points parallel to the $b$-axis while the former is canted to the
$a$-axis to a certain degree. This magnetic structure is quite similar as that
in olivine-type Fe2SiO4 and Fe2SiS4Ehrenberg and Fuess (1993); Baron  _et al._
(1999). Therefore similar physical interpretation of this magnetic structure
considering the indirect superexchange interactions between Fe cations via Se
atoms could be given. As marked in Fig. 3(a), the Fe2-Se-Fe2 angle is more
closer to 180∘ and gives rise to a strong AFM interaction, while the angle of
Fe1-Se-Fe1 is much smaller which may largely reduce the AFM interactionBaron
_et al._ (1999). A competition between crystal-field anisotropy and AFM
exchange via spin-orbit coupling results in the canting of Fe1
momentsEhrenberg and Fuess (1993).

Thirdly, at 3.5 K, the above magnetic model with $\mathbf{q_{1}}$=(0,0,0) can
only partially fit the neutron diffraction patterns and could not explain new
magnetic peaks appear below Tm2 (Fig. S1). The most prominent one is indexed
as (1,0.5,1) whose temperature dependent intensities are shown in Fig. 3(b).
The new magnetic peaks developed below Tm2 can be well defined by a new
propagation vector $\mathbf{q_{2}}$=(0,0.5,0). We found that only a combined
model which includes both the magnetic structure with $\mathbf{q_{1}}$ and
that with $\mathbf{q_{2}}$ could achieve a good fit to the data collected at
3.5 K. The refined two magnetic structures with different wave-vector are
illustrated in Fig. 3(c). One can see that the one with $\mathbf{q_{1}}$ is
quite similar as that at 80 K, except for slightly increased ordered moment
and different canting angle for Fe1 moment. The magnetic structure with
$\mathbf{q_{2}}$ is also non-collinear and more complex, the magnetic unit
cell is doubled along $b$-axis. The moments of Fe2 have two different values
and one of them is almost negligible (0.03 $\mu_{B}$).

The magnetic order at 3.5 K have two different propagation wave-vectors. These
two modulations may either reside in different domains independently, known as
the multi-domain state, or coexist in a single domain in the form of a
superimposed double-$\mathbf{q}$ state. These two states may be
indistinguishable in diffraction experiments performed on powder samples. The
multi-$\mathbf{q}$ magnetic order is not a common case in magnetic material.
Some famous examples include double-$\mathbf{q}$ spin-density wave in Fe-based
superconductor Sr1-xNaxFe2As2 with tetragonal crystal symmetryAllred  _et al._
(2016) and triple-$\mathbf{q}$ magnetic state in Na2Co2TeO6 with hexgonal
crystal symmetryChen  _et al._ (2021). For Fe2SiSe4, whether there is a
double-domain state or a double-$\mathbf{q}$ state, this state in a
orthorhombic crystal symmetry is particularly rarely observed.

To explore the underlying mechanism for the complex magnetic structure of
Fe2SiSe4 at 3.5 K, we investigate its exchange couplings based on a simple
Heisenberg model:

$H=\sum_{i<j}J_{ij}S_{i}S_{j}$ (1)

Our DFT calculations reveal the values of exchange coupling constant $J_{i}$
between Fe ions at different sites (illustrated in Fig. 3(c)) and the results
are listed in Table I. We can see that all of the $J_{i}$ values are positive,
which means each side of the sawtooth chains are AFM coupling. The intra-chain
exchange coupling $J_{1}$, $J_{2}$, $J_{6}$, and $J_{7}$ are connected by an
exchange path Fe-Si-Se-Fe with the distance d($J_{1}$) = 3.78 $\mathring{A}$.
The inter-chain exchange coupling $J_{3}$, $J_{4}$, $J_{5}$ are instead
realized through exchange path Fe-Se-Se-Fe with the distance d($J_{4}$) = 4.89
$\mathring{A}$. Furthermore, the exchange coupling values of inter-chain are
almost twice as large as those of intra-chain. This is probably due to the
strong hybridization of Se 4$p$ and Fe 3$d$ orbits near the Fermi level. A
very large Fe-Fe distance d($J_{8}$) = 7.17 $\mathring{A}$ results in a much
weaker magnetic interaction. Therefore, the competition of AFM interactions on
different sides of sawtooth chain combined with the inter-chain exchange
interaction may be responsible for a complex noncollinear double-$\mathbf{q}$
magnetic order in Fe2SiSe4. Based on this model, we employ a superimposed
double-$\mathbf{q}$ magnetic structure to constrain the direction of the
magnetic moment. The calculated spin configuration at low temperature of
Fe2SiSe4 is shown in Table S9 and reach qualitative consistent with the
superimposed double-$\mathbf{q}$ magnetic structure determined
experimentally(Table S8).

Table 1: Calculated magnetic exchange coupling constants between different Fe sites, positive value means AFM coupling. Constant | $J_{1}$ | $J_{2}$ | $J_{3}$ | $J_{4}$ | $J_{5}$ | $J_{6}$ | $J_{7}$ | $J_{8}$
---|---|---|---|---|---|---|---|---
Value (meV) | 8 | 9 | 15 | 12 | 15 | 9 | 8 | 0.3
Figure 4: Isothermal magnetizations under magnetic field applied along $a$\-
,$b$\- and $c$-axis at selected temperatures are shown respectively. Figure 5:
(a) Magnetic hysteresis under $H\parallel c$ for Fe2SiSe4 at selected
temperatures. (b) The temperature dependent ac susceptibilities measured under
an oscillated AC field of 5.0 Oe applied along $c$-axis.

###  III.3 Isothermal magnetization, AC susceptibility and heat capacity

Further measurements on Fe2SiSe4 single crystal were performed to get insights
on its magnetic state. As shown in Fig. 4, the isothermal magnetizations under
$H\parallel a$ and $H\parallel b$ are linear with increasing field for most
temperatures which are consistent with an AFM state. Nevertheless, for
$H\parallel c$, nonlinear M(H) curve appears below Tm2=50 K. Especially for 30
K and 40 K data, they firstly show a ferromagnetic-like fast magnetization to
about 0.13$\mu_{B}$/Fe at low field then increases linear with field. This
fast magnetization seems to gradually disappear below Tm3=25 K. The above
observations suggest the magnetic transition at Tm2 is a ferrimagnetic one.
The magnetic hysteresis behavior at 35 K in Fig. 5(a) and the peak in the
imaginary part of AC susceptibility $\chi^{\prime\prime}$ at Tm2 (Fig. 5(b))
further support this conclusion.

In previous section, the double-$\mathbf{q}$ magnetic structure below Tm2
identified from neutron diffraction is actually antiferromagnetic-like without
net moment component in any direction. This seems to contradict the
ferrimagnetic behavior observed in magnetization measurements. The explanation
might be, the neutron experiment is carried out in zero field which means that
the ground state without magnetic field may actually have no net
magnetization. In addition, the net magnetization 0.13$\mu_{B}$/Fe observed in
Fig. 3(c) is quite small. This small ferrimagnetic component is not easily
detected due to the neutron instrumental resolution. We also need to point
out, the moment of Fe has no component along $c$-axis at 80 K but Fe2 actually
has a $c$-axis component of 1.14 $\mu_{B}$ at 3.5 K which could be seen from
Table S7 and Table S8. This indicates the magnetic anisotropy has notable
change below Tm2.

Another question that needs to be answered is the nature of magnetic
transition at Tm3=25 K. The experimentally measured specific heat of Fe2SiSe4
is shown in Fig. 6. From the raw data Cp(T), the magnetic transitions at Tm1
and Tm2 can be identified by the jumps, while there is no detectable feature
at Tm3. Then if using an Einstein model-based curve to account for the phonon
part of the specific heat, we can obtain the magnetic specific heat Cmag shown
by the blue curve in Fig. 6. The Cmag(T) shows a shoulder at Tm3=25 K. We
found this should-like transition also exists in isostructural Fe2SiO4 at
similar temperature (20 K) and its origin has been uncovered as a Schottky
anomaly arising from the spin-orbit manifold of the Fe2+ ionAronson  _et al._
(2007). To be specific, for Fe2SiO4, Aronson $et$ $al$. found that the lowest
cystal-field splitting t2g level is further split into five states due to
spin-orbit coupling. These spin-orbit energy levels are determined through
inelastic neutron experiments and have dominate influence on the low
temperature physical properties. Then the shoulder in specific heat can be
quantitatively calculated and simulated by a Schottky anomaly arising from the
spin-orbit excitationsAronson  _et al._ (2007). This explanation may also
apply to the anomalous magnetic transition at Tm3 for Fe2SiSe4. As shown in
the inset of Fig. 6, above Tm1, the magnetic entropy of Fe2SiSe4 is close to
the limiting value $2Rln5$ which is calculated by assuming the spin-orbit
manifold is fully occupied and there are $2S+1=5$ equally occupied states for
each Fe2+. This gives further support to the above physical picture.

Figure 6: Temperature dependent specific heat of Fe2SiSe4. The magnetic
specific heat Cmag is derived by subtracting phonon contributions through an
Einstein model-based curve fit on the data above 110 K. Then the calculated
magnetic entropy is shown in the inset.

The thermal population in different spin-orbital states may drastically change
across Tm3, which lead to rich magnetic behaviors in Fe2SiSe4. As demonstrated
in Fig. 4 and Fig. 5, the magnetization switches from AFM to ferromagnetic
behavior then back to AFM behavior as deceasing temperature. In addition,
magnetic field induced spin-flop transitions are also observed for M(H) curve
at 15 K and 2 K under $H\parallel c$. Combined with the semiconductor nature
of this compound which will be discussed in the next section, Fe2SiSe4 may
find applications in optoelectronics and magnetic devices.

###  III.4 Band structure

Figure 7: Calculated band structure (a) and total and projected density of
states (b) of Fe2SiSe4.

The sawtooth lattice is known to have flat-bands similar as kagome lattice.
The band structure of Fe2SiSe4 obtained by DFT calculations is shown in Fig.
7. The result reveals Fe2SiSe4 is a semiconductor with indirect band gap of
0.7 eV. The gap value is close to that determined from experiment. As shown in
section III of Supplemental MaterialSup , the band gap is determined to be
0.40 eV through fitting the resistance curve and 0.66 eV by the diffusion
reflectance spectroscopy (DRS). Considering the errors of different methods,
0.66 eV might be a more accurate value of the band gap of Fe2SiSe4.

From Fig. 7(a), Fe2SiSe4 is found to have flat bands from the G-X
crystallographic direction, similar as other sawtooth material such as
Mn2SiSe4 and Fe2GeSe4Gudelli  _et al._ (2015); Nhalil  _et al._ (2019). The
corresponding total and projected density of states (DOS) for Fe2SiSe4 are
shown in Fig. 7(b). The states near the Fermi level mainly come from the
contributions from Fe 3$d$ states and Se 4$p$ states. Strong hybridization of
Se 4$p$ and Fe 3$d$ orbits could be observed below the Fermi level.

For kagome metals where the flat-band is quite close to the Fermi level, novel
properties including emergent ferromagnetism, anomalous Hall effect,
superconductivity and topological phases have been extensively reported in
recent yearsLin  _et al._ (2018); Wang  _et al._ (2018); Nayak  _et al._
(2016); Ortiz  _et al._ (2020); Yin  _et al._ (2019); Kang  _et al._ (2020);
Li  _et al._ (2021). We find the band gap of Fe2SiSe4 is not large, there are
possibilities that one could use pressure or chemical doping to tune it into a
metal. Then it would be worthwhile to check whether a quasi-flat-band near
Fermi level would enable realization of other versatile quantum phenomena.

##  IV Conclusions

In summary, we present a comprehensive study on the sawtooth lattice
chalcogenide olivines Fe2SiSe4. This material has shown intriguing magnetic
properties. Three magnetic transitions are identified as an AFM one at 110 K,
a ferrimagnetic one at 50 K and the last one at 25 K possibly due to the spin-
orbit excitation. We determined the magnetic structures at different
temperatures through neutron diffraction and discover a non-collinear
double-$\mathbf{q}$ magnetic order below 50 K. DFT calculations suggest this
complex magnetic structure may be due to the competition of AFM interactions
on different sides of sawtooth chain combined with the inter-chain exchange
interaction. Through band structure calculation and spectral experiment,
Fe2SiSe4 is identified as a magnetic semiconductor with indirect band gap of
0.66 eV and quasi-flat-band. We propose that Fe2SiSe4 may provide a new
material playground for further researches on magnetic devices and flat-band
effect through chemical doping.

## Acknowledgement

This work was supported by the National Key Research and Development Program
of China (No. 2018YFA0306001, 2022YFA1402802), the National Natural Science
Foundation of China (No. 12074426, No. 12004426, No. 11974432, No. 92165204),
Shenzhen International Quantum Academy (Grant No. SIQA202102), National Key
Scientific Instrument and Equipment Development Project of NSFC (No.
11227906), the Fundamental Research Funds for the Central Universities, and
the Research Funds of Renmin University of China (Grants No. 22XNKJ40), NSAF
(Grant No. U2030106).

## References

  * Hagemann  _et al._ (2000) I. S. Hagemann, P. G. Khalifah, A. P. Ramirez, and R. J. Cava, Phys. Rev. B 62, R771 (2000).
  * Lau  _et al._ (2006) G. C. Lau, B. G. Ueland, R. S. Freitas, M. L. Dahlberg, P. Schiffer, and R. J. Cava, Phys. Rev. B 73, 012413 (2006).
  * Huber and Altman (2010) S. D. Huber and E. Altman, Phys. Rev. B 82, 184502 (2010).
  * Gudelli  _et al._ (2015) V. K. Gudelli, V. Kanchana, and G. Vaitheeswaran, Journal of Physics: Condensed Matter 28, 025502 (2015).
  * Blundell and Núñez-Regueiro (2003) S. A. Blundell and M. D. Núñez-Regueiro, The European Physical Journal B - Condensed Matter and Complex Systems 31, 453 (2003).
  * Yang  _et al._ (2020) Y. Yang, S.-J. Ran, X. Chen, Z.-Z. Sun, S.-S. Gong, Z. Wang, and G. Su, Phys. Rev. B 101, 045133 (2020).
  * Yang  _et al._ (2021) H. Yang, H. Nakano, and H. Katsura, Phys. Rev. Res. 3, 023210 (2021).
  * Gnezdilov  _et al._ (2019) V. P. Gnezdilov, Y. G. Pashkevich, V. S. Kurnosov, O. V. Zhuravlev, D. Wulferding, P. Lemmens, D. Menzel, E. S. Kozlyakova, A. Y. Akhrorov, E. S. Kuznetsova, P. S. Berdonosov, V. A. Dolgikh, O. S. Volkova, and A. N. Vasiliev, Phys. Rev. B 99, 064413 (2019).
  * Lin  _et al._ (2018) Z. Lin, J.-H. Choi, Q. Zhang, W. Qin, S. Yi, P. Wang, L. Li, Y. Wang, H. Zhang, Z. Sun, L. Wei, S. Zhang, T. Guo, Q. Lu, J.-H. Cho, C. Zeng, and Z. Zhang, Phys. Rev. Lett. 121, 096401 (2018).
  * Wang  _et al._ (2018) Q. Wang, Y. Xu, R. Lou, Z. Liu, M. Li, Y. Huang, D. Shen, H. Weng, S. Wang, and H. Lei, Nature Communications 9, 3681 (2018).
  * Ortiz  _et al._ (2020) B. R. Ortiz, S. M. L. Teicher, Y. Hu, J. L. Zuo, P. M. Sarte, E. C. Schueller, A. M. M. Abeykoon, M. J. Krogstad, S. Rosenkranz, R. Osborn, R. Seshadri, L. Balents, J. He, and S. D. Wilson, Phys. Rev. Lett. 125, 247002 (2020).
  * Nie  _et al._ (2022) L. Nie, K. Sun, W. Ma, D. Song, L. Zheng, Z. Liang, P. Wu, F. Yu, J. Li, M. Shan, D. Zhao, S. Li, B. Kang, Z. Wu, Y. Zhou, K. Liu, Z. Xiang, J. Ying, Z. Wang, T. Wu, and X. Chen, Nature 604, 59 (2022).
  * Yin  _et al._ (2021) Q. Yin, Z. Tu, C. Gong, Y. Fu, S. Yan, and H. Lei, Chinese Physics Letters 38, 037403 (2021).
  * Ehrenberg and Fuess (1993) H. Ehrenberg and H. Fuess, Journal of Physics: Condensed Matter 5, 3663 (1993).
  * Baron  _et al._ (1999) V. Baron, Ã. Amcoff, and T. Ericsson, Journal of Magnetism and Magnetic Materials 195, 81 (1999).
  * Nhalil  _et al._ (2019) H. Nhalil, R. Baral, B. O. Khamala, A. Cosio, S. R. Singamaneni, M. Fitta, D. Antonio, K. Gofryk, R. R. Zope, T. Baruah, B. Saparov, and H. S. Nair, Phys. Rev. B 99, 184434 (2019).
  * Ohgushi and Ueda (2005) K. Ohgushi and Y. Ueda, Phys. Rev. Lett. 95, 217202 (2005).
  * White  _et al._ (2012) J. S. White, T. Honda, K. Kimura, T. Kimura, C. Niedermayer, O. Zaharko, A. Poole, B. Roessli, and M. Kenzelmann, Phys. Rev. Lett. 108, 077204 (2012).
  * Aronson  _et al._ (2007) M. Aronson, L. Stixrude, M. K. Davis, W. Gannon, and K. Ahilan, American Mineralogist 92, 481 (2007).
  * Cheng  _et al._ (2016) P. Cheng, H. Zhang, W. Bao, A. Schneidewind, P. Link, A. Grünwald, R. Georgii, L. Hao, and Y. Liu, Nucl. Instrum. Methods Phys. Res., Sect. A 821, 17 (2016).
  * Rodríguez-Carvajal (1993) J. Rodríguez-Carvajal, Phys. B 192, 55 (1993).
  * Rietveld (1969) H. M. Rietveld, J. Appl. Crystallogr. 2, 65 (1969).
  * Kresse and Hafner (1993) G. Kresse and J. Hafner, Phys. Rev. B 47, 558 (1993).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Phys. Rev. B 54, 11169 (1996).
  * Perdew  _et al._ (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
  * Blöchl (1994) P. E. Blöchl, Phys. Rev. B 50, 17953 (1994).
  * Dudarev  _et al._ (1998) S. L. Dudarev, G. A. Botton, S. Y. Savrasov, C. J. Humphreys, and A. P. Sutton, Phys. Rev. B 57, 1505 (1998).
  * Xiang  _et al._ (2011) H. J. Xiang, E. J. Kan, S.-H. Wei, M.-H. Whangbo, and X. G. Gong, Phys. Rev. B 84, 224429 (2011).
  * Wang  _et al._ (2009) X. F. Wang, T. Wu, G. Wu, H. Chen, Y. L. Xie, J. J. Ying, Y. J. Yan, R. H. Liu, and X. H. Chen, Phys. Rev. Lett. 102, 117005 (2009).
  * Yan  _et al._ (2008) J.-Q. Yan, A. Kreyssig, S. Nandi, N. Ni, S. L. Bud’ko, A. Kracher, R. J. McQueeney, R. W. McCallum, T. A. Lograsso, A. I. Goldman, and P. C. Canfield, Phys. Rev. B 78, 024516 (2008).
  * Li  _et al._ (2018) X.-G. Li, J.-M. Sheng, C.-K. Tian, Y.-Y. Wang, T.-L. Xia, L. Wang, F. Ye, W. Tian, J.-C. Wang, J.-J. Liu, H.-X. Zhang, W. Bao, and P. Cheng, Europhysics Letters 122, 67006 (2018).
  * Zhang  _et al._ (2009) G. M. Zhang, Y. H. Su, Z. Y. Lu, Z. Y. Weng, D. H. Lee, and T. Xiang, Europhysics Letters 86, 37006 (2009).
  * (33) See Supplemental Material at [URL will be inserted by publisher] for details about structural and composition analysis of NaCr(Te1-xSex)2 through single-crystal x-ray diffraction and EDS. .
  * Allred  _et al._ (2016) J. M. Allred, K. M. Taddei, D. E. Bugaris, M. J. Krogstad, S. H. Lapidus, D. Y. Chung, H. Claus, M. G. Kanatzidis, D. E. Brown, J. Kang, R. M. Fernandes, I. Eremin, S. Rosenkranz, O. Chmaissem, and R. Osborn, Nature Physics 12, 493 (2016).
  * Chen  _et al._ (2021) W. Chen, X. Li, Z. Hu, Z. Hu, L. Yue, R. Sutarto, F. He, K. Iida, K. Kamazawa, W. Yu, X. Lin, and Y. Li, Phys. Rev. B 103, L180404 (2021).
  * Nayak  _et al._ (2016) A. K. Nayak, J. E. Fischer, Y. Sun, B. Yan, J. Karel, A. C. Komarek, C. Shekhar, N. Kumar, W. Schnelle, J. Kubler, C. Felser, and S. S. P. Parkin, Science Advances 2, e1501870 (2016).
  * Yin  _et al._ (2019) J.-X. Yin, S. S. Zhang, G. Chang, Q. Wang, S. S. Tsirkin, Z. Guguchia, B. Lian, H. Zhou, K. Jiang, I. Belopolski, N. Shumiya, D. Multer, M. Litskevich, T. A. Cochran, H. Lin, Z. Wang, T. Neupert, S. Jia, H. Lei, and M. Z. Hasan, Nature Physics 15, 443 (2019).
  * Kang  _et al._ (2020) M. Kang, L. Ye, S. Fang, J.-S. You, A. Levitan, M. Han, J. I. Facio, C. Jozwiak, A. Bostwick, E. Rotenberg, M. K. Chan, R. D. McDonald, D. Graf, K. Kaznatcheev, E. Vescovo, D. C. Bell, E. Kaxiras, J. van den Brink, M. Richter, M. Prasad Ghimire, J. G. Checkelsky, and R. Comin, Nature Materials 19, 163 (2020).
  * Li  _et al._ (2021) M. Li, Q. Wang, G. Wang, Z. Yuan, W. Song, R. Lou, Z. Liu, Y. Huang, Z. Liu, H. Lei, Z. Yin, and S. Wang, Nature Communications 12, 3129 (2021).