# The formation of hard VHE spectra from GRB afterglow via Two-Zone
Synchrotron Self-Compton Emission

Dmitry Khangulyan Graduate School of Artificial Intelligence and Science,
Rikkyo University, Nishi-Ikebukuro 3-34-1, Toshima-ku, Tokyo 171-8501, Japan
Andrew M. Taylor DESY, D-15738 Zeuthen, Germany  Felix Aharonian Dublin
Institute for Advanced Studies, School of Cosmic Physics, 31 Fitzwilliam
Place, Dublin 2, Ireland  Max-Planck-Institut für Kernphysik, Saupfercheckweg
1, 69117 Heidelberg, Germany

###### Abstract

Electron Compton scattering of target photons into the gamma-ray energy band
(inverse Compton scattering –IC–) is commonly expected to dominate the very
high energy spectra in gamma-ray bursts especially during the afterglow phase.
For sufficiently large center-of-mass energies in these collisions, the effect
of the electron recoil starts reducing the scattering cross section (the
Klein-Nishina regime). The IC spectra generated in the Klein-Nishina regime is
softer and has a smaller flux level compared to the synchrotron spectra
produced by the same electrons. The detection of afterglow emission from
nearby GRB190829A in the very high energy (VHE) domain with H.E.S.S. has
revealed an unexpected feature: the slope of the VHE spectrum matches well the
slope of the X-ray spectra, despite expectations that for the IC production
process, the impact of the Klein-Nishina effect should be strong. The multi-
wavelength spectral energy distribution appears to be inconsistent with
predictions of one-zone synchrotron-self-Compton models. We study the possible
impact of two-zone configuration on the properties of IC emission when the
magnetic field strength differs considerably between the two zones.
Synchrotron photons from the strong magnetic field zone provide the dominant
target for cooling of the electrons in the weak magnetic field zone, which
results in a formation of hard electron distribution and consequently of a
hard IC emission. We show that the two-zone model can provide a good
description of the X-ray XRT and VHE H.E.S.S. data.

Non-thermal radiation sources(1119) — Gamma-ray transient sources(1853) —
Gamma-ray bursts(629) — Gamma-ray astronomy(628) — Particle astrophysics(96) —
X-ray sources(1822)

##  1 Introduction

The very high energy (VHE; $>100\,\mathrm{GeV}$) emission detected from gamma-
ray burst (GRB) afterglows with H.E.S.S. and MAGIC (Abdalla et al., 2019;
MAGIC Collaboration et al., 2019a, b; H. E. S. S. Collaboration et al., 2021)
is considered by many to have inverse Compton (IC) origin (see, e.g, Zhang,
2019). The emission component produced by relativistic protons is expected to
have a significantly lower flux, due to the very low radiative efficiency of
hadronic interactions (see, e.g., Abdalla et al., 2019). If the VHE emission
is produced by relativistic electrons, then because of the so-called
synchrotron burn-off limit (Guilbert et al., 1983) the synchrotron component
is expected to reach the VHE regime only if the bulk Lorentz factor is very
high, $\Gamma\geq 10^{3}$. Such high bulk Lorentz factors are excluded during
the afterglow phase by energy conservation arguments (e.g., related to self-
similar solution for relativistic blast wave obtained by Blandford & McKee,
1976) making IC scattering the most feasible radiation mechanism for the VHE
GRB emission during the afterglow period. However, the hard intrinsic spectral
slope inferred from observations by H.E.S.S. of GRB190829A afterglow cannot be
easily reproduced with standard IC models (see, e.g., H. E. S. S.
Collaboration et al., 2021). This leaves one of two possibilities: (i) invoke
alternative radiation mechanisms, or (ii) develop a more sophisticated IC
scenario to provide a better description of the observational data.

Synchrotron radiation is a very efficient radiative emission mechanism of
electrons during the afterglow phase of GRBs. If the synchrotron component
extends into the VHE domain, it can reproduce the flux level and spectral
slope revealed with H.E.S.S. from GRB190829A afterglow (H. E. S. S.
Collaboration et al., 2021). While the conservation of energy, used to
constrain the bulk Lorentz factor, is a robust argument, the burn-off energy
limit can be avoided in certain non-standard scenarios. For example, if the
strength of the accelerating electric field, ${\cal E}$, exceeds the strength
of the magnetic field, $B$ (in a plasma such configurations require non-ideal
magnetohydrodynamics) then synchrotron emission can extend beyond the burn-off
limit by the factor of ${\cal E}/B$. Alternatively, in highly turbulent
magnetic fields magnetobremsstrahlung radiation can extend beyond the burn-off
limit (Kelner et al., 2013). If the correlation length of the magnetic field
is large compared to the photon formation length, $m_{e}c^{2}/e\bar{B}$ (here
$m_{e}$ and $e$ are electron mass and charge, respectively; $c$ and $\bar{B}$
are light speed and averaged magnetic field), then the radiation is generated
in the synchrotron regime, resulting in the burn-off limit for the synchrotron
maximum energy (for a detailed consideration, see, e.g., in Kelner et al.,
2013; Derishev & Aharonian, 2019). However, if the correlation length is short
compared to the photon formation length, then the electrons instead emit in
the jitter regime, and the emission peaks at higher energy compared to the
synchrotron case, alleviating the limit from the burn-off limit (Kelner et
al., 2013). Finally, the electron synchrotron spectrum can extend beyond the
burn-off limit in two-zone systems, where the physical conditions at the
acceleration site and in the radiation production region differ substantially
(Kumar et al., 2012; Khangulyan et al., 2021). In conclusion, there are
several ways of expanding the energy spectrum of magneto-bremsstrahlung to
high or even very high energies. However, the feasibility of these scenarios
depends on the implementation of many factors and requires extreme
assumptions.

In contrast, IC scattering is a natural and very effective channel of VHE
gamma-ray production. Although the recent observations of VHE gamma rays
during the GRB afterglows challenge the simple one-zone IC model, more
sophisticated scenarios cannot be excluded. In this paper, we study the
spectral properties of gamma rays in the two-zone IC model in which the
production region of the target (synchrotron) photons and the IC gamma-ray
emitter are separated. One can propose several possible realizations for such
a two-zone setup. For example, one may expect quite different conditions at
the forward and reverse shocks, which propagate through the circumburst medium
(CBM) and the jet, respectively. If the emission from the reverse shock
appears to be important at certain frequencies then a two-zone description for
GRB afterglow emission should be considered (see, e.g., Dichiara et al., 2022;
Salafia et al., 2022). Alternatively, the shock region itself can be quite
complex potentially providing quite different physical conditions for particle
acceleration and radiation. Indeed, simulations suggest that downstream shock
material, the dominant emission site during the afterglow phase, is expected
to be highly inhomogeneous, an aspect usually neglected in GRB afterglow
emission modelling. Below we consider the impact of a strongly inhomogeneous
magnetic field on the properties of IC emission. We show that under reasonable
assumptions, even a two-zone synchrotron self-Compton (SSC) scenario can
provide a considerably improved description of the broadband spectra reported
from GRB190829A.

##  2 Standard One-Zone SSC Scenario

The standard GRB afterglow emission framework postulates that this emission is
generated via the synchrotron and IC channels, with synchrotron radiation
providing the dominant target for IC scattering – the so called SSC scenario.
The analysis of the spectral energy distribution (SED) in SSC models is
straightforward if the IC emission is generated in the Thomson regime (see,
e.g. Sari & Esin, 2001), as in this case the energy loss rate, $\dot{E}$, has
a simple form $\dot{E}\propto E^{2}$ (here $E$ is electron energy). In this
regime, a power-law injection of non-thermal electrons, $q\propto E^{-\alpha}$
(here $\alpha$ is the injection index, for conventional acceleration
mechanisms one typically assumes $\alpha\approx 2$), leads to the formation of
a broken-power-law distribution of radiating electrons. The synchrotron and IC
(Thomson) components generated by these electrons also reflect this broken-
power-law shape, with the IC component dominating at higher energies. The
subsequent broadband SED produced is double-humped, with the relative
emissivity of the synchrotron and IC components being determined by
phenomenological parameters (typically, by the radiation efficiency, i.e., by
the fraction of energy radiated away). The photon index of the synchrotron
spectrum, produced by electrons with energies above the cooling break, is
$\gamma_{\mathrm{s}}=(\alpha+2)/2$, provided that $\alpha>1$. In the single
zone SSC scenario, the corresponding IC spectrum has the same photon index, if
generated in the Thomson regime.

Figure 1: A sketch that illustrates the formation of the particle spectrum in
the case of dominant synchrotron losses and dominant IC losses. The part of
the spectrum formed in the fast cooling regime is shown.  Figure 2: A sketch
that illustrates the formation of the SED in the case of dominant synchrotron
losses and dominant IC losses.

Typically, during the afterglow phase the (synchrotron) X-ray spectrum is
observed to be hard, with a photon index $\sim 2$. Thus the photons detected
in the X-ray band provide a non-negligible target for IC scattering. In the
plasma co-moving frame, the energy of the electron, $E$, generating the VHE
emission, detected at energy111Note that we prime the quantities in the
progenitor frame, and we neglect the cosmological redshift effect.
$\varepsilon_{\mathrm{vhe}}^{\prime}$, satisfy the condition:
$E>\varepsilon_{\mathrm{vhe}}^{\prime}/\Gamma$. If electrons of this energy
up-scatter photons from a component detected by the observer at energy
$\varepsilon_{\mathrm{x}}^{\prime}$, then the typical product of the target
photon and electron energies, which determines the scattering regime, is

$\frac{E\varepsilon_{\mathrm{x}}}{m_{e}^{2}c^{4}}>\frac{\varepsilon_{\mathrm{vhe}}^{\prime}\varepsilon_{\mathrm{x}}^{\prime}}{m_{e}^{2}c^{4}\Gamma^{2}}\approx
4\quantity(\frac{\Gamma}{10})^{-2}\left(\frac{\varepsilon_{\mathrm{vhe}}^{\prime}}{0.1\,\mathrm{TeV}}\right)\left(\frac{\varepsilon_{\mathrm{x}}^{\prime}}{1\,\mathrm{keV}}\right)\,.$
(1)

Here $m_{e}$ and $c$ are the electron mass and speed of light, respectively.
Unless the bulk Lorentz factor is high, $\Gamma\geq 10^{2}$, the electrons
that produce the VHE emission up-scatter a considerable part of the photon
targets in the Klein-Nishina regime. The study of the VHE properties of GRB
afterglows should therefore be conducted with models that account for the
change of the IC cross-section in the relativistic regime.

The influence of the Klein-Nishina regime on the SED is two-fold, as one must
account for both the change of the emission and energy loss rates (see, e.g.,
Derishev et al., 2003; Nakar et al., 2009). In the fast cooling regime, the
particle spectrum, $\differential{N}=n\differential{E}$, is determined by the
injection spectrum, $q$, and by the cooling time $\tau=E/|\dot{E}|$:

$n(E)=\frac{\tau(E)}{E}\int\limits_{E}^{\infty}\differential{\tilde{E}}\,q\quantity(\tilde{E})\,.$
(2)

If the injection is a power-law $q\propto E^{-\alpha}$, then the particle
spectrum is

$n(E)\propto\tau(E)E^{-\alpha}\,.$ (3)

(Note that here we assume that the injection spectrum is sufficiently steep so
as to ensure the integral is dominated by the low energy limit).

If the synchrotron losses dominate over the Compton losses (more specifically
if the energy density of the magnetic field is larger than the energy density
of the target photons) then $\tau(E)\propto E^{-1}$, and a power-law spectral
injection also yields a power-law distribution of particles: $n(E)\propto
E^{-(\alpha+1)}$ (see Fig. 1 for a sketch of the cooled particle spectrum).
Subsequently, a power-law synchrotron component is produced with photon index
$\gamma_{\mathrm{s}}$.

The inverse Compton of radiation has the same power-law photon index as long
as the scattering takes place in the Thompson regime. In the Klein-Nishina
regime, the IC slope should (asymptotically, i.e., ignoring the logarithmic
term) approach $\gamma_{\mathrm{kn}}\approx(\alpha+2)$ (provided that the
emitting electrons obey a power-law energy distribution with index $\alpha+1$,
Blumenthal & Gould, 1970). Thus, since the slope of IC component generated in
the Thomson regime matches that of the synchrotron radiation,
$\gamma_{\mathrm{s}}$, the Klein-Nishina effect causes a spectral softening by
$\Delta\gamma\approx\gamma_{\mathrm{kn}}-\gamma_{\mathrm{s}}\approx(\alpha+2)/2$.
For example, if $\alpha\approx 2$ then the spectral slope changes from
$\gamma_{\mathrm{s}}\approx 2$ to $\gamma_{\mathrm{kn}}\approx 4$, and the
spectral softening is $\Delta\gamma\approx 2$. A schematic of the SED is shown
in Fig. 2. One should note that for a broad target photon distribution, the
transition to the Klein-Nishina regime is spread over a broad energy range and
can have a rather complex character.

The situation changes dramatically when the energy density of target photons
is larger than the energy density of the magnetic field. In this case, the
impact of the Klein-Nishina effect on the formation of the electron spectrum
becomes a dominant factor. The radiative cooling time $\tau(E)$ can be
approximated by a broken power-law function: for sufficiently low electron
energies, the IC interaction proceeds in the Thomson regime, thus
$\tau(E)\propto E^{-1}$. At higher energies, the IC interactions occur in the
Klein-Nishina regime where the energy loss rate is energy-independent, thus
$\tau(E)\propto E$. Finally at even higher energies, denoted $E_{*}$, the
synchrotron losses (as their rate increases with particle energy) begin to
dominate over the IC energy losses, and the original energy dependence of the
cooling time is recovered: $\tau(E)\propto E^{-1}$. As follows from Eq. (3),
for a power-law injection spectrum, the particle spectrum formed in the fast
cooling regime should also be a double-broken-power-law (with the power-law
index changing as $\alpha+1\rightarrow\alpha-1\rightarrow\alpha+1$: see Fig.
1). The $E^{-(\alpha+1)}$ part of the spectrum formed under dominant (Thomson
regime) IC losses changes to, $\propto E^{-(\alpha-1)}$, formed under the
dominant IC (Klein-Nishina regime) losses. Finally, above $E_{*}$, the
spectrum softens back to $E^{-(\alpha+1)}$. We note, however, that the
transition to the Klein-Nishina regime proceeds smoothly, therefore the
spectrum does not follow precisely the schematic shape explained above. For
example, as can be seen from Fig. 4, the IC cooling time in the transition
regime is better approximated as a constant, $\tau\approx\mathrm{const}$.
Therefore, the corresponding transformation of the electron spectrum is better
approximated as $\alpha+1\rightarrow\alpha\rightarrow\alpha+1$ (note that this
power law index is indicated in the bottom panel of Fig. 4 with a black guide
line).

As for the synchrotron radiation, electrons cooled by IC in the Thomson regime
produce a spectrum with photon index $\gamma_{\mathrm{s}}$; at higher energies
the hardening of the electron spectrum due to the dominant Klein-Nishina
energy losses results in a hard synchrotron spectrum with photon index in the
range between $\gamma_{\mathrm{s}}$ and
$\gamma_{\mathrm{s,kn}}\approx\alpha/2$ ($\gamma_{\mathrm{s,kn}}$ is the
limiting value achieved under IC cooling in the deep Klein-Nishina regime: see
Fig. 2). In the transition region with an approximately constant IC cooling
time, the slope of the synchrotron spectrum is approximately $(\alpha+1)/2$,
as indicated by the black guide lines in Figs. 5 and 6. Finally, the emission
produced by electrons with energies exceeding $E_{*}$ has the standard
synchrotron slope $\gamma_{\mathrm{s}}$. As the synchrotron and IC energy loss
rates for particles with $E_{*}$ are equal, the narrow-band luminosity of the
synchrotron and IC components produced by particles with $E_{*}$ are (almost)
equal.

The spectral shape of the IC component is different to that of the synchrotron
spectrum. The component generated in the Thomson regime has a spectral index
of $\gamma_{\mathrm{s}}$. At higher energies the impact of the Klein-Nishina
effect on the particle spectrum is partially compensated by the reduction of
the cross section. For example, in the limiting regime, a spectrum $\propto
E^{-(\alpha-1)}$ generates in the Klein-Nishina regime a $E^{-\alpha}$ IC
spectrum. For $\alpha\approx 2$ a Thomson spectrum with photon index
$(\alpha+2)/2$ transits smoothly into the Klein-Nishina spectrum with photon
index $\alpha$. However, in the region of transition to the Klein-Nishina
regime, this asymptotic photon index might be quite a coarse approximation.
Moreover, above $E_{*}$ the synchrotron losses dominate, thus the Klein-
Nishina spectrum eventually softens to $\alpha+2$ above $E_{*}$. Note that in
the Klein-Nishina regime almost all the electron energy is transferred to the
up-scattered photon, so the photon energy in the co-moving frame is equal to
that of the incident electron energy, $\varepsilon_{\mathrm{ic}}\approx
E_{*}$.

Observations of GRB190829A with H.E.S.S. revealed that VHE component,
corrected for the extragalactic background light (EBL) attenuation, is best
described as a single power-law spectrum extending up to $3\,\mathrm{TeV}$
with a hard photon index of $\gamma_{\mathrm{vhe}}=2.07\pm 0.09$ (H. E. S. S.
Collaboration et al., 2021). Strikingly, this slope matches well the slope of
the X-ray spectrum measured with Swift-XRT (e.g.,
$\gamma_{\mathrm{xrt}}=2.03\pm 0.06$ during the first night H. E. S. S.
Collaboration et al., 2021). Also, the Swift-XRT and H.E.S.S. observations
revealed that the fluxes in the X-ray and VHE bands appeared to be similar
(potentially a natural feature of pair loading feedback, see Derishev & Piran,
2016, 2019, for detail).

In the VHE band the influence of the Klein-Nishina effect should be
noticeable. However, this spectral effect was not observed in the H.E.S.S.
measurements. In the framework of the simple one-zone analysis introduced
above, the slope and flux level match implies that the cooling of TeV emitting
electrons proceeds in the Klein-Nishina regime, and that the X-ray synchrotron
is produced by particles with energy exceeding $E_{*}$. As the hard VHE
spectrum extends up to $3\,\mathrm{TeV}$, then
$E_{*}>0.3\quantity(\Gamma/10)^{-1}\,\mathrm{TeV}$. The synchrotron emission
produced by the high-energy electrons is detected by the observer at

$\varepsilon_{\mathrm{syn}}^{\prime}>60\quantity(\frac{\Gamma}{10})^{-1}\frac{B}{1\,\mathrm{G}}\,\mathrm{keV}\,.$
(4)

This estimate shows that a very low magnetic field of $\sim\,\mathrm{mG}$
level is required by the VHE measurements. Such a low magnetic field, however,
is incompatible with the required radiation efficiency of the production
region given the adiabatic cooling time is $\tau_{\mathrm{ad}}\sim
t_{\mathrm{tr}}^{\prime}\Gamma$, where $t_{\mathrm{tr}}^{\prime}$ is time
since the GRB trigger (as measured by a distant observer at rest in the
progenitor reference frame). The broad-band SED obtained with Swift-XRT and
H.E.S.S. therefore cannot be reproduced in the framework of the standard one-
zone SSC scenario (see also Huang et al., 2022). To resolve the spectral issue
in SSC scenario one needs either: (1) assume that there is an important low-
energy target photon field, probably of external origin; or (2) consider a
two-zone scenario.

The former scenario requires the presence of an external target that provides
a target of an energy density comparable to that of the magnetic field in the
plasma co-moving frame:

$w_{\mathrm{ext}}\sim 4\times
10^{-2}\quantity(\frac{B}{1\,\mathrm{G}})^{2}\,\mathrm{erg\,cm^{-3}}\,.$ (5)

If the photons are isotropic in the progenitor frame, then we obtain
$w_{\mathrm{ext}}^{\prime}\sim 4\times
10^{-4}\quantity({10B}/{(\Gamma\,\mathrm{G}}))^{2}\,\mathrm{erg\,cm^{-3}}$.
The VHE emission detected from GRB190829A lasted for almost $\Delta t=50$ h
(H. E. S. S. Collaboration et al., 2021), and the forward shock covered a
distance of $\Delta R^{\prime}\sim\Gamma^{2}\Delta tc\sim
10^{17}\quantity(\Gamma/10)^{2}\,\mathrm{cm}$. The luminosity of the photon
field should therefore be

$L_{\mathrm{ext}}^{\prime}\sim 4\pi\Delta
R^{\prime}{}^{2}w_{\mathrm{ext}}^{\prime}c\sim
10^{42}\quantity(\frac{B}{1\,\mathrm{G}})^{2}\,\mathrm{erg\,s^{-1}}\,.$ (6)

If the magnetic field is weak, $B\ll 1\,\mathrm{G}$, then an external photon
field of reasonable luminosity can provide a sufficiently dense external
photon field (see, e.g., Zhang et al., 2021), however external IC scenarios
with an equivalent Gauss-strength magnetic field cannot be realized.

##  3 Two-Zone SSC Emission Scenario

###  3.1 Physical justification

Figure 3: Examples of three different geometries that allow the scenario
realization. Scenario (i): two distinct regions with typical sizes of $r_{1}$
and $r_{2}$ separated by a distance $r_{0}$; scenario (ii): two converging
shells of radius $r_{1}$ and $r_{0}$; scenario (iii): a large number of
compact regions (of typical size $r_{1}$) with strong magnetic field embedded
within a larger zone of size $r_{0}$.

We consider the emission region consisting of two zones: the first zone with a
strong magnetic field, $B_{1}$, and the second zone with a weak magnetic
field, $B_{2}$, with $B_{1}\gg B_{2}$. Should particles themselves also easily
mix between the two zones, then one would not expect a significant difference
between the energy distributions of particles in these zones. We here,
however, assume that the particle exchange between the zones is inefficient,
and thus two distinct particle distributions, $n_{1}$ and $n_{2}$, are formed
in the two zones.

The target photons, however, travel freely between the two zones. The specific
realization of the scenario, in particular the shapes and relative location of
the zones, determines the actual distribution of target photons in the zones.
Let us qualitatively consider several possible realizations of the two-zone
scenario: (i) two distinct regions with typical sizes of $r_{1}$ and $r_{2}$
separated by a distance $r_{0}$; (ii) two converging shells of radius $r_{1}$
and $r_{0}$; (iii) $N$ compact regions (of typical size $r_{1}$) with strong
magnetic field embedded within a larger zone of size $r_{0}$. These three
possibilities are shown in Fig. 3. Although less apparent, scenario (iii) is
two-zone in the sense that the physical conditions and processes are the same
in the compact regions, and differ substantially from those in the larger
zone.

The synchrotron luminosity of each of the zones is $L_{1}$ and $L_{2}$,
respectively. In scenario (iii) we define $L_{1}$ as the total luminosity of
$N$ regions of enhanced B-field. We consider a situation $L_{1}\gg L_{2}$.
Thus, when considering the processes in the first zone, we can ignore the
photons supplied by the second zone. The energy density of the locally
generated photons in the first zone is

$w_{1\rightarrow 1}\sim\frac{L_{1}}{r_{1}^{2}Nc}\,,$ (7)

where $N=1$ for scenarios (i) and (ii). Equation (7) ignores a numerical
factor, which depends on the production region geometry and the distribution
of emitting particles. For example, in the case of a spherical homogeneous
production region, the volume average energy density of target photons is
given by Eq. 7 with a factor $9/(16\pi)$ (for detail, see in Atoyan &
Aharonian, 1996). We note that such factors do not affect our conclusions, we
therefore safely ignore them.

In the second zone one needs to account for the contribution of locally
generated photons:

$w_{2\rightarrow 2}^{(i)}\sim\frac{L_{2}}{r_{2}^{2}c}\mbox{\quad
and\quad}w_{2\rightarrow 2}^{(ii)/(iii)}\sim\frac{L_{2}}{r_{0}^{2}c}$ (8)

and the photons supplied from the first zone, $w_{1\rightarrow 2}$. For the
each of the above defined geometries one obtains

$w_{1\rightarrow 2}\sim\frac{L_{1}}{r_{0}^{2}c}\,.$ (9)

The suggested scenario is realized if the photon field produced in the first
zone (being locally a subdominant) provides the dominant target for the
particle cooling in the second zone:

$w_{1\rightarrow 1}\ll\frac{B_{1}^{2}}{8\pi}\mbox{\quad
and\quad}w_{1\rightarrow 2}\gg\frac{B_{2}^{2}}{8\pi}\,.$ (10)

The photon field in the second zone is diluted compared to the first zone:
$w_{1\rightarrow 1}>w_{1\rightarrow 2}$, thus the scenario requires that
$B_{1}\gg B_{2}$. The difference of the magnetic fields determines the
dilution of the photon field, $\kappa=w_{1\rightarrow 2}/w_{1\rightarrow 1}$,
that allows the scenario realization (i.e, the conditions given by Eq. 10).

The possible ratio of the magnetic fields should be determined by the physical
arguments unique to each specific realization of the scenario. However, from
the general point of view, it is obvious that if the photon field is
significantly diluted in the second zone, $\kappa\ll 1$, the required
difference between the magnetic field strength becomes larger, making the
realization of the scenario less feasible (although not excluded). For
example, a strong dilution might be expected in scenario (i) provided that
$r_{0}\gg r_{1}$. In contrast, in scenario (ii) the dilution of the photon
field in the second zone is small, by a factor of $\sim 2$, provided that two
shells are of comparable radius, $r_{1}\approx r_{0}$. Similarly, in scenario
(iii) one obtains

$\frac{w_{1\rightarrow 2}}{w_{1\rightarrow
1}}\sim\frac{r_{1}^{2}N}{r_{0}^{2}}=\frac{fr_{0}}{r_{1}}\,,$ (11)

where $f$ is filling factor. If the above ratio is not small (i.e.,
${w_{1\rightarrow 2}}/{w_{1\rightarrow 1}}\gtrsim 1$) then the photon field is
nearly homogeneous in the entire production region, i.e., $\kappa\approx 1$.

For the sake of simplicity we will consider a single common photon target
being present in the two zones. In the first place, this seems to be a
perfectly suitable choice for scenarios (ii) and (iii) if $r_{1}\approx r_{0}$
and $fr_{0}/r_{1}\gtrsim 1$, respectively. Even if these conditions are not
fulfilled, the model calculations should reproduce correctly the part of SED
formed in the fast cooling regime (provided that IC losses dominate over the
synchrotron cooling in the second zone: $w_{1\rightarrow 2}\gg
B_{2}^{2}/(8\pi)$).

Although the scenario can be realized also in scenario (i), if the magnetic
field in the second zone is sufficiently weak to remain subdominant compared
to the significantly diluted photon field provided from the first zone,
scenarios (ii) and (iii) seem to be less demanding. In particular, these
geometries can be formed during the afterglow phase of GRBs. The shells
assumed in scenario (ii) may correspond to the reverse and forward shocks.
Also an onion-like structure may be formed in the inner part of the forward
shock downstream region, where the competing processes of magnetic field
amplification and decay may lead to the formation of a layer with an enhanced
magnetic field. If the magnetic field amplification in the downstream proceeds
in a highly non-homogeneous manner, then instead of a shell-like structure one
should expect rather a large number of magnetized blobs in the production
region, i.e., scenario (iii). Although scenarios (ii) and (iii) are
characterized by quite similar geometries, the angular distribution of the
target photons in the second zone may be quite different in these two cases.
While in scenario (iii), the target photons are nearly isotropic, scenario
(ii) features a substantial anisotropy of the target photons in the second
zone (as depicted in Fig. 3). As the emitting particles are isotropized in the
plasma frame, this photon anisotropy should not have any impact on the cooling
process. However, one may need to account for anisotropic IC cross-section
(see, e.g., Aharonian & Atoyan, 1981) for accurate computation of the IC
spectra. For example, if the emission generated in the direction of the
observer is predominately produced by scattering target photons at small
scattering angles, then the IC spectra appear to be harder compared to the
spectra computed with angle-averaged IC cross-section (see, e.g., Khangulyan
et al., 2008).

Because of the Doppler boosting effect, the observer can detect the emission
coming from a patch of the shell with a typical size of $R^{\prime}/\Gamma$,
where $R^{\prime}\sim t^{\prime}_{\mathrm{tr}}\Gamma^{2}c$. Thus, one obtains
the patch size as $t^{\prime}_{\mathrm{tr}}\Gamma c\gg 10^{15}\,\mathrm{cm}$
(provided that $t^{\prime}_{\mathrm{tr}}>1\,\mathrm{h}$ for the afterglow
period). The realization of scenario (iii) requires that the size of the blobs
is small, $r_{1}\ll t^{\prime}_{\mathrm{tr}}\Gamma c$. Verification of this
condition from the first principles may require detailed plasma simulations,
which are beyond the scope of this study. As in the case of GRB afterglow, the
GeV emission seems to belong to the same component as the synchrotron, we may
therefore speculate that the acceleration in the blobs are limited by the
synchrotron cooling and the acceleration process is efficient,
$\eta_{\mathrm{acc}}\sim 1$ (here $\eta_{\mathrm{acc}}$ is the acceleration
efficiency). Thus, the size of the blobs should be sufficiently large to
confine particles with energy

$E\approx
60\left(\frac{B}{1\,\mathrm{G}}\right)^{\nicefrac{{-1}}{{2}}}\,\mathrm{TeV}\,.$
(12)

The corresponding gyro radius,

$R_{\cal{G}}\approx 2\times
10^{11}\left(\frac{B}{1\,\mathrm{G}}\right)^{\nicefrac{{-3}}{{2}}}\,\mathrm{cm}\,,$
(13)

is significantly smaller than the patch size. This likely implies that there
are no fundamental constrains from the plasma physics forbidding the scenario
realization.

###  3.2 Mathematical setup

For each of the zones we consider the injection-cooling equation:

$\partialderivative{n_{1}}{t}+\partialderivative{\dot{E}_{1}n_{1}}{E}=q_{1}(E)-\frac{n_{1}}{\tau_{1\rightarrow
2}}+\frac{n_{2}}{\tau_{2\rightarrow 1}}\,,$ (14)
$\partialderivative{n_{2}}{t}+\partialderivative{\dot{E}_{2}n_{2}}{E}=q_{2}(E)+\frac{n_{1}}{\tau_{1\rightarrow
2}}-\frac{n_{2}}{\tau_{2\rightarrow 1}}\,.$ (15)

Here $q_{i}$ is the injection term; $\dot{E}_{i}$ is the energy loss rate in
each zone; and $\tau_{i\rightarrow j}$ is the probability of particle escape
from zone $i$ to zone $j$. To illustrate the possible impact of the two-zone
setup on the IC spectrum we assume that in the energy range of interest
$\tau_{i\rightarrow j}\gg\tau_{1},\tau_{2}$, where $\tau_{i}=E/|\dot{E}_{i}|$
is radiative cooling time in zone $i$. Since we are interested in the high-
energy part of the spectrum, which is formed in the fast cooling regime, we
consider the following simplified equations:

$\partialderivative{\dot{E}_{i}n_{i}}{E}=q_{i}(E)\,.$ (16)

The magnetic field strengths differ significantly in each zone, thus we do not
adopt a universal cutoff energy in the injection spectrum, but instead find
different injection rates, $A_{i}$, and cutoff energies, $E_{{\rm cut},i}$,
within the two zones. We, however, assume that within both zones the injection
function has a common power-law spectral index:

$q_{i}=A_{i}E^{-\alpha}\exp[-\frac{E}{E_{\mathrm{cut}}{}_{,i}}]\,,$ (17)

where the cutoff energy is found through the balance of acceleration and loss
time-scales:

$\frac{\eta_{\mathrm{acc}}E_{\mathrm{cut}}{}_{,i}}{eB_{i}c}=-\frac{E_{\mathrm{cut}}{}_{,i}}{\dot{E}_{i}(E_{\mathrm{cut}}{}_{,i})}\,.$
(18)

Here $e$ is electron charge. Whilst the acceleration parameter
$\eta_{\mathrm{acc}}$ is assumed to be the same in both zones, the energy
losses are computed independently for each zone.

We assume that dimensionless parameters $\kappa_{1}$ and $\kappa_{2}$
($\kappa_{1}+\kappa_{2}=1$) define the fraction of the total energy injected
into zone one and two, respectively:

$\frac{1}{\kappa_{i}}\int\limits_{E_{\mathrm{min}}}^{\infty}\differential{E}Eq_{i}(E)=L_{0}\,,$
(19)

where $L_{0}$ is the total power injected in the production region. The
minimum energy, $E_{\mathrm{min}}$, we set to a value of $\Gamma m_{e}c^{2}$
(here $m_{e}$ is electron mass and $c$ is the speed of light). While Eq. (18)
defines the cutoff energy, Eq. (19) determines the normalization coefficients
in each zone, $A_{i}$.

As in the high-energy regime, the synchrotron and IC losses are expected to
provide the dominant energy loss channels, we therefore only take account of
these two energy loss mechanisms. Since the photon field is common between the
zones, the difference between the energy loss rate in the zones is due to the
different synchrotron losses within each zone:

$\dot{E}_{i}=\dot{E}_{\mathrm{syn}}{}_{,i}+\dot{E}_{\mathrm{ic}}\,.$ (20)

The synchrotron energy losses in zone $i$ are determined by the magnetic field
strength

$\dot{E}_{\mathrm{syn}}{}_{,i}=-\frac{16\pi}{3}\frac{e^{4}E^{2}}{m_{e}^{4}c^{7}}\quantity(\frac{2}{3}\frac{B_{i}^{2}}{8\pi})\,.$
(21)

Note that the equation above is averaged over pitch angle. The synchrotron
cooling time is

$\tau_{\mathrm{syn}}{}_{,i}=\frac{E}{|\dot{E}_{\mathrm{syn}}{}_{,i}|}\approx
400\quantity(\frac{B_{i}}{1\,\mathrm{G}})^{-2}\quantity(\frac{E}{1\,\mathrm{TeV}})^{-1}\,\mathrm{s}\,.$
(22)

The IC losses are determined by the energy distribution and number density of
target (synchrotron) photons. As photons can freely cross the zone boundaries
we assume that the photon distribution is the same throughout the entire
production region, i.e., it includes the contributions from both zones. We
compute the synchrotron emission using the particle distribution in each zone
and the corresponding magnetic field:

$\derivative{N_{\mathrm{syn}}{}_{,i}}{\varepsilon\differential{t}}\approx\int\limits_{m_{e}c^{2}}^{\infty}\differential{E}n_{i}K_{\mathrm{syn}}{}_{,\varepsilon}\quantity(E,B_{i})\,,$
(23)

where $\varepsilon$ is the target (synchrotron) photon energy. For the
synchrotron integral kernel, $K_{\mathrm{syn}}{}_{,\varepsilon}$, we use a
simple analytic approximation for the pitch angle averaged synchrotron
spectrum (for detail see in Aharonian et al., 2010). Finally, we compute the
energy distribution of the target photons as

$\derivative{N_{\mathrm{syn}}}{\varepsilon\differential{V}}\approx\frac{1}{R^{2}c}\quantity(\derivative{N_{\mathrm{syn}}{}_{,1}}{\varepsilon\differential{t}}+\derivative{N_{\mathrm{syn}}{}_{,2}}{\varepsilon\differential{t}})\,.$
(24)

Here $R$ is size of the production region.

The rate of IC scattering is determined by the angle averaged scattering cross
section (for detail see in Jones, 1968):

$\derivative{\nu_{\mathrm{ic}}}{\varepsilon_{\gamma}\differential{\varepsilon}}=\frac{8\pi
cr_{0}^{2}}{bE}\derivative{N_{\mathrm{syn}}}{V\differential{\varepsilon}}\times$
$\left[1+\frac{z^{2}}{2(1-z)}+\frac{z}{b(1-z)}-\frac{2z^{2}}{b^{2}(1-z)^{2}}-\right.$
$\left.\frac{z^{3}}{2b(1-z)^{2}}-\frac{2z}{b(1-z)}\,\ln\frac{b(1-z)}{z}\right]\,.$
(25)

Here $r_{0}=e^{2}/m_{e}c^{2}$ is the electron classical radius; the Klein-
Nishina parameter is given by $b=4\varepsilon E/(m_{e}^{2}c^{4})$; and $z$ is
the ratio of the up-scattered photon to electron energy,
$z=\varepsilon_{\gamma}/E$. The IC energy loss rate depends on the energy
distribution of target photons as

$\dot{E}_{\mathrm{ic}}\approx\int\limits_{0}^{\infty}\differential{\varepsilon}\int\limits_{\varepsilon_{\mathrm{min}}{}_{,\gamma}}^{\varepsilon_{\mathrm{max}}{}_{,\gamma}}\differential{\varepsilon_{\gamma}}(\varepsilon-\varepsilon_{\gamma})\derivative{\nu_{\mathrm{ic}}}{\varepsilon_{\gamma}\differential{\varepsilon}}\,,$
(26)

where the maximum/minimum energy of up-scattered gamma-ray,
$\varepsilon_{\mathrm{max/min}}{}_{,\gamma}$, is determined by kinematic
constraints. If electrons up-scatter low-energy target photons (i.e., the
Klein-Nishina parameter is small, $b\ll 1$), then the IC energy loss rate
depends only on the energy density of the target photons, $w_{\mathrm{ph}}$:

$\dot{E}_{\mathrm{T}}{}_{,i}=-\frac{32\pi}{9}\frac{e^{4}E^{2}}{m_{e}^{4}c^{7}}w_{\mathrm{ph}}\,,$
(27)

analogous to the corresponding angle averaged energy loss rate in a magnetic
field given in Eq. (21).

###  3.3 Model calculations

For the model calculations, magnetic field values of $B_{1}=1\,\mathrm{G}$ and
$B_{2}=10^{-3}\,\mathrm{G}$ are assumed. The injection power is set to $\sim
10^{39}\,\mathrm{erg\,s^{-1}}$, and for the size of the production region we
consider a value close to $10^{16}\,\mathrm{cm}$. If one considers this size
in the context of a GRB afterglow, one should compare it to the forward shock
radius, which depends on the time passed since the trigger,
$t_{\mathrm{tr}}^{\prime}$:

$R\sim\Gamma^{2}t_{\mathrm{tr}}^{\prime}c\sim 3\times
10^{16}\quantity(\frac{\Gamma}{10})^{2}\frac{t_{\mathrm{tr}}^{\prime}}{3\,\mathrm{h}}\,\mathrm{cm}\,.$
(28)

The typical energy density of the target photons in the production region is

$w_{\mathrm{ph}}\sim 4\times
10^{-5}\kappa_{1}\eta_{\mathrm{rad}}\quantity(\frac{R}{3\times
10^{16}\,\mathrm{cm}})^{-2}\,\mathrm{erg\,cm^{-3}}\,,$ (29)

where $\eta_{\mathrm{rad}}$ is the radiation efficiency in zone 1 (in what
follows we ignore this factor, setting $\eta_{\mathrm{rad}}=1$, for the sake
of simplicity). This energy density corresponds to an equivalent magnetic
field strength of

$B_{\mathrm{eq}}\sim 3\times
10^{-2}\kappa_{1}^{\nicefrac{{1}}{{2}}}\quantity(\frac{R}{3\times
10^{16}\,\mathrm{cm}})^{-1}\,\mathrm{G}\,.$ (30)

This photon field is the dominant target in zone 2, whereas it is negligible
in zone 1. The corresponding cooling time scales are shown in Fig. 4 (top
panel). Whilst at high energies (approaching $1\,\mathrm{PeV}$), the Klein-
Nishina losses approach their asymptotic energy-dependence,
$\tau_{\mathrm{kn}}\propto E$, for the parameter set considered, the particles
cool in the transition regime with $\tau\propto\mathrm{const}$. Thus the
spectrum formed is not as hard as expected from our earlier qualitative
analysis.

The effect of the onset of Klein Nishina cooling on the electron spectrum is
shown in Fig. 4 (bottom panel), where the energy distribution of electrons in
both zones are shown. For the calculations here we adopted the following
parameters: linear size $R=10^{16}\,\mathrm{cm}$; total power of acceleration
of non-thermal particles $L_{0}=10^{39}\mathrm{erg\,s^{-1}}$, which is
distributed between the zones with $\kappa_{1}=0.90$ and $\kappa_{2}=0.10$;
the injection index $\alpha=2.2$ (the “main case”). Finally, the acceleration
efficiency was set to $\eta_{\mathrm{acc}}=10^{2}$, for which the cutoff
energy in zone 1 is determined to be:

$E_{\mathrm{cut}}{}_{,1}\approx
6~{}\left(\frac{\eta_{\mathrm{acc}}}{10^{2}}\right)^{\nicefrac{{-1}}{{2}}}\quantity(\frac{B_{1}}{1\,\mathrm{G}})^{\nicefrac{{-1}}{{2}}}\,\mathrm{TeV}\,.$
(31)

For this acceleration efficiency the cutoff energy in zone 2 is at $\approx
200\,\mathrm{TeV}$, which is close to the energy at which the synchrotron
losses dominate over the IC losses, $E_{*}\approx 20\,\mathrm{TeV}$, thus the
influence of the high energy cutoff becomes prominent at energies just above
the Klein-Nishina hardening energy scale.

The energy dependence of the electron distribution is directly reflected in
the synchrotron spectrum from zone 2. As can be seen from Fig. 5, this
component is subdominant to the luminous synchrotron component from zone 1.
The photon index of the hardest part of the spectrum is $(\alpha+1)/2\approx
1.5$, which is considerably softer than the limiting photon index of
$\gamma_{\mathrm{s,kn}}(=\alpha/2)$. This is caused by the smooth broad
transition to the Klein-Nishina regime. While the broad transition from the
Thomson to Klein-Nishina regimes causes the electron distribution to be not as
hard as naively expected, the IC component appears to be somewhat harder than
in the limiting case. As can be seen in Fig. 5, a power-law component extends
from a few $\,\mathrm{GeV}$ to beyond $10\,\mathrm{TeV}$ with a photon index
of $\approx(\alpha+1)/2$. Note that for our calculations we set $\alpha=2.2$,
and the production region bulk Lorentz factor was assumed to be $\Gamma=10$.

Figure 4: Top panel: Synchrotron, IC cooling time together with the
acceleration time. Bottom panel: Electron distribution in two zones. Black
guide lines indicate power-law approximations. Figure 5: Spectral energy
distribution of synchrotron and IC emission from two zones. Black guide lines
show the power-law approximations.

To illustrate the influence of the model parameters, we performed calculations
for a range of different parameter sets. The results of these calculations are
shown in Fig. 6. For the “case A” we adopted a different value for the
injection index: $\alpha=2$ instead of $\alpha=2.2$ used in the “main case”.
For the “case B” we adopted a different value for the acceleration efficiency:
$\eta_{\mathrm{acc}}=10^{4}$ instead of $\eta_{\mathrm{acc}}=10^{2}$ used in
the “main case”. The adopted model parameter values are summarized in Table 2.

Figure 6: Spectral energy distribution of synchrotron and IC emission from two
zones. Black guide lines show the power-law approximations. Top panel: Case A;
Bottom panel: Case B. Thin lines in the top panel correspond to a case when
the electron distribution in the first zone features a cooling break at
$E\approx 10\,\mathrm{GeV}$.

Low-energy target photons can be an important role in the formation of a hard
VHE spectrum in the case of the conventional one-zone SSC models. To
demonstrate the relatively small influence of low-energy target photons in the
framework of our considered two-zone approach, in the top panel of Fig. 6 we
also plot the SED obtained under the same conditions assuming that the
particle spectrum in the first zone features a cooling break at $E\approx
10\,\mathrm{GeV}$. The corresponding spectra are shown with thin lines. Under
this assumption, the synchrotron spectrum from the first zone features the
cooling breaks, as expected. The IC spectrum is more strongly suppressed: one
sees here the impact of both the cooling break and reduced target photon
density. The reduction of the IC loss rate leads to a considerable enhancement
of the synchrotron emission from the second zone (note that this component
still remains subdominant). The IC spectra from the second zone shows,
however, only minor changes, noticeably only close to the high- and low-energy
cutoffs regions. This quite weak influence of the target photon spectrum on
the spectral properties of the IC component from the second zone is caused by
the fact that the IC losses determine the particle spectrum, as we assume that
the emission is generated in the fast cooling regime. Therefore, the electron
spectrum adjusts to the rate of the dominant losses, and the spectral
properties of the IC component are largely determined by the injection
spectrum.

##  4 Discussion and Conclusion

The need for studying energy losses in the inhomogeneous emission region
downstream can be easily realized by considering the evolution of the magnetic
field from the upstream to downstream regions. Based on the hydrodynamics of
the forward shock propagating through the CBM, one can obtain the following
estimate for the downstream magnetic field strength:

$B\sim 3\times 10^{2}\frac{\Gamma}{10}\frac{B_{\mathrm{cbm}}}{10\,\mathrm{\mu
G}}\,\mathrm{\mu G}\,.$ (32)

This estimate depends on the typical strength of the CBM magnetic field,
$B_{\mathrm{cbm}}$, and accounts for the transformation of this field to the
forward shock rest frame, and for the increase of the field strength at a
weakly magnetized relativistic shock due to the shock compression.

The magnetic field given by Eq. (32) appears significantly below the Gauss-
level required for the afterglow radiation. Therefore, one needs to assume an
efficient magnetic field amplification process, which can increase the energy
density of the magnetic field to the level comparable to the plasma energy
density in the downstream:

$w\sim n_{\mathrm{cbm}}m_{p}\Gamma^{2}\approx
0.15\quantity(\frac{n_{\mathrm{cbm}}}{1\,\mathrm{cm^{3}}})\quantity(\frac{\Gamma}{10})^{2}\,\mathrm{erg\,cm^{-3}}\,,$
(33)

where $n_{\mathrm{cmb}}$ is CBM density. This estimates shows that the
magnetic field in the downstream can be amplified up to a strength of

$B_{\mathrm{eq}}=\sqrt{8\pi w}\sim
2\quantity(\frac{n_{\mathrm{cbm}}}{1\,\mathrm{cm^{3}}})^{\nicefrac{{1}}{{2}}}\quantity(\frac{\Gamma}{10})\,\mathrm{G}\,.$
(34)

Gauss-strength magnetic fields in the afterglow production region are also
favored on theoretical grounds by afterglow emission modeling. If the magnetic
field is indeed amplified by a factor of $\sim 10^{3}$, it is natural to
further assume that this amplification is inhomogeneous throughout the volume
resulting in a magnetic field configuration with strong spatial fluctuations.
For example, magnetic field amplification by turbulent dynamo shows that the
magnetic energy is predominantly localized in small blobs (Zhang et al.,
2009). Moreover, this may be a general effect: the field amplification
predominately operates on small scale fields (Kazantsev, 1968).

The highly inhomogeneous structure of the downstream region can have important
implications for the properties of the non-thermal emission generated. In
particular, such a structure in the production region can significantly alter
the synchrotron radiation emission, with clumps of highly amplified magnetic
field leading to the synchrotron emission extending significantly beyond the
one-zone synchrotron burn-off limit (Khangulyan et al., 2021). This scenario
requires that particles are accelerated in a region of weak magnetic field,
and subsequently penetrate into a second zone of amplified magnetic field,
where they rapidly cool producing VHE synchrotron radiation. The requirement
of effective particle exchange between the two zones of strong and weak
magnetic field is an important element of this scenario.

It should be noted, however, that efficient particle exchange between the
zones is a significant assumption. Processes exist, which can hinder particle
exchange between the two zones. For example, if the change of the magnetic
field strength is relatively smooth, the magnetic adiabatic invariant prevents
particles from the zone of weak magnetic field reaching a strong magnetic
field zone (see the discussion in Khangulyan et al., 2021). The particle
escape from the zone of strong to weak magnetic field is not forbidden by the
magnetic adiabatic invariant, but it seems feasible that one can neglect this
process. Because of the much higher rate of synchrotron losses in the strong
magnetic field zone, the total number of particles in this zone is naturally
significantly reduced to that in the weak magnetic field zone, particularly
for the highest energy particles with energies close to the maximum energy.

Figure 7: Spectral energy distribution of synchrotron and IC emission from two
zones. Filled zones show X-ray (XRT) and VHE intrinsic (H.E.S.S.) spectra (for
detail see in H. E. S. S. Collaboration et al., 2021).

On the other hand, synchrotron photons can freely travel between the two
zones. The photon exchange between the zones have two major effects: (i)
altering particle energy losses, and (ii) change the properties of IC
emission. We suggest a simple model that allows one to study these two
effects. We find that for feasible model parameters, IC scattering dominates
the cooling process in the zone of weak magnetic field. Due to the Klein-
Nishina effect, the particle spectrum formed in the fast cooling regime
appears to be significantly harder than the spectrum formed for the case when
synchrotron losses dominate. While the synchrotron emission from this zone may
appear completely sub-luminous with respect to the synchrotron emission
generated in the strong magnetic field zone, the IC component from the weak
magnetic field zone would be expected to dominate.

The second signature of the hard particle spectrum is expected in the IC
component generated by these particles. This spectrum appears to be hard, with
a photon index coinciding with the value expected for the synchrotron/Thomson
spectra, $(\alpha+2)/2$ (where $\alpha$ is the injection index). The IC
spectrum therefore appears to have the same slope as the dominant synchrotron
emission. The relative flux through these two channels is determined by the
phenomenological parameters, $\kappa_{i}$, which determine the ratio of the
acceleration powers in the two zones. Our simulations presented in Fig. 7 show
that for an acceleration spectrum with a spectral slope of $\alpha=2.1$, which
allows the slope of the X-ray spectra for GRB190829A to be reproduced, the
X-ray and VHE flux ratio seen in GRB190829A implies $\kappa_{1}=0.7$ and
$\kappa_{2}=0.3$ (see column “GRB190829A” in Table 2). This suggests that
acceleration processes of comparable power operate in the both zones. However,
the acceleration in the zone of stronger magnetic field is somewhat more
efficient. The obtained hard VHE IC spectrum extends beyond $10\,\mathrm{TeV}$
for a modest bulk Lorentz factor of $\Gamma=10$. This implies that a hard
multi-TeV IC component can be generated also during the late afterglow phases,
when the forward shock transits into the mildly relativistic regime. During
the prompt or early afterglow phases, when the bulk Lorentz factor can be
significantly larger, $\Gamma\geq 100$, the intrinsically hard IC component
can extend up to the ultra high energy domain ($\geq 100\,\mathrm{TeV}$).
However, we note that the extragalactic EBL attenuation is severe already in
the VHE domain for even the most local GRB redshift values.

Table 1: Used notations Parameter | Notation | Comment
---|---|---
Production region size | $R$ | Similar to the forward shock radius
Magnetic field | $B$ |  $B_{1,2}$ in zone 1 or 2; $B_{\mathrm{cbm}}$ in the shock upstream
Bulk Lorentz factor | $\Gamma$ | We do not distinguish the bulk Lorentz factors of emitting plasma or forward shock
Injection power | $L_{0}$ | Eq. (19) is in the co-moving frame, but note that $L_{0}$ is a Lorentz invariant
Normalization factors | $\kappa_{1,2}$ | Power distribution between two zone; $\kappa_{1}+\kappa_{2}=1$
Energy density | $w$ | e.g., $w_{\mathrm{ext}}$ is energy density of external photon fields
Electron energy | $E$ | –
Electron energy density | $n$ | i.e., $\differential{N}=n\differential{E}$
Injection index | $\alpha$ | see, e.g., Eq. 17
Injection rate | $q$ | i.e., $\differential{N}=q\differential{E}\differential{t}$
Acceleration efficiency | $\eta_{\mathrm{acc}}$ | see, e.g., Eq. (18)
Cooling time | $\tau$ | Synchrotron, IC, and adiabatic cooling or escape
Energy loss rate | $\dot{E}$ | Synchrotron, IC, and adiabatic
Photon energy | $\varepsilon$ | For photons produced through the synchrotron or IC channels
Photon index | $\gamma$ | In particular, $\gamma_{\mathrm{s}}$, $\gamma_{\mathrm{kn}}$
Radiation efficiency | $\eta_{\mathrm{rad}}$ | We assume $\eta_{\mathrm{rad}}=1$
Circumburst medium density | $n_{\mathrm{cmb}}$ | i.e. $\differential{N}_{\mathrm{cmb}}=n_{\mathrm{cmb}}\differential{V}$
Table 2: Parameter values used for model calculations Parameter | Notation | units | Main case | Case A | Case B | GRB190829A
---|---|---|---|---|---|---
|  |  | in Figs. 4&5 | in Fig. 6 | in Fig. 6 | in Fig. 7
Size | $R$ | $\,\mathrm{cm}$ | $10^{16}$ | $10^{16}$ | $10^{16}$ | $10^{16}$
Bulk Lorentz factor | $\Gamma$ | – | $10$ | $10$ | $10$ | $10$
Luminosity | $L_{0}$ | $\mathrm{erg\,s^{-1}}$ | $10^{39}$ | $10^{39}$ | $10^{39}$ | $2.2\times 10^{38}$
Strong magnetic field | $B_{1}$ | $\,\mathrm{G}$ | $1$ | $1$ | $1$ | $1$
Weak magnetic field | $B_{2}$ | $\,\mathrm{G}$ | $10^{-3}$ | $10^{-3}$ | $10^{-3}$ | $10^{-3}$
Zone 1 power fraction | $\kappa_{1}$ | – | $0.9$ | $0.9$ | $0.9$ | $0.7$
Zone 2 power fraction | $\kappa_{2}$ | – | $0.1$ | $0.1$ | $0.1$ | $0.3$
Injection slope | $\alpha$ | – | 2.2 | 2 | 2.2 | 2.1
Acceleration efficiency | $\eta_{\mathrm{acc}}$ | – | $10^{2}$ | $10^{2}$ | $10^{4}$ | $10^{2}$
Radiation efficiency | $\eta_{\mathrm{rad}}$ | – | $1$ | $1$ | $1$ | $1$
Distance | $D$ | $\,\mathrm{Mpc}$ | – | – | – | $400$

DK acknowledges support by the RSF grant No. 21-12-00416.

## References

  * Abdalla et al. (2019) Abdalla, H., Adam, R., Aharonian, F., et al. 2019, Nature, 575, 464, doi: 10.1038/s41586-019-1743-9
  * Aharonian et al. (2010) Aharonian, F. A., Kelner, S. R., & Prosekin, A. Y. 2010, Phys. Rev. D, 82, 043002, doi: 10.1103/PhysRevD.82.043002
  * Aharonian & Atoyan (1981) Aharonian, F. A., & Atoyan, A. M. 1981, Ap&SS, 79, 321, doi: 10.1007/BF00649428
  * Atoyan & Aharonian (1996) Atoyan, A. M., & Aharonian, F. A. 1996, MNRAS, 278, 525, doi: 10.1093/mnras/278.2.525
  * Blandford & McKee (1976) Blandford, R. D., & McKee, C. F. 1976, Physics of Fluids, 19, 1130, doi: 10.1063/1.861619
  * Blumenthal & Gould (1970) Blumenthal, G. R., & Gould, R. J. 1970, Reviews of Modern Physics, 42, 237, doi: 10.1103/RevModPhys.42.237
  * Derishev & Aharonian (2019) Derishev, E., & Aharonian, F. 2019, ApJ, 887, 181, doi: 10.3847/1538-4357/ab536a
  * Derishev & Piran (2019) Derishev, E., & Piran, T. 2019, ApJ, 880, L27, doi: 10.3847/2041-8213/ab2d8a
  * Derishev et al. (2003) Derishev, E. V., Kocharovsky, V. V., Kocharovsky, V. V., & Mészáros, P. 2003, in American Institute of Physics Conference Series, Vol. 662, Gamma-Ray Burst and Afterglow Astronomy 2001: A Workshop Celebrating the First Year of the HETE Mission, ed. G. R. Ricker & R. K. Vanderspek, 292–294, doi: 10.1063/1.1579362
  * Derishev & Piran (2016) Derishev, E. V., & Piran, T. 2016, MNRAS, 460, 2036, doi: 10.1093/mnras/stw1175
  * Dichiara et al. (2022) Dichiara, S., Troja, E., Lipunov, V., et al. 2022, MNRAS, 512, 2337, doi: 10.1093/mnras/stac454
  * Guilbert et al. (1983) Guilbert, P. W., Fabian, A. C., & Rees, M. J. 1983, MNRAS, 205, 593, doi: 10.1093/mnras/205.3.593
  * H. E. S. S. Collaboration et al. (2021) H. E. S. S. Collaboration, Abdalla, H., Aharonian, F., et al. 2021, Science, 372, 1081, doi: 10.1126/science.abe8560
  * Huang et al. (2022) Huang, Z.-Q., Kirk, J. G., Giacinti, G., & Reville, B. 2022, ApJ, 925, 182, doi: 10.3847/1538-4357/ac3f38
  * Jones (1968) Jones, F. C. 1968, Physical Review, 167, 1159, doi: 10.1103/PhysRev.167.1159
  * Kazantsev (1968) Kazantsev, A. P. 1968, Soviet Journal of Experimental and Theoretical Physics, 26, 1031
  * Kelner et al. (2013) Kelner, S. R., Aharonian, F. A., & Khangulyan, D. 2013, ApJ, 774, 61, doi: 10.1088/0004-637X/774/1/61
  * Khangulyan et al. (2008) Khangulyan, D., Aharonian, F., & Bosch-Ramon, V. 2008, MNRAS, 383, 467, doi: 10.1111/j.1365-2966.2007.12572.x
  * Khangulyan et al. (2021) Khangulyan, D., Aharonian, F., Romoli, C., & Taylor, A. 2021, ApJ, 914, 76, doi: 10.3847/1538-4357/abfcbf
  * Kumar et al. (2012) Kumar, P., Hernández, R. A., Bošnjak, Ž., & Barniol Duran, R. 2012, MNRAS, 427, L40, doi: 10.1111/j.1745-3933.2012.01341.x
  * MAGIC Collaboration et al. (2019a) MAGIC Collaboration, Acciari, V. A., Ansoldi, S., et al. 2019a, Nature, 575, 455, doi: 10.1038/s41586-019-1750-x
  * MAGIC Collaboration et al. (2019b) —. 2019b, Nature, 575, 459, doi: 10.1038/s41586-019-1754-6
  * Nakar et al. (2009) Nakar, E., Ando, S., & Sari, R. 2009, ApJ, 703, 675, doi: 10.1088/0004-637X/703/1/675
  * Salafia et al. (2022) Salafia, O. S., Ravasio, M. E., Yang, J., et al. 2022, ApJ, 931, L19, doi: 10.3847/2041-8213/ac6c28
  * Sari & Esin (2001) Sari, R., & Esin, A. A. 2001, ApJ, 548, 787, doi: 10.1086/319003
  * Zdziarski (1989) Zdziarski, A. A. 1989, ApJ, 342, 1108, doi: 10.1086/167668
  * Zhang et al. (2009) Zhang, W., MacFadyen, A., & Wang, P. 2009, ApJ, 692, L40, doi: 10.1088/0004-637X/692/1/L40
  * Zhang (2019) Zhang, B. 2019, Nature, 575, 448, doi: 10.1038/d41586-019-03503-6
  * Zhang et al. (2021) Zhang, B. T., Murase, K., Veres, P., & Mészáros, P. 2021, ApJ, 920, 55, doi: 10.3847/1538-4357/ac0cfc


# Experimental Impact Analysis of Cyberattacks in Power Systems using Digital
Real-Time Testbeds

Kalinath Katuri1, Ioannis Zografopoulos2, Ha Thi Nguyen1, Charalambos
Konstantinou2  1Eversource Energy Center, Department of Electrical and
Computer Engineering, University of Connecticut, Storrs, CT, USA
2CEMSE Division, King Abdullah University of Science and Technology (KAUST)
E-mail: {kalinath.katuri<EMAIL_ADDRESS>{ioannis.zografopoulos,
<EMAIL_ADDRESS>
###### Abstract

Smart grid advancements and the increased integration of digital devices have
transformed the existing power grid into a cyber-physical energy system. This
reshaping of the current power system can make it vulnerable to cyberattacks,
which could cause irreversible damage to the energy infrastructure resulting
in the loss of power, equipment damage, etc. Constant threats emphasize the
importance of cybersecurity investigations. At the same time, developing
cyber-physical system (CPS) simulation testbeds is crucial for vulnerability
assessment and the implementation and validation of security solutions. In
this paper, two separate real-time CPS testbeds are developed based on the
availability of local research facilities for impact analysis of denial-of-
service (DoS) attacks on microgrids. The two configurations are implemented
using two different digital real-time simulator systems, one using the real-
time digital simulator (RTDS) with a hardware-in-the-loop (HIL) setup and the
other one using OPAL-RT with ExataCPS to emulate the cyber-layer
infrastructure. Both testbeds demonstrate the impact of DoS attacks on
microgrid control and protection operation.

###### Index Terms:

Cybersecurity, cyber-physical co-simulation, denial-of-service, digital real-
time testbeds.

##  I Introduction

The integration of computing and communication resources in power systems,
with the goal of developing a robust smart grid infrastructure, is
transforming power transmission, distribution, and utilization. The quest for
secure and resilient energy power systems (EPS) has promoted the deployment of
distributed and decentralized grid architectures. According to the 2020 U.S.
Department of Energy smart grid system report [1], smart grids are the
application of digital and cyber infrastructure to the physical system in
order to perform the functions of sensing, communications, control, computing,
and data and information management to assist planning and operations.
Reliable communication is required to leverage the benefits of smart grids and
the support services provided by distributed energy resources (DER). As a
result, the cyber layer of EPS and information and communication technologies
(ICT) are becoming an integral part of the smart grid paradigm.

As the power grid transitions to a cyber-physical system, it becomes highly
dependent on digital communication systems, increasing the risk of
cyberattacks. Some of the most recent examples of cyberattacks on the power
grid include the one in Ukraine on December 23, 2015, which resulted in the
loss of power to more than 225,000 customers and affected more than 50
substations [2], and the attack on the nuclear power plant in India during
2019 [3]. The greater threat to cyber-physical energy systems (CPES) and the
prohibitive cost and operational risks of conducting cybersecurity
investigations on actual power systems motivate the need for CPES testbeds
[4]. Such testbed setups can capture the spatial and temporal characteristics
of distributed EPS, while co-simulation studies can comprehensively evaluate
the impacts of cyber contingencies (e.g., communication delays) on the
physical system’s dynamic performance [5].

Research teams around the world have developed diverse CPES testbed topologies
to mirror real-world applications and enable cybersecurity research [6, 7]. In
[8], a real-time cyber-physical testbed with Real Time Digital Simulator
(RTDS) and OPNET is developed to investigate the impact of man-in-the-middle
(MITM) attacks on the 11-Bus system. In [9], a testbed was developed with both
real-time and non-real-time capabilities. The authors use RTDS and
PowerFactory for the physical system components, i.e., RTDS is communicated to
the cyber layer through IEC 61850 and DNP3, and PowerFactory is connected to a
virtual substation through an OPC interface. Researchers at Sandia National
Laboratories have developed a virtual control system testbed that incorporates
various commercially available components for the cyberattack impact analysis
on the energy critical infrastructure [10]. A comprehensive survey of
industrial control system testbeds and their capabilities, along with the
essential data required to perform cybersecurity studies, is furnished in
[11].

The security investigations performed on cyber-physical testbeds include, i)
vulnerability research, ii) mitigation research, iii) security validation, and
iv) impact analysis [9, 12]. In [13], the impact of coordinated powerbot
attacks on the power system is studied, while [14] and [15] demonstrate the
impact of denial-of-service (DoS) attacks on power grids. The impact of cross-
layer firmware attacks on microgrid (MG) inverter controllers is discussed in
[16], along with methodologies to design host-based defenses to overcome such
incidents. In [17], cyberattacks against the IEC 61850 protocol (in
photovoltaic (PV) installations) are studied. In [18] and [19], the authors
present the impacts of communication-based attacks targeting CPES protection
systems. The results of the aforementioned research underline the importance
of CPES security studies. Equally important is the design and constant
improvement of unified cyber-physical simulation and/or hardware platforms
where we can investigate different threats, varying system topologies, and
operation conditions.

This paper presents two different real-time simulation platforms for CPES to
evaluate the impact of diverse cyberattack scenarios on EPS. The work
illustrates the investigation of cybersecurity studies on different platforms
based on the availability of local research facilities. Specifically, our
contributions are as follows:

  * •

Impact analysis of cyberattacks targeting distribution systems is conducted
with one of the testbeds developed by University of Connecticut (UCONN)’s
Eversource energy center using RTDS Novacor chasis and a hardware relay.

  * •

The consequences of communication delays – due to DoS conditions – are
illustrated using toolkits provided by _OPAL-RT_ and _ExataCPS_ in King
Abdullah University of Science and Technology (KAUST)’s CPES testbed.

  * •

Development of a hardware-in-the-loop (HIL) implementation to investigate the
overcurrent conditions of a PV system and a battery under the impacts of
malicious cyberattacks degrading critical operations.

The rest of the paper is organized as follows. Section II provides an overview
of the testbed setups at UCONN and KAUST and their capabilities. The impact of
trip delays and setpoint changes is demonstrated via HIL simulations in
Section III. The impact of attacks targeting the grid’s cyber infrastructure
and their effects on generation and frequency are presented in Section IV.
Finally, Section V concludes the paper by highlighting the importance of CPES
testbeds for future secure power systems.

##  II Experimental Laboratory setup

This section outlines the implementation of two different experimental
laboratory setups for the impact analysis of MGs under cyberattacks. These
testbeds are built leveraging real-time simulation platforms.

###  II-A HIL Testbed with RTDS at UCONN

The CPES testbed is developed on the RTDS platform. For the impact analysis, a
controller HIL (CHIL)-based testbed is preferred not only to evaluate the
logic or settings of a protection system but also to assess closed-loop
controller performance. The schematic and lab setup of HIL developed for the
relay under test are shown in Fig. 1.

|
---|---
a) Schematic | b) Lab setup
Figure 1: Hardware-in-the-loop (HIL) testbed with RTDS.

The SEL 451 relay hardware in this work is used as a non-directional over-
current relay. The analog output card (GTAO) of the RTDS is used to retrieve
the voltage and current signals from the simulation. And these signals are
connected directly to the relay at a peak-to-peak voltage of 10 V. The relay
output is connected to the digital input (GTDI) card of RTDS (through a 5V DC
supply). Two cyberattack scenarios – DoS and MITM attacks – are simulated, and
their impacts on the MG and the relay performance are evaluated in real-time.

###  II-B Co-simulation Testbed with OPAL-RT at KAUST

This real-time co-simulation testbed is composed of two distinct layers, as
shown in Fig. 2. The setup investigates the impact that different durations of
DoS attacks could have on the physical portion of the EPS.

Figure 2: Cyber-physical co-simulation setup with OPAL-RT.

####  II-B1 Physical Model

On the physical layer, we model a MG in which we distinguish two power
sources; namely, a diesel generator (Genset) and a battery energy storage
system (BESS) with nameplate capacities of 1 MW and 100 kW, respectively.
Apart from the generation buses, we have also included three load buses, a
non-critical 300 kW lumped residential load, a 200 kW critical load (e.g.,
hospital), and a commercial load with a power demand of 700 kW. Both
generation and load components are interfaced through the point-of-common-
coupling (PCC), where the MG central controller is also located, with the main
grid. The aforementioned physical layer is simulated using the _OPAL-RT_ real-
time simulator and the _eMegasim_ power electronics toolkit.

####  II-B2 Cyber Model

On the cyber layer, for each of the buses (generation or load) residing on the
physical layer, a corresponding cyber layer node is emulated. We use the
IEEE-1815 Distributed Network Protocol 3 (_DNP3_) for the communication
between the cyber layer nodes. Furthermore, it can be seen from Fig. 2, that
the cyber layer nodes are connected in a radial architecture where the MG
central controller at the PCC serves as the _DNP3_ master, and the remaining
nodes constitute the _DNP3_ outstations. For the modeling and simulation of
the cyber layer, we leverage the _EXataCPS_ toolkit that allows the emulation
of communication networks and can be seamlessly integrated within the _OPAL-
RT_ environment.

Figure 3: Single line diagram of the microgrid network simulated in RTDS.

##  III RTDS-HIL Testbed: Cases & Results

The MG simulated in this work consists of a 20 kV distribution network with
two types of DERs, an aggregated representation of rooftop PV of 1 MW, and a
large-scale BESS of 1MW capacity. The single-line diagram (SLD) of the MG
implemented in the RTDS platform is shown in Fig. 3. Through a 10 MVA, 20/0.4
kV transformer, a 0.48 kV Genset is connected to the network. A total load of
4 MW is connected on both the 20 kV and 0.48 kV levels. In this work, the
attacker is considered to have gained access to the communication network, and
the test cases are the possible scenarios of DoS and MITM attacks on the
network.

The impacts of cyberattacks are studied by imitating the attack’s physical
behavior on the RTDS network. Three different test cases (under attack) are
considered. The first two cases imitate the behavior of a DoS attack during a
three-phase short-circuit fault at F1 and F2 locations, and the third case
imitates the behavior of a MITM attack on the PV system. The physical impact
of all three cases is studied with the help of the SEL 451 relay used in the
CHIL experiment, connected at locations R1, R2, and R3, respectively [20].
Relays are connected through a current transformer as shown in Fig. 3. For
these test cases, three-phase short-circuit faults are simulated due to their
severe fault currents.

###  III-A Battery Trip Delay Impact on Microgrid

Through this test case, the impact of DoS attacks on a protection trip during
a three-phase short-circuit fault at location F1 on Bus-7 is studied. During
the fault, due to the low short-circuit fault path, the battery will
contribute to the fault current, resulting in the tripping of the relay
connected to the high voltage (HV) side of the transformer at the R1 location,
which is referred to as “sympathetic tripping”.

Figure 4: Battery test case.

To study the impact of a DoS attack during this fault, it is assumed that the
trip command from the relay to the breaker is delayed, which causes a delay in
isolating the battery from the faulted network. The impact of this delayed
breaker operation is simulated by interfering with the relay output from SEL
451 in the CHIL simulation. The results are plotted in Fig. 4. A $2$ $secs$
delay is assumed in the breaker operation.

The comparison of the relay trip is shown in Fig. 4 with normal operation
denoted using black and under cyberattack in red. The current contribution
from the BESS during the fault is momentary, as the fault is temporary in
nature. It is observed that, without any delay, the voltage at PCC is
recovered at a slower rate after the fault is cleared. However, this does not
trigger any voltage protection, and the system regains stability after the
fault. In the $2$ $secs$ delay case, even though the frequency stabilizes to
its nominal value after a few cycles, the inverter voltage fails to recover to
the rated 1.0 pu range, resulting in the activation of the voltage protection
on the BESS side.

###  III-B Delay in Relay Operation for Feeder Faults

This test case is used to examine the effects of operations being delayed as a
result of a DoS attack during a fault in the middle feeder. The main
overcurrent feeder protection relay is connected at R2 and operates in this
capacity. It is expected that a coordinated DoS attack can cause the trip
feedback to the relay at R2 to be delayed. The delayed response of the relay,
and consequently the circuit breaker under attack, can cause disturbances in
the voltage and frequencies of the diesel generator, hence affecting other
connected equipment in the MG. If the fault persists, this delay in relay
operation would result in unrecoverable damage to the system, depending on its
duration. The prolonged unstable voltage and frequency profiles may also
result in the tripping of inverter-based resources connected to the system,
which could potentially lead to MG collapse.

Figure 5: Genset response to the delay in fault on feeder case.

Fig. 5 shows the results for the $2$ $secs$ delayed impact (red) compared to
the no delay (black) response for the fault in Bus-4 at F2. The red curves
show that the frequency of the Genset exhibits a large oscillatory behavior,
and the system voltage on Bus-2 is observed to be abnormal for a considerable
period. These results not only show the severe nature of the impact on the
Genset, but also give an understanding of the damage to voltage-sensitive
equipment connected to the network. Fig. 6 depicts the operation of voltage
protection for the PV-connected bus. Fig. 7 shows the impact of delayed
breaker operation in the system. Both the BESS and the PV systems have tripped
on voltage protection as a result of the voltages recorded at Bus-7 and Bus-3
falling below the set value of voltage protection for UV2 ($<$ 0.45 pu) for
longer than the set period ($0.16$ $secs$). The voltage values used are
compliant with the recommendations provided in IEEE 1547:2018.

Figure 6: PV response to fault on middle feeder. Figure 7: BESS response to
fault on middle feeder.

###  III-C PQ Setpoint Change for PV System

In the secondary control of the PV plant, active (P) and reactive (Q) power is
controlled by the utility set points. Such set points enable the effective use
of green energy available at any time of the day. Through this case study, the
impact of MITM attacks on PV systems is studied. When an attacker gains access
to the control network, they can either deny the service or change the control
set points of the PV plant without the knowledge of the utility. This
deliberate action would result in a sudden increase in reverse power flow to
the distribution grid. This sudden change in reference setpoints of P and Q
would result in an increase in the current flow from the PV to the network,
which can trigger the overcurrent relay connected at R3 at the PCC point for
the PV plant in Fig. 3. Furthermore, it may also activate internal voltage and
frequency protections depending on the system state during the change.

Figure 8: PV response to PQ set point change.

The results in Fig. 8 show the system conditions when the P and Q outputs of
the PV system increase. The effect is demonstrated not only by the higher
current output but also depicted in the voltage profile. The current change
activates the overcurrent relay, subsequently resulting in breaker operation,
whereas the voltage profile results in voltage protection pickup. The spikes
in P and Q profiles are related to the disconnection of the PV system during
the breaker operation.

##  IV OPAL-RT Testbed: Cases & Results

In the following case study, we leverage a real-time co-simulation testbed to
examine the cross-layer impact of DoS cyberattacks. Specifically, high-
fidelity models are considered for both the cyber layer (accounting for
communication delays) and the physical layers. Then, using the real-time co-
simulation environment (Fig. 2), we demonstrate how the impact of different
DoS conditions (encountered on the cyber layer) can propagate to the actual
physical grid compromising its stable operation and affecting its steady-state
stability.

We evaluate different DoS attack examples and their impact on the EPS
frequency and power measurements. DoS cyberattacks, impeding the timely
delivery of critical control signals (e.g., islanding or load-shedding
commands), can introduce abnormal system states and lead to severe system
maloperation. In our simulations, the DoS severity is indicated by the
duration of the introduced time-delay attacks, which can, in turn, have minor
(e.g., uneconomic operation if power flow data are delayed) or major (e.g.,
generation instabilities if control commands do not respect the constraints of
the EPS) impacts on system objectives. It should be noted that although the
discussed attacks (Section IV-A) are performed solely on the cyber layer (by
delaying communications), their impact on grid properties (frequency,
generated power, etc.) can be severe. In the rest of the section, we present
simulation results illustrating the impact of DoS attacks on the isolated MG
setup. The specifics about the MG’s physical and cyber modeling are detailed
in Sections II-B1 and II-B2.

(a)

(b)

Figure 9: Frequency and power measurements of the MG during DoS attacks.
LABEL:sub@fig:freq_all Impact of $2$, $5$, and $15$ $secs$ load-shedding delay
on frequency, LABEL:sub@fig:pow_all Impact on the power generation profiles.

###  IV-A Simulation Results

In this part, we present the potential impact of cyber-layer-based attacks
(introducing communication delays) on MG operations. Fig. 9 illustrates the
frequency fluctuations and power sharing between the generator and BESS during
three different DoS attack scenarios. In more detail, the MG system (Fig. 2)
is islanded from the main grid by tripping the circuit breaker located at the
PCC at $t=10$ $secs$. The inherent MG generation capacity, i.e., power sourced
from the Genset and the BESS, cannot meet the aggregated load demand, as can
happen in weak MG architectures. To avoid under-frequency conditions,
defensive load-shedding should be performed, prioritizing critical and
commercial loads (as described in Section II-B1) and curtailing the power
demand introduced by residential loads.

Immediately after the MG isolation ($t=10$ $secs$), the MG controller issues a
load-shedding command to maintain grid frequency stability. However, at the
same time, an adversary performs a DoS attack that “floods” the network’s
communication capacity. As a result, the DNP3 messages issued by the MG master
controller are delivered with arbitrary delays to the outstation devices
(e.g., load controllers). In our case, the load-shedding command for the
residential load is delayed. Fig. 9a illustrates the impact of such time DoS
delays on the MG frequency if the load-shedding mechanism is delayed by $2,5,$
and $15~{}secs$, respectively.

Even in the less severe DoS attack case, where the load-shedding command is
delayed for $2~{}secs$, we can observe that the MG frequency drops below the
nominal operational range, e.g., $59.5$ Hz [21]. In the second scenario, where
the load-shedding command is delayed by $5~{}secs$, a significant frequency
decrease occurs too, i.e., the frequency drops below $56$ Hz, in which case
corrective mechanisms should be immediately enforced to protect the grid
assets from this abnormal frequency condition. It should be noted that
regardless of the insufficient MG power generation, in both of the
aforementioned scenarios, the MG is able to maintain its stable operation,
although the Genset has to momentarily operate above its rated capacity for
the duration of the DoS.

The MG resilience to DoS-based cyberattacks is compromised in the last
scenario. Fig. 9a and Fig. 9b illustrate the impact on frequency and generator
stability during an attack that delays the load-shedding command delivery for
$15~{}secs$. The results in both graphs explicitly demonstrate the potential
impacts that cyber contingencies, e.g., communication delays, could impose on
the MG dynamics. In more detail, in Fig. 9a, we illustrate that the MG
frequency drops below $55$ Hz if the command is delayed by more than
$5~{}secs$, which could prompt further load-shedding to maintain frequency
stability. Furthermore, as shown in Fig. 9b, the Genset is operated above its
nominal capacity for more than $10~{}secs$, eventually rendering the whole MG
system unstable. At this point, the MG system – currently operating
autonomously – would experience a blackout and would not be able to support
the $200kW$ critical load.

##  V Conclusion

The risk of cyberattacks on power grids and energy infrastructure has been
continuously increasing. The advancement of smart grid infrastructure poses a
greater risk of cyberattacks due to unknown vulnerabilities. The growing
importance of cybersecurity research has motivated researchers to develop
unique testbeds focused on different aspects of cybersecurity research.

In this work, two different testbeds are studied using real-time simulations
for the impact analysis of cyberattacks on MGs. The two testbeds show that
cybersecurity research can be conducted on different platforms based on the
availability of local facilities. The first setup implements a HIL experiment
using the RTDS simulator, while the second focuses on a cyber-physical setup
using the _OPAL-RT_ simulator. Both configurations investigate DoS attack
impacts in real-time environments, and the results from both testbeds show
that cyberattacks can cause damage and even cascading failures or power system
collapse. More research on CPES security is essential to provide insights for
system operators and propose solutions to prevent, detect, and mitigate
impacts. In future work, different types of cyberattacks, such as MITM and
replay attacks, will be studied on both testbeds to evaluate the impacts of
these attacks on power grids and provide potential mitigation strategies.

## References

  * [1] Department of Energy (DOE), “2020 smart grid system report.” [Online]. Available: https://www.energy.gov/oe/articles/2020-smart-grid-system-report-download
  * [2] D. E. Whitehead, K. Owens, D. Gammel, and J. Smith, “Ukraine cyber-induced power outage: Analysis and practical mitigation strategies,” in _70th Annual Conf. for Protective Relay Engineers (CPRE)_ , 2017.
  * [3] M.-G. P. K. Mallick and V. Retd, “Cyber attack on kudankulam nuclear power plant a wake up call.” [Online]. Available: https://www.vifindia.org/sites/default/files/cyber-attack-on-kudankulam-nuclear-power-plant.pdf
  * [4] I. Zografopoulos, N. D. Hatziargyriou, and C. Konstantinou, “Distributed energy resources cybersecurity outlook: Vulnerabilities, attacks, impacts, and mitigations,” 2022. [Online]. Available: https://arxiv.org/abs/2205.11171
  * [5] I. Zografopoulos, J. Ospina, X. Liu, and C. Konstantinou, “Cyber-physical energy systems security: Threat modeling, risk assessment, resources, metrics, and case studies,” _IEEE Access_ , vol. 9, pp. 29 775–29 818, 2021.
  * [6] M. H. Cintuglu, O. A. Mohammed, K. Akkaya, and A. S. Uluagac, “A survey on smart grid cyber-physical system testbeds,” _IEEE Communications Surveys & Tutorials_, vol. 19, no. 1, pp. 446–464, 2016.
  * [7] A. Keliris, C. Konstantinou, N. G. Tsoutsos, R. Baiad, and M. Maniatakos, “Enabling multi-layer cyber-security assessment of industrial control systems through hardware-in-the-loop testbeds,” in _21st ASP-DAC_ , 2016, pp. 511–518.
  * [8] B. Chen, K. L. Butler-Purry, A. Goulart, and D. Kundur, “Implementing a real-time cyber-physical system test bed in RTDS and OPNET,” in _2014 North American Power Symposium (NAPS)_ , 2014, pp. 1–6.
  * [9] A. Hahn, A. Ashok, S. Sridhar, and M. Govindarasu, “Cyber-physical security testbeds: Architecture, application, and evaluation for smart grid,” _IEEE Transactions on Smart Grid_ , vol. 4, no. 2, pp. 847–855, 2013.
  * [10] G. Pollock, W. Atkins, M. Schwartz, R. Chavez, J. Urrea, N. Pattengale, M. McDonald, R. Cassidy, R. Halbgewachs, B. Richardson, and J. Mulder, “Modeling and simulation for cyber-physical system security research, development and applications,” 01 2010.
  * [11] M. Conti, D. Donadel, and F. Turrin, “A survey on industrial control system testbeds and datasets for security research,” _IEEE Communications Surveys & Tutorials_, vol. 23, no. 4, pp. 2248–2294, 2021.
  * [12] I. Zografopoulos, C. Konstantinou, N. G. Tsoutsos, D. Zhu, and R. Broadwater, “Security assessment and impact analysis of cyberattacks in integrated t&d power systems,” in _Proceedings of the 9th Workshop on Modeling and Simulation of Cyber-Physical Energy Systems_ , 2021.
  * [13] S. Soltan, P. Mittal, and H. V. Poor, “Blackiot: Iot botnet of high wattage devices can disrupt the power grid,” in _Proceedings of the 27th USENIX Conference on Security Symposium_ , 2018, p. 15–32.
  * [14] R. Kalluri, L. Mahendra, R. S. Kumar, and G. G. Prasad, “Simulation and impact analysis of denial-of-service attacks on power scada,” in _2016 National Power Systems Conference (NPSC)_ , 2016, pp. 1–5.
  * [15] C. Deng, “Distributed resilient control for cyber-physical systems under denial-of-service attacks,” in _2019 23rd International Conference on Mechatronics Technology (ICMT)_ , 2019, pp. 1–5.
  * [16] I. Zografopoulos, A. P. Kuruvila, K. Basu, and C. Konstantinou, “Time series-based detection and impact analysis of firmware attacks in microgrids,” _Energy Reports_ , vol. 8, pp. 11 221–11 234, 2022.
  * [17] B. Kang, P. Maynard, K. McLaughlin, S. Sezer, F. Andrén, C. Seitl, F. Kupzog, and T. Strasser, “Investigating cyber-physical attacks against iec 61850 photovoltaic inverter installations,” in _IEEE 20th Conf. on Emerging Technologies & Factory Automation (ETFA)_, 2015, pp. 1–8.
  * [18] A. A. Jahromi, A. Kemmeugne, D. Kundur, and A. Haddadi, “Cyber-physical attacks targeting communication-assisted protection schemes,” _IEEE Transactions on Power Systems_ , vol. 35, no. 1, pp. 440–450, 2019.
  * [19] S. M. S. Hussain, M. A. Aftab, S. M. Farooq, I. Ali, T. S. Ustun, and C. Konstantinou, “An effective security scheme for attacks on sample value messages in iec 61850 automated substations,” _IEEE Open Access Journal of Power and Energy_ , vol. 10, pp. 304–315, 2023.
  * [20] K. Katuri, H. T. Nguyen, T. T. Hoang, Z. Wang, and E. Anagnostou, “Lea-based hardware in the loop validation for der integration impact on distribution protection,” in _2022 International Conference on Renewable Energies and Smart Technologies (REST)_ , vol. I, 2022, pp. 1–5.
  * [21] North American Electric Reliability Corporation (NERC), “2019 Frequency Response Annual Analysis,” https://www.nerc.com/, 2019.