# Weibel-dominated quasi-perpendicular shock: hybrid simulations and in-situ
observations

J.A. Kropotina,1, A.A. Petrukovich2, O.M. Chugunova2, A.M. Bykov1
1Ioffe Institute
2Space Research Institute, Russian Academy of Sciences, Moscow, Russia
E-mail<EMAIL_ADDRESS>
(Accepted 2023 July 4. Received 2023 July 4; in original form 2023 April 27)

###### Abstract

We directly compare hybrid kinetic simulations and in-situ observations of a
high Mach number high-beta shock in the Solar wind. We launch virtual probes
to demonstrate that the model quantitatively reproduces the observations. The
observed wave properties are caused by the ion Weibel instability in the shock
foot. Parameters of reflected ions in the shock foot are extracted from
simulations, and their coordinate dependencies are linearly approximated.
These approximations could be used in analytical models. Due to strong
magnetic variations at ramp the reflected ions density can be locally very
high (nearly that of the incoming flow), which makes favourable conditions for
the instability.

###### keywords:

Physical Data and Processes: instabilities, Physical Data and Processes: shock
waves, Software: simulations, Physical Data and Processes: plasmas

††pubyear: 2023††pagerange: Weibel-dominated quasi-perpendicular shock: hybrid
simulations and in-situ observations–Weibel-dominated quasi-perpendicular
shock: hybrid simulations and in-situ observations

##  1 Introduction

Collisionless shocks propagating through low magnetized plasma appear in
various astrophysical objects. The synchrotron radiation of gamma-ray burst
afterglows is likely associated with energetic particles and magnetic fields,
produced by the Weibel instability mediated shocks driven by relativistic
outflows in the low-magnetized circumburst medium (see e.g. Medvedev & Loeb,
1999; Milosavljević & Nakar, 2006; Lemoine et al., 2019). The spectacular
merger and accretion events in the clusters of galaxies are accompanied by the
observed large scale shocks which are propagating through a hot intercluster
plasma (see e.g. Churazov et al., 2021, 2023; Markevitch & Vikhlinin, 2007;
Bykov et al., 2019). The shock Alfv‘en Mach numbers and the ratio of thermal
to magnetic pressure (i.e. the plasma parameter $\beta$) can be large in the
intercluster medium. Also, the magnetic field in the cold expanding supernova
ejecta is expected to be extremely low, if one assumes that it is the field of
progenitor star scaled according to the magnetic flux conservation law
(Ellison et al., 2005; Telezhinsky et al., 2012). Hence, the reverse shocks
observed in supernova remnants (see Raymond, 2018, for a review) are likely to
be unmagnetized.

Numerical models and laser plasma experiments (Huntington et al., 2015; Park
et al., 2015; Marcowith et al., 2016) indicate that in the environments
initially negligible magnetic fields can substantially grow due to the ion
Weibel instability (IWI) which appears during the interaction of incoming and
reflected flows (Chang et al., 1990; Burgess et al., 2016). The electron
Weibel instability arises first and quickly thermalizes electrons. After that
the much slower IWI instability comes into play and forms a collisionless
shock with a strong electromagnetic turbulence. As a result the magnetic
energy density can reach about 1–2 % of the upstream bulk kinetic energy
density (Kato & Takabe, 2008). These magnetic fields not only shape
collisionless shocks but are also favorable for magnetic reconnection and
subsequent particle acceleration (Bohdan et al., 2020). Thus low-magnetized
Weibel-mediated shocks might be a site of electrons pre-acceleration and
injection into the first order Fermi acceleration. Meanwhile strong magnetic
fields near shock transition increase the momentum which particles need to
enter the Fermi process. So the summary impact of the IWI on the Fermi
acceleration efficiency is still an open question. Understanding the
microstructure and properties of Weibel-mediated shocks is required to solve
this problem.

Weibel-mediated shocks have been extensively studied by means of kinetic
simulations (Kato & Takabe, 2008, 2010; Spitkovsky, 2008; Bohdan et al.,
2020). Particle-in-cell (PIC) codes used for the simulations are highly
resource-intensive because they operate on electron scales and must resolve
the Debye length to avoid a nonphysical heating. As far as collisionless
shocks form on much greater ion scales, some tricks are routinely used to
artificially bring electron scales closer to the ion ones and reduce the
computation time. Those tricks include reducing of the proton-to-electron mass
ratio and increasing the upstream temperature, both increasing the ratio of
the Debye radius to the ion inertial length. The upstream flow velocity (in
the shock front reference frame) $V_{sh}$ must be increased correspondingly to
keep the sonic Mach number. For this reason PIC simulations usually deal with
relativistic or subrelativistic shocks (at least $V_{sh}\sim 0.1c$, where $c$
is the velocity of light).

In the nonrelativisic case the magnetisation parameter can be estimated as
$\sigma\equiv[B_{0}^{2}/8\pi]/[n_{0}(m_{i}+m_{e})V^{2}/2]$, where $B_{0}$,
$n_{0}$ and $V$ are a far upstream magnetic field, number density and flow
velocity; $m_{i}$ and $m_{e}$ are proton and electron masses. It was proposed
in Kato & Takabe (2008) that $\sigma$ must be lower than $10^{-4}$ for a shock
to be Weibel-dominated. This makes doubtful the presence of Weibel-dominated
shocks in older than 1000 yr supernova remnants. Meanwhile in the laser plasma
experiment of Park et al. (2015) Weibel structures were found in shock with
$\sigma\sim 10^{-3}$. In the solar wind $\sigma$ are even higher (at least
$\sim 10^{-2}$). Moreover, in Burgess et al. (2016) Weibel structures were
obtained in hybrid simulations of low Mach number low-beta shocks with sigma
about 0.1. This points out that shocks can be Weibel-dominated in moderately
magnetized regime. However, the authors pointed out that hybrid codes have
some limitations, and thus their results should be confronted with
observations or PIC simulations. They also proposed that the Cluster and MMS
spacecraft are capable of resolving Weibel structures.

Near-Earth observations with spacecraft provide unique possibility to study
collisionless shock structure and dynamics in-situ, measuring electric and
magnetic field, as well as electron and ion distribution functions. However,
low-magnetized (high-$\beta$) conditions are not frequent in solar wind
plasmas. Among a very few early direct observations of the low magnetized
shocks in the solar wind (e.g. Formisano et al., 1975; Winterhalter &
Kivelson, 1988), ISEE 1 and ISEE 2 spacecraft revealed important details of
the high $\beta$ terrestrial bow shock structure (Farris et al., 1992). The
large separation between the spacecraft (over 2500 km) allowed simultaneous
upstream and downstream solar wind plasma measurements for a period of about 8
min. Large amplitude magnetic field and density fluctuations were measured and
a hot dense field-aligned ion beam escaping from the downstream region of the
shock was detected. The beam was associated by Farris et al. (1992) with short
periodic magnetic holes detected in front of the bow shock. Recent studies of
high $\beta$ shocks with MMS, Cluster and Geotail spacecrafts were reviewed by
Petrukovich & Chugunova (2021). Often the transition region of a high-beta
shock contained quasi-periodic linearly polarised pulsations, most probably
associated with the IWI (Sundberg et al., 2017; Petrukovich et al., 2019;
Petrukovich & Chugunova, 2021). But the relation between observed quasi-
periodic oscillations and shock structure was never studied in details.

In this paper we build a hybrid kinetic model of a nonrelativistic high-beta
shock observed by MMS and directly compare it with the observations. In our
model we study the growth of the magnetic variance in the foot region and find
it consistent with the predictions of the kinetic linear theory for the IWI.
Also we launch a virtual probe to study the nature of the observed
quasiperiodic oscillations. We conclude that the observed nonrelativistic
($V\sim 400$ km/s) quasiperpendicular shock is formed due to the IWI and has a
typical structure with normal-aligned filaments of density and magnetic field.
Weibel structures are non-propagating in the plasma reference frame, but they
are convected supersonically along the shock surface. This happens because the
mean flow velocity along shock surface is substantial in the foot region
occupied by reflected ions. Hence waves minima and maxima come across the
slowly moving spacecraft and lead to the observed pulsations.

The quantitative agreement of our hybrid kinetic model and in-situ
observations, as well as qualitative agreement with Burgess et al. (2016)
allows to verify that hybrid codes are capable of reproducing Weibel-dominated
shocks. We also determine the properties of the reflected ions beam
responsible for the development of the instability.

The paper is organized as follows: in Section 2 we provide the linear theory
of the IWI; in Section 3 we describe an observed event in the Solar wind; in
Section 4 we introduce our kinetic numerical model and in Section 5 we discuss
the simulated shock properties and compare them with the observed ones. The
discussion and conclusions are given in Sections 6 and 7 respectively.

##  2 Theory

The transverse electromagnetic Weibel (1959) instability is widely discussed
for a long time in the modeling of collective processes in plasma with
anisotropic particle distributions, both in the laboratory plasma
installations (e.g. Morse & Nielson, 1971; Davidson et al., 1972, 2004) and in
the space environment (see e.g. Balogh & Treumann, 2013; Bykov & Treumann,
2011; Sironi et al., 2015; Marcowith et al., 2016; Pelletier et al., 2017;
Takabe, 2023). The ion beam Weibel instability in a cold unmagnetized cross-
field ion beam moving relative to the static cold magnetized electrons was
considered by Chang et al. (1990). Besides the well-known modified two-stream
and lower-hybrid drift instabilities they found a purely growing
electromagnetic mode which they called the IWI. Their approach was generalized
in Park et al. (2015); Burgess et al. (2016) for the case of two opposite cold
unmagnetized cross-field ion beams. In the center of mass reference frame the
growth rate is given by

$\Gamma^{2}=\frac{k^{2}n_{c}n_{b}(V_{c}-V_{b})^{2}}{(nc+n_{b})^{2}(1+k^{2}c^{2}/\omega_{pi}^{2})},$
(1)

where $V$ and $n$ are a flow velocity and number density of both ion
populations, $k$ is a wavenumber, $\omega_{pi}$ is the ion plasma frequency.
Here the subscript $c$ denotes the denser core and the subscript $b$ — the
fainter beam (note, however, that the expression is symmetric, so the
subscripts can be exchanged). The growth rate is independent on $B_{0}$ and
becomes asymptotic to $|V_{c}-V_{b}|\sqrt{n_{c}n_{b}}/(n_{c}+n_{b})$ for
$k\gg\omega_{pi}/c$. The wavevector is perpendicular to the beams.

Kato & Takabe (2010) studied the IWI kinetically taking the parameters from
their PIC simulation. In case when a magnetic field and a wavevector are along
$z$, and both beams lie in the $x-y$ plane the dispersion equation reads as

$\det\Lambda=0,$ (2)

where

$\displaystyle\Lambda_{xx}=1-\left(\frac{kc}{\omega}\right)^{2}+\frac{1}{2}\left(\frac{\omega_{pe}}{\omega}\right)^{2}\xi_{0}[Z(\xi_{1})+Z(\xi_{-1})]+$
$\displaystyle+\sum_{s}\left(\alpha_{s}+2\left(\frac{V_{x,s}}{V_{T,s}}\right)^{2}(1+\alpha_{s})\right)\left(\frac{\omega_{ps}}{\omega}\right)^{2},$
(3)
$\displaystyle\Lambda_{yy}=1-\left(\frac{kc}{\omega}\right)^{2}+\frac{1}{2}\left(\frac{\omega_{pe}}{\omega}\right)^{2}\xi_{0}[Z(\xi_{1})+Z(\xi_{-1})]+$
$\displaystyle+\sum_{s}\left(\alpha_{s}+2\left(\frac{V_{y,s}}{V_{T,s}}\right)^{2}(1+\alpha_{s})\right)\left(\frac{\omega_{ps}}{\omega}\right)^{2},$
(4)
$\displaystyle\Lambda_{zz}=1+2\left(\frac{\omega_{pe}}{kV_{T,e}}\right)^{2}[1+\xi_{0}Z(\xi_{0})]+$
$\displaystyle+2\sum_{s}\left(\frac{\omega_{ps}}{kV_{T,s}}\right)^{2}(1+\alpha_{s}),$
(5)
$\displaystyle\Lambda_{xy}=\frac{i}{2}\left(\frac{\omega_{pe}}{\omega}\right)^{2}\xi_{0}[Z(\xi_{1})-Z(\xi_{-1})]+$
$\displaystyle+2\sum_{s}\left(\frac{\omega_{ps}}{\omega}\right)^{2}\frac{V_{x,s}}{V_{T,s}}\frac{V_{y,s}}{V_{T,s}}(1+\alpha_{s}),$
(6)
$\displaystyle\Lambda_{yx}=-\frac{i}{2}\left(\frac{\omega_{pe}}{\omega}\right)^{2}\xi_{0}[Z(\xi_{1})-Z(\xi_{-1})]+$
$\displaystyle+2\sum_{s}\left(\frac{\omega_{ps}}{\omega}\right)^{2}\frac{V_{x,s}}{V_{T,s}}\frac{V_{y,s}}{V_{T,s}}(1+\alpha_{s}),$
(7)
$\displaystyle\Lambda_{xz}=\Lambda_{zx}=2\sum_{s}\left(\frac{\omega_{ps}}{\omega}\right)^{2}\frac{V_{x,s}}{V_{T,s}}\frac{\omega}{kV_{T,s}}(1+\alpha_{s}),$
(8)
$\displaystyle\Lambda_{yz}=\Lambda_{zy}=2\sum_{s}\left(\frac{\omega_{ps}}{\omega}\right)^{2}\frac{V_{y,s}}{V_{T,s}}\frac{\omega}{kV_{T,s}}(1+\alpha_{s}),$
(9)
$\displaystyle\xi_{n}=\frac{\omega-n\Omega_{e}}{kV_{T,e}},\quad\alpha_{s}=\left(\frac{\omega}{kV_{T,s}}\right)Z\left(\frac{\omega}{kV_{T,s}}\right),$
(10) $\displaystyle\Omega_{e}=-\frac{eB}{mc},\quad\omega_{ps}=\sqrt{\frac{4\pi
n_{s}e_{s}^{2}}{m_{s}}},\quad V_{T,s}=\sqrt{\frac{2k_{B}T_{s}}{m_{s}}},$ (11)
$\displaystyle
Z(\xi)\equiv\pi^{-\frac{1}{2}}\int_{-\infty}^{\infty}\frac{e^{-z^{2}}}{z-\xi}dz.$
(12)

Here $\omega$ is a complex frequency, and the summation is only over ion
sorts, i. e. for $s=b,c$ (beam, core).

This dispersion equation will be solved numerically in section 5.5 with
parameters taken from our simulations. The solution includes a purely growing
mode which corresponds to the kinetic IWI. The increment is typically much
lower than (1).

##  3 Observations

For the analysis we used measurements of NASA Magnetospheric Multiscale (MMS)
project from magnetic field (FGM) Russell et al. (2016) and plasma (FPI)
Pollock et al. (2016) experiments. In order to directly compare simulations
and observations we chose the bow shock crossing by MMS spacecraft on November
25, 2017 (see Petrukovich & Chugunova (2021)). This is a high-beta strong
collisionless shock in a region with an ambient magnetic field $B_{0}$ as low
as $0.9$ nT. The Alfvén Mach number in the shock rest frame is $M_{a}\approx
60$ and the shock inclination angle is $\theta\approx 65^{\circ}$. with a
total ion number density $n_{i}\approx 9$ cm-3. The protons’ temperature is
$T_{p}\approx 1.1$ eV, and the electrons’ temperature is $T_{e}\approx 13.4$
eV.

Shocks with such parameters are rare in solar wind, only about 30 well
documented cases for $\beta>30$ were found in the observations by modern
spacecraft (Petrukovich & Chugunova, 2021). About half of these cases have
rather rich internal structure with extended variations, similar to the event
presented here, while the other have the appearance closer to a more standard
shock (a single magnetic field and density jump). It should be noted that this
difference in appearance is not due to the angle between shock normal and
upstream magnetic field (parallel shocks are known to have more extended
variations than perpendicular ones). Most of considered shocks have this angle
larger than 45∘.

Though the shock crossing takes some minutes (Fig. 1), the physical width of
the shock layer is only about one proton cyclotron radius in the (very low)
upstream magnetic field. The spacecraft gradually crosses the shock from
downstream to upstream and observes relatively stable picture of periodically
($\sim$15 s) emerging activations, gradually thermalizing the solar wind ion
flow. Each activation, in turn, consists of high-amplitude magnetic variations
with a period about 1 s, coupled with pulses of a downstream-like plasma flow.
Sometimes these density and magnetic field peaks are higher than the
downstream averaged plasma density and field values. Between the pulses more
upstream-like flow is observed with a substantial fraction of reflected and
accelerated ions.

Available observations with four closely separated spacecraft allow to
determine the wavelength of 1-sec pulsations of about 150 km as well as the
propagation velocity and direction. These waves have linear polarisation and
are almost standing in the plasma rest frame, consistent with the expectation
for the Weibel mode. Later on we compare these values with those obtained in
simulations.

Figure 1: MMS observations of the bow shock crossing at November 25, 2017. (a)
Ion omnidirectional spectrogram; (b) ion number density; (c) magnetic field
magnitude; (d) $B_{y}$ magnetic field; (e) wavelet dynamic spectrum of
magnetic field $B_{y}$.

##  4 Simulations

We modeled a shock with parameters taken from the observations by means of the
hybrid code “Maximus” (Kropotina et al., 2019, 2021). We used a 3d cartesian
grid sized $L_{x}\times L_{y}\times L_{z}=2500\times 150\times 150$ cells,
each cell $0.1l_{i}^{3}$, where $l_{i}$ is the proton inertial length. The
shock was launched via the rigid piston method, when a super-Alfvénic flow
with a bulk velocity $V_{x}=-45V_{a}$ was reflected from a conductive wall at
$x=0$. This resulted in a formation of a shock front moving in the positive
$x$ direction with $V_{f}\approx 15.3V_{a}$. Thus in the shock front frame
$M_{a}=60.3$. The initial magnetic field lay in the $x-z$ plane at an angle
$\theta=65^{\circ}$ to the shock normal. Average values for Helium content (4
% He(+2) by number) and temperature (He(+2) is four times hotter than protons)
were used in the model. Thus the total mass density was $\rho_{0}\approx
1.7\cdot 10^{-23}$ g/cm3. Plasma parameters for all particle sorts were
$\beta_{p}\equiv 8\pi 0.96n_{i}T_{p}/B_{0}^{2}=4.8$, $\beta_{He}=8\pi
0.04n_{i}4T_{p}/B_{0}^{2}=0.8$, $\beta_{e}=8\pi
1.04n_{i}T_{e}/B_{0}^{2}=62.4$, and total
$\beta=\beta_{p}+\beta_{He}+\beta_{e}=68$. Electrons were treated as
neutralizing massless fluid with the adiabatic equation of state and standard
adiabatic index $\Gamma=5/3$.

It should be noted that hybrid codes cannot capture electron kinetics, thus
the highest-frequency modes might be modeled incorrectly. However, our model
is highly resource-intensive even within the hybrid approach. Meanwhile the
same simulation box size seems to be unreachable in frames of full PIC
modeling, especially with realistic electron-to-proton mass ratios (which in
turn might affect the results). For this reason we chose the hybrid approach.
The comparison with observations will validate our method at least in the
sense of reproducing wave directions, amplitudes, spectra, and polarisation.

In the hybrid code all quantities are normalized, i.e. a magnetic field and a
mass density are measured in $B_{0}$ and $\rho_{0}$, lengths are measured in
$l_{i}=c/\omega_{pi}\approx 2.3\cdot 10^{7}n_{0}^{-0.5}$ cm, times — in the
inverse proton gyrofrequencies $\Omega^{-1}=m_{p}c/eB_{0}\approx 11.6$ s,
velocities — in the Alfvén velocities $V_{a}=B_{0}/\sqrt{4\pi\rho_{0}}\approx
5.8$ km/s, temperatures are given in energy units $m_{p}V_{a}^{2}\approx 0.4$
eV. To make units self-consistent we took as $n_{0}$ the number density of a
pure proton plasma with the same $\rho_{0}$ (see Matthews (1994)). Thus for
the proton-helium plasma $n_{0}$ was slightly greater than the electron number
density $n_{e}$. All simulation parameters are listed in Table 1.

$B_{0}$ | 0.9 nT | $T_{p}$ | 1.1 eV
---|---|---|---
$n_{i}$ | 9 cm-3 | $T_{He}$ | 4.4 eV
$n_{0}$ | 10 cm-3 | $T_{e}$ | 13.4 eV
$l_{i}$ | 68 km | $L_{x}$ | 2500 cells
$\Omega_{ci}^{-1}$ | 11.6 s | $L_{y}$ | 150 cells
$V_{a}$ | 5.8 km / s | $L_{z}$ | 150 cells
$M_{a}$ | 60.3 | cell size | 0.1 $l_{i}^{3}$
$\theta$ | $65^{\circ}$ |  |
Table 1: Simulation parameters

##  5 Results

###  5.1 Shock structure

Figure 2: Shock dynamics. From top to bottom: cross-section averaged
transverse magnetic field $B_{y}$, magnetic field magnitude and density. An
approximate front position is marked by dashed gray lines. Black arrows
represent a trajectory of a virtual probe. The insets represent the closer
view of the front.

The shock temporal evolution is color-coded in Fig. 2. The front is formed at
$t\approx 30$ s ($\sim 3\Omega_{ci}^{-1}$) and propagates with a nearly
uniform velocity $V_{f}\approx 15.3V_{a}\sim 89$ km/s (its trajectory is shown
by dashed gray lines in all panels). To mimic the MMS shock crossing four
virtual probes were launched. They were located in vertexes of a right
tetrahedron with an edge equal to 0.3 $l_{i}\approx 22$ km. These virtual
spacecraft moved along the shock normal from the downstream to the upstream
with $V_{p,x}=16.3V_{a}$ ($\approx 1V_{a}$ in the front rest frame) and
measured magnetic field and plasma parameters on their way. Their trajectories
are shown in Fig. 2 by a black arrow (the distance between the probes is not
resolved). The resulting temporal profiles are discussed in section 5.2. We
checked that transverse probe movement with $V_{p,y}=V_{p,z}=0.3V_{a}$ didn’t
introduce any differences.

Figure 3: Shock structure. Top panel: $x-V_{x}$ and $x-V_{y}$ protons phase
spaces averaged over $y$ and $z$. Other panels (from top to bottom): velocity,
magnetic field and number density maps in two projections. $V_{x}$ is given in
the front rest frame. Colored rectangles mark downstream and upstream zones
where spectral analyses were made. Magnetic field lines are superimposed in
black. Figure 4: Flow velocity, magnetic field and density along the shock
normal $y=z=0$ at $t=10\Omega_{ci}^{-1}$. Selected regions are the same as in
Fig. 3 Figure 5: Overview of the Probe 1 shock crossing. From top to bottom:
protons phase space $x-E$; ion number density; magnetic fields; wavelet
dynamic spectrum of $B_{y}$.

The shock structure at $t=10\Omega^{-1}\approx 2$ min after initialisation is
shown in Figures 3 and 4. The shock front has a complicated filamentary
structure, and its surface is highly corrugated (rippled). Ion phase spaces in
the top panels of Fig. 3 show that the shock foot extends over several
thousands kilometers. Specularly reflected ions gyrate in the upstream
magnetic field and induce a cross-field current along $x$ and $y$. This is
auspicious for the IWI. A detailed inspection of $V_{x}$ and $n$ maps revealed
that high density filaments correspond to higher negative velocities along
normal. This indicates that density enhancements in the foot are not due
spatial variations of the number of shock-reflected ions. On the contrary
density variations appear further in the foot and are convected towards the
shock front due to the plasma bulk flow.

In the shock coplanarity plane (Fig. 3, right column) thin structures
(filaments) are seen in the foot region. They make a small angle with a shock
normal and grow rapidly towards the front. The overall picture resemble those
seen in simulations of Kato & Takabe (2008, 2010), who concluded that such
structures appear due to the IWI. Structures seen in both planes are also very
alike those discussed in Burgess et al. (2016). They identify narrow Weibel
filaments with width close to $2l_{i}$ in the coplanarity plane and somewhat
wider oblique ’tongues’ in the perpendicular plane which they call ’the AIC-
like ripples’.

Magnetic fluctuations associated with these filaments are very strong near the
front. In the downstream (yellow rectangle) and in the close foot region (cyan
rectangle) their amplitude reaches $10B_{0}$, and between them, at ramp, it is
twice as large. Thus the magnetic energy density reaches about 2 % of the bulk
kinetic energy. Number density peaks reach about $10\rho_{0}$.

###  5.2 Virtual probes

Instantaneous shock profiles along normal (Fig. 4) are much smoother than
those seen by MMS (Fig. 1). Hence we cannot directly compare observations with
simulations if we ignore relative motion of the shock and spacecraft. In Fig.
5 we show what observes one of the virtual probes, which starts at
$t=3.5\Omega^{-1}$ at $x=37.5l_{i}$. The picture is very alike the
observations (see Fig. 1) and severely differs from the instantaneous profiles
(Fig. 4). The reason is that the plasma moves across the probes much faster
than the probes move across the shock. Moreover, the plasma transverse
movements are highly oscillatory. We found that along probes’ ways $V_{y}$
varies between $-30V_{a}$ and $15V_{a}$, and $V_{z}$ — between $-15V_{a}$ and
$5V_{a}$.

To study how the observed oscillations change while the probe is moving across
shock we performed the Morlet wavelet transform of $B_{y}$ projection and
found that the picture qualitatively resembles MMS observations (Cf. the
bottoms panels of Fig. 1 and Fig. 5). The main difference between the model
and observations is the absence in the model of the clear-cut bunching of 1-Hz
oscillations in the $\sim$15 sec ’packets’, though some enhancements are
observed in the frequency spectra at about 0.1 Hz (bottom panel of Fig. 5).

It should be noted that our virtual spacecraft crossed the shock much faster
than the real ones (two minutes vs five). It was done in order to show the
whole transition with reasonable computational efforts. However to check the
impact of probes velocity we also launched several probes starting ar
different points and moving with $V_{p,x}=15.5V_{a}$ and $V_{p,x}=15.7V_{a}$
(i. e. $0.2V_{a}$ and $0.4V_{a}$ relatively to the front).

Fig. 6 shows the results of three probes starting at $x=38l_{i}$, $41l_{i}$
and $42.5l_{i}$ and moving with $V_{p,x}=15.7V_{a}$. They should cross the
shock in $\sim 5$ minutes, just as the real spacecraft. Black lines show the
probes measurements, and red ones correspond to the smoothed data.

The modeled data seems to be more noisy than due to limited number of
particles per cell. Also the highest frequencies might be affected by the grid
resolution and the lack of electrons kinetics. Hence the modeled curves do not
ideally reproduce the observed ones. However rather prominent wave packets
appear in the case of slower probes, which stay longer in each region. More
upstream-like regions with lower density and magnetic fields alternate with
more downstream-like ones. The wave packets are less clearly separated from
each other than in observations. This probably indicates that the observed
shock is more variable. One of the reasons might be that the longest waves are
restricted by simulation box sizes.

Overall, it is possible that prolonged shock crossings like the presented one
are observed due to extremely low proper shock speed (of the order of km/s).
The detailed discussion of this issue is beyond the scope of this paper. Here
we concentrate on comparison of the observed and simulated structures.

Figure 6: Measurements of probes moving with $0.4V_{a}$ relatively to the
shock front. Left column – a probe starting at $x=38l_{i}$; middle column – at
$x=41l_{i}$; right column – at $x=42.5l_{i}$. Red lines show the same data
smoothed by a Savitzky-Golay filter.

We compare the plasma wave properties observed at the near-Earth shock, with
those in the recordings of the virtual probes with $V_{p,x}=16.3V_{a}$ (Table
2). Four probes in observations and simulations allow to determine not only
the temporal sequence of measured parameters, but also determine spatial
gradient (hence, wavevector) on a scale of separation. In both cases magnetic
oscillations were linearly polarised (the maximum variance eigenvalue is at
least 4-5 times larger, than the medium variance one) and the wavevector was
close to the local magnetic field ($\theta_{kB}$ < 40o). Dominating frequency
(measured as frequency of the peak in the spectrum) in observations is larger:
1.25 Hz vs 0.38 and 0.8 Hz (two equivalent peaks are present in the simulation
interval).

The doppler shift in frequency is in both cases close to the observed
frequency, which means that the waves are standing in the rest frame and are
purely convected with the plasma flow. Observed and simulated wavelengths are
about 100-200 km. Concluding, we consider the properties of magnetic
oscillations in observation and simulation as very similar in their principal
characteristics.

Table 2: Wave analysis data for P#2 and P#3. Parameter | observation | simulation
---|---|---
| 23:40:15–23:40:21.5 UT | 97.1–102.7 s
observed frequency, Hz | 1.25 | 0.37, 0.8
eigenvalues | 4.72, 5.78, 46.2 | 11.9 22.47 63.99
observed wave speed, km/s | 173 | 85
wavelength, km | 138 | 230, 106
Doppler shift, Hz | 1.3 | 0.37, 0.81
$\theta_{kB}$ | 35o | 37o

There are oscillations with frequencies between $10^{-1}$ and 1 Hz with
quasiperiodic enhancements and frequency growth towards upstream.

###  5.3 Shock dynamics

Supercritical collisionless shocks are known to be quasi-stationary. The
transition appears due to partial reflection of the incoming ions and reforms
quasi-periodically. Shock reformation is a topic of great interest, actively
investigated by means of numerical models and in-situ observations (see, e. g.
Turner et al., 2021; Yang et al., 2020; Johlander et al., 2022). There are two
mechanisms of this process: (I) the accumulation of reflected ions in a foot
until their density becomes comparable to that at ramp, and (II) the front
interaction with waves convected by the upstream flow (Marcowith et al.,
2016). The insets of Fig. 2 show some signatures of the first type
reformation: the front velocity and the cross-section averaged magnetic field
at the overshoot slightly vary with time. However, field variations are
relatively weak, and a density profile is nearly stationary. So the
“classical” picture of shock reformation did not reveal in this case. The more
thorough investigation of this problem is beyond the scope of this paper.

Meanwhile waves generated by the IWI and convected by the flow substantially
contribute to the shock nonstationarity as well. To demonstrate this we made a
real-time movie of the probe with $V_{p,x}=15.7V_{a}$ recordings together with
its movement through the shock. The movie is available online in the
supplementary materials. Fig. 7 shows one frame of this video. In the upper
row color maps of $B_{y}$, $n$ and $V_{x}$ are given. The velocity is in the
front rest frame. The probe position is marked by a white triangle, and its
recordings are shown in the bottom row. The red line corresponds to the data
measured until the current moment, and the blue one – to the future
recordings.

From the color maps we can see that the front is highly corrugated. In the
movie all these structures move both towards the shock and across it, leading
to a lively structure and oscillatory probe measurements. We chose the moment
when the probe is near a density peak, which corresponds to a higher negative
$V_{x}$. Such regions appear where transverse magnetic fields are low and the
upstream plasma easily penetrate downstream. These “paths of least resistance”
are surrounded by regions with higher $B_{\perp}$, where hotter reflected ions
lead to pressure increase. This pressure compresses colder “paths of least
resistance” up to nearly downstream density. So thin dense filaments appear,
clearly visible in the upper middle panel of Fig. 7.

Figure 7: A snapshot of the shock transition with a virtual probe marked by a
white triangle. Upper row: color maps of $B_{y}$, $n$ and $V_{x}$, lower row:
the probe recordings of the same quantities. The current moment is at the
conjunction of red (past) and blue (future) lines. The corresponding movie is
available online.

###  5.4 Spectral analyses

Figure 8: Panel (d): $k_{x}-k_{z}$ maps of spectral power in $B_{x}$
fluctuations in the downstream region (inside the yellow rectangle in Fig. 3);
panel (e): the same in the close foot region (cyan rectangle in Fig. 3);
panels (a) and (b): 1d spectral power of $B_{x}(x)$ at $z=0$ in the
corresponding regions; panels (c) and (e): 1d spectral power of $B_{x}(z)$ at
the left edges of the corresponding regions. Figure 9: The same as in Fig. 8,
but for $k_{x}-k_{y}$

To better understand the structure of the shock we made spectral analyses of
magnetic field fluctuations in regions marked by yellow and cyan rectangles in
Figures 3 and 4. Figures 8 and 9 show $B_{x}$ spectral power density in
$k_{x}-k_{z}$ and $k_{x}-k_{y}$ planes in the close downstream region marked
by a yellow rectangle, and in the close upstream region marked by a cyan
rectangle. The spectral power density is $|\tilde{B}_{x}|^{2}$, where
$\tilde{B}_{x}$ is a discrete Fourier transform of $B_{x}$. From 1d and 2d
spectra it can be seen that in the close upstream region $k_{z}\approx
3.5l_{i}^{-1}$ and $k_{x}\approx k_{y}\approx 1l_{i}^{-1}$, which corresponds
to $\lambda\sim 2l_{i}\sim 140$ km. In the close downstream wavelengths are
larger. We also checked that in the close upstream $B_{x}$ fluctuations are
stronger than those of $B_{y,z}$, while in the close downstream the spectral
power in all three magnetic field projections is comparable.

The upstream wavevector direction and waves polarisation in the coplanarity
plane $x-z$ agree with those expected for the IWI, i. e. k is along the mean
magnetic field and, hence, is nearly perpendicular to the cross-field beam in
the $x-y$ plane (see top panels of Fig. 3).

###  5.5 Growth rate analyses

To demonstrate that the shock transition is governed by the IWI we directly
compare the growth rate and spectral properties with linear predictions. We
expect a zero frequency mode, thus if we move towards the shock front with the
flow, the waves amplitude grows as

$b_{0}\exp\left(\int\Gamma(x(t))dt\right),$ (13)

where $b_{0}$ is its initial value and $\Gamma$ is an increment.

To check this we studied the evolution of magnetic field fluctuations
amplitude towards shock front. We calculated it as a standard deviation over a
transverse slice embedded in the upstream flow. The result is shown in Fig.
10, a.

We integrated the $x-V_{x}$ and $x-V_{y}$ phase spaces (see the top panels of
Fig. 3) to estimate velocity, density and temperature of the core and beam. We
arbitrarily placed a boundary between them at $V_{x}=-15V_{a}\approx-90$ km/s
and $V_{y}=-5V_{a}\approx-30$ km/s. We also checked that the result is nearly
the same for boundaries at $V_{x}=0$ and $V_{y}=-10V_{a}$. Core and beam
densities, as well as their thermal and flow velocities in the center of mass
rest frame are shown in the panels (b) and (c) of Fig. 10.

Knowing physical parameters at each point we could find hydrodynamic and
kinetic increments and maximal wavenumbers from (1) and (2) respectively. It
should be noted that the kinetic approach suggests treating protons and helium
ions separately. However, we considered unmagnetized ions, so only their
plasma frequencies and thermal velocities are important. The latter are equal
in our model because ions temperatures are mass-proportional. The ratio of
He(+2) and proton plasma frequencies depends only on their number densities,
as if they were both protons. So we solved (2) considering a pure proton
plasma.

Knowing $\Gamma$ and taking initial $b_{0}$ from simulations we directly
compared simulated and theoretical growth using (13). Theoretical amplitudes
are shown in the panel (a) of Fig. 10 by orange and blue curves, and the
predicted wavenumber is superposed on the actual spectrum in the panel (e).
The curve color codes the corresponding increment. We also checked that the
real frequency found from the dispersion equation was zero.

The linear kinetic theory can also predict waves polarisation. In respect that
$\Lambda_{ij}E_{j}=0$, and
$c\mathbf{k}\times\mathbf{E}=\omega\delta\mathbf{B}$, where $\mathbf{E}$ is an
electric field and $\delta\mathbf{B}$ is a magnetic field variation, and
$\mathbf{k}$ is along $z$, we find that

$(\Lambda_{yx}\Lambda_{zz}-\Lambda_{yz}\Lambda_{zx})\delta
B_{y}=(\Lambda_{zy}\Lambda_{yz}-\Lambda_{yy}\Lambda_{zz})\delta B_{x}.$

In panel (d) of Fig. 10 we compare
$|\Lambda_{yx}\Lambda_{zz}-\Lambda_{yz}\Lambda_{zx}|^{2}\langle\delta
B_{y}^{2}\rangle$ and
$|\Lambda_{zy}\Lambda_{yz}-\Lambda_{yy}\Lambda_{zz}|^{2}\langle\delta
B_{x}^{2}\rangle$, where the magnetic variance is taken from simulations, and
$\Lambda_{ij}$ is calculated from the beam properties. The curves do not
perfectly coincide, but they resemble each other even in a highly nonlinear
regime.

It can be seen that the hydrodynamic increment is far too large, but the
simulated growth rate is reasonably explained by the kinetic linear theory
until the wave amplitude approaches about $0.1B_{0}$. After that the system
gradually enters a nonlinear regime, and the predicted growth rate outplays
the actual one. Note also that (2) was obtained for a uniform medium, and the
actual increment may differ due to strong gradients. The predicted wavenumbers
are slightly higher than the simulated spectral maxima, but the simulated
spectrum is rather broad, so the agreement is satisfactory. The polarisation
properties of the IWI are also well reproduced. So we can conclude, that the
IWI governs the shock transition.

Figure 10: Growth rate analyses: (a) the simulated growth of magnetic
fluctuations amplitude towards shock front compared with the predictions of
hydrodynamic and kinetic linear theories; (b) beam and core number densities;
(c) beam and core thermal velocities and flow velocities along $y$ in the
center of mass rest frame; (d) polarisation analyses (see text); (e) color-
coded $B_{x}$ fourier spectrum with overlaid linear prediction of $k_{max}$
(the curve color codes the maximal predicted growth rate).

It should be noted that the investigated shock has a sound Mach number $M_{s}$
as low as 7. However, linear analyses in Nishigai & Amano (2021) indicated
that a shock must have both $M_{a}$ and $M_{s}$ as great as $\sim 20-40$ to be
Weibel-dominated. The authors argued that the instability behaves Weibel-like
if the growth rate is much greater than the ion cyclotron frequency. For an
Alfven Mach number of 60 and a sound Mach number of 7, Fig. 3 of Nishigai &
Amano (2021) predicts the growth rate which is comparable and slightly larger
than the ion cyclotron frequency. In our simulations it locally reaches $\sim
10$.

In Nishigai & Amano (2021) the reflected ions in the foot were parameterised
as a ring distribution with a number density about $0.2n_{0}$, the radius of
the ring equal to the upstream flow velocity, and thermal spread equal to the
one of the upstream flow. From Fig. 10 it can be seen that the actual
quantities strongly vary along the shock normal. Closer to the shock the
reflected ions density nearly reaches that of the incoming flow. In the
regions where the density is about $0.2n_{0}$ the beam velocity is greater
than the upstream one. That is why the IWI growth rate exceeds that predicted
by Nishigai & Amano (2021).

Simple analytical models like that in Nishigai & Amano (2021) are a powerfool
tool to scan a wide range of parameters with minimal computational efforts. So
it is useful to precise them with parameterisation of ions distributions based
on numerical models. To make the first step in this direction we approximated
the reflected ion density and flow velocity in the simulated shock foot by
simple linear functions of coordinate:

$n_{b}/n_{0}=2.0(1.0+(x-x_{sh})/R_{g}),$
$|V_{b,y}-V_{c,y}|/V_{a}=0.66M_{a}(1.0+(x-x_{sh})/R_{g}),$

where $R_{g}$ is the effective particle gyroradius in the foot. It appeared to
be equal to $0.4M_{a}l_{i}$ in our case because the transverse magnetic field
is greater than that far upstream. A thermal velocity of the beam varied only
slightly and was close to $25V_{a}$ (an order higher than that of the core).
The beam velocity along $x$ could not be approximated linearly, but it quickly
reached a relatively stable value $|V_{b,x}-V_{c,x}|/V_{a}\approx M_{a}$. The
difference between $V_{b,x}$ and $V_{b,y}$ is due to the shock drift
acceleration by a motional electric field along $y$ axis Sagdeev (1966). As a
result $V_{b,y}$ eventually exceeds $V_{b,x}$, and $B_{x}$ variation becomes
stronger than that of $B_{y}$.

##  6 discussion

Near Earth spacecraft plasma observations afford unique possibility to sample
in-situ such important astrophysical phenomena as collisionless shocks. A rich
variety of shock structures was discovered, depending on basic plasma
constants and geometry (Mach number, plasma $\beta$, magnetic field direction,
etc). However these experiments are essentially limited by number of
spacecraft simultaneously available — measurements can be performed only in
few points, while the spatial structure at large is only inferred. On the
other hand, numerical modelling affords the full access to spatio-temporal
structure of the shock transition. To make calculations to be completed in
realistic time, simplifications to physical models are usually introduced,
which may question the applicability of results.

In our work, we were able to prove that rather typically observed high
$\beta$, high-Ma shock structure with the developed high-amplitude magnetic
fluctuations is well reproduced with our hybrid model with Helium and
adiabatic electrons. Consistency is found in the general appearance of the
shock transition (Fig.1 and Fig.5), as well as in quantitative characteristics
of the dominating plasma wave mode (Tab. 2).

It is shown that the temporal profiles of the shock crossing depend
substantially on the relative velocity of the shock front and the spacecraft
probes. Slowly flying probes (as mostly in space experiment) are able to
detect the strong temporal variability of the shock front, while high-speed
motion results in rather simple almost instantaneous profile cuts (Fig. 4). It
is not always possible to determine the spacecraft-shock relative velocity in
orbit and possibility of such strong dependence of observations on relative
motion should be taken into account.

Of course, observed differences of the shock structure (see discussions in
Petrukovich et al. (2019); Petrukovich & Chugunova (2021) might be due to some
differences in shock parameters such as magnetic field angle and Mach number.
Some our simulation runs, not shown here, reveal significant variance of shock
structure across parameter range and model details, even if all cases are
high-$\beta$ shocks. This parametric dependence of shock structure is left for
the future studies.

Yet another advantage offered by simulations, is the ability to access the 3D
spatial structure of the transition region in full details. The cuts of
simulation box (like Fig. 3) reveal the complicated breathing filament
structure with varying scale in different directions. These filaments move
rapidly along the shock front and create the magnetic and plasma variability
observed by the probes. The amplitude of these variability is very large,
magnetic amplitudes are order of magnitude larger than the background magnetic
field. Such variability might provide sites of magnetic reconnection and
particle acceleration, though in our case we have not seen it neither in the
simulations nor in the observations, probably due to high $\beta$ value.

Plasma properties also principally change across filaments: more sheath-type
thermalized and high-density streams interchange with more upstream-type with
low density and high percentage of the reflected ions. Close to the ramp the
reflected ions density nearly reaches that of the incoming flow. The detailed
physics of such complicated shock transition remains to be studied with point-
by-point comparison of observations and simulations. It is important for such
a study, as it was stated above, that our numerical model is closely
compatible with observation in all comparable properties.

Finally, our results represent one more proof that high-$\beta$ shock
transition is dominated with the Weibel-like plasma wave mode. We determine
polarisation as well as dispersion characteristics, which coincide in
observations and modeling. The complicated spatial structure detected,
suggests that such mode needs to be considered in the deeply non-linear
regime, practically shaping the process of plasma flow thermalisation. Knowing
parameters of magnetic variations allows to analyse variants of shock-related
particle acceleration and diffusion at such astrophysical objects. To improve
analytical models of shock transition we extracted the parameters of the
reflected ions distribution from our simulation and found that the beam
density and flow velocity could be well approximated by linear functions.

##  7 Conclusions

We demonstrated that hybrid kinetic models can quantitatively reproduce
observed properties of strong Weibel-dominated high-beta quasiperpendicular
shocks. Hybrid models are much less resource-intensive than PIC ones and do
not need high upstream temperatures and subrelativistic flow velocities. This
gives a possibility to study a large field of shock parameters and find the
conditions when shocks become Weibel-dominated. Strong magnetic variations at
ramps of such shocks could prevent particles injection into the first order
Fermi acceleration process. On the other hand such variations might cause
magnetic reconnection which in turn produces nonthermal particles. So the net
impact of the IWI on particle acceleration is still to be determined.

We also extracted from simulations distributions of reflected ions in the
shock foot. This allows to improve existing analytical models of such shocks.

## Acknowledgements

JK and AB acknowledge the Russian Science Fund grant 21-72-20020, which
supported the plasma numerical modeling presented here. Some of the modeling
was performed at the Joint Supercomputer Center JSCC RAS and at the “Tornado”
subsystem of the St. Petersburg Polytechnic University supercomputing center.
AP and OC acknowledge the Russian Science Fund grant 19-12-00313, which
supported the observation analysis and comparison with simulations. Authors
are grateful to NASA MMS project team for excellent space project and
observations. We are very grateful to the reviewer Dr. Takanobu Amano, whose
fruitful suggestions greatly improved this paper.

## Data Availability

MMS spacecraft data are open at the NASA CDAWeb data archive
https://cdaweb.gsfc.nasa.gov/.

## References

  * Balogh & Treumann (2013) Balogh A., Treumann R. A., 2013, Physics of Collisionless Shocks: Space Plasma Shock Waves. ISSI Scientific Report Series, Springer, New York, doi:10.1007/978-1-4614-6099-2, https://cds.cern.ch/record/1517899
  * Bohdan et al. (2020) Bohdan A., Pohl M., Niemiec J., Vafin S., Matsumoto Y., Amano T., Hoshino M., 2020, ApJ, 893, 6
  * Burgess et al. (2016) Burgess D., Hellinger P., Gingell I., Trávníček P. M., 2016, Journal of Plasma Physics, 82, 905820401
  * Bykov & Treumann (2011) Bykov A. M., Treumann R. A., 2011, A&ARv, 19, 42
  * Bykov et al. (2019) Bykov A. M., Vazza F., Kropotina J. A., Levenfish K. P., Paerels F. B. S., 2019, Space Sci. Rev., 215, 14
  * Chang et al. (1990) Chang C. L., Wong H. K., Wu C. S., 1990, Phys. Rev. Lett., 65, 1104
  * Churazov et al. (2021) Churazov E., Khabibullin I., Lyskova N., Sunyaev R., Bykov A. M., 2021, A&A, 651, A41
  * Churazov et al. (2023) Churazov E., Khabibullin I., Bykov A. M., Lyskova N., Sunyaev R., 2023, A&A, 670, A156
  * Davidson et al. (1972) Davidson R. C., Hammer D. A., Haber I., Wagner C. E., 1972, Physics of Fluids, 15, 317
  * Davidson et al. (2004) Davidson R. C., Kaganovich I., Qin H., Startsev E. A., Welch D. R., Rose D. V., Uhm H. S., 2004, Physical Review Accelerators and Beams, 7, 114801
  * Ellison et al. (2005) Ellison D. C., Decourchelle A., Ballet J., 2005, A&A, 429, 569
  * Farris et al. (1992) Farris M. H., Russell C. T., Thomsen M. F., Gosling J. T., 1992, J. Geophys. Res., 97, 19121
  * Formisano et al. (1975) Formisano V., Russell C. T., Means J. D., Greenstadt E. W., Scarf F. L., Neugebauter M., 1975, J. Geophys. Res., 80, 2013
  * Huntington et al. (2015) Huntington C. M., et al., 2015, Nature Physics, 11, 173
  * Johlander et al. (2022) Johlander A., et al., 2022, Geophysical Research Letters, 49, e2021GL096335
  * Kato & Takabe (2008) Kato T. N., Takabe H., 2008, ApJL, 681, L93
  * Kato & Takabe (2010) Kato T. N., Takabe H., 2010, ApJ, 721, 828
  * Kropotina et al. (2019) Kropotina J., Bykov A., Krassilchtchikov A., Levenfish K., 2019, in Voevodin V., Sobolev S., eds, Supercomputing. Springer International Publishing, Cham, pp 242–253
  * Kropotina et al. (2021) Kropotina J. A., Webster L., Artemyev A. V., Bykov A. M., Vainchtein D. L., Vasko I. Y., 2021, ApJ, 913, 142
  * Lemoine et al. (2019) Lemoine M., Gremillet L., Pelletier G., Vanthieghem A., 2019, Phys. Rev. Lett., 123, 035101
  * Marcowith et al. (2016) Marcowith A., et al., 2016, Reports on Progress in Physics, 79, 046901
  * Markevitch & Vikhlinin (2007) Markevitch M., Vikhlinin A., 2007, Physics reports, 443, 1
  * Matthews (1994) Matthews A. P., 1994, J. Comput. Phys., 112, 102
  * Medvedev & Loeb (1999) Medvedev M. V., Loeb A., 1999, ApJ, 526, 697
  * Milosavljević & Nakar (2006) Milosavljević M., Nakar E., 2006, ApJ, 641, 978
  * Morse & Nielson (1971) Morse R. L., Nielson C. W., 1971, Physics of Fluids, 14, 830
  * Nishigai & Amano (2021) Nishigai T., Amano T., 2021, Physics of Plasmas, 28, 072903
  * Park et al. (2015) Park H. S., et al., 2015, Physics of Plasmas, 22, 056311
  * Pelletier et al. (2017) Pelletier G., Bykov A., Ellison D., Lemoine M., 2017, Space Sci. Rev., 207, 319
  * Petrukovich & Chugunova (2021) Petrukovich A. A., Chugunova O. M., 2021, Journal of Geophysical Research (Space Physics), 126, e29004
  * Petrukovich et al. (2019) Petrukovich A. A., Chugunova O. M., Shustov P. I., 2019, Annales Geophysicae, 37, 877
  * Pollock et al. (2016) Pollock C., et al., 2016, Space Science Reviews, 199, 331
  * Raymond (2018) Raymond J. C., 2018, Space Sci. Rev., 214, 28
  * Russell et al. (2016) Russell C. T., et al., 2016, Space Science Reviews, 199, 189
  * Sagdeev (1966) Sagdeev R. Z., 1966, Rev Plasma Phys, 4, 23
  * Sironi et al. (2015) Sironi L., Keshet U., Lemoine M., 2015, Space Sci. Rev., 191, 519
  * Spitkovsky (2008) Spitkovsky A., 2008, ApJ, 673, L39
  * Sundberg et al. (2017) Sundberg T., Burgess D., Scholer M., Masters A., Sulaiman A. H., 2017, ApJ, 836, L4
  * Takabe (2023) Takabe H., 2023, Physics of Plasmas, 30, 030901
  * Telezhinsky et al. (2012) Telezhinsky I., Dwarkadas V. V., Pohl M., 2012, Astroparticle Physics, 35, 300
  * Turner et al. (2021) Turner D. L., et al., 2021, The Astrophysical Journal Letters, 911, L31
  * Weibel (1959) Weibel E. S., 1959, Phys. Rev. Lett., 2, 83
  * Winterhalter & Kivelson (1988) Winterhalter D., Kivelson M. G., 1988, Geophys. Res. Lett., 15, 1161
  * Yang et al. (2020) Yang Z., et al., 2020, The Astrophysical Journal Letters, 901, L6


# The Bellman equation and optimal local flipping strategies for kinetic Ising
models

F. Caravelli  Theoretical Division (T4), Los Alamos National Laboratory, Los
Alamos, New Mexico 87545, USA

###### Abstract

There is a deep connection between thermodynamics, information and work
extraction. Ever since the birth of thermodynamics, various types of Maxwell
demons have been introduced in order to deepen our understanding of the second
law. Thanks to them it has been shown that there is a deep connection between
thermodynamics and information, and between information and work in a thermal
system. In this paper we study the problem of energy extraction from a
thermodynamic system satisfying detailed balance, from an agent with perfect
information, e.g. that has an optimal strategy, given by the solution of the
Bellman equation, in the context of Ising models. We call these agents
kobolds, in contrast to Maxwell’s demons which do not necessarily need to
satisfy detailed balance. This in stark contrast with typical Monte Carlo
algorithms, which choose an action at random at each time step. It is thus
natural to compare the behavior of these kobolds to a Metropolis algorithm.
For various Ising models, we study numerically and analytically the properties
of the optimal strategies, showing that there is a transition in the behavior
of the kobold as a function of the parameter characterizing its strategy.

Figure 1: Representation of the game: the kobolds (a minor demon) observes the
phase space of the Ising model and decides to attempt to flip a spin according
to an optimal strategy in order to gain energy, attempting to locally lower
the energy. However, detailed balance is enforced, and the spin flip can be
rejected by the requirement of time reversibility.

##  I Introduction

In later years there has been a lot of interest in the connection between
thermodynamics, information and the role that thought experiments involving
beings such as Maxwell demons play revmodphys ; price ; pekola ; kosloff ;
brunner ; plesch . An example of a Maxwell demon maxwell  in a thermodynamic
system is a being able to observe the speed of particles approaching a gate
separating two gases, and open it if only the particle speed is such that the
temperature of one of the gases can be raised, or acts in order for the
pressure of one can be raised (temperature vs pressure demons). These thought
experiments are useful to understand if the second law of thermodynamics can
be locally violated. This line of thought has led later many researchers to
actually propose various versions of the Maxwell demon, including notoriously
Szilard szilard ; lubkin , Brillouin brillouin , Landauer landauer  and
Bennett bennett . A typical Maxwell demon is able, quintessentially but with
some restriction in certain cases, to do anything on the system. Depending on
the point of view one takes, these thought experiments often can say something
both about the system and the demon itself. If one assumes that thermodynamics
is valid, then one can uses these arguments to infer what demons cannot do
skordos . On the other hand, if one assumes that the ability of the demons are
valid, how thermodynamics can be violated. Either way, these thought
experiments present a valid challenge for physicists, forcing them to think
carefully about the nature of thermodynamic laws. As in the case of Landauer,
and later Bekenstein and Hawking bekenstein1 ; bekenstein2 ; hawking , these
gedankexperiments with demons can lead to new discoveries about extreme
regimes, ranging from nanoscale devices to black holes. A resolution of the
violation of the second law can be obtained by assigning information to
molecules or atoms, and introducing logical irreversibility on the operation
of the operation performed by the demon, leading to an entropy increase due to
information erasure. If however no information is erased and the operation is
reversible, the second law is not violated. In fact, physical implementations
of a Maxwell demon’s have been proposed lloyd  and realized phot , confirming
such picture. These ideas can also be extended from the classical to the
quantum realm, proving their extreme generality lande ; peres ; zurek , and
has also important applications to quantum thermodynamics delrio ; horodecki
and non-equilibrium statistical mechanics jarzynski .

Maxwell demons are often very far from real beings, in the sense that they
stretch and bend slightly the roles of the games, and time is not as
important. Innovation in technology is however a constant struggle between
satisfying physical laws while gaining as much as possible from the system of
interest, whether this is energy or computing power, or something else. In
this spirit, we consider here more mundane types of Maxwell demons, those
which play by physical laws but with optimal strategies in complex landscapes.
An optimal strategy is one which, given any location of the agent in a state
space of a system, performs an action which maximizes a certain function.
These optimal strategies however require perfect “information”, e.g. the demon
knows everything about the phase space of the system beforehand, and in
particular has full phase space observability. We will use, for this purpose,
the complex landscape of an Ising model. We note that this is not information
as it is commonly stated in the theory of communication, but in the game
theoretic sense. This simply means that the Kobold has perfect knowledge of
energy landscape in the phase space, but has to satisfy detailed balance when
performing a move; it thus follows a Markov decision process.

In particular, this paper is concerned with the following problem. Consider a
thermodynamic demon which knows everything about the dynamics of an Ising
model, and desires to reach the ground state of the system by acting locally
on a certain spin, while at fixed temperature. We choose the demon to act
locally to closely represent the behavior of a physical device, or a Monte
Carlo algorithm. For this reason, we introduce a kinetics for the Ising system
which preserves detailed balance and thus microscopic reversibility. The
question we then ask is: how hard would it be to act optimally, and to
navigate through the phase space of the model in order to reach one of the
states with minimum energy, while satisfying the detailed balance? Such entity
is not quite as powerful as a Maxwell demon, as it still has to obey the laws
of detailed balance and can only act locally on a single spin. Thus it
resembles more of a kobold of german folklore, or goblin in english culture
rather than a full fledged demon. However, the kobold has a perfect knowledge
of the state space of the model (perfect information). A cartoon
representation of the game is shown in Fig. 1. In addition, the kobold values
time, and has an associated parameter with which discounting occurs;
discounting is associated to a loss of energy as a function of time, and is a
parameter that can be changed by the kobold. While microscopic reversibility
is preserved at any intermediate step, there is a particular state at which
the kobold decides to not act anymore, and saves the energy stored. In this
sense, eventually microscopic reversibility is broken by the presence of a “no
action” by the kobold.

To answer this question, we need to specify some further details. First, we
assume that we have a statistical physics system, described by spins variables
$s_{i}=\pm 1$, and an energy written in the form

$H=\sum_{ij}J_{ij}s_{i}s_{j}.$ (1)

By navigating, we mean that we are able to observe our current phase space
state $\phi=[s_{1},\cdots,s_{n}]$, and that we need to act locally via a map
$a:\phi\rightarrow\phi$, where the action $a$ can act only on one spin or none
in order to gain energy from the system.

This setup is not different from the typical Markov Chain Monte Carlo, used to
thermalize or anneal a statistical model. For instance, the Metropolis-
Hastings metropolis ; hastings  or Glauber glauber  algorithm for the time
evolution of a spin model: in the case of the latter the transition rates are
determined by a spin flip. The model dynamics is given by the spin flip
probability

$W_{ij}=p(\phi_{i}\rightarrow\phi_{j})=\frac{e^{-\beta\Delta
U(\phi_{i},\phi_{j})}}{1+e^{-\beta\Delta U(\phi_{i},\phi_{j})}}$ (2)

implying that the system satisfies detailed balance, and can thus approach
thermodynamic equilibrium onsager . Then, we have an external clock with a
discrete time which labels the operations and discounting.

Let us note that Glauber dynamics preserves the entropy of the model. This is
because it the probabilty distributions and transition rate satisfy
$P_{i}W_{ij}=P_{j}W_{ji}$, from which we obtain that the Schnakenberg entropy
production schnakenberg  is zero, e.g.

$\Pi(t)=\frac{k}{2}\sum_{ij}(P_{i}(t)W_{ij}-P_{j}(t)W_{ji})\log\frac{P_{i}(t)W_{ij}}{P_{j}(t)W_{ji}}=0.$
(3)

Thus, the change in free energy of the model reduces to the change in
potential energy of the system. Interestingly, this also means that a kobold
acts as a reversible computer frank . Another point of view is that the system
represents the battery, and the kobold is an operation by which we wish the
extract energy from the system without wasting it in entropy production.

One question one might ask is why is a kobold different from a standard Monte
Carlo “demon”. We could in fact consider a system out of equilibrium and let
it thermalize, end extract the energy as a result of thermalization. However,
if $\bar{E}$ is the expected thermalization energy, if $E(0)<\bar{E}$ the
demon would actually have to supply energy to the system. A kobold, on the
other hand, has an optimal strategy for any initial condition, and can act on
its discounting factor to try to maximize it. In doing so, it can also reduce
fluctuations over the gained energy.

The paper is organized as follows. We will first discuss the theoretical
underpinning of kobolds, envisaged as a stochastic Bellman equation for the
optimal strategy bellman1 ; bellman2 ; reinfl . We will derive an analytical
formulation for the Bellman equation restricted to the case of Glauber
dynamics, and show that this can be written in terms of projector operators on
the phase space. We will then provide numerical results both for the optimal
strategies and discounted gained energy. In particular, we will show that
depending on the discounting factor, kobold can take a greedy or a wise
strategy approach; a greedy strategy is one in which the local spin flip
action is such that it always minimizes the energy, while a wise strategy
allows for local increases in energy. In this respect, a greedy approach can
be compared to a low temperature Metropolis algorithm, while a wise approach
to a higher temperature one. However, the analogy ends there, as the kobold
tries independently to maximize the energy extracted from the system. This is
the reason why we compare the behavior of such an agent to a “Metropolis”
agent.

Conclusions follow.

##  II Optimal strategies and the Bellman equation

###  II.1 The discounted reward

We assume the following rules for the abilities and limitations of the kobold.
First, the kobold knows everything about the system: the couplings, a complete
picture of the phase space, and thus also the location of the ground state. In
the language of economics and game theory, the kobold has perfect information
(rather than imperfect, which would occur if the kobold only knows a certain
part of the phase space). However, it has some restrictions, most importantly
it has to respect constraints such as the detailed balance. This implies that
while it might decide to attempt a spin flip to gain energy, whether this
occurs or not is determined by the acceptance rule at the temperature of
system. The second restriction is that the kobold can only operate locally in
time and space, meaning that it has to perform actions sequentially, and one
spin at the time, but any spin of choice. The third restriction is that the
kobold does not get to choose the initial state from which it starts to
operate on the system. What the kobold gains by playing this game is energy.
If at each time step the kobold gains $\Delta U_{t}$ of energy via a certain
action (spin flips), after $T$ steps the kobold will have gained

$R_{\gamma}=\sum_{t=1}^{\tau}R_{\gamma}(t)=\sum_{t=1}^{\tau}\gamma^{t}\Delta
U(\phi_{t},\phi_{t-1});$ (4)

in the equation above $\gamma\in(0,1]$ is a discounting factor: the kobold
knows the rules of finance, and knows that energy today is better than energy
tomorrow.

In practice, if $\gamma<\gamma^{\prime}<1$, this means that even if the
actions are such that the energy is positive at each time step one has
$R_{\gamma}<R_{\gamma^{\prime}}$. When $\gamma=1$, we simply have

$R_{1}=\sum_{t=1}^{\tau}R_{1}(t)=\sum_{t=1}^{\tau}(U_{t}-U_{t-1})=U_{\tau}-U_{0}.$
(5)

In the equation above, and in the context of the Bellman equation, the
parameter $\gamma$ characterizes the kobold strategy, as it defines how
important is time to the kobold.

###  II.2 Bellman equation

In general, a kobold agent wants to maximize the reward over time (possibly an
infinite horizon) via local actions $a_{\phi}=\pi(\phi)$, meaning that by
observing the local state of the system, the kobold can attempt to flip one
spin at time $t$ with the intention of gaining energy. Thus, the parameter
$\gamma$ represents how important is the speed at which the kobold attempts to
reach the ground state. It is important to stress that the form of the
discounted energy of eqn. (4) can be obtained by assumptions on the time
invariance of the optimal strategy, and is thus natural.

In the following, we use a notation in which $\pi(\phi)$ determines both an
action and a state, as the two in this context are the same. If the action
$\pi(\phi)$ is to flip the k-th spin of $\phi$, then this determines a new
state $\phi^{\prime}$. Of course, if the kobold reaches the ground state, it
will want to perform no action. This implies that if the system contains $N$
spins, there can be $N+1$ possible actions. We say attempt because of course
the acceptance probability of such action is the determined by the system’s
temperature, and thus if the currently the kobold is in state $\phi$, the
action $\pi(\phi)$ of the kobold can lead to a state $\phi^{\prime}$ with
probability $P_{T}(\phi;\pi(\phi),\phi^{\prime})$, with $T$ being the
temperature. The average energy reward is then given by

$\displaystyle\langle R_{\gamma}(t)\rangle\equiv\langle
U\rangle_{t,\pi(\phi)}=\sum_{\phi^{\prime}}P_{T}(\phi;\pi(\phi),\phi^{\prime})\gamma^{t}\Delta
U(\phi_{t},\phi_{t-1}).$

For the kobold, an ideal world would be such that
$P_{T}(\phi;\pi(\phi),\phi^{\prime})=\delta_{\phi,\pi(\phi)}$. The setup is
now such that the kobold is playing a (reversible) Markov Decision Process in
discrete time. A Monte Carlo algorithm plays typically the same game, with the
difference that the action $\pi(\phi)$ is random in some form, e.g. a random
single (Metropolis or Glauber) or multiple (Kawasaki) spin flip eli . The
kobold instead plays like an economist: uses dynamic programming and the
notion of discounting, knows everything about the system and applies an
optimal strategy $\pi^{*}(\phi)$. Under these assumptions, an infinite-horizon
decision problem takes the form of a maximization, e.g.

$\displaystyle
V^{*}_{T,\gamma}(\phi)=\text{max}_{a_{0},\cdots,a_{\tau}}\langle
R_{\gamma}(a_{0},\cdots,a_{\tau})\rangle,$ (6)

where $\phi$ is the starting state and $a_{i}$ the actions taken at each time
step. The optimal reward $V_{T,\gamma}^{*}(\phi)$ is the maximum (discounted)
reward that the kobold can obtain starting from a certain state in phase
space. For $\gamma=1$, this corresponds to the maximum energy the kobold can
extract. We then see the reason of such gedankexperiment. The kobold
represents the best possible algorithm designed to reach the ground state of
the system, and then also represents the best possible line of action an
algorithm can take at fixed temperature. Here, the assumption is that the
system is not annealed, probabilities are time independent and thus the kobold
can act with an infinite time horizon. For this type of problems, the optimal
solution $\pi^{*}(\phi)$, and the optimal reward $V_{T,\gamma}^{*}(\phi)$ can
be obtained by solving the stochastic Bellman equation. Above, $\pi^{*}(\phi)$
is the best possible action that the kobold takes if it finds itself in state
$\phi$. The Bellman equation is given by the linear relationship

$V_{T,\gamma}^{*}(\phi)=\sum_{\phi^{\prime}}P_{T}(\phi;\pi^{*}(\phi),\phi^{\prime})\Big{(}-\Delta
U(\phi,\phi^{\prime})+\gamma V^{*}_{T,\gamma}(\phi^{\prime})\Big{)}.$ (7)

The equation above is written implicitly: to solve it one would have already
to have $\pi^{*}$, which we do not. There are many ways to solve it however,
and we use here the strategy iteration method which starts from a random
initial state $V_{0}(\phi)$, and then uses the iteration

$\displaystyle V_{T,\gamma;k+1}(\phi)$ $\displaystyle=$ $\displaystyle
max_{a}\sum_{\phi^{\prime}}P_{T}(\phi;a,\phi^{\prime})\Big{(}-\Delta
U(\phi,\phi^{\prime})$ $\displaystyle\hskip 99.58464pt+\gamma
V_{T,\gamma;k}(\phi^{\prime})\Big{)},$ $\displaystyle\pi^{*}(\phi)$
$\displaystyle=$ $\displaystyle arg\
max_{a}\sum_{\phi^{\prime}}P_{T}(\phi;a,\phi^{\prime})\Big{(}-\Delta
U(\phi,\phi^{\prime})$ $\displaystyle\hskip 108.12054pt+\gamma
V^{*}_{T}(\phi^{\prime})\Big{)}.$

which gives, at convergence, the optimal strategy $\pi^{*}(\phi)$.

Note that here we immediately face a computational problem. While the number
of actions is $N+1$, the number of states $\phi$ scales exponentially with the
size of the system, and thus unlike the kobold we will have to use small
systems by Monte Carlo standards. Despite such curse of dimensionality, the
Bellman solution for optimal sequential Markov decision processes is regarded
as the most feasible one reinfl . Nonetheless, this will be sufficient to
obtain a picture of the complexity of Ising models from the point of view of
the kobold, as we will discuss in a moment.

###  II.3 Application in kinetic statistical mechanics.

Let us provide an immediate comment on why such technique is useful. Since we
are essentially finding the functions $V^{*}$ and $\pi^{*}$ for every state of
the model, if one is interested only in the energy of the model, exhaustive
search (or brute force) works much better than solving the Bellman equation.
However, Bellman’s optimal strategy and discounted energy give a more complete
picture of how the state space of an Ising model is tangled, and how one could
unentangle it and navigate through the states, or make less blind moves in
Monte Carlo algorithms.

First, we note that eqn. (7) can be written in an explicit form using the
Glauber transition of eqn. (2), given a certain action $\pi$ on the phase
space:

$\displaystyle
P_{T}(\phi;\pi(\phi),\phi^{\prime})=\frac{\delta_{\phi\phi^{\prime}}+\delta_{\phi^{\prime}\pi(\phi)}e^{-\beta\Delta
U(\phi,\pi(\phi))}}{1+e^{-\beta\Delta U(\phi,\pi(\phi))}}$ (8)

from which it follows that we can write the Bellman equation in the form

$\sum_{\phi^{\prime}}\mathcal{O}(\phi,\phi^{\prime})V_{T,\gamma}^{*}(\phi^{\prime})=-\frac{\Delta
U\big{(}\phi,\pi^{*}(\phi)\big{)}}{1+e^{-\beta\Delta U(\phi,\pi^{*}(\phi))}}$
(9)

where

$\displaystyle\mathcal{O}(\phi,\phi^{\prime})$ $\displaystyle=$
$\displaystyle\delta_{\phi\phi^{\prime}}\frac{1+e^{-\beta\Delta
U(\phi,\pi^{*}(\phi))}-\gamma}{1+e^{-\beta\Delta U(\phi,\pi^{*}(\phi))}}$ (10)
$\displaystyle-$ $\displaystyle
P_{\pi^{*}}\big{(}\phi,\phi^{\prime}\big{)}\frac{\gamma e^{-\beta\Delta
U(\phi,\pi^{*}(\phi))}}{1+e^{-\beta\Delta U(\phi,\pi^{*}(\phi))}}$ (11)

Above, $P_{\pi^{*}}\big{(}\phi,\phi^{\prime}\big{)}$ is a state transition
matrix, e.g.
$P_{\pi^{*}}\big{(}\phi,\phi^{\prime}\big{)}=\delta_{\phi^{\prime},\pi^{*}(\phi)}$.
After a brief calculation, it follows that we can write the exact solution for
$V^{*}$ in the form:

$\displaystyle V^{*}_{T,\gamma}(\phi)=\sum_{\phi^{\prime}}(I-\gamma
DP_{\pi^{*}})^{-1}_{\phi,\phi^{\prime}}\tilde{U}(\phi^{\prime}).$ (12)

where

$\displaystyle\tilde{U}(\phi)$ $\displaystyle=$ $\displaystyle-\Delta
U(\phi,\pi^{*}(\phi))R(\phi),\ \ \
D(\phi,\phi^{\prime})=\delta_{\phi,\phi^{\prime}}R(\phi)$ (13)
$\displaystyle\text{and where }R(\phi)=\frac{e^{-\beta\Delta
U(\phi,\pi^{*}(\phi))}}{1-\gamma+e^{-\beta\Delta U(\phi,\pi^{*}(\phi))}}.\ \ \
\ \ \ $

Assuming a random initial condition, the average maximum discounted energy
obtained from the kobold is given by

$\langle
V_{T,\gamma}\rangle=\frac{1}{2^{N}}\sum_{\phi}V^{*}_{T,\gamma}(\phi),$ (14)

where $N$ is the number of spins. This is the average (discounted) energy that
an optimal player (a kobold) can achieve by acting locally on the system, and
against the temperature, by reaching the ground state. In general, the Bellman
equation selects one possible action for each phase space state, but there
might be some other stochastic policies, which we do not consider here, such
that $\pi$ is a stochastic function as well. This helps however in evaluating
the amount of information that a policy contains. If we have $M$ possible
actions per phase state, and $K$ states in phase space, then we have $M^{K}$
possible policies. We can then evaluate the information in the policy using
the entropy (in a base 2) $S_{\pi}=K\log_{2}M$, which is an estimation of the
amount of information associated to a kobold, including an optimal one.

###  II.4 Example: classical Zener Hamiltonian.

Let us consider first a simple application of the equations above. We take as
a model the classical Zener Hamiltonian for an Ising spin in an external
field, which is given by $H=\mu_{0}hs$, with $s=\pm 1$ representing a single
spin. In this case we have only two states, and the optimal strategy is easy
to guess: it is a spin flip if $sign(s)=-sign(h)$, and do nothing otherwise.
Let us assume $h>0$. Then,

$\displaystyle\begin{pmatrix}V_{T,\gamma}^{*}(+)\\\
V_{T,\gamma}^{*}(-)\end{pmatrix}$ $\displaystyle=$
$\displaystyle\begin{pmatrix}\frac{2-\gamma}{2}-\frac{\gamma}{1+\gamma}&0\\\
\frac{\gamma
e^{\beta\mu_{0}h}}{1+e^{\beta\mu_{0}h}}&\frac{1+e^{\beta\mu_{0}h}-\gamma}{1+e^{\beta\mu_{0}h}}\end{pmatrix}^{-1}\begin{pmatrix}0\\\
\mu_{0}h\end{pmatrix}$ $\displaystyle=$
$\displaystyle\frac{1}{d}\begin{pmatrix}\frac{1+e^{\beta\mu_{0}h}-\gamma}{1+e^{\beta\mu_{0}h}}&-\frac{\gamma
e^{\beta\mu_{0}h}}{1+e^{\beta\mu_{0}h}}\\\
0&\frac{2-\gamma-\gamma^{2}}{2+2\gamma}\end{pmatrix}\begin{pmatrix}0\\\
\mu_{0}h\end{pmatrix}$ $\displaystyle=$
$\displaystyle\frac{\mu_{0}h}{1+e^{\beta\mu_{0}h}-\gamma}\begin{pmatrix}-(\gamma
e^{\beta\mu_{0}h})\frac{2+2\gamma}{2-\gamma-\gamma^{2}}\\\
1+e^{\beta\mu_{0}h}\end{pmatrix}$

where
$d=\frac{2-\gamma-\gamma^{2}}{2+2\gamma}\frac{1+e^{\beta\mu_{0}h}-\gamma}{1+e^{\beta\mu_{0}h}}$.
It follows that the average energy the kobold can gain is

$\displaystyle\langle
V^{*}_{T,\gamma}\rangle=\frac{\mu_{0}h}{2(1+e^{\beta\mu_{0}h}-\gamma)}\big{(}1+e^{\beta\mu_{0}h}(1-2\frac{\gamma+\gamma^{2}}{2-\gamma-\gamma^{2}})\big{)}$

It follows that with a proper discounting strategy, e.g.
$\gamma\in[0,\frac{1}{2}(\sqrt{3}-1)]$, the Kobold can always extract energy
from the system, at any temperature.

###  II.5 Properties of $P_{\pi^{*}}$

In this section we discuss the properties of $P_{\pi^{*}}$, which are useful
to understand the behavior of a kobold in the case of an Ising system via
local spin flipping.

In fact, as it turns out, the matrix $P_{\pi^{*}}(\phi,\phi^{\prime})$ is a
projector operator, as its spectrum is only composed of $0$’s and $1$’s.
Intuitively, the proof follows from the fact that we have a directed graph of
outdegree equal to one, with an absorbing state, which is what we discuss
below.

The important quantity in the Bellman equation for kinetic Ising models is the
transition matrix $P_{\pi^{*}}$, which could be in principle a permutation
matrix. However, as it turns out, the optimal strategy is always such that
$P_{\pi^{*}}$ is a projector operator, e.g. $P_{\pi^{*}}^{2}=P_{\pi^{*}}$, as
we will prove in a moment. The kobold’s strategy is to lower the energy, by
transitioning between states and eventually to one of the absorbing states.
Then, the basins of attraction of a particular absorbing state $G$ as a sink
in a directed tree, similar to Fig. 5, which follows from the fact that the
outdegree of every node is always one. If we have $G$ absorbing states, then
the matrix $P_{\pi^{*}}$ can be written in block diagonal form. Let us call
these sub-blocks $P^{G}_{\pi*}$; then, simply one has
$P_{\pi^{*}}=\bigoplus_{G}P^{G}_{\pi*}$ and then we can focus on each sub-
block. Now, we can label the nodes such that if $i$ is the numerical value of
a certain node and $i^{\prime}$ down the tree, then $l(i)>l(i^{\prime})$. Such
labeling is always possible because we have a directed acyclic graph, and one
example is shown in Fig. 5. The fact that it is acyclic follows from the fact
that we have an absorbing state for a graph with outdegree one. In fact,
assume by absurd that one has a cycle in such a graph with outdegree one. If
it is a cycle, there cannot be an absorbing state, and thus there must be at
least one node with degree three and with outdegree two, since one directed
edge goes inside the cycle and the other in the direction of the absorbing
state. However, this is incompatible with the fact that the outdegree must
always be one. With the labeling $l$ we see that these energy state
transitions correspond to an upper triangular matrix elements of
$P^{G}_{\pi^{*}}$. Since the absorbing state is the only state which has null
action (the kobold will want to remain in that state), this means the
absorbing state state is the only element with a one on the diagonal. Each
sub-block can always be written, via the state relabeling, as

$\displaystyle P^{G}_{\pi^{*}}=\begin{pmatrix}0&\times&\times&\times&\times\\\
\cline{1-1}\cr\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&0&\times&\times&\times\\\
\cline{2-2}\cr&\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&0&\times&\times\\\
\cline{3-3}\cr&\makebox(0.0,0.0)[]{\text{\huge
0}}&\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&0&\times\\\
\cline{4-4}\cr&&&\lx@intercol\hfil\hfil\lx@intercol\vrule\lx@intercol&1\\\
\cline{5-5}\cr\end{pmatrix}$ (15)

from which it follows, since the matrix is triangular, that if such sub-block
is $D$ dimensional, then $D-1$ eigenvalues are $0$, and only one is $1$. This
is enough to show it is is a projector operator, and that it satisfies the
condition $(P_{\pi^{*}}^{G})^{2}=P_{\pi^{*}}^{G}$. Then, the spectrum of
$P_{\pi^{*}}$ is simply determined by the number of absorbing states. In fact,
the number of $1$s correspond to the number of absorbing states of the system,
and the rest of the spectrum contains only zeros. Since every diagonal sub-
block is a projector operator, so is $P_{\pi^{*}}$.

This result is useful for the following reason. First, in the limit
$\gamma\rightarrow 1$, we have $D\rightarrow I$ and $R(\phi)\rightarrow 1$. We
can write explicitly the inverse, in the neighborhood of $\gamma=1^{-}$ and
using the property that $P^{2}_{\pi^{*}}=P_{\pi^{*}}$, as

$\displaystyle(I-\gamma
P_{\pi^{*}})^{-1}=I+\frac{\gamma}{1-\gamma}P_{\pi^{*}},$ (16)

and then the discounted energy values are simply given by

$\displaystyle V^{*}_{T,\gamma\approx 1}(\phi)$ $\displaystyle=$
$\displaystyle-\sum_{\phi^{\prime}}\Big{(}I+\frac{\gamma}{1-\gamma}P_{\pi^{*}}\Big{)}_{\phi,\phi^{\prime}}\Delta
U(\phi^{\prime},\pi^{*}(\phi^{\prime}))$

We can thus use the exact inverse to regularize the limit $\gamma\rightarrow
1$, which is otherwise ill-defined in the general case. We can use the
following assumption. Typically, $\partial_{\gamma}\pi^{*}=0$ almost
everywhere. Then, assuming that the optimal strategy can analytically extended
from $\gamma=1-\epsilon$ to $\gamma=1^{-}$, i.e. that it is constant, we can
obtain a discounted effective energy of the form

$\displaystyle\tilde{V}^{*}_{T,\gamma}(\phi)$ $\displaystyle=$
$\displaystyle\lim_{\gamma\rightarrow 1^{-}}(1-\gamma)V^{*}_{T,\gamma}(\phi)$
$\displaystyle=$
$\displaystyle\sum_{\phi^{\prime}}(P_{\pi^{*}})_{\phi,\phi^{\prime}}\Delta
U(\phi^{\prime},\pi^{*}(\phi^{\prime})),$

which is a regularized $\gamma=1$ limit for the optimal strategy. In this
regime, all absorbing states are ground states.

In fact, we have been careful in calling $G$ the absorbing state and not the
ground state. If $\gamma=1$, of course all absorbing states are ground states,
but at $\gamma<1$ this is not guaranteed, although likely for $\gamma\approx
1$. Intuitively, this is because the kobold might find more rewarding to stop
at a certain state rather than attempting to reach the ground state, given
that there is a cost in how long the game takes, and in the intermediate steps
that might actually increase rather than lowering the energy.

##  III Numerical results

We consider, for the purpose of this paper, four types of Ising models, of
which three ferromagnetic and one frustrated. The testbed of our analysis are
the Ising models given below

$H=\begin{cases}\frac{1}{2}\sum_{i=1}^{N-1}s_{i}s_{i+1}&\ \ \ \ \ \text{Ising
1D},\\\ \frac{1}{2}\sum_{<ij>=1}^{N_{1}N_{2}}s_{i}s_{j}&\ \ \ \ \ \text{Ising
2D},\\\ \frac{1}{2N}\sum_{ij=1}^{N}s_{i}s_{j}&\ \ \ \ \ \text{Curie-
Wei}\ss,\\\ \frac{1}{2\sqrt{N}}\sum_{ij=1}^{N}J_{ij}s_{i}s_{j}&\ \ \ \ \
\text{Spin Glass}\ \ \ J_{ij}=\pm 1.\end{cases}$ (18)

where in the latter case we consider $P(J_{ij}=\pm)=1/2$. Since solving the
Bellman equation requires solving iteratively a vectorial equation of the size
of the phase space, we are forced to study relatively small systems, with $12$
spins.

Figure 2: Top: maximum discounted energy $V^{*}_{T,\gamma}$ as a function of
the states $\phi$ for $N=12$, ordered from high to low, and $\langle
V^{*}_{T,\gamma}\rangle$ shown in the inset. We see the difference between the
spin glass, Ising 2D and Ising 1D and the Curie-Weiß model. Since
$\gamma=0.95$, the low energy states are associated to ground states of the
model. Bottom: The profile of the function $\langle V_{T,\gamma}^{*}\rangle$
as a function of $T$ and $\gamma$. While the dependence on $T$ is negligible
for these values of $T$, the dependence on $\gamma$ is strong, due to the fact
that high discounting rates imply less value to the energy obtained in the
future.

Figure 3:  Change in kobold strategy as a function of temperature and
discounting rate. A greedy strategy implies only lowering the energy at every
step, while a wise one can also accept intermediate energy increases.

Figure 4: Kobold strategy tree (top are absorbing states) as a function of the
$\gamma$ for $T=3$ for the Ising 2D and spin glass cases. Red curves are low
discounting, while blue are high discounting rates. We see the emergence of
two extra absorbing states on top of the Ising 2D ground states.

We iterate eqn. (II.2) until
$\frac{1}{2^{N}}\|V_{k+1}-V_{k+1}\|^{2}<\epsilon$, with $\epsilon=10^{-5}$,
read the strategy out and sort the spin states with descending discounted
potential $V^{*}(\phi)_{T,\gamma}$. The plots are shown in Fig. 2 for
$\gamma=0.95$. In the inset of Fig. 2, we plot $\langle
V_{T,\gamma}^{*}\rangle$ for each model at $\gamma=0.95$. The quantity
$\langle V_{T,\gamma}^{*}\rangle$ is the average utility (discounted energy)
that the kobold can obtain from the Ising model, assuming that we let him
start from any spin state at at infinite time. Such average shows, in
arbitrary units of energy and assuming all being equal, which model lead to
more gains via a Monte Carlo method, withing the assumption of only local
moves. Intuitively, we see that the spin glass lead to more gains than the
others, with the Ising 2D being harder than Ising 1D model, and the mean field
ferromagnetic model being the easiest in comparison. Such hierarchy seemingly
makes sense from a computational perspective, as we would expect the spin
glass to have a longer way down to the ground state. However, we note that the
two key parameters we need analyze are $\langle V_{T,\gamma}^{*}\rangle$ as a
function of temperature, of course, and as a function of the parameter
$\gamma$. If $\gamma\approx 1$, time for the kobold is not an issue. At strong
discounting, however, longer chains to reach the ground state imply losses.
This can be seen in Fig. 2 (bottom), in which we see that for all models at
stronger discounting ($\gamma\rightarrow 0$) the effective value is reduced.

An interesting question is how the optimal strategy of the kobold changes as a
function of $T$ and $\gamma$. We explore this in two ways. First, we analyze
whether the only allowed strategies are those that go down in energy, which we
call “greedy” kobold strategies ($S=+1$). On the other hand, if the kobold
uses strategies such that at a certain point it actually increases momentarily
the energy, these are ”wise” ($S=-1$), as they allow the kobold to reach a
lower state faster. We find that for the 1D Ising model and the Curie-Weiß
model the strategy is always greedy. However, we plot $S$ as a function of
$\gamma$ and $T$ in Fig. 3. We see that for the 2D Ising model, the parameter
$S$ does not depend on $T$, but it depends on $\gamma$ and for $\gamma>0.8$,
the strategy becomes wise, while it is greedy for $\gamma<0.8$. For the spin
glass case, instead, we find a dependence on the temperature as well, with a
transition between a greedy and wise strategy at approximately $\gamma=0.62$
for $T=3$, and $\gamma\approx 0.55$ for $T=0.1$.

The maximum number of steps required (on average) can be estimated by the
depth of the optimal policy tree. For the cases of the Ising 2D and Spin
Glasses, these are shown in Fig. 4, showing the difference between a high and
low discounting rates for $T=3$. As we see from the change of the trees, shown
in different colors, the kobolds uses a different strategy at different
discounting rates. Interestingly, at high discounting rates the kobolds
develops two extra absorbing states, on top of the ground states of the model
for the Ising 2D case.

We can see the difference between the kobold strategy and the Monte Carlo
thermalization process in Fig. 6 for the case of the Ising 2D and Spin Glass.
The random strategy fluctuates strongly in energy, while the kobold’s strategy
reaches the ground state in a few steps.

Figure 5: The graph representation of the an optimal strategy $\pi^{*}$
restricted to the basin of attraction of the absorbing state. Since the graph
is directed and acyclic, there is a node labeling $l(i)$ such that
$l(i)>l(i^{\prime})$ if $i^{\prime}$ can be reached by a directed path from
$i$.

Figure 6: Example of Markov Chain Monte Carlo for the Ising and Spin Glass
examples, for $N=12$. We compare the dynamics to the Ising ground state of the
associated model.

##  IV Conclusions

The present paper introduced optimal flipping policies for kinetic Ising
models, using the stochastic Bellman equation as a prototypical Maxwell demon
eager to extract the energy from a thermal system. As we have shown in this
paper, these demons have less power than a typical Maxwell demon, as they can
only operate locally on a single spin, and have to satisfy the detailed
balance. They thus parallel the typical strategy of a Monte Carlo algorithm,
with the different that instead than thermalizing the model, their intent is
lowering the energy given the Markov chain transition probabilities. As we
have shown, their strategy and approach strongly depends on how fast they want
to extract energy from the system, changing from a greedy to a wise approach
when discounting is high or low respectively. Unfortunately, the Bellman
equation still is plagued by the curse of dimensionality, and only small
systems could be analyzed. Nonetheless, we have shown that these present
interesting strategy changes also for small systems.

This approach can also be interpreted as the optimal sequential strategy to
optimize Ising models in an uncertain but time-invariant environment. As a way
to overcome the curse of dimensionality, we note that this is exactly the same
problematic that reinforcement learning aims to tackle, which is the natural
extension of this work rl . In fact, in reinforcement learning agents have
“imperfect information”, and do not necessarily know fully the state space,
and in our case the energy landscape. In this sense, a natural extension of
our work is the one in which a kobold is “trained” via a reinforcement
learning algorithm, learning the strategy by testing the results of actions on
the phase space, and updating iteratively the strategy. This will be the focus
of future works.

###### Acknowledgements.

This work was carried out under the auspices of the NNSA of the U.S. DoE at
LANL under Contract No. DE-AC52-06NA25396, and in particular grant PRD20190195
from the LDRD. This paper was inspired by “Meister Eder und sein Pumuckl”.

## References

  * (1) K. Muruyama, F. Nori, V. Vedral, Rev. Mod. Phys. 81, 1 (2009)
  * (2) G.N. Price, S.T Bannerman, K. Viering, E. Narevicius, M. G. Raizen, Phys. Rev. Lett., 100, 093004 (2008).
  * (3) J.V. Koski, V.F. Maisi,T. Sagawa, and J.P. Pekola, Phys. Rev. Lett., 113, 030601 (2014).
  * (4) R. Kosloff, L. Amikam, Annu. Rev. Phys. Chem., 65, 365 (2014).
  * (5) N. Brunner, N. Linden, S. Popescu, P. Skrzypczyk, Phys. Rev. E, 85, 051117 (2012)
  * (6) M. Plesch, O. Dahlsten, J. Goold and Vlatko Vedral, Sci. Rep., 4, 6995 (2014).
  * (7) J.C. Maxwell, Theory of Heat (Longman, London, 1871).
  * (8) L. Szilard, Z. Phys. 53, 840, (1929).
  * (9) E. Lubkin, Int. J. Theor. Phys. 26, 523 (1987)
  * (10) L. Brillouin, J. Appl. Phys. 22, 334 (1951)
  * (11) R. Landauer, IBM J. Res. Dev. 5, 183, (1961)
  * (12) C.H. Bennett, Int. J. Theor. Phys. 21, 905, (1982)
  * (13) P. A. Skordos, W. H. Zurek, Am. J. Phys. 60, 876 (1992)
  * (14) J. D. Bekenstein, Phys. Rev. D 7, 2333 (1973)
  * (15) J. D. Bekenstein, Phys. Rev. D 9, 3292 (1974)
  * (16) S. W. Hawking, Nature (London) 248, 30 (1974)
  * (17) S. Lloyd, Phys. Rev. A, 56, 3374 (1997)
  * (18) M. D. Vidrighin, O. Dahlsten, M. Barbieri, M. S. Kim, V. Vedral, I. A. Walmsley, Phys. Rev. Lett. 116, 050401 (2015)
  * (19) A. Lande, Phys. Rev. 87, 267 (1952)
  * (20) A. Peres, Phys. Rev. Lett. 63, 1114 (1989)
  * (21) W. H. Zurek, Phys. Rev. A 67, 012320 (2003)
  * (22) M. Horodecki and J. Oppenheim, Nat. Comm., 4, 2059 (2013)
  * (23) L. Del Rio et al, Nature 474, 61–63 (2011)
  * (24) C. Jarzynski, Eur. Phys. J. B, 64, 331340 (2008).
  * (25) R. J. Glauber, J. Math. Phys. 4, 294 (1963);
  * (26) N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, J. of Chem. Phys. 21, pp 1087–1092 (1953)
  * (27) W.K. Hastings, In: Biometrika 57, pp 97–109 (1970)
  * (28) L. Onsager, Phys. Rev. 37, 405–426; II 38, 2265–2279 (1931)
  * (29) J. Schnakenberg, Rev. Mod. Phys. 48, 571 (1976)
  * (30) M.P. Frank, In: Phillips, I., Hafizur, R. (eds.) RC 2017. LNCS, vol. 10301, pp. 19–34, Springer, Cham (2017)
  * (31) R. E. Bellman, Dynamic Programming, Princeton University Press. Dover paperback edition (2003)
  * (32) R.E. Bellman, Proc. Nat. Acad. Sci. U S A. 38(8): 716–719 (1952)
  * (33) R. S. Sutton, A. G. Barto, Reinforcement Learning: An Introduction, Second edition, MIT University Press, Boston (2015)
  * (34) P. L. Krapivsky,S. Redner, E. Ben-Naim, A kinetic view of statistical physics, Cambridge University Press, Boston (2010) ISBN:9780511780516
  * (35) A. Barto, R. S. Sutton, Reinforcement Learning: An Introduction, Adaptive Computation and Machine Learning series (2016)