# Acceleration of Solar Energetic Particles by the shock of Interplanetary
Coronal Mass Ejection

Shanwlee Sow Mondal Astronomy and Astrophysics Division, Physical Research
Laboratory, Ahmedabad 380009, India  Indian Institute of Technology,
Gandhinagar, Gujarat 382355, India  Aveek Sarkar  Astronomy and Astrophysics
Division, Physical Research Laboratory, Ahmedabad 380009, India  Bhargav
Vaidya  Discipline of Astronomy, Astrophysics and Space Engineering, Indian
Institute of Technology Indore, Simrol, Indore 453552, India  Andrea Mignone
Dipartimento di Fisica, Universita di Torino, via P. Giuria 1, I-10125 Torino,
Italy

###### Abstract

Interplanetary Coronal Mass Ejection (ICME) shocks are known to accelerate
particles and contribute significantly to Solar Energetic Particle (SEP)
events. We have performed Magnetohydrodynamic-Particle in Cell (MHD-PIC)
simulations of ICME shocks to understand the acceleration mechanism. These
shocks vary in Alfvénic Mach numbers as well as in magnetic field orientations
(parallel & quasi-perpendicular). We find that Diffusive Shock Acceleration
(DSA) plays a significant role in accelerating particles in a parallel ICME
shock. In contrast, Shock Drift Acceleration (SDA) plays a pivotal role in a
quasi-perpendicular shock. High-Mach shocks are seen to accelerate particles
more efficiently. Our simulations suggest that background turbulence and local
particle velocity distribution around the shock can indirectly hint at the
acceleration mechanism. Our results also point towards a few possible in situ
observations that could validate our understanding of the topic.

methods: numerical, Sun: corona, Sun: coronal mass ejections (CMEs), Sun:
heliosphere, Sun: particle emission, Sun: UV radiation, Sun: CME, Sun: shock

##  1 Introduction

Solar Energetic Particles (SEPs) are non-thermal particles emanating from the
Sun. Their energy can be several orders of magnitude higher than the ambient
solar wind particles. These particles can be electrons, protons, or heavy
ions. Usually, SEPs are detected soon after solar flares or when the
heliospheric particle detectors encounter Interplanetary Coronal Mass
Ejections (ICMEs). They have been detected at different heliospheric distances
by multiple in situ spacecrafts. SEPs can achieve energies up to several
hundred MeV (Reames, 1999) and sometimes even to GeV (Ryan et al., 2000).
Intense SEP events may become hazardous for humans in space and can also
damage the electronics of space instruments (Feynman & Gabriel, 2000).

These energetic events are considered to originate via two different
mechanisms – particle acceleration in CME-driven shock waves (Gradual SEP) and
magnetic reconnection (Impulsive SEP) (Reames, 1999). Impulsive SEP events
occur for a shorter duration ($\leq$ one day), with an occurrence frequency of
about 1000 per year during solar maxima (Reames, 2002). These are of low
intensities. Gradual SEP events continue for a longer duration (several days).
They are relatively rare (a few tens per year)(Reames, 2002) and produce
proton flux several orders of magnitude larger than impulsive SEP events
(Klein & Dalla, 2017). It is not easy to distinguish between these two
processes in practice as they are intertwined and may co-occur. Magnetic
reconnection may happen at the current sheet located behind the CME and at
regions where the CME interacts with the ambient solar wind or fragmented
current sheets of the post-eruptive environment. Shock acceleration may occur
in the region where the reconnection jet interacts with the ambient coronal
material and produces an outflow termination shock and shock waves driven by
the super magnetosonic motion of the CME (Klein & Dalla, 2017).

Multiple in situ observations have demonstrated (e.g., Sheeley et al. (1983);
Gopalswamy (2003); Kahler (2004)) the existence of SEPs during the passage of
ICMEs. It is believed that the solar wind particles interacting with the
shocks ahead of the ICMEs may get energized to such high energies (Lee & Fisk,
1982).

SEPs’ peak intensities are found to be correlated with the speed and other
physical parameters of the associated CMEs (Kouloumvakos et al., 2019; Desai &
Giacalone, 2016; Kahler, 2001; Kahler & Vourlidas, 2005, 2013). The strength
of ICME driven shock waves is usually determined by the Alfvénic Mach number
($M_{A}=u/v_{A}$), where $v_{A}=B/\sqrt{4\pi\rho}$ is the local solar wind
Alfvén speed, and $u$ is the upstream solar wind speed in the shock rest
frame. Commonly observed ICME-driven shocks possess $M_{A}\approx$ 2-4
(Berdichevsky et al., 2000; Oh et al., 2007). However, for stronger ICMEs,
$M_{A}$ can go beyond five (Lugaz et al., 2015). A powerful shock was recorded
(Russell et al., 2013) on 2012 July 23, by the STEREO-A spacecraft, where the
Alfvénic Mach number was determined to be $\approx$ 21 (Riley et al., 2016).
Proton energy spectra obtained from the particle detectors onboard STEREO-A
have shown that particles are energized up to 100 MeV (Russell et al., 2013).

Like other astrophysical shocks, ICME shocks are also collisionless in nature.
The particle mean-free path of the system is much larger than the shock length
scale (say, the width of the shock) or gyroradii of the particles. Primarily
two mechanisms are thought to be responsible for particle acceleration in
these kinds of systems – Diffusive Shock Acceleration (DSA) (Fermi, 1954;
Krymskii, 1977; Bell, 1978; Blandford & Ostriker, 1978) and Shock Drift
Acceleration (SDA) (Armstrong et al., 1985; Decker, 1988). Depending on the
upstream magnetic field’s orientation with respect to the shock normal, the
respective mechanism dominates. For example, in the case of quasi-parallel
shocks, where the ambient magnetic field makes an angle in the range
$0^{\circ}<\theta<45^{\circ}$ with the shock normal, DSA works as the dominant
mechanism in accelerating the charged particles. Particles get scattered by
magnetic turbulence and undergo repetitive head-on interactions with upstream
and downstream plasma across the shock, thereby gaining energy. DSA is
considered to be the most efficient mechanism responsible for the origin of
highest-energy particles, mostly observed at quasi-parallel shocks. For strong
($M_{A}\gg 1$) non-relativistic quasi-parallel shocks, the energy spectra of
DSA maintains a power-law, $f(E)\propto E^{-1.5}$ (Caprioli & Spitkovsky,
2014a). DSA is considered to be the promising theory behind the genesis of
gradual solar energetic particle (SEP) events.

On the other hand, in quasi-perpendicular shocks, where the angle between the
upstream magnetic field and the shock normal lies in the range
$45^{\circ}<\theta<90^{\circ}$, SDA plays a pivotal role in particle
acceleration. In this case particles gain energy by drifting along the shock
in the direction of the induced electric field (Decker, 1983; Reames, 2012).

Particle acceleration is common in other astrophysical shocks such as in
supernova remnants (Helder et al., 2012) or in extragalactic jets (Romero et
al., 2017). These shocks are considered to be the primary sources of galactic
cosmic rays (CRs). Many efforts have already been made to understand the
origin and transport of high energy particles (CRs) — their complex non-linear
interactions with the background thermal plasma and the ambient magnetic field
(Fermi, 1949, 1954; Blandford & Ostriker, 1978).

The most obvious framework to study the evolution of CRs numerically around a
shock would be Particle in Cell (PIC) method. In the PIC, electrons and ions
are treated as particles, therefore suitable to study plasma kinetics in the
presence of a background electromagnetic field (Spitkovsky, 2005). However,
one needs to resolve the electron skin depth in the PIC code, which makes the
simultaneous study of shock evolution and particle acceleration highly
expensive.

Recently the field has been enriched with the development of hybrid-PIC codes
where the electrons are considered as part of the background thermal fluid,
and ions are treated as kinetic particles (Lipatov, 2002; Gargaté et al.,
2007). Such simulations need to resolve the ion-skin depth of the medium,
making the simulations computationally expensive. With the availability of
powerful machines, numerous works have been done where simulations were
performed (Caprioli & Spitkovsky, 2013, 2014a, 2014b, 2014c) with relatively
large domains and for long duration, but still, only the initial stages of
shock evolution and particle acceleration have been captured. However, the
transition of charged particles from non-relativistic to relativistic regimes
and the fraction of ions participating in the process would only be known when
we can study the long term evolution of collisionless shocks. The MHD-PIC
approach is more relevant for such purposes (Bai et al., 2015; van Marle et
al., 2018; Mignone et al., 2018). This method can describe the non-linear
interaction between thermal (ions+electrons) and non-thermal energetic
particles. The method conveniently ignores microscopic plasma scales but
captures the kinetic effect of ions by resolving the particles’ gyroradii.

The MHD formulations carry out the evolution of the thermal plasma consisting
of both electrons and ions, and the dynamics of the particles (non-thermal
ions) are studied using the PIC method. All the non-thermal particles are
extracted from the thermal fluid itself. After extracting a certain fraction
of ions from the background plasma, the Lorentz force experienced by these
charged particles in the background fluid’s electric and magnetic fields is
calculated. This force now acts as the feedback force that the particles exert
on the background fluid and thus modifies its evolution. Recent usage of
adaptive mesh refinement on the background thermal plasma (van Marle et al.,
2018) has made the MHD-PIC more computationally affordable method.

Numerical study of particle acceleration in heliospheric shocks was initiated
using the Monte-Carlo simulation (Baring et al., 1997). They showed that the
particle energy spectra obtained from their simulation are in good agreement
with the in situ energy spectra obtained from the interplanetary shock,
suggesting an active DSA mechanism in the quasi-parallel heliospheric shock.

An extensive study of Particle acceleration by CME shocks has been done using
the Particle Acceleration and Transport in Heliosphere (PATH) code (Zank et
al., 2000). This model includes spherically symmetric Solar wind, in which a
CME-driven shock wave propagates. Particles are injected at the shock using an
injection model where the number of injected particles is a small fraction of
the thermal Solar wind particles. Injected particles get accelerated through
DSA. Some of these particles also escape far upstream through diffusion. The
outcome of such modeling includes the temporal evolution of energetic particle
spectra and intensity profile at all the spatial locations (both upstream and
downstream), and the determination of the particle injection energy and the
maximum energy of particles accelerated at the shock.

The same PATH code was later modified to include shocks with arbitrary
strengths (Rice et al., 2003), giving a better estimation of the maximum
particle energy for each shock. The transport of energetic particles (Li et
al., 2003) escaping the shock, simulating the acceleration of heavy ions (Li
et al., 2005b) also have been studied using the same. It has also been used
(Zank et al., 2004, 2006) to study particle acceleration at a perpendicular
interplanetary shock. Modeling individual SEP events using the PATH code have
been pursued by several authors (Li et al., 2005a; Zank et al., 2007;
Verkhoglyadova et al., 2009, 2010). Li et al. (2012) have studied the particle
acceleration in oblique shock. Their results suggest that close to the Sun
($r<$ 0.1 AU), quasi-parallel shocks are better particle accelerators than
quasi-perpendicular shocks.

A 2-D extension of the original PATH code (Zank et al., 2000) was developed by
(Hu et al., 2017). This new model can take care of the evolution of the
background solar wind and propagate the CME shock in a 2-D domain. It also
calculates the particle acceleration at the shock along with their diffusion
and convection in the upstream and downstream regions. This 2-D model can also
capture the longitudinal distribution of the CME.

Using the hybrid-PIC code, dHybrid Gargaté et al. (2007) have studied the
particle acceleration in ICME shocks. This kind of study can be categorized
more as a local shock simulation. Instead of taking the whole longitudinal
distribution of the shock, a particular region is considered for the
simulation, and acceleration of non-thermal particles is studied. The
simulation could successfully capture the detailed spatial and temporal
information of the electromagnetic waves generated due to counter-streaming
ions in the shock upstream. However, the hybrid method’s high demand for
computational resources restricts them to simulate only up to 15 sec of the
particle acceleration. The MHD-PIC method is a natural choice to deal with
such difficulties.

The present work utilizes the MHD-PIC version of the PLUTO code (Mignone et
al., 2018). Unlike the hybrid-PIC code, through this approach, one can track
the shock for several minutes and evolve the energy spectrum of non-thermal
particles. We simulate ICME shocks with various strengths (with different
$M_{A}$) and magnetic configurations (parallel shock as well as perpendicular
shock). The evolution of non-thermal protons is studied until their energy
spectrum gets saturated. The aforementioned simulations explain the physical
processes behind gradual SEP events.

The rest of the article is organized as follows. We describe the overall
simulation setup in section 2. Parallel and quasi-perpendicular shocks with
high Mach numbers are elaborated in section 3. Results of low Mach shocks are
detailed in section 4. Finally, we summarize our results in section 5.

##  2 Simulation setup

A complete derivation of the MHD-PIC formalism is given in Bai et al. (2015).
The MHD-PIC version of the PLUTO code is also explained in Mignone et al.
(2018).

We have simulated ICME shocks of different Alfvénic Mach numbers and analysed
particle acceleration in each of them. ICME-shocks observed at 1 AU vary over
a wide range of Alfvénic Mach numbers. A comparative study is performed to
understand the role of the shocks’ strength in the particle energisation
process. In all our simulations density, temperature and magnetic field scales
are taken to be $n_{0}=10$ cm-3, $T_{0}=5.6\times 10^{4}$ K and $B_{0}=2.817$
nT. Our length, time and velocity scales are chosen in terms of ion skin depth
($c/\omega_{pi}=72$ km), inverse of ion cyclotron frequency
($\Omega^{-1}=3.69$ sec) and Alfvén velocity ($v_{A}=19.5$ km sec-1). Here,
$c$ is the speed of light and $\omega_{pi}$ is the ion plasma frequency. To
emulate the solar wind (SW) plasma at $1$ AU, simulation domains are initially
filled with a plasma of uniform density $n_{bg}=0.5n_{0}$ ($5$ cm-3), and
temperature $9T_{0}$ ($5\times 10^{5}$ K) with uniform magnetic field strength
of $B_{bg}=3.55B_{0}$ ($10$ nT) whose orientation is different for parallel
and quasi-perpendicular cases. Each simulation domain consists of a 2D
rectangular box of length $L_{x}=1.5\times 10^{5}(c/\omega_{pi})$ and
$L_{y}=5\times 10^{3}(c/\omega_{pi})$ with uniform resolution $\Delta x=\Delta
y=10(c/\omega_{pi})$.

###  2.1 Shock generation

Generally, in the inner heliosphere, ICMEs travel faster than the ambient
solar wind, with both reaching supersonic and super Alfvénic speeds. A shock
wave gets developed ahead of the ICME. In between the shock and the ICME, a
turbulent region gets matured known as the sheath region (Kilpua et al.,
2017).

In numerical simulations such as the ones considered in this work, shocks are
typically formed by using the piston method (Gargaté et al., 2014; Bai et al.,
2015). A leftward super Alfvénic flow driven continuously from the right
boundary gives rise to a rightward propagating shock after getting reflected
from the left conductive wall. The shock’s strength and speed are controlled
by adjusting the velocity of the upstream flow. The top and bottom boundaries
maintain the periodic boundary condition.

The method described above can simulate an ICME shock quite nicely. Using this
method, we capture the shock and part of the sheath region of an ICME in the
downstream reference frame. With other parameters chosen to emulate solar wind
parameters at 1 AU, the inflow velocity is identified with the relative
velocity between the SW and the CME. For our purpose we have chosen shocks
with strong compression ratio,
$r=\frac{\rho_{d}}{\rho_{u}}\approx\frac{\gamma+1}{\gamma-1}\approx 4$, where
$\gamma=5/3$ is the adiabatic index and $\rho_{u}$ & $\rho_{d}$ represent the
upstream and downstream densities, respectively.

The shock converts the kinetic energy of the upstream flow to the thermal
energy of the downstream plasma. As a result the downstream temperature is
increased to $T_{d}=v_{u}^{2}/(r-1)\sim v_{u}^{2}/3$ (in code units), where
$v_{u}$ is the upstream velocity in the downstream reference frame. For a
strong shock the upstream temperature ($T_{u}$) is independent of the
downstream temperature (Bai et al., 2015).

###  2.2 Particle Injection

Unlike the Hybrid-PIC, the present MHD-PIC is not capable of generating non-
thermal particles self-consistently from the thermal plasma. Therefore a
strategy has to be adopted to inject protons that participate in the
acceleration process. Even though electrons do also get energized in reality,
limited computational power prohibits one from simulating their dynamics.
Therefore we assume all electrons are thermalized and do not participate at
the energization process. Only protons are capable of participating in the
non-thermal process. In this article, proton and particle both represent the
same entity. The participating protons are injected throughout the simulation.
Equivalent mass, momentum, and energy of protons are subtracted from the
background to ensure conservation.

The amount of injected protons is chosen in accordance with Bai et al. (2015),
inspired by Caprioli & Spitkovsky (2014a). The fraction of the injected proton
mass to the mass swept away by the shock, $\eta$ is chosen to be $3\times
10^{-3}$ in all our production runs. We have checked that if we increase this
fraction by a significant amount, the injected protons can disrupt the shock
through backreaction. On the other hand, a too-small value of $\eta$ may not
contribute to the backreaction, and the proton acceleration process may become
inefficient. The location of the shock front is tracked every time using the
PLUTO shock tracking algorithm (Mignone et al., 2018). Once the shock is
detected, the amount of particles we inject is a fixed fraction, $\eta$ of the
mass swept by the shock as mentioned above. The number of injected particles
is equivalent to 10 particles per cell at background fluid density $n_{0}$.

Since the background thermal plasma and injected particles are treated
differently in the MHD-PIC code, an ad hoc injection prescription needs to be
adopted for proton injection. An injection recipe can be adopted by
calibrating the MHD-PIC simulation with a more self-consistent hybrid-PIC
simulation. Nevertheless, in the absence of such results, so far monoenergetic
injection recipes have been used in MHD-PIC (Bai et al., 2015; Mignone et al.,
2018; van Marle et al., 2018). However, to inject particles with a wide range
of energies, we have primarily used a Maxwellian distribution with downstream
characteristic temperature.

The following sections are dedicated to the numerical setup and results
obtained from simulations of particle acceleration in different shocks.

##  3 Simulation of ICME shock with $M_{A}\approx 19$

###  3.1 Numerical setup

Keeping in mind the detection of a strong ICME shock (Russell et al., 2013;
Riley et al., 2016) at $1$ AU with Alfvénic Mach $\approx 19$ in the
downstream reference frame, we initially simulate an ICME-shock with the same
Alfvénic Mach number. An initial background magnetic field of strength
$B_{bg}=3.55B_{0}$ is applied along the $x$ direction, which is parallel to
the shock normal. Number of particles injected at the shock front depends on
$\eta$ times the amount of mass swept by the shock, with $\eta=3\times
10^{-3}$ as mentioned in section 2.2. In a system like this, an instability is
expected to develop due to the current generated by the energetic protons
escaping the shock front thereby perturbing the upstream magnetic field (Bell,
2004). The present resolution is sufficient to resolve the most unstable mode
of such instability. Near the shock front the wavelength of this mode is
approximately $3\pi(B_{bg}/B_{0})/\eta M_{A}=587(c/\omega_{pi})$ (Bai et al.,
2015). Considering the SW and ICME speed to be $400$ km s-1 and $2250$ km s-1,
respectively, we take their relative speed to be $-1850$ km s-1. Inflow speed
of $-94.87v_{A}$ ($=-1850$ km s-1) produces a rightward propagating shock of
Alfvénic Mach number $\approx 19$. The detailed generation mechanism of the
shock is described in section 2.1. Once the shock is formed, particles are
injected following the prescription discussed in section 2.2.

Figure 1: Density evolution of the parallel shock ($M_{A}\approx 19$), moving
towards right. Only a fraction of the simulation domain is depicted here. (An
animation of this figure is available here.) Figure 2: Fluctuating mass
density ($\delta\rho/\rho_{0}$) power spectrum during the near-saturation
phase of the simulation (averaged over time $t=1000$ to $1170$ $\Omega^{-1}$
). Here $\rho_{0}$ is the initial upstream density. The spectrum is derived at
the location $330(c/\omega_{pi})$ ahead of the shock front (Figure 1), over a
width $2400(c/\omega_{pi})$. Figure 3: Transverse magnetic field evolution
around the shock front. At every time snap, the shock location is re-centered
at $x=0$. (An animation of this figure is available here.) Figure 4:
Dimensionless transverse magnetic energy power spectrum at $t=1170$
$\Omega^{-1}$ is calculated in a region of width $\approx 2400(c/\omega_{pi})$
just ahead of the shock. Blue and magenta dashed lines indicate the wave modes
in resonance with the particles having energies $E_{sh}$ and $300E_{sh}$
respectively. Theoretical prediction of the fastest growing mode of the Bell’s
instability is marked with the green dashed line, calculated with the
transversely averaged current at the shock front. The black dashed line
indicates the Kolmogorov power-law slope.

###  3.2 Shock propagation

Density snapshots of the propagating shock are shown in Figure 1. It is to be
noted that only a small part of the whole simulated domain is shown here. In
the beginning, injected protons move freely in the upstream region. The motion
of the energetic particles across the shock generates a current $\bm{J}$,
primarily parallel to the ambient upstream magnetic field. The current
initially perturbs the mean upstream magnetic field to generate Alfvén waves.
Once the current gets stronger due to the high flux of energetic particles, it
generates an instability known as Bell instability after its founder (Bell,
2004). That, in turn, produces a perpendicular ($\delta\bm{B}\sim B_{\perp}$)
magnetic field in the upstream (Bell, 2005; Reville & Bell, 2012). This
magnetic field generates a force ($-\bm{J}\times\delta\bm{B}$), which acts on
the local plasma and evacuates the region. As a result, low-density cavities
are produced. Suprathermal particles occupy these cavities. Once the maximum
current carrying particles’ gyroradii become comparable to the cavities’ size,
the growth of the instability ceases (Caprioli & Spitkovsky, 2013). Freely
moving protons in the region get scattered by the density cavities.

To estimate the density cavities’ size we consider the average density power
spectrum over a short time window ($t=1000$ to $1170~{}\Omega^{-1}$), towards
the end of the simulation. The density power spectrum (Figure 2) of the shock
upstream plasma is derived from a region $330(c/\omega_{pi})$ ahead of the
shock front, having width $2400(c/\omega_{pi})$. Sufficient cavities have
formed at this location. The peak of the power spectrum, which is at
$55(c/\omega_{pi})$ indicates the dominant length scale of the cavities. The
magnetic field at the same location averaged over the horizontal and
transverse directions is found to be $\approx 3.97B_{0}$. Considering the
typical particle velocity to be $\approx 232v_{A}$ (local particle velocity
distribution peaks around this velocity), we find the particle-gyroradius to
be $\approx 58(c/\omega_{pi})$. The similarity in gyroradius and cavity size
indicates the instability is saturated.

Figure 5: Transversely averaged $x$-component of the current density
($J_{ion,x}$) in normalized unit. The dashed line indicates the location of
the shock and corresponding current is considered to calculate the theoretical
value of the most unstable Bell mode (Figure 4).

####  3.2.1 Role of instabilities in shock upstream

Figure 3 shows the magnetic field evolution close to the shock front. For
every time snap, the shock front is re-centered at $x=0$. It is seen that in
the shock precursor, along with the density, the magnetic field also becomes
turbulent. The power spectrum of the perpendicular magnetic field, depicted in
Figure 4, describes the nature of the magnetic turbulence. The maximum wave
number ($k_{max}$) is determined by the simulation resolution, whereas the
minimum wave number ($k_{min}$) is set by the largest length scale of the
chosen region. Alfvén waves, initiated by the energetic particles in the shock
upstream is natural to expect (Bell, 2004). These waves, resonating with the
gyrofrequency of the streaming energetic particles, generate an instability in
the system known as resonant streaming instability (Bell, 2004; Amato & Blasi,
2009). The mode of such instability achieves maximum growth rate at a
wavenumber $k$ satisfying $kr_{g}(E,B_{bg})=1$, where $E$ is the energy of the
particle gyrating with gyroradius $r_{g}$ in the background magnetic field of
strength $B_{bg}$ (Caprioli & Spitkovsky, 2014b). We find that in the region
the longest (for particles with $E=E_{sh}$) and shortest (particles with
$E=300E_{sh}$) wavenumbers of the resonant mode to be $3.74\times
10^{-2}(\frac{c}{\omega_{pi}})^{-1}$ and $2.13\times
10^{-3}(\frac{c}{\omega_{pi}})^{-1}$, respectively. While these two modes
encompass the peak of the power spectrum (Figure 4), wave mode corresponding
to the spectral peak resonates with the gyrofrequency of the particles having
energy $E=3E_{sh}$. However, the non-resonant hybrid (NRH) instability
introduced by Bell (2004) also may exist in this atmosphere. The theoretical
estimation of the fastest growing mode of this instability is found to be
$k_{Bell}=\frac{J_{ion,x}}{2B_{bg}}$, where $J_{ion,x}$ corresponds to the
transversely averaged horizontal component of the current, generated due to
the motion of the non-thermal protons, shown in Figure 5. Considering the
upstream shock vicinity value of $J_{ion,x}$ and the initial background
magnetic field $B_{bg}$, we estimate the fastest growing mode of the Bell
instability to be $k_{Bell}=1.59\times 10^{-2}(\frac{c}{\omega_{pi}})^{-1}$,
which nearly coincides with the wave mode corresponding to the peak of the
transverse magnetic power spectrum. Among the two modes the non-resonant mode
dominates when $(\dfrac{J_{ion,x}r_{g}}{B_{bg}})\gg 1$ (Bell, 2004; Amato &
Blasi, 2009). Using the local values of $J_{ion,x}$, $r_{g}$ and $B_{bg}$ of
the immediate shock upstream region, one finds the value of this non-
dimensional quantity to be $\approx 1.5$, which favors the NRH instability in
the region. Following the peak, the magnetic field power spectrum falls off
with the characteristic Kolmogorov slope ($k^{-5/3}$).

Figure 6 shows the magnetic field map in a region close to the shock front.
Instabilities described above generate magnetic fluctuations enhancing the
local magnetic field by as much as twice the initial value at the shock
precursor. Following Rankine–Hugoniot shock condition, it is expected that the
transverse magnetic components would get further amplified by a factor of 4.
In other words, the already pre-amplified magnetic field of the upstream would
be amplified further up to 8 times with respect to the background due to shock
compression, and propagate into the downstream region. However, the total
magnetic field contour in Figure 6 shows an amplification by a factor of up to
15 at some locations of downstream. Clearly, a different magnetic field
amplification mechanism is at play here. The initial shock surface gets
perturbed by the upstream density disturbance and gets corrugated. At the
corrugated surface, because of the non-aligned pressure and density gradient,
vorticity gets introduced and drives the Richtmeyer Mechkov Instability (RMI)
(Brouillette, 2002). The turbulence introduced by the RMI helps to enhance the
magnetic field further by stretching the field lines (Sano et al., 2012),
indicating the possible activation of turbulent dynamo behind the shock. The
role of turbulent dynamo in the shock downstream has also been reported by
earlier authors (e.g., Giacalone & Jokipii (2007); Mizuno et al. (2011)).

####  3.2.2 Magnetic field enhancement in the shock downstream

For further investigation in the turbulence saturation phase, the kinetic and
magnetic energy power spectrum of the region is plotted in Figure 7. In the
downstream region, $388(c/\omega_{pi})$ away from the shock front, a strip of
width $800(c/\omega_{pi})$ is considered to calculate the power spectrum. For
the kinetic energy spectrum the parameter $\rho^{1/2}[\bm{u}-\bar{\bm{u}}]$ is
employed (Podesta et al., 2007), whereas for the magnetic energy spectrum
$\bm{B}-\bm{\bar{B}}$ is considered. Here $\rho$, $\bm{u}$, $\bm{B}$,
$\bar{\bm{u}}$ and $\bar{\bm{B}}$ are mass density, velocity, magnetic field,
mean velocity and mean magnetic field of the region respectively. To calculate
the mean fields ($\bar{\bm{u}}$,$\bar{\bm{B}}$), a time window of about
$347~{}\Omega^{-1}$ is considered. Higher magnetic energy towards the small
scale is an indicator of possible small scale active dynamo in the region.

Figure 6: Total magnetic field contour close to the shock front at
$t=1170~{}\Omega^{-1}$

.

Figure 7: Dimensionless kinetic and magnetic energy power spectrum at
$t=1170~{}\Omega^{-1}$ in a region of width $8000(c/\omega_{pi})$ in the
downstream region $400(c/\omega_{pi})$ away from the shock front.

####  3.2.3 Particle energization

Figure 8: 2D energy spectrum, showing the particle energy distribution as a
function of position $x$. Here $E_{sh}$ is the shock energy given by
$E_{sh}=v_{u}^{2}/2$ where $v_{u}$ is the upstream plasma speed in the
downstream reference frame. Location of the shock front is indicated by the
white dashed line. High energy particles outrunning the shock are evident from
this 2D spectrum confirming the commencement of particle acceleration. Figure
9: Time evolution of the maximum energy of the particles (normalized with
respect to $E_{sh}$) present in the entire domain. The curve plateaus at large
times suggesting near-saturation. Figure 10: Evolution of downstream particle-
energy spectra ($E^{3/2}$ compensated), extracted from a region of width
$1000(c/\omega_{pi})$, just behind the shock. Flattening of the spectra
ensures DSA is an efficient mechanism in particle acceleration.

Turbulence in the shock upstream and downstream can scatter particles across
the shock and energize them. As mentioned earlier, we inject Maxwellian
distribution of particles with characteristic downstream temperature. A 2D
particle energy spectrum at the evolved state of the simulation is shown in
Figure 8. This figure depicts the particles’ energy distribution around the
shock. It clearly shows the existence of escaped high-energy particles in the
shock upstream.

Figure 9 displays the maximum energy attained by the protons over time. We ran
our simulation until maximum energy curve reaches its plateau. Figure 10
depicts the evolution of 1D energy spectrum ($E^{3/2}$ compensated) of the
particles. These spectra are extracted from a region of width
$1000(c/\omega_{pi})$ just behind the shock. During the initial phase (e.g.,
$50~{}\Omega^{-1}$), the spectrum almost follows the injected Maxwellian
distribution, eventually particles gain energy over time. During the near-
saturation phase, particles are seen to have gained energy up to $E=500E_{sh}$
($\approx 10$ MeV) in the downstream region just behind the shock. The major
part of the evolved spectrum follows $-3/2$ slope, indicating the active
presence of the DSA mechanism.

Figure 11: Saturated particle energy spectra for three different transverse
domain sizes. Particle energization is seen to increase with the domain size.
Figure 12: Saturated particle energy spectra for mono-energetic (Einj = 10Esh)
and maxwellian injection.

The maximum energy of the SEPs during CME shocks, however, has been seen to
have reached 100 MeV (Russell et al., 2013). A possible reason behind not
being able to produce such highly energetic protons in our simulation could be
the size of the simulation box, which restricts the longer wave modes. Figure
11 compares the energy spectra for three different simulations with varying
transverse box sizes, 250, 500 & 5000 $(c/\omega_{pi})$, during their near-
saturation phase. The tail of protons’ energy spectrum extends towards the
higher energy with an increase in transverse box size, thereby indicating the
role of the same in constraining the maximum energy gained by the protons.

In reality, solar wind demonstrates characteristic fluctuations. The ambient
solar wind of our present simulation is initiated without any turbulent
fluctuation. We envisage that introduction of such turbulent fluctuation in
the initial atmosphere may play an important role in the particle energization
process.

We have also compared our Maxwellian injection recipe with the mono-energetic
injection recipe, where $E_{inj}=10E_{sh}$, typically employed by previous
investigators (Bai et al., 2015; Mignone et al., 2018). In section 2.2, we
have already mentioned that in our present simulation, we choose to have a
Maxwellian injection with downstream plasma temperature. For this case we had
carefully chosen the combination of $\eta=3\times 10^{-3}$ and the number of
particles per cell to be 10, so that the particles do not disrupt the shock.
On the other hand, the mono-energetic injection with the above combination of
$\eta$ and the number of particles per cell disrupts the shock indicating
over-injection. We, therefore, have chosen to bring down $\eta$ and also the
number of particles per cell for the mono-energetic injection case to $2\times
10^{-3}$ and $4$, respectively. A comparison between the two energy spectra in
the case of two injections are shown in Figure 12. In the case of Maxwellian
injection, particles are seen to have gained energy up to $500E_{sh}$
(equivalent to 10 MeV), whereas for mono-energetic injection recipe particles
have gained energy up to $1000E_{sh}$ (equivalent to 20 MeV). The possible
reason behind this discrepancy is that the particles have wide energy
distribution in Maxwellian injection. Lower energy particles in the
distribution have low probability of participating in the energization
process. On the other hand, in mono-energetic injection, all particles being
already suprathermal, the entire population may participate in the
acceleration process and thus may end up with a higher energy tail in the
energy spectrum.

###  3.3 Quasi-perpendicular shock

It is well known that the particle acceleration mechanism in quasi-
perpendicular shocks is much different than the parallel shocks. We,
therefore, simulate a quasi-perpendicular shock with the same Alfvénic Mach
number ($\approx 19$). This time, the initial magnetic field makes an angle
$\psi=75^{\circ}$ with the shock normal direction, making the shock quasi-
perpendicular. The corresponding magnetic field components are now
$B_{x}=B_{bg}\cos(\psi)$ and $B_{y}=B_{bg}\sin(\psi)$ with $B_{bg}=3.55B_{0}$
as mentioned in section 3.1. Initially particles are injected following the
same prescription mentioned in section 2.2. Figure 13 shows the particle
energy spectrum at the end of the simulation. It shows that the particles only
get energized up to $40-45E_{sh}$.

Earlier hybrid-PIC (Caprioli & Spitkovsky, 2014a) simulation has shown a
decrease in acceleration efficiency for quasi-perpendicular shock. Naturally,
this suggests the injection of less number of suprathermal particles in the
present case. We ran our quasi-perpendicular simulation for two more
combinations of $\eta$ and the number of protons per cell, keeping the mass
density of individual protons same as our previous injection. The two sets
include $\eta=1.5\times 10^{-3}$ and $6\times 10^{-4}$ with $5$ and $2$
protons per cell, respectively. The result is compared with our previous
simulation where we considered $\eta=3\times 10^{-3}$ and 10 protons per cell.
Our new choices ensure less amount of suprathermal particles in the domain.
Though the resultant energy spectra from all three combinations vary in terms
of the total number of protons, all of them show a very similar energization
(Figure 14), indicating the fact that the population of injected suprathermal
protons does not affect the overall energization process in quasi-
perpendicular shock. The following analyses are performed on the simulation
done with the injection prescribed in Section 2.2.

To understand the energization process, we need to understand which velocity
component gains energy. Figure 15 depicts 2-D particle velocity distribution
of the downstream region. Both the velocity distribution functions show
anisotropic nature. Figure 15 (a) shows that the velocity distribution
$V_{x}-V_{y}$ is more skewed along an axis making $\approx 4^{\circ}$ with the
$V_{x}$ axis, while $V_{y}-V_{z}$ is more skewed along the $V_{z}$ axis. This
indicates particles are getting accelerated on a plane making an angle
$\approx 4^{\circ}$ with the X-Z plane.

We further investigated the mean magnetic field direction of the downstream
region of interest (ROI) during the simulation’s saturation phase. Here we
remind the reader that this is a 2.5D simulation. Therefore, even though
computation is performed only on a 2D plane, all three components of a vector
field are evolved in the simulation. As our simulation is performed in the X-Y
plane, we assume that the local magnetic field $\bm{B}$ makes an angle
$\theta$ with the Y-axis. To accommodate the 3rd component of the magnetic
field, we also assume $\phi$ to be the angle between $\bm{B}$’s projection
onto the X-Z plane and the Y-Z plane. We derive the value of $\theta$ & $\phi$
values at every grid of the ROI and the respective histograms are plotted to
see the overall orientation of the field. Figures 16 (a) & (b) demonstrate
that in the ROI $\phi\approx 90^{\circ}$ is the dominant angle, which implies
the mean $\bm{B}$ lies on the X-Y plane. Figure 16 (c) & (d) demonstrates that
in the same downstream region $\theta\approx 4^{\circ}$ is the dominant angle.
Overall Figure 16 concludes that the mean $\bm{B}$ lies in the X-Y plane and
makes an angle $\approx 4^{\circ}$ with the Y-axis.

After knowing the accelerated components of the particles’ velocity and the
orientation of the mean magnetic field in the ROI, it is now understood that
particles are predominantly getting accelerated in the plane perpendicular to
the mean magnetic field (i.e., the plane making an angle $\approx 4^{\circ}$
with the X-Z plane). The convective electric field
($\bm{E}=-\bm{u}\times\bm{B}$) is the reason behind this. The field acts on
the plane perpendicular to the mean magnetic field and accelerates protons on
the same plane. Concurrently, the local grad B force, developed across the
shock front due to shock compression, makes the particles drift on the plane
of the shock. Additionally, particles’ velocity components perpendicular to
the mean magnetic field, namely the $V_{\perp}$ component, increases in
magnitude. The parallel component $V_{\parallel}$ still maintains the injected
distribution, resulting in an anisotropic 2-D velocity distribution behind the
shock. The mechanism is popularly known as SDA and plays a major role in
accelerating protons in our quasi-perpendicular shock simulation.

We should point out that although the Alfvénic Mach of both the parallel and
quasi-perpendicular shocks are same in our above mentioned simulations, the
energy gained by the protons in the quasi-perpendicular shock case is
significantly less. This can be explained with Figure 17. Only a few high
energy protons are seen to have escaped the shock front and move into the
upstream region. After getting accelerated by SDA, protons are mostly advected
into the downstream region within a few gyrocycles, resulting in the
truncation of further energization process. The few high energy protons that
managed to escape the shock front are insufficient to generate strong
fluctuations in the upstream region. Nearly non-fluctuating density
($\delta\rho/\rho_{0}\approx 0$) in the upstream implies that there is no
active participation of DSA in our simulation to accelerate protons. However,
in reality, inherent fluctuations of the SW can excite DSA and may accelerate
particles further. This will be investigated in our future work.

Figure 13: Injected Maxwellian and saturated energy spectra for quasi-
perpendicular shock with Mach number $\approx$ 19. Figure 14: Energy spectra
derived from the simulation of quasi-perpendicular shock with $M_{A}\approx
19$, by changing injection fraction ($\eta$) and particles per cell (ppc).

Figure 15: 2D velocity distribution function of the quasi-perpendicular shock
with $M_{A}\approx$ 19\. Particles were collected from just behind the shock.
The region is marked in Figure 16

.

Figure 16: The orientation of the magnetic field just behind the shock. The
particular region of interest is marked with vertical dashed lines (a) Contour
plot of the angles ($\phi$) between the Y-Z plane and the projection of
$\bm{B}$ onto the X-Z plane. (b) A histogram depicting the dominant $\phi$
angle of the region of interest, showing that the mean magnetic field of the
region lies on the X-Y plane. (c) Contour plot of the angle between $\bm{B}$
and the Y axis. (d) Histogram plot of the region of interest to show the
dominant $\theta$ angle showing the mean field makes an angle $\approx
4^{\circ}$ with the Y axis. Particle velocity distribution function of the
same region is plotted in Figure 15

.

Figure 17: 2D energy spectrum, showing the particle energy distribution as a
function of position $x$. Here $E_{sh}$ is the shock energy given by
$E_{sh}=v_{u}^{2}/2$ where $v_{u}$ is the upstream plasma speed in the
downstream reference frame. Location of the shock front is indicated by the
white dashed line. As evident from this plot, high energy particles outrunning
the shock are very less in number and hence cannot seed strong fluctuations in
the upstream.

##  4 ICME simulation for low mach shock

ICME-shocks with low Alfvénic Mach numbers (Berdichevsky et al., 2000; Oh et
al., 2007; Lugaz et al., 2015) are more common. Keeping that in mind, we also
have simulated ICME-shocks with low Alfvénic Mach $\approx 10$ and $4$.
Generation of shock and particle injection recipes are identical to the cases
discussed in section 3. The only parameter that is changed in order to produce
shocks with low Machs is the relative velocity between the SW and ICME plasma.
For Alfvénic Mach $\approx 10$, we have injected a SW flow with relative speed
$-50.05\,v_{A}$ $(=-976$ km s${}^{-1})$, whereas for Alfvénic Mach $\approx 4$
the relative speed is $-20\,v_{A}$ $(=-390$ km s${}^{-1})$. Both these shocks
are parallel in nature.

Figure 18 exhibits the particle energy spectrum for the $M_{A}\approx 10$
shock. Particles are seen to get energised up to $300E_{sh}$ in this
simulation. The flat part of the spectrum directs towards the DSA mechanism
for particle acceleration. On the other hand, Figure 19, which depicts
particle energy spectrum for the shock with Alfvénic Mach $\approx 4$, shows
no indication of particle acceleration.

Therefore, these sets of numerical simulations follow our intuitive guess that
Alfvénic Mach number is a good indicator of shocks’ capability to accelerate
particles. Low Mach shocks may not work as good accelerators.

Figure 18: Evolution of downstream spectra for a parallel shock with $M_{A}$
$\approx 10$, calculated over a region of width $1000(c/\omega_{pi})$ just
behind the shock. Figure 19: Evolution of downstream spectra for a parallel
shock with $M_{A}$ $\approx 4$, calculated over a region of width
$1000(c/\omega_{pi})$ just behind the shock.

##  5 Summary and discussion

This paper has investigated the particle acceleration process in ICME shocks.
We used the MHD-PIC module of the PLUTO code (Mignone et al., 2018). The
neutral background plasma is treated as a magnetohydrodynamic fluid, and non-
thermal proton dynamics is taken care of using a PIC code. The reactions of
particle dynamics onto the MHD fluid and vice versa are appropriately managed.
Compared to earlier heliospheric shock models (Zank et al., 2000), this model
is more local in nature. Unlike the former, it can only simulate a small
region across the shock. But, particles evolve self consistently on the
background thermal plasma. This hybrid model is more similar to Gargaté et al.
(2014)’s model, but computationally much less expensive. Such a hybrid
simulation enables us to follow the particle dynamics and fluid evolution
until the particle energy spectrum nearly saturates. These simulations empower
us to predict various in situ observations and help us plan future missions.

We started our investigation with a CME shock having an extremely high
Alfvénic Mach number ($M_{A}\approx 19$). Even though such a strong shock is
rare in the heliosphere, one was evidenced lately (on 2012 July 23). The
particle acceleration process has been probed initially for a parallel shock,
where the shock propagation direction exactly matches with the ambient solar
wind magnetic field. This being an idealistic case, we also have studied a
quasi-perpendicular shock, where the ambient magnetic field forms a
$75^{\circ}$ angle with the propagation direction. CME shocks with lower Mach
numbers (e.g., $M_{A}\approx$ 10 & 4) have also been investigated in the
present study. Particles are injected at the shock front in all our
simulations. During injection, particles’ distribution follows a Maxwellian
with characteristic downstream temperature (i.e., energy $E_{sh}$). In the
following, we outline our results as well as suggest ways to further explore
the observation data that may confirm various physical processes occuring
around such CME shocks.

Instability and dynamo in the parallel shock – The first simulation, which
delineates a plane parallel shock with high Alfvénic Mach ($M_{A}\approx 19$),
demonstrates that such a shock can energize protons up to $500E_{sh}(\approx
10MeV)$. The slope of the particle energy spectrum (Figure 10) indicates the
role of DSA in the particle acceleration process. We also note that once
particles get energized, many of them leave the shock front and move to the
far upstream region (Figure 8).

The upstream of the shock is seen to become turbulent with the growth of
density cavities. Such cavities are due to the development of non-resonant
Bell instability in the upstream region. As the instability grows, the sizes
of these cavities increase. Once the local ion-gyroradius becomes comparable
to the cavity size, the instability ceases to grow (Figure 2). We have further
investigated the existence of resonant streaming instability in the domain. We
find that the longest and shortest resonant modes encompass the peak of the
magnetic power spectrum. However, the fastest growing mode of non-resonant
Bell-instability turns out to be responsible for the peak of the power
spectrum of the region’s transverse magnetic field (Figure 4).

The turbulent magnetic field of the shock downstream region is seen to have
enhanced more than is expected from the shock compression, hinting at the
presence of a dynamo mechanism in the region. Due to non-aligned pressure and
density gradient at the corrugated shock front, instability like Richtmyer-
Meshkov instability may kick in and develop turbulence. Such turbulence
stretches the field lines and thereby enhances the local magnetic field.
Kinetic or magnetic diffusivity are not explicitly specified in our
simulation. Nevertheless, when we compare the kinetic and magnetic energy
spectra (Figure 7) the dominant magnetic energy towards the small scale
($k>0.07~{}(c/\omega_{pi})^{-1}$) indicates the development of small scale
dynamo in the region.

Quasi-perpendicular shock demonstrates velocity anisotropy – We have also
performed simulation of quasi-perpendicular shock with a high Alfvénic Mach
number ($M_{A}\approx 19$). The particle energy spectrum shown in Figure 13
indicates energization up to $40-45E_{sh}$ which is much smaller than the
particle energisation in the parallel shock with the same Alfvénic Mach. In
quasi-perpendicular shock, the particle energization is expected due to Shock
Drift Acceleration process (SDA), where the convective electric field
$\bm{E}=-\bm{u}\times\bm{B}$ plays a vital role in accelerating the particles.
Activation of SDA in the present case is evident from the particles’ velocity
distribution as shown in Figure 15. The distribution shows the development of
anisotropy on the plane perpendicular to the local mean magnetic field.

Our simulations demonstrate that, in quasi-perpendicular shocks (Figure 17),
high energetic particles rarely escape to the upstream and therefore do not
develop strong turbulence in the region, as it happens in the parallel shock.
This also prevents the particles from participating in DSA and accelerating
further. However, once again, the SW turbulence may play a crucial role in
accelerating these particles, which is not incorporated in the present
simulations.

Low Mach shocks are less effective to accelerate particles – As expected, our
simulations of low-Mach number shocks show no effective particle acceleration.
We have simulated parallel shocks with Alfvénic Mach, $M_{A}\approx$ 10 and 4.
Though the shock with $M_{A}\approx$ 10 could energize particles upto
$300E_{sh}$, the shock with $M_{A}\approx$ 4 does not show evidence of
particle acceleration.

It is worth mentioning that our simulations were unable to produce particles
with very high energy (e.g., up to 100 MeV, as in the case of real
observations Russell et al. (2013)). Multiple reasons can be responsible for
this. While traveling from the Sun to 1 AU, the CME shocks encounter variable
backgrounds over a long duration of time. This may help particles to get
energized continuously. On the other hand, our localized simulations are
provided with a background atmosphere of 1 AU, and therefore can only energize
particles to whatever energies is possible over this background. Furthermore,
our present simulations start with a uniform solar wind background. The
inherent fluctuations of the real solar wind may energize particles to a
greater extent, which needs to be investigated in the future. The restricted
input energy distribution of the particles may also be responsible for the
limited energization process. In the simulations, we inject particles with a
Maxwellian distribution having a characteristic temperature of the downstream.
However, these may not be the only particles getting induced at the shock
front in reality. Some amount of much higher energetic particles that are
already energized due to magnetic reconnection may also get induced at the
shock front for further processing and may give rise to the observed high
energy tails. Such a possibility has been proposed earlier (Klein & Trottet,
2001). As others (e.g., Bai et al. (2015); Mignone et al. (2018)), we also
have attempted to inject particles with $10E_{sh}$. The resultant energy
distribution gets extended towards the higher end but only by a few MeVs as
compared to the Maxwellian injection (Figure 12). We have also shown that the
limited size of the simulation domain may restrict the growth of more extended
modes in the domain and regulate energy growth beyond a threshold.

Understanding particle acceleration in shocks is a major topic in overall
astrophysics. Shocks of different strengths and magnetic configurations are
ubiquitous in the heliosphere. Heliospheric shocks can be quasi-parallel or
quasi-perpendicular, depending on the relative angle between their propagation
direction and ambient upstream magnetic field. Therefore, shocks that are
potential accelerators (relatively high $M_{A}$) should be able to demonstrate
the characteristics of both DSA and SDA mechanisms. The issue of particle
acceleration, therefore, can be explored directly by using heliospheric in
situ data. In the following paragraph, we highlight a few observables that can
confirm our current understanding of the subject.

In principle, SEP energy spectra should be able to demonstrate the
characteristic -3/2 slope if DSA is the dominant mechanism. On the other hand,
the 2D velocity distribution of the particles should be able to demonstrate
anisotropy if SDA plays a pivotal role in accelerating particles. The
existence of non-resonant Bell instability can be verified by measuring the
size of the upstream density cavities or by analyzing the upstream magnetic
turbulence. It will also be interesting to see if the downstream magnetic
field gets enhanced due to turbulent dynamo.

It is possible to perform local numerical simulations of wide ranges of ICME
shocks and compare the simulation outputs with real observations. ICME shocks
with varied strength at different heliospheric distances can be compared with
the simulation results. We can also corroborate the effect of different solar
wind background plasma on ICME shocks. We envisage data from the Integrated
Science Investigation of the Sun (IS$\odot$IS) and Electromagnetic Fields
Investigation (FIELD) experiments on board the Parker Solar Probe (PSP), the
Energetic Particle Detector (EPD), Solar Wind Plasma Analyser (SWA) and
Magnetometer (MAG) on board Solar Orbiter or the Aditya Solarwind Particle
EXperiment (ASPEX) on board Indian Space Research Organisation’s Aditya-L1
mission will be able to probe in this direction.

Computations were carried out on the Physical Research Laboratory’s VIKRAM
cluster. We thank the anonymous referee for her/his insightful comments, which
helped us improve the manuscript. BV would like to thank the support provided
by the Max Planck Partner Group established at IIT Indore.

## References

  * Amato & Blasi (2009) Amato, E., & Blasi, P. 2009, MNRAS, 392, 1591, doi: 10.1111/j.1365-2966.2008.14200.x
  * Armstrong et al. (1985) Armstrong, T. P., Pesses, M. E., & Decker, R. B. 1985, Washington DC American Geophysical Union Geophysical Monograph Series, 35, 271, doi: 10.1029/GM035p0271
  * Bai et al. (2015) Bai, X.-N., Caprioli, D., Sironi, L., & Spitkovsky, A. 2015, ApJ, 809, 55, doi: 10.1088/0004-637X/809/1/55
  * Baring et al. (1997) Baring, M. G., W. Ogilvie, C. Ellison, & Forsyth, R. J. 1997, ApJ, 476, 889, doi: 10.1086/303645
  * Bell (1978) Bell, A. R. 1978, MNRAS, 182, 147, doi: 10.1093/mnras/182.2.147
  * Bell (2004) —. 2004, MNRAS, 353, 550, doi: 10.1111/j.1365-2966.2004.08097.x
  * Bell (2005) —. 2005, MNRAS, 358, 181, doi: 10.1111/j.1365-2966.2005.08774.x
  * Berdichevsky et al. (2000) Berdichevsky, D. B., Szabo, A., Lepping, R. P., Viñas, A. F., & Mariani, F. 2000, Journal of Geophysical Research: Space Physics, 105, 27289, doi: 10.1029/1999JA000367
  * Blandford & Ostriker (1978) Blandford, R. D., & Ostriker, J. P. 1978, ApJ, 221, L29, doi: 10.1086/182658
  * Brouillette (2002) Brouillette, M. 2002, Annual Review of Fluid Mechanics, 34, 445, doi: 10.1146/annurev.fluid.34.090101.162238
  * Caprioli & Spitkovsky (2013) Caprioli, D., & Spitkovsky, A. 2013, ApJ, 765, L20, doi: 10.1088/2041-8205/765/1/L20
  * Caprioli & Spitkovsky (2014a) Caprioli, D., & Spitkovsky, A. 2014a, The Astrophysical Journal, 783, 91, doi: 10.1088/0004-637x/783/2/91
  * Caprioli & Spitkovsky (2014b) —. 2014b, The Astrophysical Journal, 794, 46, doi: 10.1088/0004-637x/794/1/46
  * Caprioli & Spitkovsky (2014c) —. 2014c, The Astrophysical Journal, 794, 47, doi: 10.1088/0004-637x/794/1/47
  * Decker (1983) Decker, R. B. 1983, J. Geophys. Res., 88, 9959, doi: 10.1029/JA088iA12p09959
  * Decker (1988) —. 1988, Space Sci. Rev., 48, 195, doi: 10.1007/BF00226009
  * Desai & Giacalone (2016) Desai, M., & Giacalone, J. 2016, Living Reviews in Solar Physics, 13, 3, doi: 10.1007/s41116-016-0002-5
  * Fermi (1949) Fermi, E. 1949, Physical Review, 75, 1169, doi: 10.1103/PhysRev.75.1169
  * Fermi (1954) —. 1954, ApJ, 119, 1, doi: 10.1086/145789
  * Feynman & Gabriel (2000) Feynman, J., & Gabriel, S. B. 2000, J. Geophys. Res., 105, 10543, doi: 10.1029/1999JA000141
  * Gargaté et al. (2007) Gargaté, L., Bingham, R., Fonseca, R. A., & Silva, L. O. 2007, Computer Physics Communications, 176, 419, doi: 10.1016/j.cpc.2006.11.013
  * Gargaté et al. (2014) Gargaté, L., Fonseca, R. A., Silva, L. O., Bamford, R. A., & Bingham, R. 2014, The Astrophysical Journal, 792, 9, doi: 10.1088/0004-637x/792/1/9
  * Giacalone & Jokipii (2007) Giacalone, J., & Jokipii, J. R. 2007, ApJ, 663, L41, doi: 10.1086/519994
  * Gopalswamy (2003) Gopalswamy, N. 2003, Advances in Space Research, 31, 869, doi: 10.1016/S0273-1177(02)00888-8
  * Helder et al. (2012) Helder, E. A., Vink, J., Bykov, A. M., et al. 2012, Space Sci. Rev., 173, 369, doi: 10.1007/s11214-012-9919-8
  * Hu et al. (2017) Hu, J., Li, G., Ao, X., Zank, G. P., & Verkhoglyadova, O. 2017, Journal of Geophysical Research (Space Physics), 122, 10,938, doi: 10.1002/2017JA024077
  * Kahler (2001) Kahler, S. W. 2001, Journal of Geophysical Research: Space Physics, 106, 20947, doi: 10.1029/2000JA002231
  * Kahler (2004) —. 2004, The Astrophysical Journal, 603, 330, doi: 10.1086/381358
  * Kahler & Vourlidas (2005) Kahler, S. W., & Vourlidas, A. 2005, Journal of Geophysical Research: Space Physics, 110, doi: 10.1029/2005JA011073
  * Kahler & Vourlidas (2013) —. 2013, The Astrophysical Journal, 769, 143, doi: 10.1088/0004-637x/769/2/143
  * Kilpua et al. (2017) Kilpua, E., Koskinen, H. E. J., & Pulkkinen, T. I. 2017, Living Reviews in Solar Physics, 14, 5, doi: 10.1007/s41116-017-0009-6
  * Klein & Dalla (2017) Klein, K.-L., & Dalla, S. 2017, Space Sci. Rev., 212, 1107, doi: 10.1007/s11214-017-0382-4
  * Klein & Trottet (2001) Klein, K.-L., & Trottet, G. 2001, Space Sci. Rev., 95, 215
  * Kouloumvakos et al. (2019) Kouloumvakos, A., Rouillard, A. P., Wu, Y., et al. 2019, ApJ, 876, 80, doi: 10.3847/1538-4357/ab15d7
  * Krymskii (1977) Krymskii, G. F. 1977, Akademiia Nauk SSSR Doklady, 234, 1306
  * Lee & Fisk (1982) Lee, M. A., & Fisk, L. A. 1982, Space Sci. Rev., 32, 205, doi: 10.1007/BF00225185
  * Li et al. (2012) Li, G., Shalchi, A., Ao, X., Zank, G., & Verkhoglyadova, O. P. 2012, Advances in Space Research, 49, 1067, doi: 10.1016/j.asr.2011.12.027
  * Li et al. (2005a) Li, G., Zank, G. P., Desai, M. I., Mason, G. M., & Rice, W. K. M. 2005a, Washington DC American Geophysical Union Geophysical Monograph Series, 156, 51, doi: 10.1029/156GM07
  * Li et al. (2003) Li, G., Zank, G. P., & Rice, W. K. M. 2003, Journal of Geophysical Research (Space Physics), 108, 1082, doi: 10.1029/2002JA009666
  * Li et al. (2005b) —. 2005b, Journal of Geophysical Research (Space Physics), 110, A06104, doi: 10.1029/2004JA010600
  * Lipatov (2002) Lipatov, A. S. 2002, The hybrid multiscale simulation technology: an introduction with application to astrophysical and laboratory plasmas
  * Lugaz et al. (2015) Lugaz, N., Farrugia, C. J., Smith, C. W., & Paulson, K. 2015, Journal of Geophysical Research (Space Physics), 120, 2409, doi: 10.1002/2014JA020848
  * Mignone et al. (2018) Mignone, A., Bodo, G., Vaidya, B., & Mattia, G. 2018, ApJ, 859, 13, doi: 10.3847/1538-4357/aabccd
  * Mizuno et al. (2011) Mizuno, Y., Pohl, M., Niemiec, J., et al. 2011, ApJ, 726, 62, doi: 10.1088/0004-637X/726/2/62
  * Oh et al. (2007) Oh, S. Y., Yi, Y., & Kim, Y. H. 2007, Sol. Phys., 245, 391, doi: 10.1007/s11207-007-9042-2
  * Podesta et al. (2007) Podesta, J. J., Roberts, D. A., & Goldstein, M. L. 2007, ApJ, 664, 543, doi: 10.1086/519211
  * Reames (1999) Reames, D. V. 1999, Space Sci. Rev., 90, 413, doi: 10.1023/A:1005105831781
  * Reames (2002) —. 2002, ApJ, 571, L63, doi: 10.1086/341149
  * Reames (2012) —. 2012, ApJ, 757, 93, doi: 10.1088/0004-637X/757/1/93
  * Reville & Bell (2012) Reville, B., & Bell, A. R. 2012, MNRAS, 419, 2433, doi: 10.1111/j.1365-2966.2011.19892.x
  * Rice et al. (2003) Rice, W. K. M., Zank, G. P., & Li, G. 2003, Journal of Geophysical Research (Space Physics), 108, 1369, doi: 10.1029/2002JA009756
  * Riley et al. (2016) Riley, P., Caplan, R. M., Giacalone, J., Lario, D., & Liu, Y. 2016, ApJ, 819, 57, doi: 10.3847/0004-637X/819/1/57
  * Romero et al. (2017) Romero, G. E., Boettcher, M., Markoff, S., & Tavecchio, F. 2017, Space Sci. Rev., 207, 5, doi: 10.1007/s11214-016-0328-2
  * Russell et al. (2013) Russell, C. T., Mewaldt, R. A., Luhmann, J. G., et al. 2013, The Astrophysical Journal, 770, 38, doi: 10.1088/0004-637x/770/1/38
  * Russell et al. (2013) Russell, C. T., Mewaldt, R. A., Luhmann, J. G., et al. 2013, ApJ, 770, 38, doi: 10.1088/0004-637X/770/1/38
  * Ryan et al. (2000) Ryan, J. M., Lockwood, J. A., & Debrunner, H. 2000, Space Sci. Rev., 93, 35, doi: 10.1023/A:1026580008909
  * Sano et al. (2012) Sano, T., Nishihara, K., Matsuoka, C., & Inoue, T. 2012, ApJ, 758, 126, doi: 10.1088/0004-637X/758/2/126
  * Sheeley et al. (1983) Sheeley, N. R., J., Howard, R. A., Koomen, M. J., et al. 1983, in NASA Conference Publication, Vol. 228, NASA Conference Publication, 0.693
  * Spitkovsky (2005) Spitkovsky, A. 2005, in American Institute of Physics Conference Series, Vol. 801, Astrophysical Sources of High Energy Particles and Radiation, ed. T. Bulik, B. Rudak, & G. Madejski, 345–350, doi: 10.1063/1.2141897
  * van Marle et al. (2018) van Marle, A. J., Casse, F., & Marcowith, A. 2018, MNRAS, 473, 3394, doi: 10.1093/mnras/stx2509
  * Verkhoglyadova et al. (2009) Verkhoglyadova, O. P., Li, G., Zank, G. P., Hu, Q., & Mewaldt, R. A. 2009, ApJ, 693, 894, doi: 10.1088/0004-637X/693/1/894
  * Verkhoglyadova et al. (2010) Verkhoglyadova, O. P., Li, G., Zank, G. P., et al. 2010, Journal of Geophysical Research (Space Physics), 115, A12103, doi: 10.1029/2010JA015615
  * Zank et al. (2006) Zank, G. P., Li, G., Florinski, V., et al. 2006, Journal of Geophysical Research (Space Physics), 111, A06108, doi: 10.1029/2005JA011524
  * Zank et al. (2004) —. 2004, Journal of Geophysical Research (Space Physics), 109, A04107, doi: 10.1029/2003JA010301
  * Zank et al. (2007) Zank, G. P., Li, G., & Verkhoglyadova, O. 2007, Space Sci. Rev., 130, 255, doi: 10.1007/s11214-007-9214-2
  * Zank et al. (2000) Zank, G. P., Rice, W. K. M., & Wu, C. C. 2000, J. Geophys. Res., 105, 25079, doi: 10.1029/1999JA000455


# Access Control Using Spatially Invariant Permutation of Feature Maps for
Semantic Segmentation Models

Hiroki Ito, MaungMaung AprilPyone, and Hitoshi Kiya Tokyo Metropolitan
University, Japan

###### Abstract

In this paper, we propose an access control method that uses the spatially
invariant permutation of feature maps with a secret key for protecting
semantic segmentation models. Segmentation models are trained and tested by
permuting selected feature maps with a secret key. The proposed method allows
rightful users with the correct key not only to access a model to full
capacity but also to degrade the performance for unauthorized users.
Conventional access control methods have focused only on image classification
tasks, and these methods have never been applied to semantic segmentation
tasks. In an experiment, the protected models were demonstrated to allow
rightful users to obtain almost the same performance as that of non-protected
models but also to be robust against access by unauthorized users without a
key. In addition, a conventional method with block-wise transformations was
also verified to have degraded performance under semantic segmentation models.

##  1 Introduction

Deep neural networks (DNNs) have led to major breakthroughs in computer vision
for a wide range of applications. Convolutional neural networks (CNNs) are a
type of DNN. Current commercial applications for image recognition, object
detection, and semantic segmentation are primarily powered by CNNs [1, 2].
Therefore, CNNs have become the de facto standard for visual recognition
systems for many different applications. However, training successful CNNs
requires three ingredients: a huge amount of data, GPU-accelerated computing
resources, and efficient algorithms, and this is not a trivial task.
Therefore, trained CNNs have great business value. Considering the expenses
necessary for the expertise, money, and time taken to train a CNN model, a
model should be regarded as a kind of intellectual property (IP).

There are two aspects of IP protection for DNN models: ownership verification
and access control. The former focuses on identifying the ownership of models,
and the latter addresses protecting the functionality of DNN models from
unauthorized access. Ownership verification methods are inspired by digital
watermarking, and they embed watermarks into DNN models so that the embedded
watermarks can be used to verify the ownership of the models in question [3,
4, 5, 6, 7, 8, 9, 10, 11].

Although the above watermarking methods can facilitate the identification of
the ownership of models, in reality, a stolen model can be exploited in many
different ways. For example, an attacker can use a model for their own benefit
without arousing suspicion, or a stolen model can be used for model inversion
attacks [12] and adversarial attacks [13]. Therefore, it is crucial to
investigate mechanisms for protecting DNN models from unauthorized access and
misuse. In this paper, we focus on protecting a model from misuse when it has
been stolen (i.e., access control).

A method for protecting a model against unauthorized access was inspired by
adversarial examples [13, 14, 15] and image encryption [16, 17, 18, 19, 20,
21], and it was proposed to utilize secret perturbation to control the access
of a model [22]. Another study introduced a secret key for protecting a model
[23, 24]. The secret key-based protection method [24] uses a key-based
transformation that was originally used by an adversarial defense in [25],
which was in turn inspired by perceptual image encryption methods [20, 26, 21,
27, 28, 29, 30, 31]. This model protection method utilizes a secret key in
such a way that a stolen model cannot be used to its full capacity without a
correct secret key.

However, these methods were evaluated only on image classification tasks, and
it is not known how well they perform on other advanced tasks. Therefore, in
this paper, we consider protecting semantic segmentation models from
unauthorized access for the first time, and we propose a novel access control
method that uses the spatially invariant permutation of feature maps on the
basis of a secret key. The proposed method allows rightful users with the
correct key to access a model to full capacity and degrade the performance for
unauthorized users. In experiments on semantic segmentation, we evaluated the
access control performance of models trained by using feature map permutation.
The results show that the protected models provided almost the same
segmentation performance as that of non-protected models against authorized
access, while the segmentation accuracy seriously dropped when an incorrect
key was given. Furthermore, the conventional method proposed for image
classification tasks [24] was demonstrated to have degraded performance under
semantic segmentation tasks.

##  2 Proposed Method

###  2.1 Overview

Figure 1 illustrates an overview of access control for protecting semantic
segmentation models from unauthorized access.

Figure 1: Overview of access control

A protected model is prepared by training a network with secret key $K$.
Authorized users input test images into the protected model with correct key
$K$, in which the model provides almost the same segmentation map as that
predicted by using a non-protected model. In contrast, when unauthorized users
who do not know key $K$ input test images into a protected model without any
key or with incorrect (estimated) key $K^{\prime}$, the model provides a
degraded map.

###  2.2 Semantic Segmentation

The goal of semantic segmentation is to understand what is in an image at the
pixel level. Figure 2 shows an example of semantic segmentation. The
segmentation model predicts a segmentation map from an input image, where each
pixel in the segmentation map represents a class label.

Figure 2: Example of semantic segmentation

The mean intersection-over-union (mean IoU) [32, 33] is used as a metric for
evaluating the segmentation performance. The mean IoU value of an image is
calculated by averaging the IoU of each class, which is defined by

$\mathrm{IoU}=\frac{TP}{TP+FP+FN},$ (1)

where $TP$, $FP$, and $FN$ mean true positive, false positive, and false
negative values calculated from a predicted segmentation map and ground truth
one, respectively. In addition, the metric ranges from zero to one, where a
value of one means that a predicted segmentation map is the same as that of
the ground truth, and a value of zero indicates they have no overlap.

###  2.3 Training Model with Key $K$

To protect semantic segmentation models, we train models by randomly permuting
feature maps with secret key $K$. The permutation is applied to feature maps
in a network as in Fig. 3. In the figure, a fully convolutional network (FCN)
[33] with a ResNet-50 [34] backbone is illustrated as an example, although the
permutation is not limited to the FCN. There are six feature maps in Fig. 3,
and a number of feature maps from the six maps are chosen to be permuted prior
to permutation. In this paper, a feature map $x$ with a dimension of ($c\times
h\times w$), where $c$ is the number of channels, $h$ is the height, and $w$
is the width of the feature map, is transformed with key $K$ at each iteration
for training a model. There are two steps in the process of transforming a
feature map as below (see Fig. 4).

  1. 1.

Generate a random vector with a size of $c$ such that

$[\alpha_{1},.,\alpha_{i},.,\alpha_{i^{\prime}},\ldots,\alpha_{c}],\alpha_{i}\in\\{1,\ldots,c\\},$
(2)

where $\alpha_{i}\neq\alpha_{i^{\prime}}$ if $i\neq i^{\prime}$.

  2. 2.

Replace all elements of $x$, $x(i,j,k)$, $i\in\\{1,\ldots,c\\}$,
$j\in\\{1,\ldots,h\\}$, and $k\in\\{1,\ldots,w\\}$ with $x(a_{i},j,k)$ to
produce permuted feature map $x^{\prime}$, where an element of $x^{\prime}$,
$x^{\prime}(i,j,k)$ is equal to $x(a_{i},j,k)$.

If multiple feature maps are chosen to be permuted, the above steps are
applied to each feature map.

The above feature map permutation is spatially invariant as illustrated in
Fig. 4. This spatially invariant property is important when the permutation is
applied to applications that are required to output images like semantic
segmentation. In contrast, an example of spatially variant permutation is
given in Fig. 5. The conventional protection method for image classification
[24] is carried out by using a spatially variant permutation method, so it is
not available for protecting other models such as semantic segmentation models
as described later. In this paper, a model protection method for semantic
segmentation is discussed for the first time.

Figure 3: Semantic segmentation model (FCN with ResNet-50 backbone) with
feature map permutation

Figure 4: Spatially invariant feature map permutation

Figure 5: Spatially variant feature map permutation

###  2.4 Applying Queries to Model

As shown in Fig. 1, authorized users have key $K$, and key $K$ is also used
for semantic segmentation. In the proposed method, a query image is applied to
a model trained with $K$, and the model is protected by permuting feature maps
with key $K$ as well as for training the model. Protected models are expected
to satisfy the following requirements.

  1. 1.

Providing almost the same performance as that of using unprotected models to
authorized users.

  2. 2.

Degrading performance for unauthorized users even when they estimate key $K$.

##  3 Experiments and Results

The effectiveness of the proposed method was evaluated in terms of
segmentation performance and robustness against unauthorized uses.

###  3.1 Experimental Setup

We used a FCN with a ResNet-50 backbone (see Fig. 3) for semantic
segmentation. Segmentation models were trained by using the PASCAL visual
object classes dataset released in 2012 [35] for semantic segmentation, where
the dataset has a training set with 1464 pairs (i.e., images and corresponding
ground truths) and a validation set with 1449 pairs. We split the training set
into training and development sets; 1318 samples were used for training, and
146 samples were used for development when training models. The whole
validation set was used for testing the models. The performance of the trained
models was calculated by averaging the mean IoU (see Eq. (1)) values of
predictions. The conventional method using block-wise transformations [24]
requires a fixed input image size. Therefore, to compare the conventional
method with the proposed one, all input images and ground truths were resized
to $256\times 256$. In addition, standard data-augmentation methods, i.e.,
random resized crop and horizontal flip, were performed in the training.

All networks were trained for 100 epochs by using the stochastic gradient
descent (SGD) optimizer with a weight decay of 0.005 and a momentum of 0.9.
The learning rate was initially set to 0.1 and scheduled by cosine annealing
with warm restarts [36]. The minimum learning rate was 0.0001, and the
learning rate was restarted every ten epochs. The batch size was 256. We used
cross-entropy loss to calculate loss. After the training, we selected the
model that provided the lowest loss value under the validation.

###  3.2 Performance Evaluation under Correct Key $K$

In this experiment, one feature map was chosen from six feature maps in the
network, and segmentation models were then trained by permuting the chosen
feature map with key $K$. The trained models were evaluated for authorized
users with $K$. “Correct ($K$)” in Table 1 shows the result under this
condition, where the model trained by permuting feature map 1 is referred to
as Model-1 as an example, and “Baseline” denotes that the model was trained
and tested without any feature map permutation.

Table 1: Segmentation accuracy (mean IoU) of protected models.
Best accuracies are shown in bold. Selected feature map | Correct ($K$) | No-perm | Incorrect ($K^{\prime}$)
---|---|---|---
1 (Model-1) | 0.5290 | 0.1678 | 0.0356
2 (Model-2) | 0.4975 | 0.3735 | 0.0004
3 (Model-3) | 0.4520 | 0.3799 | 0.0001
4 (Model-4) | 0.5446 | 0.4040 | 0.0013
5 (Model-5) | 0.5759 | 0.5328 | 0.0000
6 (Model-6) | 0.5888 | 0.1595 | 0.0000
Baseline | 0.5966

From Correct ($K$) in Table 1, several models such as Model-6 and Model-5 had
a high segmentation accuracy, which was almost the same as that of the
baseline, although a couple of models had a slightly degraded accuracy. Figure
6 also shows an example of prediction results. From this figure, the proposed
model was demonstrated to maintain prediction results similar to the baseline.

Input | Ground Truth | Baseline | Proposed
---|---|---|---
Correct | Plain | Incorrect
|  |  |  |  |
|  | 0.8532 | 0.8570 | 0.0000 | 0.0000
|  |  |  |  |
|  | 0.4691 | 0.5062 | 0.1145 | 0.0000
|  |  |  |  |
|  | 0.7912 | 0.3679 | 0.0630 | 0.0000
|  |  |  |  |
|  | 0.6231 | 0.9279 | 0.6126 | 0.0000
Figure 6: Example of prediction results. Mean IoU values are given under
predictions.

###  3.3 Robustness against Unauthorized Access

We assume that unauthorized users have no key $K$ and that they know both the
method for protecting models and the permuted feature maps. To evaluate
robustness against unauthorized access, we also evaluated the six protected
models under two key conditions: No-perm and Incorrect, as shown in Table 1.
“No-perm” indicates that protected models were tested without any feature map
permutation. “Incorrect” denotes that protected models were tested by
permuting a feature map used in training with incorrect (randomly generated)
key $K^{\prime}$.

Table 1 shows the results under these conditions, where the results for
Incorrect were averaged over 100 incorrect keys. From the table, the protected
models provided a low segmentation accuracy, so the proposed models were
robust against these attacks. An example of prediction results is also
illustrated in Fig. 6. The robustness of the proposed models can be visually
confirmed from the figure as well.

###  3.4 Comparison with State-of-the-art Method

A method for model protection was proposed in [24]. In the method, input
images are transformed by using three block-wise transformations with a secret
key: pixel shuffling (SHF), negative/positive transformation (NP), and format-
preserving Feistel-based encryption (FFX) [37]. Although this conventional
method can achieve high performance in image classification models, it has
never been applied to other models such as semantic segmentation ones. To be
compared with the proposed method, it was also applied to semantic
segmentation models.

We trained segmentation models with different block sizes and tested the
models with three key conditions, i.e., Correct, Plain, and Incorrect, where
“Plain” used plain images as input ones, and “Incorrect” used images encrypted
by using an incorrect key. As shown in Table 2, the performance of all
transformations heavily decreased compared with the proposed method.

Table 2: Segmentation accuracy (mean IoU) of transformations of conventional method Block size | SHF | NP | FFX
---|---|---|---
Correct | Plain | Incorrect | Correct | Plain | Incorrect | Correct | Plain | Incorrect
2 | 0.5447 | 0.5112 | 0.4800 | 0.5143 | 0.4370 | 0.0529 | 0.4224 | 0.0047 | 0.0000
4 | 0.4731 | 0.4536 | 0.3671 | 0.4706 | 0.3359 | 0.1505 | 0.3823 | 0.0157 | 0.0012
8 | 0.3599 | 0.3691 | 0.2527 | 0.4094 | 0.3838 | 0.0771 | 0.3613 | 0.0001 | 0.0012
16 | 0.2214 | 0.1994 | 0.1150 | 0.3439 | 0.2114 | 0.0832 | 0.2611 | 0.0007 | 0.0079
32 | 0.1062 | 0.0598 | 0.0398 | 0.2928 | 0.0499 | 0.1239 | 0.1995 | 0.0152 | 0.0770
Model-6 (Proposed) | 0.5888 (Correct) | 0.1595 (Plain) | 0.0000 (Incorrect)
Baseline (non-protected) | 0.5966

In particular, when using a large block size, the accuracy for the correct key
was low. In contrast, the accuracy for the plain and incorrect keys was high,
so the performance was confirmed to be poor for protecting semantic
segmentation models.

Semantic segmentation models are required to output visual information as an
image, so transformations applied to images or feature maps for training and
testing models have to be spatially invariant, but the conventional block-wise
transformations are not spatially invariant. That is why the performance of
the block-wise transformations was poor for semantic segmentation models.

##  4 Conclusion

We proposed an access control method that uses the spatially invariant
permutation of feature maps for protecting semantic segmentation models for
the first time. Semantic segmentation models are required to output visual
information as an image, so transformations for model protection have to be
spatially invariant, but conventional transformations are not spatially
invariant. The proposed method allows us not only to obtain a high
segmentation accuracy but also for there to be robustness against various
attacks by unauthorized users. In experiments, the effectiveness of the
proposed method was verified in terms of segmentation performance and
robustness against unauthorized access. In contrast, the conventional
protection method with block-wise transformations that was proposed for image
classification models was demonstrated to not be applicable to segmentation
models. As future work, we plan to evaluate the robustness against more
diverse attacks such as key estimation attacks.

## Acknowledgment

This study was partially supported by JSPS KAKENHI (Grant Number JP21H01327)
and Support Center for Advanced Telecommunications Technology Research,
Foundation (SCAT).

## References

  * [1] A. Krizhevsky, I. Sutskever, and G. E. Hinton, “Imagenet classification with deep convolutional neural networks,” _Commun. ACM_ , vol. 60, no. 6, pp. 84–90, May 2017.
  * [2] X. Liu, Z. Deng, and Y. Yang, “Recent progress in semantic image segmentation,” _Artif. Intell. Rev._ , vol. 52, no. 2, pp. 1089–1106, Aug. 2019.
  * [3] M. Xue, J. Wang, and W. Liu, “DNN intellectual property protection: Taxonomy, attacks and evaluations (invited paper),” in _Proceedings of the 2021 on Great Lakes Symposium on VLSI_ , Jun. 2021, pp. 455–460.
  * [4] Y. Uchida, Y. Nagai, S. Sakazawa, and S. Satoh, “Embedding watermarks into deep neural networks,” in _Proceedings of the 2017 ACM on International Conference on Multimedia Retrieval_ , 2017, pp. 269–277.
  * [5] J. Zhang, Z. Gu, J. Jang, H. Wu, M. P. Stoecklin, H. Huang, and I. Molloy, “Protecting intellectual property of deep neural networks with watermarking,” in _Proceedings of the 2018 on Asia Conference on Computer and Communications Security_ , 2018, pp. 159–172.
  * [6] B. Darvish Rouhani, H. Chen, and F. Koushanfar, “Deepsigns: An end-to-end watermarking framework for ownership protection of deep neural networks,” in _Proceedings of the Twenty-Fourth International Conference on ASPLOS_ , 2019, pp. 485–497.
  * [7] E. Le Merrer, P. Pérez, and G. Trédan, “Adversarial frontier stitching for remote neural network watermarking,” _Neural Computing and Applications_ , vol. 32, no. 13, pp. 9233–9244, Jul. 2020.
  * [8] L. Fan, K. W. Ng, and C. S. Chan, “Rethinking deep neural network ownership verification: Embedding passports to defeat ambiguity attacks,” in _Advances in Neural Information Processing Systems_ , vol. 32, 2019, pp. 4716–4725.
  * [9] S. Sakazawa, E. Myodo, K. Tasaka, and H. Yanagihara, “Visual decoding of hidden watermark in trained deep neural network,” in _2019 IEEE Conference on Multimedia Information Processing and Retrieval_ , 2019, pp. 371–374.
  * [10] H. Chen, B. D. Rouhani, C. Fu, J. Zhao, and F. Koushanfar, “Deepmarks: A secure fingerprinting framework for digital rights management of deep learning models,” in _Proceedings of the 2019 on International Conference on Multimedia Retrieval_ , 2019, pp. 105–113.
  * [11] Y. Adi, C. Baum, M. Cisse, B. Pinkas, and J. Keshet, “Turning your weakness into a strength: Watermarking deep neural networks by backdooring,” in _27th USENIX Security Symposium (USENIX Security 18)_ , Aug. 2018, pp. 1615–1631.
  * [12] M. Fredrikson, S. Jha, and T. Ristenpart, “Model inversion attacks that exploit confidence information and basic countermeasures,” in _Proceedings of the 22nd ACM SIGSAC Conference on Computer and Communications Security_ , 2015, pp. 1322–1333.
  * [13] C. Szegedy, W. Zaremba, I. Sutskever, J. Bruna, D. Erhan, I. J. Goodfellow, and R. Fergus, “Intriguing properties of neural networks,” in _2nd International Conference on Learning Representations_ , 2014.
  * [14] I. J. Goodfellow, J. Shlens, and C. Szegedy, “Explaining and harnessing adversarial examples,” in _3rd International Conference on Learning Representations_ , 2015.
  * [15] N. Papernot, P. McDaniel, I. Goodfellow, S. Jha, Z. B. Celik, and A. Swami, “Practical black-box attacks against machine learning,” in _Proceedings of the 2017 ACM on Asia Conference on Computer and Communications Security_ , 2017, pp. 506–519.
  * [16] M. Tanaka, “Learnable image encryption,” in _2018 IEEE International Conference on Consumer Electronics-Taiwan (ICCE-TW)_ , 2018, pp. 1–2.
  * [17] K. Madono, M. Tanaka, M. Onishi, and T. Ogawa, “Block-wise scrambled image recognition using adaptation network,” in _Artificial Intelligence of Things (AIoT), Workshop on AAAI conference Artificial Intellignece, (AAAI-WS)_ , 2020.
  * [18] W. Sirichotedumrong and H. Kiya, “A gan-based image transformation scheme for privacy-preserving deep neural networks,” in _28th European Signal Processing Conference, EUSIPCO_ , 2020, pp. 745–749.
  * [19] H. Ito, Y. Kinoshita, M. Aprilpyone, and H. Kiya, “Image to perturbation: An image transformation network for generating visually protected images for privacy-preserving deep neural networks,” _IEEE Access_ , vol. 9, pp. 64 629–64 638, 2021.
  * [20] W. Sirichotedumrong, Y. Kinoshita, and H. Kiya, “Pixel-based image encryption without key management for privacy-preserving deep neural networks,” _IEEE Access_ , vol. 7, pp. 177 844–177 855, 2019.
  * [21] W. Sirichotedumrong, T. Maekawa, Y. Kinoshita, and H. Kiya, “Privacy-preserving deep neural networks with pixel-based image encryption considering data augmentation in the encrypted domain,” in _2019 IEEE International Conference on Image Processing_ , 2019, pp. 674–678.
  * [22] M. Chen and M. Wu, “Protect your deep neural networks from piracy,” in _2018 IEEE International Workshop on Information Forensics and Security_ , 2018, pp. 1–7.
  * [23] M. AprilPyone and H. Kiya, “Training dnn model with secret key for model protection,” in _2020 IEEE 9th Global Conference on Consumer Electronics_ , 2020, pp. 818–821.
  * [24] ——, “A protection method of trained CNN model with a secret key from unauthorized access,” _APSIPA Transactions on Signal and Information Processing_ , vol. 10, p. e10, 2021.
  * [25] ——, “Block-wise image transformation with secret key for adversarially robust defense,” _IEEE Transactions on Information Forensics and Security_ , vol. 16, pp. 2709–2723, 2021.
  * [26] T. Chuman, W. Sirichotedumrong, and H. Kiya, “Encryption-then-compression systems using grayscale-based image encryption for jpeg images,” _IEEE Transactions on Information Forensics and Security_ , vol. 14, no. 6, pp. 1515–1525, 2019.
  * [27] W. Sirichotedumrong and H. Kiya, “Grayscale-based block scrambling image encryption using ycbcr color space for encryption-then-compression systems,” _APSIPA Transactions on Signal and Information Processing_ , vol. 8, p. e7, 2019.
  * [28] K. Kurihara, S. Imaizumi, S. Shiota, and H. Kiya, “An encryption-then-compression system for lossless image compression standards,” _IEICE Transactions on Information and Systems_ , vol. E100.D, no. 1, pp. 52–56, 2017.
  * [29] T. Chuman, K. Kurihara, and H. Kiya, “Security evaluation for block scrambling-based etc systems against extended jigsaw puzzle solver attacks,” in _2017 IEEE International Conference on Multimedia and Expo_ , 2017, pp. 229–234.
  * [30] A. P. M. Maung and H. Kiya, “Training DNN model with secret key for model protection,” in _2020 IEEE 9th Global Conference on Consumer Electronics_ , 2020, pp. 818–821.
  * [31] M. Swanson, M. Kobayashi, and A. Tewfik, “Multimedia data-embedding and watermarking technologies,” _Proceedings of the IEEE_ , vol. 86, no. 6, pp. 1064–1087, 1998.
  * [32] H. Rezatofighi, N. Tsoi, J. Gwak, A. Sadeghian, I. Reid, and S. Savarese, “Generalized intersection over union: A metric and a loss for bounding box regression,” in _2019 IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , 2019, pp. 658–666.
  * [33] J. Long, E. Shelhamer, and T. Darrell, “Fully convolutional networks for semantic segmentation,” in _2015 IEEE Conference on Computer Vision and Pattern Recognition_ , 2015, pp. 3431–3440.
  * [34] K. He, X. Zhang, S. Ren, and J. Sun, “Deep residual learning for image recognition,” in _2016 IEEE Conference on Computer Vision and Pattern Recognition_ , 2016, pp. 770–778.
  * [35] M. Everingham, S. M. A. Eslami, L. Van Gool, C. K. I. Williams, J. Winn, and A. Zisserman, “The pascal visual object classes challenge: A retrospective,” _International Journal of Computer Vision_ , vol. 111, no. 1, pp. 98–136, Jan. 2015.
  * [36] I. Loshchilov and F. Hutter, “SGDR: stochastic gradient descent with warm restarts,” in _5th International Conference on Learning Representations_ , 2017.
  * [37] M. Bellare, P. Rogaway, and T. Spies, “Addendum to “the ffx mode of operation for format-preserving encryption” a parameter collection for enciphering strings of arbitrary radix and length,” in _Draft 1.0_ , 2010. [Online]. Available: https://csrc.nist.rip/groups/ST/toolkit/BCM/documents/proposedmodes/ffx/ffx-spec2.pdf