# Mathematical aspects of the Digital Annealer’s simulated annealing algorithm

Bruno Hideki Fukushima-Kimura  Faculty of Science, Hokkaido University, Kita
10, Nishi 8, Kita-ku, Sapporo, Hokkaido 060-0810, Japan
<EMAIL_ADDRESS>,  Noe Kawamoto  Graduate School of Science,
Hokkaido University, Kita 10, Nishi 8, Kita-ku, Sapporo, Hokkaido 060-0810,
Japan<EMAIL_ADDRESS>,  Eitaro Noda  Graduate School of
Science, Hokkaido University, Kita 10, Nishi 8, Kita-ku, Sapporo, Hokkaido
060-0810, Japan<EMAIL_ADDRESS>and  Akira Sakai  Faculty of
Science, Hokkaido University, Kita 10, Nishi 8, Kita-ku, Sapporo, Hokkaido
060-0810, Japan<EMAIL_ADDRESS>
###### Abstract.

The Digital Annealer is a CMOS hardware designed by Fujitsu Laboratories for
high-speed solving of Quadratic Unconstrained Binary Optimization (QUBO)
problems that could be difficult to solve by means of existing general-purpose
computers. In this paper, we present a mathematical description of the first-
generation Digital Annealer’s Algorithm from the Markov chain theory
perspective, establish a relationship between its stationary distribution with
the Gibbs-Boltzmann distribution, and provide a necessary and sufficient
condition on its cooling schedule that ensures asymptotic convergence to the
ground states.

##  1\. Introduction

Problems that involve, in a certain sense, an optimal choice of a
configuration or a set of parameters among a large number of possibilities
commonly arise in several problems of practical and theoretical interest,
especially in various fields related to engineering, computer science,
artificial intelligence, machine learning, logistics, and very-large-scale
integration (VLSI) [34]. The study of the so-called combinatorial optimization
problems [19, 29] has been revealed to be relevant primarily because of their
recurrent appearance in real-world issues and has been accompanied by a large
number of proposals for their solution.

The Metropolis algorithm [25] is a widespread Markov chain Monte Carlo (MCMC)
method for simulating the evolution of a solid in contact with a heat
reservoir to thermal equilibrium. Almost thirty years after its introduction,
Kirkpatrick et al. [18] and Černy [6] independently recognized the connection
between combinatorial optimization problems and statistical mechanics
motivated by their investigations in proposing a Monte Carlo algorithm
approach for integrated circuit design problems and the traveling salesman
problem. Their approach was essentially based upon the analogy coming from
condensed matter physics of the thermal annealing of solids for obtaining low-
energy states.

Kirkpatrick et al. [18] showed that some combinatorial optimization problems
of practical interest could be unambiguously translated (see [21] for more
examples) into the minimization problem of the Hamiltonian of an Ising model
in such a way that each optimal solution for the original problem corresponds
to a ground state of the Hamiltonian, and vice versa. The simulated annealing
algorithm (SA) is a method whose control parameter is regarded as the
temperature of a physical system and consists of iterations of the Metropolis
algorithm initially evaluated at very high temperatures and subsequently at
slowly decreasing values of the temperature in such a way as to guarantee the
convergence to a candidate for a minimal energy state. Still, under the
analogy with physical systems, it is well-known that if the cooling is too
rapid, the system cannot achieve thermal equilibrium for each temperature
value, which may result in a configuration with defects in the form of high-
energy, metastable, locally optimal structures. A treatment of this topic from
a theoretical point of view based on the Markov chain theory is present in [1,
12, 15, 22, 26, 32], where specifically in [15] was derived a necessary and
sufficient condition on the cooling speed that guarantees asymptotic
convergence of the SA to the ground states.

Achieving rigorous results in order to obtain a deep understanding of
algorithms proposed to solve combinatorial optimization problems is becoming
highly desirable, especially due to the increasing demand for methods for
efficiently solving large-scale problems. In [20], the authors explored the
mathematical foundations of dynamical system algorithms to support their
applications in coherent Ising machines (CIM) [33, 36], Kerr-nonlinear
parametric oscillators [13], and the simulated bifurcation (SB) algorithm [14]
in the search for Hamiltonian minimizers. In [9], the authors provided a
theoretical treatment of a simulated annealing algorithm based on
probabilistic cellular automata called Stochastic Cellular Automata Annealing
(SCA) [35] motivated by the potential to make use of its parallelizability to
provide solutions within a shorter amount of time as compared to SA. The
authors showed that an appropriate cooling rate in the form
$T_{k}\propto\frac{1}{\log(k)}$ is sufficient to guarantee the asymptotic
convergence to ground states. Experimental results in [9, 10, 11, 17] showed
that the SCA performed better than SA in searching for ground states in most
of the considered scenarios, while its variation $\varepsilon$-SCA (or RPA)
surpassed the performance of SCA and SA in all scenarios; however, there is
still no mathematical foundation that supports the convergence of the
$\varepsilon$-SCA.

The so-called Digital Annealer [2, 23, 24, 28, 31], developed by Fujitsu
Laboratories, is a CMOS hardware designed to solve fully connected Quadratic
Unconstrained Binary Optimization (QUBO) problems whose performance has been
compared with other state-of-the-art dedicated solvers in terms of solution
quality (success probabilities) and scalability (time to solution). For an up-
to-date, extensive review of several dedicated hardware solvers for the Ising
model, including Fujitsu’s Digital Annealer benchmarks against multiple well-
known Ising machines, we highly recommend [27]. The first-generation Digital
Annealer’s Algorithm [2], or simply DA, is a physics-inspired algorithm
(motivated by SA) employed in the search for ground states of a Hamiltonian by
relying on a parallel-trial scheme and an escape mechanism (called a _dynamic
offset_) as an attempt to increase spin-flip acceptance probabilities and
avoid local minima. In [2], the DA was benchmarked against SA, PT (PT + ICM)
[37, 38], and PTDA [2], considering sparse and fully connected Ising models:
two-dimensional (with periodic boundary condition) and Sherrington-Kirkpatrick
(SK) spin-glasses with bimodal and Gaussian disorder. For sparse models, the
PT + ICM exhibited the lowest time to solution (TTS) with a clear scaling
advantage; moreover, for sparse bimodal models, the parallel-trial scheme
incorporated by the DA was not enough to decrease the TTS and provide better
scaling than SA. On the other hand, the performance for fully connected models
was significantly better, where the analysis of TTS for DA against SA, PT, and
PTDA revealed a significant and consistent speedup of at least two orders of
magnitude, and the DA achieved higher success rates and a small scaling
advantage over SA. In [17], the authors benchmarked the DA against SA, SCA,
and $\varepsilon$-SCA, considering extensive series of models with different
densities for the spin-spin interactions $J_{i,j}\in\\{-1,0,1\\}$. The authors
obtained higher success probabilities for DA against SA and SCA in all
instances, while $\varepsilon$-SCA outperformed DA except in the cases where
the non-null interactions were dense and anti-ferromagnetic. In [11], DA
clearly outperformed SA, SCA, and $\varepsilon$-SCA in obtaining solutions for
the traveling salesman problem with higher success probability.

Although further iterations of the first-generation DA also rely on replica
exchange methods [28], in this paper, we restrict ourselves to the
mathematical description of the DA in the particular scenario where only the
parallel search technique is regarded. Even though we restrict ourselves to
such a particular case, this simple setting is mathematically interesting from
the Markov chain theory point-of-view and has resulted in significant
practical improvement over SA and other novel simulated annealing algorithms
in several instances [11, 17, 23, 27]. We show that the equilibrium
distribution corresponding to the finite fixed temperature Markov chain
assumes a non-trivial form and does not necessarily coincide with the Gibbs-
Boltzmann distribution. Moreover, based on the techniques developed by Hajek
[15], we obtain a necessary and sufficient condition on the cooling schedule
that guarantees asymptotic convergence of the algorithm to a ground state.

##  2\. Digital Annealer’s Algorithm

This paper’s main problem of interest concerns finding points of global
minima, the so-called ground states, of a Hamiltonian function for a given
Ising model on a graph. Let $G=(V,E)$ be an undirected graph with a finite
vertex set, containing no multi-edges and no self-loops, and let us consider
the set $\Omega=\\{-1,+1\\}^{V}$ of all spin configurations
$\boldsymbol{\sigma}=(\sigma_{x})_{x\in V}$ whose spin values are either $+1$
or $-1$. Given a family of spin-spin coupling constants $(J_{x,y})_{x,y\in V}$
which is symmetric (that is, $J_{x,y}=J_{y,x}$ for every pair $x,y$ in $V$)
and satisfies $J_{x,y}=0$ whenever $\\{x,y\\}\notin E$, and a family of local
external fields $(h_{x})_{x\in V}$, let the Hamiltonian of the Ising model be
the real-valued function $H$ defined on $\Omega$ given by

(2.1) $\displaystyle H(\boldsymbol{\sigma})=-\frac{1}{2}\sum_{x,y\in
V}J_{x,y}\sigma_{x}\sigma_{y}-\sum_{x\in V}h_{x}\sigma_{x}$

for each spin configuration $\boldsymbol{\sigma}=(\sigma_{x})_{x\in V}$. Then,
the set of ground states of $H$ is defined by

(2.2)
$\text{GS}=\\{\boldsymbol{\sigma}\in\Omega:H(\boldsymbol{\sigma})=\min_{\boldsymbol{\tau}\in\Omega}H(\boldsymbol{\tau})\\}.$

Before we proceed, let us introduce some notation. Given a configuration
$\boldsymbol{\sigma}$ and a vertex $x$, let us denote by
$\boldsymbol{\sigma}^{x}$ the spin configuration in $\Omega$ which coincides
with $\boldsymbol{\sigma}$ except for its spin value at $x$, more precisely,
we define

(2.3) $(\boldsymbol{\sigma}^{x})_{y}=\begin{cases}-\sigma_{x}&\text{if
$y=x$,}\\\ \sigma_{y}&\text{otherwise.}\end{cases}$

We define the energetic cost $E_{x}(\boldsymbol{\sigma})$ of flipping the spin
value of $\boldsymbol{\sigma}$ at $x$ by letting
$E_{x}(\boldsymbol{\sigma})=H(\boldsymbol{\sigma}^{x})-H(\boldsymbol{\sigma})$,
and denote by $E_{x}(\boldsymbol{\sigma})^{+}$ its positive part given by
$E_{x}(\boldsymbol{\sigma})^{+}=\max\\{E_{x}(\boldsymbol{\sigma}),0\\}$. Note
that $E_{x}(\boldsymbol{\sigma})$ can be expressed in terms of the cavity
field $\tilde{h}_{x}(\boldsymbol{\sigma})$, whose expression given by
$\tilde{h}_{x}(\boldsymbol{\sigma})=\sum_{y\in V}J_{x,y}\sigma_{y}+h_{x}$, in
the form

(2.4)
$E_{x}(\boldsymbol{\sigma})=2\tilde{h}_{x}(\boldsymbol{\sigma})\sigma_{x}.$

(a) SA

(b) DA

Figure 1. Comparison between the spin updates performed by SA and DA. We
illustrate the spin update rules for these algorithms by decomposing them into
two steps: the “Random select” phase and the “Flip?” phase. To help the reader
visualize their operation, we represent the spin values $+1$ and $-1$ in red
and blue, respectively, which, after being evaluated in the “Flip?” phase, are
represented in transparent colors in case their new values are the same as in
the original configuration.

The Digital Annealer’s Algorithm works as follows. Given a state
$X_{k-1}=\boldsymbol{\sigma}$ at time $k-1$, we propose a parallel-trial where
each vertex $x$ is assigned as eligible to flip its spin value according to
the Metropolis criterion

(2.5)
$\exp(-\beta_{k}E_{x}(\boldsymbol{\sigma})^{+})=\min\\{\exp(-\beta_{k}E_{x}(\boldsymbol{\sigma})),1\\}$

corresponding to temperature $1/{\beta_{k}}$. If the set $S$ of all eligible
vertices contains at least one element, then a vertex $x$ is chosen uniformly
at random from $S$, and we place $X_{k}=\boldsymbol{\sigma}^{x}$; otherwise,
nothing is done, and we consider $X_{k}=\boldsymbol{\sigma}$. Note that,
differently from [23, 24], we are disregarding the escape mechanism so as to
guarantee the Markov property holds. Although the SA and DA can flip at most
one spin per update, they differ greatly in how they update spin
configurations; see Figure 1. In Figure 1(a), given a spin configuration, we
select one spin uniformly at random and decide whether we flip it or not
according to Metropolis criterion (2.5). In the end, the spin value is updated
while the other spins remain unchanged. In Figure 1(b), differently from the
SA, we first apply the “Flip?” phase, where we propose parallel spin-flips
according to the Metropolis criterion (2.5). Note that the only spins
evaluated in the “Random select” phase are those represented in opaque colors,
which correspond to the elements of the set $S$ of spins that accepted spin-
flips, and we select only one of them, uniformly at random, to replace its
corresponding value in the initial configuration.

According to practical simulations [2, 11, 17, 24, 31], it is expected that,
by properly decreasing the temperature toward zero at each step, the updated
configuration approaches a ground state. Therefore, a good approximation for a
ground state can be achieved, provided we run the algorithm for a sufficiently
large number of steps. In this paper, we prove Theorem 2.1 to provide a
mathematical foundation for these observations.

In our framework, the Digital Annealer’s Algorithm transition matrix
$P_{\beta}^{\text{DA}}$ at inverse temperature $\beta\in(0,\infty)$ is defined
by

(2.6)
$P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\tau})=\begin{cases}\displaystyle\sum_{\begin{subarray}{c}S\subseteq
V\\\ S\ni x\end{subarray}}\frac{1}{|S|}\displaystyle\prod_{y\in S}e^{-\beta
E_{y}(\boldsymbol{\sigma})^{+}}\displaystyle\prod_{y\in V\backslash
S}(1-e^{-\beta E_{y}(\boldsymbol{\sigma})^{+}})&\text{if
$\boldsymbol{\tau}=\boldsymbol{\sigma}^{x}$ for some $x$ in $V$,}\\\\[20.0pt]
\displaystyle\prod_{y\in V}(1-e^{-\beta
E_{y}(\boldsymbol{\sigma})^{+}})&\text{if
$\boldsymbol{\tau}=\boldsymbol{\sigma}$, and }\\\\[20.0pt]
0&\text{otherwise}.\end{cases}$

First, let us note that the spin-flip probability
$P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})$ can be
decomposed in the form

(2.7)
$P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=\frac{1}{|V|}e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}}+\sum_{\begin{subarray}{c}S\subseteq V\\\ S\ni
x\end{subarray}}\left(\frac{1}{|S|}-\frac{1}{|V|}\right)\displaystyle\prod_{y\in
S}e^{-\beta E_{y}(\boldsymbol{\sigma})^{+}}\displaystyle\prod_{y\in
V\backslash S}(1-e^{-\beta E_{y}(\boldsymbol{\sigma})^{+}}),$

where the first term in the right-hand side of the equation (2.7) coincides
with the standard single site spin-flip Metropolis dynamics transition
probability from $\boldsymbol{\sigma}$ to $\boldsymbol{\sigma}^{x}$.
Therefore, as discussed in [2], the DA, in fact, provides us with a larger
spin-flip rate compared to the Metropolis dynamics, which is generally
attributed as a decisive property that can lead to greater effectiveness over
SA.

The main result of this paper, Theorem 2.1 stated below, provides us with a
necessary and sufficient condition on the speed of convergence of the
temperature $\beta_{k}^{-1}$ to zero so that the algorithm converges
asymptotically to a ground state. Its proof is based on the techniques
developed in [15], and it is presented in Section 4.

######  Theorem 2.1.

Let $(\beta_{k})_{k\geq 1}$ be a cooling schedule, that is, a nondecreasing
sequence of positive real numbers such that
$\lim_{k\to\infty}\beta_{k}=+\infty$, and let $(X_{k})_{k\geq 0}$ be the
discrete-time inhomogeneous Markov chain satisfying

(2.8)
$\mathbb{P}(X_{k}=\boldsymbol{\sigma}_{k}|X_{k-1}=\boldsymbol{\sigma}_{k-1},\dots,X_{0}=\boldsymbol{\sigma}_{0})=\mathbb{P}(X_{k}=\boldsymbol{\sigma}_{k}|X_{k-1}=\boldsymbol{\sigma}_{k-1})=P_{\beta_{k}}^{\text{DA}}(\boldsymbol{\sigma}_{k-1},\boldsymbol{\sigma}_{k})$

for every positive integer $k$ and
$\boldsymbol{\sigma}_{0},\dots,\boldsymbol{\sigma}_{k}$ in $\Omega$. Then,
there exists a constant $\gamma^{\ast}$ such that the limit

(2.9) $\lim_{k\to\infty}\mathbb{P}(X_{k}\in\text{GS})=1$

holds if and only if

(2.10) $\sum_{k=1}^{\infty}e^{-\beta_{k}\gamma^{\ast}}=+\infty.$

As we will precise more and prove in Section 4, the constant $\gamma^{\ast}$
from Theorem 2.1 coincides with the depth of the deepest local minimum that is
not a ground state. Note that, in particular, if we consider a logarithmic
cooling schedule $(\beta_{k})_{k\geq 1}$ in the form

(2.11) $\beta_{k}=\frac{1}{\gamma}\log(k)$

for some positive parameter $\gamma$, then, equation (2.9) holds if and only
if $\gamma\geq\gamma^{\ast}$. Let us note that since the necessary and
sufficient condition (2.10) is the same as the one proven for the SA in [15],
this result does not imply that the DA performs better than SA but only
guarantees its convergence. In one possible attempt to provide rigorous
statements about the superiority of one algorithm over another, it would be
necessary to weaken the notion of convergence and derive results considering
cooling schedules for finite-time simulations. Since each algorithm has its
strengths and weaknesses and its performance is problem-dependent [27], such
results are still challenging to derive and are out of the scope of the
current state of our investigations.

##  3\. Stationary distributions

The standard method used for proving the convergence of simulated annealing
algorithms relies on obtaining the properties of weak and strong ergodicity
for inhomogeneous Markov chains and on the analysis of their stationary
distributions at fixed temperatures, see [1, 3, 16, 30]. A mathematical
foundation for the convergence of SA can be found in [1, 3]. Following the
same idea, the convergence of an algorithm based on probabilistic cellular
automata, namely the SCA, was obtained in [9]. In this section, we show that,
differently from the Metropolis dynamics and the SCA, the stationary
distribution for the DA transition probability does not necessarily assume a
general formula or coincide with the Gibbs distribution. For that reason, in
Section 4 we adopt a different approach to prove convergence which does not
depend upon the knowledge of the stationary distribution.

First, let us note that $P_{\beta}^{\text{DA}}$ is an irreducible and
aperiodic transition matrix. Thus, the discrete-time homogeneous Markov chain
$(X_{k})_{k\geq 0}$ determined by $P_{\beta}^{\text{DA}}$ converges to its
(unique) stationary distribution $\pi_{\beta}^{\text{DA}}$. In this section,
we provide some examples where it is possible to determine the stationary
distribution $\pi_{\beta}^{\text{DA}}$ of $P_{\beta}^{\text{DA}}$ and show
that, for some particular cases, it differs from the Gibbs distribution
$\pi_{\beta}^{\text{G}}$ defined by

(3.1) $\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})=\frac{e^{-\beta
H(\boldsymbol{\sigma})}}{Z_{\beta}}$

for each configuration $\boldsymbol{\sigma}$ in $\Omega$, where the
normalizing factor $Z_{\beta}=\sum_{\boldsymbol{\tau}\in\Omega}e^{-\beta
H(\boldsymbol{\tau})}$ is known as the partition function.

######  Example 3.1.

Let us consider the simplest case where $V$ consists of only two vertices and
the $H$ is the Hamiltonian of a ferromagnetic Ising model without external
fields with pairwise interaction equal $J\geq 0$. Let
$\pi_{\beta}^{\text{DA}}=(\pi_{\beta}^{\text{DA}}({\uparrow\uparrow}),\pi_{\beta}^{\text{DA}}({\uparrow\downarrow}),\pi_{\beta}^{\text{DA}}({\downarrow\uparrow}),\pi_{\beta}^{\text{DA}}({\downarrow\downarrow}))$
be the stationary distribution of the transition matrix
$P_{\beta}^{\text{DA}}$, which is expressed as the $4\times 4$ matrix given by

(3.2) $P_{\beta}^{\text{DA}}=\begin{pmatrix}(1-e^{-2\beta J})^{2}&e^{-2\beta
J}(1-\frac{1}{2}e^{-2\beta J})&e^{-2\beta J}(1-\frac{1}{2}e^{-2\beta
J})&0\\\\[5.0pt] \frac{1}{2}&0&0&\frac{1}{2}\\\\[5.0pt]
\frac{1}{2}&0&0&\frac{1}{2}\\\\[5.0pt] 0&e^{-2\beta J}(1-\frac{1}{2}e^{-2\beta
J})&e^{-2\beta J}(1-\frac{1}{2}e^{-2\beta J})&(1-e^{-2\beta J})^{2}\\\
\end{pmatrix}.$

It follows that
$\pi_{\beta}^{\text{DA}}({\uparrow\uparrow})=\pi_{\beta}^{\text{DA}}({\downarrow\downarrow})$,
$\pi_{\beta}^{\text{DA}}({\uparrow\downarrow})=\pi_{\beta}^{\text{DA}}({\downarrow\uparrow})$,
and

(3.3) $\pi_{\beta}^{\text{DA}}({\uparrow\uparrow})=\frac{e^{-\beta
H(\uparrow\uparrow)}}{Z_{\beta}+2e^{-\beta J}(1-e^{-2\beta J})}$

and

(3.4) $\pi_{\beta}^{\text{DA}}({\uparrow\downarrow})=\frac{e^{-\beta
H(\uparrow\downarrow)}+e^{-\beta J}(1-e^{-2\beta J})}{Z_{\beta}+2e^{-\beta
J}(1-e^{-2\beta J})}.$

It is straightforward to show that $\pi_{\beta}^{\text{DA}}$ coincides with
$\pi_{\beta}^{\text{G}}$ if and only if $J=0$, moreover,
$\pi_{\beta}^{\text{DA}}$ converges to the uniform distribution concentrated
on the ground states $\uparrow\uparrow$ and $\downarrow\downarrow$ as $\beta$
tends to infinity.

It follows from equation (2.6) that the spin-flip probability
$P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})$ can also
be expressed as

(3.5)
$P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}},$

where $R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})$ is defined by

(3.6)
$R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=\sum_{\begin{subarray}{c}S\subseteq
V\setminus\\{x\\}\end{subarray}}\frac{1}{|S|+1}\displaystyle\prod_{y\in
S}e^{-\beta
E_{y}(\boldsymbol{\sigma})^{+}}\displaystyle\prod_{y\in(V\setminus\\{x\\})\setminus
S}(1-e^{-\beta E_{y}(\boldsymbol{\sigma})^{+}}).$

######  Proposition 3.2.

Let $\pi_{\beta}^{\text{G}}$ be the Gibbs distribution on $\Omega$ at inverse
temperature $\beta$. Then, it follows that

(3.7)
$(\pi_{\beta}^{\text{G}}P_{\beta}^{\text{DA}})(\boldsymbol{\sigma})-\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})=\sum_{x\in
V}e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}}(R(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})-R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x}))$

holds for every $\boldsymbol{\sigma}$ in $\Omega$.

###### Proof.

Let $P_{\beta}^{\text{G}}$ denote the single site spin-flip Metropolis
dynamics transition matrix at inverse temperature $\beta$ whose transition
probability from $\boldsymbol{\sigma}$ to $\boldsymbol{\sigma}^{x}$ is given
by
$P_{\beta}^{\text{G}}(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=\frac{1}{|V|}e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}}$. Let us begin by noting that identities

(3.8)
$(\pi_{\beta}^{\text{G}}P_{\beta}^{\text{DA}})(\boldsymbol{\sigma})=\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma})+\sum_{x\in
V}\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x})P_{\beta}^{\text{DA}}(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})$

and

(3.9)
$\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})=(\pi_{\beta}^{\text{G}}P_{\beta}^{\text{G}})(\boldsymbol{\sigma})=\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})P_{\beta}^{\text{G}}(\boldsymbol{\sigma},\boldsymbol{\sigma})+\sum_{x\in
V}\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x})P_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})$

imply

(3.10)
$\displaystyle(\pi_{\beta}^{\text{G}}P_{\beta}^{\text{DA}})(\boldsymbol{\sigma})-\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})=$
$\displaystyle\;\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})\left(P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma})-P_{\beta}^{\text{G}}(\boldsymbol{\sigma},\boldsymbol{\sigma})\right)$
$\displaystyle+\sum_{x\in
V}\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x})\left(P_{\beta}^{\text{DA}}(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})-P_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})\right).$

Now, by summing both sides of equation (3.5) over all vertices of $V$, we
obtain

(3.11)
$P_{\beta}^{\text{G}}(\boldsymbol{\sigma},\boldsymbol{\sigma})-P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma})=\sum_{x\in
V}\left(R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})-\frac{1}{|V|}\right)e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}},$

and, by exchanging the roles of $\boldsymbol{\sigma}$ and
$\boldsymbol{\sigma}^{x}$ in equation (3.5), we have

(3.12)
$P_{\beta}^{\text{DA}}(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})-P_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})=\left(R(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})-\frac{1}{|V|}\right)e^{-\beta
E_{x}(\boldsymbol{\sigma}^{x})^{+}}.$

Therefore, since the Gibbs distribution satisfies the detailed balance
equations for the Metropolis dynamics, i.e.,
$\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}}=\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma}^{x})e^{-\beta
E_{x}(\boldsymbol{\sigma}^{x})^{+}}$, we derive the identity

(3.13)
$\pi_{\beta}^{\text{G}}P_{\beta}^{\text{DA}}(\boldsymbol{\sigma})-\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})=\pi_{\beta}^{\text{G}}(\boldsymbol{\sigma})\left[P_{\beta}^{\text{DA}}(\boldsymbol{\sigma},\boldsymbol{\sigma})-P_{\beta}^{\text{G}}(\boldsymbol{\sigma},\boldsymbol{\sigma})+\sum_{x\in
V}\left(R(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})-\frac{1}{|V|}\right)e^{-\beta
E_{x}(\boldsymbol{\sigma})^{+}}\right],$

which combined with equation (3.11), implies equation (3.7). ∎

######  Example 3.3.

Let us consider a Hamiltonian which is free of pairwise interactions, subject
only to external fields, that is, the Hamiltonian $H$ is given in the form

(3.14) $H(\boldsymbol{\sigma})=-\sum_{x\in V}h_{x}\sigma_{x}.$

Note that, given a spin configuration $\boldsymbol{\sigma}$, the identity
$E_{y}(\boldsymbol{\sigma})=E_{y}(\boldsymbol{\sigma}^{x})$ holds whenever $x$
and $y$ are distinct vertices of $V$, therefore, we have
$R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=R(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})$
for each vertex $x$. It follows from Proposition 3.2 that
$\pi_{\beta}^{\text{G}}P_{\beta}^{\text{DA}}=\pi_{\beta}^{\text{G}}$, and
consequently, $\pi_{\beta}^{\text{DA}}=\pi_{\beta}^{\text{G}}$.

In the following proposition, we provide a generalization of Example 3.1.

######  Proposition 3.4.

Let $H$ be the Hamiltonian of a ferromagnetic Ising model without external
fields, which is given by

(3.15) $H(\boldsymbol{\sigma})=-\frac{1}{2}\sum_{x,y\in
V}J_{x,y}\sigma_{x}\sigma_{y},$

where each pairwise interaction $J_{x,y}$ is a nonnegative real number. It
follows that $\pi_{\beta}^{\text{DA}}=\pi_{\beta}^{\text{G}}$ if and only if
$J_{x,y}=0$ for all $x,y$ in $V$.

###### Proof.

Suppose that there is a pair $v,w$ of vertices in $V$ such that $J_{v,w}>0$.
Let us consider the particular case where $\boldsymbol{\sigma}$ is a ground
state of $H$, which is a configuration whose spin values are all the same.
Then, the identity

(3.16) $E_{y}(\boldsymbol{\sigma}^{x})=E_{y}(\boldsymbol{\sigma})-4J_{x,y}$

holds for any distinct points $x$ and $y$ in $V$. Corresponding to a fixed
vertex $x$ in $V$, if we are given the families $(X_{i})_{i\in V}$ and
$(Y_{i})_{i\in V}$ of independent Bernoulli random variables satisfying
$\mathbb{P}(X_{i}=1)=1-\mathbb{P}(X_{i}=0)=e^{-\beta
E_{i}(\boldsymbol{\sigma})^{+}}$ and
$\mathbb{P}(Y_{i}=1)=1-\mathbb{P}(Y_{i}=0)=e^{-\beta
E_{i}(\boldsymbol{\sigma}^{x})^{+}}$, it is straightforward to show that

(3.17)
$R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=\mathbb{E}\left[\frac{1}{\sum_{i\in
V}X_{i}}\Bigg{|}X_{x}=1\right]\geq\mathbb{E}\left[\frac{1}{\sum_{i\in
V}Y_{i}}\Bigg{|}Y_{x}=1\right]=R(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma}).$

In particular, if we assume that $x=v$, the inequality from equation (3.17) is
strict, therefore, we conclude from Proposition 3.2 that
$\pi_{\beta}^{\text{DA}}\neq\pi_{\beta}^{\text{G}}$. Conversely, suppose all
pairwise interactions are null. In that case, it follows that
$R(\boldsymbol{\sigma},\boldsymbol{\sigma}^{x})=R(\boldsymbol{\sigma}^{x},\boldsymbol{\sigma})$
for every configuration $\boldsymbol{\sigma}$ and vertex $x$, thus,
$\pi_{\beta}^{\text{G}}P_{\beta}^{\text{DA}}=\pi_{\beta}^{\text{G}}$. ∎

##  4\. Proof of Theorem 2.1

In this section, we derive Theorem 2.1 as a particular case of [15, Theorem
2]. In order to do so, let us start by formulating our problem in a more
general setting. Let us consider a nonempty finite state space $\Omega$, an
energy function $H:\Omega\to\mathbb{R}$, and a system
$\mathcal{N}=\\{N(\boldsymbol{\sigma}):\boldsymbol{\sigma}\in\Omega\\}$ of
neighbors of states. Recall that we are particularly interested in the set of
ground states (global minima) of $H$ defined as

(4.1)
$\text{GS}=\\{\boldsymbol{\sigma}\in\Omega:H(\boldsymbol{\sigma})=\min_{\boldsymbol{\tau}\in\Omega}H(\boldsymbol{\tau})\\}.$

In this setting, the property of irreducibility can be defined as follows.

######  Definition 4.1 (Irreducibility).

The pair $(\Omega,\mathcal{N})$ is said to be irreducible if for any distinct
states $\boldsymbol{\sigma}$ and $\boldsymbol{\tau}$ there is a sequence
$\boldsymbol{\sigma}_{0}=\boldsymbol{\sigma}$, $\boldsymbol{\sigma}_{1}$,
$\dots$, $\boldsymbol{\sigma}_{n}=\boldsymbol{\tau}$ in $\Omega$, where $n\geq
1$, such that $\boldsymbol{\sigma}_{k}\in N(\boldsymbol{\sigma}_{k-1})$
whenever $1\leq k\leq n$.

In addition to the property of irreducibility, it will be necessary to assume
weak reversibility. Before introducing such a notion, let us define the idea
of reachability. For any pair $\boldsymbol{\sigma},\boldsymbol{\tau}$ of
points in $\Omega$, $\boldsymbol{\tau}$ is said to be reachable from
$\boldsymbol{\sigma}$ at height $E$ if either
$\boldsymbol{\sigma}=\boldsymbol{\tau}$ and $H(\boldsymbol{\sigma})\leq E$, or
$\boldsymbol{\sigma}\neq\boldsymbol{\tau}$ and there exists a sequence of
states $\boldsymbol{\sigma}_{0}=\boldsymbol{\sigma}$,
$\boldsymbol{\sigma}_{1}$, $\dots$,
$\boldsymbol{\sigma}_{n}=\boldsymbol{\tau}$, where $n\geq 1$, such that
$\max_{0\leq i\leq n}H(\boldsymbol{\sigma}_{i})\leq E$ and
$\boldsymbol{\sigma}_{k}\in N(\boldsymbol{\sigma}_{k-1})$ for each $k$ such
that $1\leq k\leq n$.

######  Definition 4.2 (Weak reversibility).

The triple $(\Omega,H,\mathcal{N})$ is said to be weakly reversible if for any
states $\boldsymbol{\sigma}$ and $\boldsymbol{\tau}$ and any real number $E$,
$\boldsymbol{\sigma}$ is reachable from $\boldsymbol{\tau}$ at height $E$ if
and only if $\boldsymbol{\tau}$ is reachable from $\boldsymbol{\sigma}$ at
height $E$.

From now on, let us assume the properties of irreducibility and weak
reversibility. In order to state Theorem 4.3, it is necessary to introduce the
notions of a cup and local minima. A subset $C$ of $\Omega$ is called a cup if
there exists some real number $E$ such that for any $\boldsymbol{\sigma}\in
C$, the set $C$ can be expressed as the set of all states that are reachable
from $\boldsymbol{\sigma}$ at height $E$. Then, given a cup $C$, we define its
boundary by

(4.2) $\partial C=\\{\boldsymbol{\tau}\notin C:\text{$\boldsymbol{\tau}\in
N(\boldsymbol{\sigma})$ for some $\boldsymbol{\sigma}$ in $C$}\\},$

its bottom by

(4.3) $B(C)=\\{\boldsymbol{\sigma}\in
C:H(\boldsymbol{\sigma})=\min_{\boldsymbol{\tau}\in
C}H(\boldsymbol{\tau})\\},$

and its depth by

(4.4) $d(C)=\min_{\boldsymbol{\sigma}\in\partial
C}H(\boldsymbol{\sigma})-\min_{\boldsymbol{\sigma}\in
C}H(\boldsymbol{\sigma}).$

Furthermore, we say that $\boldsymbol{\sigma}$ is a local minimum if there is
no $\boldsymbol{\tau}$ satisfying
$H(\boldsymbol{\tau})<H(\boldsymbol{\sigma})$ which is reachable from
$\boldsymbol{\sigma}$ at height $H(\boldsymbol{\sigma})$. So, the depth of a
local minimum $\boldsymbol{\sigma}$ which is not a ground state is defined as
the smallest positive real number $E$ so that there exists $\boldsymbol{\tau}$
with $H(\boldsymbol{\tau})<H(\boldsymbol{\sigma})$ that is reachable from
$\boldsymbol{\sigma}$ at height $H(\boldsymbol{\sigma})+E$; otherwise, in case
$\boldsymbol{\sigma}$ is a ground state, its depth is defined as infinity. It
is straightforward to show that every local minimum of depth $d$ is contained
in the bottom of a cup of depth $d$.

Although the motivation of the problem proposed by B. Hajek in [15] was
formulated with the Metropolis dynamics in mind, his main result was proven in
a rather more general setting. Hajek considered a continuous-time process
$(W_{t})_{t\in[0,\infty)}$ defined by letting $W_{t}=X_{k}$ whenever
$U_{k}\leq t<U_{k+1}$, derived from a discrete-time homogeneous Markov process
$(U_{k},X_{k})_{k\geq 0}$ on the state space $[0,\infty)\times\Omega$ whose
transition probabilities assume the form

(4.5) $\mathbb{P}(U_{k+1}\geq
u,X_{k+1}=\boldsymbol{\tau}|\text{$(U_{i},X_{i})=(u_{i},\boldsymbol{\sigma}_{i})$
for all $0\leq i\leq
k$})=\int_{u}^{\infty}Q(\boldsymbol{\sigma}_{k},\boldsymbol{\tau},t,\lambda_{t})\Phi(dt,u_{k},\boldsymbol{\sigma}_{k}).$

In the equation above, we suppose that $(\lambda_{t})_{t\in[0,\infty)}$ is a
nonincreasing family of numbers in the interval $(0,1)$ such that
$\lim_{t\to\infty}\lambda_{t}=0$. Moreover, we assume the existence of
positive constants $a,D,c_{1}$, and $c_{2}$ such that

  * •

$\mathbb{E}[\min\\{k:U_{k}\geq
t+a\\}|U_{0}=t,X_{0}=\boldsymbol{\sigma}]\leq\frac{1}{a}$ holds for all $t\geq
0$ and $\boldsymbol{\sigma}\in\Omega$.

  * •

$\Phi(\cdot,s,\boldsymbol{\sigma})$ is a probability distribution function
such that

(4.6) $\Phi(s^{-},s,\boldsymbol{\sigma})=0$

and

(4.7) $\int_{0}^{\infty}(t-s)\Phi(dt,s,\boldsymbol{\sigma})\leq D.$

  * •

$Q(\cdot,\cdot,t,\lambda)$ is an stochastic matrix such that the conditions

(4.8) $c_{1}\lambda^{[H(\boldsymbol{\tau})-H(\boldsymbol{\sigma})]^{+}}\leq
Q(\boldsymbol{\sigma},\boldsymbol{\tau},t,\lambda)\leq
c_{2}\lambda^{[H(\boldsymbol{\tau})-H(\boldsymbol{\sigma})]^{+}}$ for
$\boldsymbol{\tau}\in N(\boldsymbol{\sigma})$

and

(4.9) $Q(\boldsymbol{\sigma},\boldsymbol{\tau},t,\lambda)=0$ for
$\boldsymbol{\tau}\notin N(\boldsymbol{\sigma})\cup\\{\boldsymbol{\sigma}\\}$

hold whenever $t\geq 0$ and $0\leq\lambda\leq 1$.

######  Theorem 4.3 ([15]).

Let us assume that $(\Omega,H,\mathcal{N})$ is irreducible and weakly
reversible, and $(W_{t})_{t\in[0,\infty)}$ is the continuous-time process
defined as above. Then, the following conditions hold.

  1. (a)

Every state $\boldsymbol{\sigma}$ which is not a local minimum satisfies

(4.10) $\lim_{t\to\infty}\mathbb{P}(W_{t}=\boldsymbol{\sigma})=0.$

  2. (b)

Given a cup $C$ of depth $d$ whose elements of its bottom are local minima of
depth $d$, then

(4.11) $\lim_{t\to\infty}\mathbb{P}(W_{t}\in B(C))=0$

if and only if

(4.12) $\int_{0}^{\infty}\lambda_{t}^{d}dt=+\infty.$

  3. (c)

If we define $\gamma^{\ast}$ as the depth of the second deepest local minimum,
then the limit

(4.13) $\lim_{t\to\infty}\mathbb{P}(W_{t}\in\text{GS})=1$

holds if and only if

(4.14) $\int_{0}^{\infty}\lambda_{t}^{\gamma^{*}}dt=+\infty.$

###### Proof of Theorem 2.1 assuming Theorem 4.3.

Let us show that the setting assumed in Theorem 2.1 fits as a particular case
of Theorem 4.3. In the following, let us consider the triple
$(\Omega,H,\mathcal{N})$, where $\Omega=\\{-1,+1\\}^{V}$ is the set of all
Ising spin configurations on the vertex set $V$ of a graph, $H$ is the
Hamiltonian given by equation (2.1), and
$\mathcal{N}=\\{N(\boldsymbol{\sigma}):\boldsymbol{\sigma}\in\Omega\\}$ is the
system of neighbors of spin configurations defined as follows. Given a
configuration $\boldsymbol{\sigma}$ in $\Omega$, let $N(\boldsymbol{\sigma})$
be the set of neighbors of $\boldsymbol{\sigma}$ defined by

(4.15) $N(\boldsymbol{\sigma})=\\{\boldsymbol{\sigma}^{x}:x\in V\\}.$

It is straightforward to show that the triple $(\Omega,H,\mathcal{N})$ is
irreducible and weakly reversible.

Let us consider the discrete-time homogeneous Markov chain
$(U_{k},X_{k})_{k\geq 0}$ on the state space $[0,\infty)\times\Omega$ whose
transition probabilities assume the form (4.5), with the transition matrix
$Q(\cdot,\cdot,t,\lambda)$, the probability distribution function
$\Phi(\cdot,s,\boldsymbol{\sigma})$ and the family of real numbers
$(\lambda_{t})_{t\in[0,\infty)}$ chosen by letting
$Q(\cdot,\cdot,t,\lambda)=P_{\log\lambda^{-1}}^{\text{DA}}$,
$\Phi(t,s,\boldsymbol{\sigma})=\mathbbm{1}_{\\{t\geq s+1\\}}$, and

$\lambda_{t}=\begin{cases}\exp(-\beta_{\lfloor t\rfloor})&\text{if $t\geq
1$,}\\\ \exp(-\beta_{1})&\text{otherwise;}\end{cases}$

where $\lfloor t\rfloor$ stands for the greatest integer less than or equal to
$t$. Note that equations (2.6) and (2.7) imply that

(4.16)
$\frac{1}{|V|}\lambda^{[H(\boldsymbol{\tau})-H(\boldsymbol{\sigma})]^{+}}\leq
Q(\boldsymbol{\sigma},\boldsymbol{\tau},t,\lambda)\leq\lambda^{[H(\boldsymbol{\tau})-H(\boldsymbol{\sigma})]^{+}}$

holds for each $\boldsymbol{\tau}$ in $N(\boldsymbol{\sigma})$, and
$Q(\boldsymbol{\sigma},\boldsymbol{\tau},t,\lambda)=0$ whenever
$\boldsymbol{\tau}$ does not belong to
$N(\boldsymbol{\sigma})\cup\\{\boldsymbol{\sigma}\\}$. Furthermore, equation
(4.5) is expressed as

(4.17) $\mathbb{P}(U_{k+1}\geq
u,X_{k+1}=\boldsymbol{\tau}|\text{$(U_{i},X_{i})=(u_{i},\boldsymbol{\sigma}_{i})$
for all $0\leq i\leq k$})=P_{\beta_{\lfloor
u_{k}+1\rfloor}}^{\text{DA}}(\boldsymbol{\sigma}_{k},\boldsymbol{\tau})\mathbbm{1}_{\\{u_{k}+1\geq
u\\}},$

which implies that

(4.18)
$\mathbb{P}(U_{k+1}=u_{k}+1,X_{k+1}=\boldsymbol{\tau}|\text{$(U_{i},X_{i})=(u_{i},\boldsymbol{\sigma}_{i})$
for all $0\leq i\leq k$})=P_{\beta_{\lfloor
u_{k}+1\rfloor}}^{\text{DA}}(\boldsymbol{\sigma}_{k},\boldsymbol{\tau}).$

In the particular case where $U_{0}=0$, it is straightforward to show that
$\mathbb{P}(U_{k}=k)=1$ for each $k\geq 0$ and the one-step transition
probabilities of $(X_{k})_{k\geq 0}$ take the form (2.8). Since the hypotheses
of Theorem 4.3 are applicable to this case, then, by using the continuous-time
stochastic process $(W_{t})_{t\in[0,\infty)}$ evaluated on integer times and
Theorem 4.3(c), the result follows. ∎

##  5\. Concluding remarks

First, let us recall that the general algorithm introduced in [2] relies on
the implementation of the so-called dynamic offset in order to encourage the
system to flip more spins and prevent it from getting trapped in a local
minimum. Thus, in order to prove statements regarding the convergence to the
ground states for such an algorithm, one should rely on a different approach
since the Markov property would no longer be valid.

On the other hand, note that in the proof of Theorem 4.3 the assumption that
$H$ is given as in equation (2.1) was not necessary, therefore, we do not have
to restrict ourselves only to QUBO problems and Theorem 2.1 can also be
applied to search for ground states of Hamiltonians other than those written
in this form. Furthermore, recall that, as discussed in Section 2, we dealt
with results regarding the asymptotic convergence to the ground states, and,
in a realistic scenario, especially when dealing with problems that involve
thousands of variables, one does not have infinite time to perform a
simulation. Relying on ideas from Freidlin and Wentzel [8], O. Catoni [4, 5,
7] obtained rigorous results regarding finite time simulation, large deviation
principles, and optimal cooling schedules for the Metropolis algorithm and
Markov chains with rare transitions. Investigations on such topics
specifically for the DA and other algorithms (such as those discussed in [9]),
including rigorous results concerning the effectiveness of one algorithm over
the other, are still in progress.

Another topic that is also worth investigating is the behavior of the
corresponding time-homogeneous Markov chain at a fixed temperature in order to
provide a complete characterization of its stationary distribution and obtain
estimates of its mixing time. As the reader can verify, providing an upper
bound for the mixing time at high temperatures such as in [9, 11], where it
was possible to show that the mixing times for certain probabilistic cellular
automata is upper bounded by a quantity proportional to $\log|V|$, is not
straightforward. Hence, there are still interesting directions to be explored.

## Declarations

This work was supported by JST CREST Grant Number PJ22180021, Japan. Data
sharing is not applicable to this article as no datasets were generated or
analyzed during the current study. The authors have no relevant financial or
non-financial interests to disclose.

## Acknowledgements

We are grateful to the following members for their continual encouragement and
stimulating discussions: Masato Motomura and Kazushi Kawamura from the Tokyo
Institute of Technology; Hiroshi Teramoto from Kansai University; Masamitsu
Aoki from the Graduate School of Mathematics at Hokkaido University. The
authors would like to thank Thiago Raszeja, Olena Karpel, Dominik Kwietniak,
Wioletta Ruszel, Cristian Spitoni, Ross Kang, and Lars Fritz for the
hospitality and for providing nice opportunities for discussion during Bruno
Hideki Fukushima Kimura’s visiting period at the AGH University of Science and
Technology and Utrecht University.

## References

  * [1] E. Aarts and J. Korst.  Simulated annealing and Boltzmann machines: A stochastic approach to combinatorial optimization and neural computing.  John Wiley & Sons, Inc., 1989.
  * [2] M. Aramon, G. Rosenberg, E. Valiante, T. Miyazawa, H. Tamura, and H.G. Katzgraber.  Physics-inspired optimization for quadratic unconstrained problems using a digital annealer.  Frontiers in Physics, 7:48, 2019.
  * [3] P. Brémaud.  Markov Chains: Gibbs Fields, Monte Carlo Simulation, and Queues.  Texts in Applied Mathematics. Springer New York, 2001.
  * [4] O. Catoni.  Sharp large deviations estimates for simulated annealing algorithms.  Annales De L Institut Henri Poincare-probabilites Et Statistiques, 27:291–383, 1991.
  * [5] O. Catoni.  Rough large deviation estimates for simulated annealing: Application to exponential schedules.  Annals of Probability, 20:1109–1146, 1992.
  * [6] V. Černý.  Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm.  Journal of Optimization Theory and Applications, 45(1):41–51, January 1985.
  * [7] C. Cot and O. Catoni.  Piecewise constant triangular cooling schedules for generalized simulated annealing algorithms.  Annals of Applied Probability, 8:375–396, 1998.
  * [8] M.I. Freidlin, J. Szucs, and A.D. Wentzell.  Random Perturbations of Dynamical Systems.  Grundlehren der mathematischen Wissenschaften. Springer New York, 2012\.
  * [9] B.H. Fukushima-Kimura, S. Handa, K. Kamakura, Y. Kamijima, K. Kawamura, and A. Sakai.  Mixing time and simulated annealing for the stochastic cellular automata.  Journal of Statistical Physics, 190(4):79, 2023.
  * [10] B.H. Fukushima-Kimura, Y. Kamijima, K. Kawamura, and A. Sakai.  Stochastic optimization via parallel dynamics: rigorous results and simulations.  Proceedings of the ISCIE International Symposium on Stochastic Systems Theory and its Applications, 2022:65–71, 2022.
  * [11] B.H. Fukushima-Kimura, Y. Kamijima, K. Kawamura, and A. Sakai.  Stochastic optimization - Glauber dynamics versus stochastic cellular automata.  Transactions of the Institute of Systems, Control and Information Engineers, 36(1):9–16, 2023.
  * [12] S.B. Gelfand and S.K. Mitter.  Analysis of simulated annealing for optimization.  In 1985 24th IEEE Conference on Decision and Control, pages 779–786, 1985.
  * [13] H. Goto.  Quantum computation based on quantum adiabatic bifurcations of Kerr-nonlinear parametric oscillators.  Journal of the Physical Society of Japan, 88(6):061015, 2019.
  * [14] H. Goto, K. Tatsumura, and A.R. Dixon.  Combinatorial optimization by simulating adiabatic bifurcations in nonlinear Hamiltonian systems.  Science Advances, 5(4):eaav2372, 2019.
  * [15] B. Hajek.  Cooling schedules for optimal annealing.  Mathematics of Operations Research, 13(2):311–329, 1988.
  * [16] D.L. Isaacson and R.W. Madsen.  Markov Chains: Theory and Applications.  Wiley Series in Probability and Statistics. Wiley, 1976.
  * [17] K. Kawamura, J. Yu, D. Okonogi, S. Jimbo, G. Inoue, A. Hyodo, A.L. García-Arias, K. Ando, B.H. Fukushima-Kimura, R. Yasudo, T. Van Chu, and M. Motomura.  Amorphica: 4-replica 512 fully connected spin 336MHz metamorphic annealer with programmable optimization strategy and compressed-spin-transfer multi-chip extension.  In 2023 IEEE International Solid-State Circuits Conference (ISSCC), pages 42–44, 2023.
  * [18] S. Kirkpatrick, C. D. Gelatt, and M. P. Vecchi.  Optimization by simulated annealing.  Science, 220(4598):671–680, 1983.
  * [19] E.L. Lawler.  Combinatorial Optimization: Networks and Matroids.  Holt, Rinehart and Winston, 1976.
  * [20] B. Liu, K. Wang, D. Xiao, and Z. Yu.  Mathematical mechanism on dynamical system algorithms of the Ising model, 2020.
  * [21] A. Lucas.  Ising formulations of many NP problems.  Frontiers in Physics, 2, 2014.
  * [22] M. Lundy and A. Mees.  Convergence of an annealing algorithm.  Mathematical Programming, 34(1):111–124, Jan 1986.
  * [23] S. Matsubara, M. Takatsu, T. Miyazawa, T. Shibasaki, Y. Watanabe, K. Takemoto, and H. Tamura.  Digital annealer for high-speed solving of combinatorial optimization problems and its applications.  In 2020 25th Asia and South Pacific Design Automation Conference (ASP-DAC), pages 667–672, 2020.
  * [24] S. Matsubara, H. Tamura, M. Takatsu, D. Yoo, B. Vatankhahghadim, H. Yamasaki, T. Miyazawa, S. Tsukamoto, Y. Watanabe, K. Takemoto, and A. Sheikholeslami.  Ising-model optimizer with parallel-trial bit-sieve engine.  In Leonard Barolli and Olivier Terzo, editors, Complex, Intelligent, and Software Intensive Systems, pages 432–438, Cham, 2018. Springer International Publishing.
  * [25] N. Metropolis, A.W. Rosenbluth, M.N. Rosenbluth, A.H. Teller, and E. Teller.  Equation of state calculations by fast computing machines.  The Journal of Chemical Physics, 21(6):1087–1092, 1953.
  * [26] D. Mitra, F. Romeo, and A. Sangiovanni-Vincentelli.  Convergence and finite-time behavior of simulated annealing.  In 1985 24th IEEE Conference on Decision and Control, pages 761–767, 1985.
  * [27] N. Mohseni, P. L. McMahon, and T. Byrnes.  Ising machines as hardware solvers of combinatorial optimization problems.  Nature Reviews Physics, 4(6):363–379, 2022.
  * [28] H. Nakayama, J. Koyama, N. Yoneoka, and T. Miyazawa.  Description: third generation digital annealer technology, 2021.
  * [29] C.H. Papadimitriou and K. Steiglitz.  Combinatorial Optimization: Algorithms and Complexity.  Prentice-Hall, Inc., USA, 1982.
  * [30] E. Seneta.  Non-negative matrices and Markov chains.  Springer Science & Business Media, 2006.
  * [31] S. Tsukamoto, M. Takatsu, S. Matsubara, and H. Tamura.  An accelerator architecture for combinatorial optimization problems.  Fujitsu Sci. Tech. J, 53(5):8–13, 2017.
  * [32] P.J. van Laarhoven and E.H. Aarts.  Simulated Annealing: Theory and Applications.  Mathematics and Its Applications. Springer Netherlands, 1987.
  * [33] Z. Wang, A. Marandi, K. Wen, R.L. Byer, and Y. Yamamoto.  Coherent Ising machine based on degenerate optical parametric oscillators.  Phys. Rev. A, 88:063853, Dec 2013.
  * [34] D.F. Wong, H.W. Leong, and H.W. Liu.  Simulated Annealing for VLSI Design.  The Springer International Series in Engineering and Computer Science. Springer US, 1988.
  * [35] K. Yamamoto, K. Kawamura, K. Ando, N. Mertig, T. Takemoto, M. Yamaoka, H. Teramoto, A. Sakai, S. Takamaeda-Yamazaki, and M. Motomura.  STATICA: A 512-spin 0.25M-weight annealing processor with an all-spin-updates-at-once architecture for combinatorial optimization with complete spin–spin interactions.  IEEE Journal of Solid-State Circuits, 56(1):165–178, 2021.
  * [36] Y. Yamamoto, K. Aihara, T. Leleu, K. Kawarabayashi, S. Kako, M. Fejer, K. Inoue, and H. Takesue.  Coherent Ising machines - optical neural networks operating at the quantum limit.  npj Quantum Information, 3(1):49, Dec 2017.
  * [37] Z. Zhu, C. Fang, and H.G. Katzgraber.  borealis-A generalized global update algorithm for Boolean optimization problems.  Optimization Letters, 14(8):2495–2514, Nov 2020.
  * [38] Z. Zhu, A.J. Ochoa, and H.G. Katzgraber.  Efficient cluster algorithm for spin glasses in any space dimension.  Phys. Rev. Lett., 115:077201, Aug 2015.


# Unidirectional modes induced by nontraditional Coriolis force in stratified
fluids

Nicolas Perez<EMAIL_ADDRESS>Pierre Delplace  Antoine Venaille
Univ Lyon, ENS de Lyon, Univ Claude Bernard, CNRS, Laboratoire de Physique
(UMR CNRS 5672), F-69342 Lyon, France

###### Abstract

Using topology, we unveil the existence of new unidirectional modes in
compressible rotating stratified fluids. We relate their emergence to the
breaking of time-reversal symmetry by rotation and vertical mirror symmetry by
stratification and gravity. We stress the role of the Coriolis force’s
nontraditional part, induced by a rotation field tangent to the surface. In
contrast with horizontally trapped equatorial waves induced by the traditional
component of Coriolis force perpendicular to the surface, we find vertically
trapped modes that are trapped along interfaces between regions with distinct
stratification properties. We show that such modes are generalized atmospheric
Lamb waves whose direction of propagation can be selected by the
nontraditional component of the Coriolis force.

Introduction – When breaking time-reversal symmetry, two-dimensional systems
may host unidirectional trapped modes. This is the case in the quantum Hall
effect Büttiker (1988), plasmas Parker  _et al._ (2020); Fu and Qin (2021),
and certain photonic crystals Wang  _et al._ (2020), when a magnetic field is
applied. Such a symmetry breaking can also be achieved by taking advantage of
rotation effects, either at a local level through self-propelled particles
Shankar  _et al._ (2020), gyroscope arrays Nash  _et al._ (2015) or curved
photonic waveguides Ozawa  _et al._ (2019), or at a global scale by a
background solid-body rotation through the effect of the Coriolis force
Delplace  _et al._ (2017); Wang  _et al._ (2015). Over planetary scales,
atmospheres and oceans are rotating quasi-two-dimensional fluids, and as such
they support the propagation of various trapped modes in the horizontal plane
Iga (1995, 2001); Godin (2012). For such flows, it is customary to retain only
the component of the Coriolis force that involves the projection of the
rotation vector on the local vertical axis, perpendicular to the fluid, which
is encoded by the Coriolis parameter $f_{t}$ (see Fig. 1). Neglecting the
horizontal contribution $f_{nt}$ of the planetary rotation is referred to as
the traditional approximation Vallis (2017). While exact for a truly two-
dimensional fluid, this traditional approximation cannot be justified from
first principles in more comprehensive three-dimensional geophysical flows.
Thereby, the consideration of nontraditional effects has been subjected to a
resurgence of interest for the last years Gerkema  _et al._ (2008); Tort  _et
al._ (2014). Here we report a striking effect due to the nontraditional
contribution of the Coriolis force: When combined with compressibility and
density stratification, it induces previously unnoticed unidirectional trapped
waves in the fluid.

To reveal the existence of such waves, we apply the machinery of topology.
Topology offers an elegant theoretical framework to unveil the existence of
robust trapped waves along boundaries or interfaces, from quantum to classical
wave systems, without having to compute the detailed spectrum of the
inhomogeneous system for a specific geometry and set of parameters. In the
geophysical realm, this approach shed new light on various well-known waves
such as equatorial waves Matsuno (1966); Delplace  _et al._ (2017), coastal
Kelvin waves Venaille and Delplace (2021) and Lamb waves Lamb (1911); Perrot
_et al._ (2019). In the equatorial case, the interface that traps
unidirectional waves is defined by the change of sign of the traditional
component of the Coriolis parameter $f_{t}$ at the equator. In contrast, we
focus here on a flow taking place in the vertical plane at the equator (Fig.
1), where this contribution vanishes, so that the only effect of rotation
occurs through the nontraditional component $f_{nt}$. This contribution being
constant in that plane, an interface thus cannot be obtained by rotation
effects alone. Instead, we show that an emergent interface defined by the
change of sign of a stratification parameter involving density gradients,
gravity and nontraditional Coriolis force, separates two regions that behave
as if they were topologically distinct, thus hosting an unidirectional trapped
mode.

Figure 1:  Geometry of the planar problem on a rotating planet. The
traditional component $f_{t}$ cancels at the equator, while the nontraditional
term $f_{nt}$ is maximal. Modes are trapped in the direction $z$ where
$S=+f_{nt}$ (or $S=-f_{nt}$) and propagate along the zonal direction $x$.

Flow model – To describe the consequence of time-reversal symmetry breaking by
nontraditional Coriolis effects, we introduce the dynamical equations for
adiabatic perturbations of a compressible fluid initially at rest in a
uniformly rotating frame at rotation rate $\mathbf{\Omega}$, say a planet’s
oceanic or atmospheric layer, or a stellar core. We consider propagation
phenomena of small length scales compared to the typical radius of the planet
or star, and define local Cartesian coordinates (Fig. 1). The medium is
assumed to be inviscid and subjected to a constant gravity
$-g\hat{\mathbf{e}}_{z}$ along the local vertical direction.

The global equilibrium state in the rotating frame has a vertical
stratification of density (and temperature, pressure) which is characterized
by the buoyancy (or Brunt-Väisälä) frequency $N(z)$ Vallis (2017), which is
real and positive as we assume a stably stratified medium. The local speed of
sound waves, denoted $c_{s}$, is also assumed to be a function of $z$ only. To
highlight most salient features induced by nontraditional effect of the
Coriolis force, we focus on the dynamics taking place in the vertical
equatorial plane $(x,z)$ (Fig. 1), where the traditional component $f_{t}$
vanishes. The dynamics is then restricted to a two-dimensional flow in the
$(x,z)$ plane.

We consider the wave dynamics governed by the primitive equations –
conservation of momentum, mass and entropy – linearized around the state of
rest. By performing a suitable rescaling of the dynamical fields with the
density profile, those equations can be recast in a Schrödinger-like matricial
form as ${\rm
i}\partial_{t}\psi=\hat{\mathcal{H}}[z,\partial_{x},\partial_{z}]\psi$, with a
four-component vector field
$\psi=\begin{pmatrix}u&w&\Theta&p\end{pmatrix}^{t}$, where
$u\hat{\mathbf{e}}_{x}+w\hat{\mathbf{e}}_{z}$ is the rescaled perturbation’s
velocity field in the normalized basis
$(\hat{\mathbf{e}}_{x},\hat{\mathbf{e}}_{z})$, $\Theta$ and $p$ are the
rescaled perturbed potential temperature and pressure respectively (see
Supplemental Material). The wave operator reads

$\begin{split}&\hat{\mathcal{H}}=\\\ &\begin{pmatrix}0&-{\rm i}f_{nt}&0&-{\rm
i}c_{s}\partial_{x}\\\ {\rm i}f_{nt}&0&-{\rm i}N&{\rm i}S-\frac{{\rm
i}c_{s}^{\prime}}{2}-{\rm i}c_{s}\partial_{z}\\\ 0&{\rm i}N&0&0\\\ -{\rm
i}c_{s}\partial_{x}&-{\rm i}S-\frac{{\rm i}c_{s}^{\prime}}{2}-{\rm
i}c_{s}\partial_{z}&0&0\end{pmatrix}\ ,\end{split}$ (1)

where $c_{s}^{\prime}$ stands for the derivative of $c_{s}$ along the vertical
coordinate $z$, and the stratification parameter $S(z)$ is defined as Perrot
_et al._ (2019):

$S(z)=\frac{c_{s}}{2g}\left(N^{2}-\frac{g^{2}}{c_{s}^{2}}\right)-\frac{c_{s}^{\prime}}{2}\
.$ (2)

This stratification parameter contains all the physical ingredients breaking
the system’s mirror symmetry in the vertical direction $z$ (gravity $g$,
stratification of density through $N$, and non-uniformity of $c_{s}$), while
time-reversal symmetry is broken only by the constant nontraditional Coriolis
parameter $f_{nt}$.

Topological properties of the bulk waves – The physical model of interest is
inhomogeneous in the $z$ direction, through $N(z)$, $S(z)$ and $c_{s}(z)$.
However, important spectral properties of the operator (1), and in particular
the existence of unidirectional trapped modes, can be inferred from the
topology of the eigenmodes $\Psi_{n}$ of the much simpler eigenvalue problem
$\omega_{n}\Psi_{n}=H\Psi_{n}$ with

$H=\begin{pmatrix}0&-{\rm i}f_{nt}&0&c_{s}k_{x}\\\ {\rm i}f_{nt}&0&-{\rm
i}N&{\rm i}S+c_{s}k_{z}\\\ 0&{\rm i}N&0&0\\\ c_{s}k_{x}&-{\rm
i}S+c_{s}k_{z}&0&0\end{pmatrix}\ .$ (3)

This dual ”bulk” problem describes plane-wave solutions
$\psi(x,z,t)=\Psi_{n}{\rm e}^{{\rm i}\left(k_{x}x+k_{z}z-\omega t\right)}$ of
(1), valid in the case of an unbounded $(x,z)$ plane, with constant parameters
$N$, $S$ and $c_{s}$. In mathematical terms, the matrix $H$ is called the
symbol of the operator $\hat{\mathcal{H}}$, and is bijectively associated to
it through the Wigner-Weyl transform Moyal (1949) (see Supplemental Material).
In the general situation (where the plane is bounded and the parameters are
functions of $z$), plane waves are not solutions of the problem, however the
symbol matrix (3) can still be formally defined as a function of $z$ (through
the parameters $N$, $S$ and $c_{s}$) and the wave vector
$\mathbf{k}=k_{x}\hat{\mathbf{e}}_{x}+k_{z}\hat{\mathbf{e}}_{z}$. For a given
triplet $(z,k_{x},k_{z})$, the symbol matrix $H$ admits four eigenvalues
$\omega_{n}$, corresponding to two low-frequency internal-gravity wavebands
($n=\pm 1$) and two high-frequency acoustic wavebands ($n=\pm 2$), both
modified by rotation (see Fig. 2). Because of the realness of the fields, each
eigenvalue $\omega$ of $H$ at $(z,k_{x},k_{z})$ is associated with the
eigenvalue $-\omega$ at $(z,-k_{x},-k_{z})$. Therefore, it is sufficient to
consider positive frequencies ($n=1,2$) in the following.

Building upon previous topological analysis in geophysics Delplace  _et al._
(2017); Perrot  _et al._ (2019); Venaille and Delplace (2021), we expect the
existence of unidirectional trapped modes in the inhomogeneous problem (1) to
be ruled by topological indices of the symbol’s eigenmodes $\Psi_{n}$. More
specifically, these topological modes are expected to arise around a typical
wavenumber $k_{x}^{\star}$, and to be trapped around a height $z^{\star}$,
provided wavebands of the symbol cross at an isolated point
$(z^{\star},k_{x}^{\star},k_{z}^{\star})$ in the $(z,k_{x},k_{z})$ space. In
that case, it is possible to assign topological charges to such degeneracy
points, given by topological indices called the Chern numbers
$\mathcal{C}_{n}$, that account for the twisting of the eigenmode bundles
around the band-crossing point in $(z,k_{x},k_{z})$ space (see Supplemental
Material). According to index theorems Nakahara (2003), we expect a number of
$\mathcal{C}_{n}$ trapped modes whose frequency branches transit towards the
waveband $n$ of the inhomogeneous problem when swapping $k_{x}$.

Figure 2:  Degeneracy points (in black) between the positive-frequency bands
($n=+1$ and $+2$) of the spectrum of (3), with $f_{nt}>0$. The pairs of
integers indicate the topological charges $(\mathcal{C}_{1},\mathcal{C}_{2})$
of the two bands involved in a given degeneracy point. The insets on the left
show their dispersion relations associated with the green, red and blue planes
in parameter space. We also show the degeneracy points (grey) between the
negative-frequency bands (not shown on the left), with
$(\mathcal{C}_{-2},\mathcal{C}_{-1})$.

The first step of our strategy then consists in identifying the degeneracy
points of the symbol’s eigenvalues. We find such crossing points between the
gravity ($n=1$) and the acoustic ($n=2$) wavebands precisely when $S=+f_{nt}$
and $S=-f_{nt}$ (Fig. 2), respectively designated with the $\pm$ sign in the
following paragraphs. Those are located at
$\left(z^{\star},k_{x}^{\star},k_{z}^{\star}\right)$, with $k_{z}^{\star}=0$,
$k_{x}^{\star}=\mp\sqrt{f_{nt}^{2}+N(z^{\star})^{2}}/c_{s}(z^{\star})$ and
$z^{\star}$ such that $S(z^{\star})=\pm f_{nt}$: their existence thus results
from the combination of mirror symmetry breaking (due to $S$) and time-
reversal symmetry breaking (due to $f_{nt}$). Since the degeneracy condition
depends on $z$ through $S$ only, one can consider the more convenient
parameter space of $(S,k_{x},k_{z})$ to compute the Chern numbers of the
symbol (see Supplemental Material).

Once the degeneracy points are found out, one needs to compute the Chern
numbers associated to the eigenstate bundles involved in the waveband
crossing. This is a standard calculation that can always be performed at least
numerically, by means of the integration of a local quantity, called the Berry
curvature Berry (1984), over a closed surface enclosing the degeneracy point
in $(S,k_{x},k_{z})$ space (see Supplemental Material). Actually, since this
degeneracy simply involves a generic conical two-band crossing point, the
possible values of the Chern numbers of the two eigenstate bundles are known
to be opposite and of absolute value $1$ (see e.g. the appendix of Venaille
and Delplace (2021)). Therefore, one can directly infer the existence of one
unidirectional wave, trapped around the interface line $z=z^{\star}$ in the
vertical $(x,z)$ plane, that is where the quantity $S(z)\mp f_{nt}$ changes
sign.

Figure 3:  Numerical solutions of the eigenvalue problem (1) projected on the
Fourier basis – such that $(\partial_{x},\partial_{t})\rightarrow{\rm
i}(k_{x},-\omega)$ –, with constant sound speed $c_{s}$, obtained with Dedalus
Burns  _et al._ (2020). ($\mathbb{a}$) Spectrum for $f_{nt}=0.1$: $S$ does not
take the value $f_{nt}$, no mode transits through the gap.
$N_{m}=\min\\{N(z)\\}$. ($\mathbb{b}$) $f_{nt}=2$: $S$ takes twice the value
of $f_{nt}$, a mode is localised at each interface $S=f_{nt}$ and transits
through the gap around $k_{x}^{\star}$. The pressure perturbations’s amplitude
$|\tilde{p}(z)|$ of the topological modes at $k_{x}=k_{0}$ are plotted with
the background stratification $S$ in green. ($\mathbb{c}$) Plot of the
velocity field of the topological modes with wavenumber $k_{0}$, with pressure
levels in background (positive values in solid lines and negative in dashed
lines).

Unidirectional modes induced by the nontraditional Coriolis component – We
then check this prediction numerically by computing the frequency spectrum of
$\hat{\mathcal{H}}$ for solutions of the form
$\psi(x,z,t)=\begin{pmatrix}\tilde{u}(z)&\tilde{w}(z)&\tilde{\Theta}(z)&\tilde{p}(z)\end{pmatrix}^{t}{\rm
e}^{{\rm i}(k_{x}x-\omega t)}$, with varying parameters in the vertical
direction $z$. We can always consider $f_{nt}>0$, since the eastward
(increasing $x$) and westward (decreasing $x$) directions can be defined with
the direction of planetary rotation. For numerical convenience, we also
consider the profile of $S(z)>0$ depicted in Fig. 3, that is periodic in $z$
and crosses the value $f_{nt}$ in both an increasing and a decreasing fashion
(the degeneracy point involved in this numerical simulation is thus the one at
$S=+f_{nt}$, but a general conclusion is given in the next paragraph,
including for the case $S=-f_{nt}$), so that we can discuss, in the same
simulation, the two opposite slopes $S^{\prime}(z^{\star})>0$ and $<0$ at an
interface where $S(z)-f_{nt}$ changes sign. In the case where $S(z)$ does not
reach the value of the nontraditional Coriolis parameter $f_{nt}$, (Fig.
3(a)), there is no mode that transits between the gravity and the acoustic
waveband, as expected from the topological analysis, since there is no
degeneracy point of the symbol and thereby no topological charge. In contrast,
if the profile of $S(z)$ is such that $S(z)-f_{nt}$ changes sign (Fig. 3(b)),
a topological charge arises in the symbol picture at
$(z^{\star},k_{x}^{\star},0)$, and accordingly, a new mode appears in the
spectrum of $\hat{\mathcal{H}}$ around $k_{x}=k_{x}^{\star}$. This mode
transits from the internal-gravity waveband to the acoustic one, so that the
number of modes gained or lost by the two wavebands when sweeping $k_{x}$ is
correctly given by the Chern numbers $\pm 1$. Such a spectral flow from one
waveband to another implies the sign of the group velocity of the
corresponding wave to be fixed, hence its unidirectionality. We have also
computed the spreading of this mode in the $z$ direction and found that it is
indeed localized where $S-f_{nt}$ changes sign (in color in Fig. 3).

The numerical analysis reveals interesting differences on the nature of the
topological trapped modes depending on the sign of the slope of $S(z)$ at the
interface where $S(z)-f_{nt}$ changes sign. The corresponding rescaled
velocity fields are shown in Fig. 3(c), and their typical vertical trapping
length around the interface is
$\mathcal{L}^{\star}=\sqrt{c_{s}(z^{\star})/|S^{\prime}(z^{\star})|}$ Perrot
_et al._ (2019): the steeper the slope of $S$ at the crossing point, the more
localised the modes.

  * •

If $S$ decreases, a first unidirectional mode appears (red in Fig. 3). For an
interface $S=+f_{nt}$ this mode propagates westward, lies around
$k_{x}^{\star}=-\sqrt{f_{nt}^{2}+N(z^{\star})^{2}}/c_{s}(z^{\star})$ and
transits from the acoustic waveband to the gravity waveband (as $k_{x}$
increases around the value $k_{x}^{\star}$) in the positive-frequency
spectrum. For $S=-f_{nt}$ it propagates eastward, lies around
$k_{x}^{\star}=\sqrt{f_{nt}^{2}+N(z^{\star})^{2}}/c_{s}(z^{\star})$ and
transits from the gravity waveband to the acoustic waveband in the positive-
frequency spectrum. In the case where $c_{s}$ is constant (e.g. in an
isothermal ideal gas) this mode is longitudinal ($w=0$), non-dispersive, with
dispersion relation $\omega=\pm c_{s}k_{x}$ if $S(z^{\star})=\mp f_{nt}$ (Fig.
3(b)). It thus shares strong similarities with the Lamb-like wave discussed in
Perrot  _et al._ (2019) and the equatorial Kelvin wave Delplace  _et al._
(2017). However its nature is different: while it is unidirectional and
confined in the vertical direction, the equatorial Kelvin wave is confined in
the meridional direction due to the gradient of $f_{t}$, the traditional
Coriolis parameter – which vanishes here –, and Lamb-like waves are not
unidirectional because their existence does not rely on time-reversal symmetry
breaking. Therefore this topological mode can be seen as the rotating
counterpart of the Lamb-like waves in stratified compressible fluids. We have
checked that this wave still exists when $c_{s}$ varies with $z$, as expected
from our topological analysis, although it becomes dispersive.

  * •

If $S$ increases, another unidirectional mode of opposite spectral flow (i.e.
transiting from the gravity waveband to the acoustic wave band if
$S(z^{\star})=f_{nt}$ and the other way round if $S(z^{\star})=-f_{nt}$)
appears around $k_{x}^{\star}$ (blue in Fig. 3). Contrary to the previous
topological mode, this one has a non-longitudinal velocity field (Fig. 3(c)).
It is reminiscent of a zero-group-velocity oscillation at frequency $N$,
modified by rotation and varying buoyancy frequency $N(z)$. As their group
velocities have the same sign, both modes propagate in the same direction,
that is westward if $S(z^{\star})=+f_{nt}$ and eastward if
$S(z^{\star})=-f_{nt}$.

According to this analysis, a topological transition occurs when the number of
crossings between the vertical profile $S(z)$ and the value $\pm f_{nt}$
changes, leading to the appearance/disappearance of unidirectional interface
waves whose frequency branches transit between the gravity and acoustic
wavebands. These modes shall not be confused with the ones described by Iga
Iga (1995, 2001), whose existence depend on the boundary conditions, neither
with those addressed by Godin Godin (2012) who describes an incompressible
wave motion solution of the compressible dynamics with a free surface and
traditional rotation, whereas compressibility, nontraditional rotation and
internal interfaces are key ingredients of the modes studied in this letter.

Discussion – Figure 3 depicts an ideal situation where $S$ takes the value
$+f_{nt}$ and not $-f_{nt}$. Westward topological modes are thus allowed to
propagate whereas the eastward ones are not. This shows that rotation selects
the direction of propagation of the Lamb-like waves studied in Perrot  _et
al._ (2019). In realistic planetary fluids, $f_{nt}$ is however expected to be
much smaller. On Earth, $f_{nt}\approx 10^{-4}$ s-1 remains negligible
compared to the typical values of $S$ in the atmosphere ($\approx-10^{-2}$
s-1), and the trapping length in the ocean is higher than its typical depth
(see Supplemental Material of Perrot  _et al._ (2019)). Therefore the Earth’s
atmosphere and the oceans are not fitted for potential observations of these
topological modes.

Nevertheless, if $S$ changes sign at a given $z_{0}$, then we expect the
topological modes’ location to be shifted even if $f_{nt}>0$ is small: if
$S^{\prime}(z_{0})<0$, the eastward Lamb-like mode should be located at
$z>z_{0}$ and the westward one at $z<z_{0}$, and conversely if
$S^{\prime}(z_{0})>0$. As far as applications are concerned, this effect of
splitting the topological modes’ location in the weak rotational regime is a
feature that is more likely to be observed, especially in the astrophysical
context. In fact, possible observations of Lamb-like waves in asteroseismology
(the study of stellar pulsations) have been previously argued, yet without
discussing the effect of rotation Perrot  _et al._ (2019); Leclerc  _et al._ .
Indeed, using the model of a polytropic star with constant adiabatic index,
$S$ is a decreasing function of the radius $r$, that varies from $+\infty$ at
the center of the star to $-\infty$ at the surface Leclerc  _et al._ . It thus
cancels at some radius $r_{0}$ and generically takes both values $\pm f_{nt}$
for any rotation rate. Therefore, in rotating stellar interiors, one can
expect a westward Lamb-like mode closer to the core ($r<r_{0}$) and an
eastward one closer to the surface ($r>r_{0}$), thus improving the chance of a
possible direct observation of the latter. These modes are well-separated if
$f_{nt}$ is high enough so that the topological modes do not completely
overlap, in other words if the sum of their respective extension in the radial
direction, $\mathcal{L}=\sqrt{c_{s}/|S^{\prime}|}_{S=\pm f_{nt}}$, is not high
compared to the distance $d$ separating the interfaces. This could be achieved
for some young stars of the upper main sequence of the Hertzprung-Russell
diagram, which can reach high rotation rates. For instance, Vega, in the Lyra
constellation, has a rotation period of about 0.7 days Takeda  _et al._
(2008), thus $f_{nt}\approx 2\,10^{-4}$ s-1 at the equator, yielding a ratio
of order 1 between $d$ and the sum of the trapping lengths.

Conclusion – Our study shows the essential role of the nontraditional Coriolis
force in the propagation properties of rotating stratified fluid waves. Thanks
to topological arguments, we have found new unidirectional waves transiting
between the low-frequency gravity and high-frequency acoustic wavebands, that
are trapped along lines determined by a balanced between rotation and vertical
stratification. The existence of these modes is robust to perturbations in the
stratification profile until the interface constraint $S=\pm f_{nt}$ is
broken. All salient features of the two-dimensional case presented in this
letter can be extended to the three-dimensional case with full Coriolis force,
which will be explained in a companion paper. The application of this work to
asteroseismology provides an promising line of research for future
investigations of topological waves in natural media, in particular in the
presence of non-Hermitian effects induced by dissipation, mean flows and
convection.

##  I Supplemental Material

Model equations and fields rescaling – We introduce a local basis $(x,y,z)$,
where $z$ is normal to the sphere, so that the gravity be
$\mathbf{g}=-g\hat{\mathbf{e}}_{z}$, $x$ follow the longitudinal direction
pointing eastward and $y$ the meridional direction pointing northward (Fig.
1). In that local basis the solid body rotation rate has two components and we
note

$2\mathbf{\Omega}=\begin{pmatrix}0\\\ f_{nt}\\\ f_{t}\\\
\end{pmatrix}=2\Omega\begin{pmatrix}0\\\ \cos{\theta}\\\ \sin{\theta}\\\
\end{pmatrix}\ ,$ (4)

where $\Omega$ is the value of the rotation rate and $\theta$ the latitude. We
are interested in small perturbations of an inviscid fluid initially at rest
in the rotating frame. Absorbing the centrifugal force into the gravity
$\mathbf{g}$, Euler’s equation reads

$\frac{\partial\mathbf{v}}{\partial
t}+(\mathbf{v}\cdot\nabla)\mathbf{v}=-\frac{1}{\rho}\nabla
P-g\hat{\mathbf{e}}_{z}-2\mathbf{\Omega}\times\mathbf{v}\ ,$ (5)

where $P$ is the pressure, $\rho$ the density and $\mathbf{v}$ the velocity
field.

We also introduce the equation of conservation of mass,

$\frac{\partial\rho}{\partial t}+\nabla\cdot(\rho\mathbf{v})=0\ ,$ (6)

and the equation of conservation of entropy (assuming adiabatic evolution of
fluid parcels),

$\frac{\partial\rho}{\partial
t}+(\mathbf{v}\cdot\nabla)\rho=c_{s}^{-2}\left[\frac{\partial P}{\partial
t}+(\mathbf{v}\cdot\nabla)P\right]\ .$ (7)

The isentropic compressibility factor of the fluid is noted
$\chi_{s}=c_{s}^{-2}$, where $c_{s}$ is the local celerity of sound waves. For
the fluid at rest, we have $\mathbf{v}\equiv\mathbf{0}$ and the fields (noted
with a $0$ index) are stationary ($\partial_{t}\equiv 0$), hence the
hydrostatic relation

$\nabla P_{0}=-g\rho_{0}\hat{\mathbf{e}}_{z}\ .$ (8)

Therefore the fluid is stratified along the $z$ direction and relation (8)
reads

$\frac{{\rm d}P_{0}(z)}{{\rm d}z}=-g\rho_{0}(z)\ .$ (9)

We also consider that $c_{s}$ is a function of $z$. We now write the
perturbation fields around the stationary state of rest as

$\mathbf{v}^{\prime}\equiv\mathbf{v},\quad\rho^{\prime}\equiv\rho-\rho_{0},\quad
P^{\prime}\equiv P-P_{0}\ ,$ (10)

and linearize equations (5), (6) and (7):

$\displaystyle\rho_{0}$
$\displaystyle\frac{\partial\mathbf{v}^{\prime}}{\partial t}=-\nabla
P^{\prime}-g\rho^{\prime}\hat{\mathbf{e}}_{z}-2\rho_{0}\mathbf{\Omega}\times\mathbf{v}^{\prime}\
,$ (11a) $\displaystyle\frac{\partial\rho^{\prime}}{\partial
t}+v_{z}^{\prime}\frac{{\rm d}\rho_{0}}{{\rm
d}z}=-\rho_{0}\nabla\cdot\mathbf{v}^{\prime}\ ,$ (11b)
$\displaystyle\frac{\partial\rho^{\prime}}{\partial
t}+v_{z}^{\prime}\frac{{\rm d}\rho_{0}}{{\rm
d}z}=c_{s}^{-2}\left(\frac{\partial P^{\prime}}{\partial
t}-g\rho_{0}v_{z}^{\prime}\right)\ .$ (11c)

It is now convenient to use the following transformations Perrot  _et al._
(2019):

$\mathbf{u}\equiv\rho_{0}^{1/2}\>\mathbf{v}^{\prime}\equiv(u,v,w),\quad\tilde{\rho}\equiv\rho_{0}^{-1/2}\>\rho^{\prime},\quad\tilde{P}\equiv\rho_{0}^{-1/2}\>P^{\prime}\
,$ (12)

so that equations (11) can be rewritten as

$\displaystyle\partial_{t}\mathbf{u}=-\nabla\tilde{P}-g\tilde{\rho}\hat{\mathbf{e}}_{z}-2\mathbf{\Omega}\times\mathbf{u}-\frac{1}{2\rho_{0}}\frac{{\rm
d}\rho_{0}}{{\rm d}z}\tilde{P}\hat{\mathbf{e}}_{z}\ ,$ (13a)
$\displaystyle\partial_{t}\left(\tilde{\rho}-c_{s}^{-2}\tilde{P}\right)=-\left(\frac{1}{\rho_{0}}\frac{{\rm
d}\rho_{0}}{{\rm d}z}+gc_{s}^{-2}\right)w\ ,$ (13b)
$\displaystyle\partial_{t}\left(\frac{\tilde{P}}{c_{s}}\right)=-c_{s}\nabla\cdot\mathbf{u}+\left(\frac{g}{c_{s}}+\frac{c_{s}}{2\rho_{0}}\frac{{\rm
d}\rho_{0}}{{\rm d}z}\right)w\ .$ (13c)

At this point we assume that the quantity

$\displaystyle-\left(\frac{1}{\rho_{0}}\frac{{\rm d}\rho_{0}}{{\rm
d}z}+gc_{s}^{-2}\right)\ ,$ (14)

that is a function of $z$, is strictly positive everywhere. This ensures the
hyperbolicity of the set of partial differential equations (13). Physically it
means that the fluid is stably stratified. Therefore we can define the
following frequencies:

$\begin{split}&N(z)\equiv\sqrt{-\left(\frac{g}{\rho_{0}}\frac{{\rm
d}\rho_{0}}{{\rm d}z}+\frac{g^{2}}{c_{s}^{2}}\right)}\\\
\text{and}\quad&S(z)=\frac{c_{s}}{2g}\left(N^{2}-\frac{g^{2}}{c_{s}^{2}}\right)-\frac{c_{s}^{\prime}}{2}\
.\end{split}$ (15)

The Brunt-Väisälä frequency $N(z)$ Vallis (2017) is the pulsation at which a
fluid particle oscillates vertically in the stably stratified fluid, at
position $z$. The derivative of $c_{s}$ with respect to $z$ is noted
$c_{s}^{\prime}$.

Finally we define the following fields for perturbations:

$p\equiv\frac{\tilde{P}}{c_{s}}\quad\text{and}\quad\Theta\equiv\frac{g}{N}(\tilde{\rho}-c_{s}^{-2}\tilde{P})\
.$ (16)

The second one corresponds to the perturbation of potential temperature. Using
definitions (15) and (16) along with equations (13), we finally get

$\displaystyle\partial_{t}\mathbf{u}=-2\mathbf{\Omega}\times\mathbf{u}-N\Theta\hat{\mathbf{e}}_{z}+\left(S-\frac{c_{s}^{\prime}}{2}\right)p\hat{\mathbf{e}}_{z}-c_{s}\nabla
p\ ,$ (17a) $\displaystyle\partial_{t}\Theta=Nw\ ,$ (17b)
$\displaystyle\partial_{t}p=-c_{s}\nabla\cdot\mathbf{u}-\left(S+\frac{c_{s}^{\prime}}{2}\right)w\
,$ (17c)

that we can write in matrix form

$\partial_{t}\psi=\hat{\mathcal{H}}\psi\ ,$ (18)

with $\psi=\begin{pmatrix}u&v&w&\Theta&p\end{pmatrix}^{t}$ and the operator

$\begin{split}&\hat{\mathcal{H}}=\\\
&\begin{pmatrix}0&f_{t}&-f_{nt}&0&-c_{s}\partial_{x}\\\
-f_{t}&0&0&0&-c_{s}\partial_{y}\\\
f_{nt}&0&0&-N&S-\frac{c_{s}^{\prime}}{2}-c_{s}\partial_{z}\\\ 0&0&N&0&0\\\
-c_{s}\partial_{x}&-c_{s}\partial_{y}&-S-\frac{c_{s}^{\prime}}{2}-c_{s}\partial_{z}&0&0\end{pmatrix}\end{split}$
(19)

Among the parameters of expression (19), $N$, $S$ and $c_{s}$ are functions of
$z$, whereas $f_{t}$ and $f_{nt}$ are functions of the meridional coordinate
$y$. In many geophysical situations, the planetary radius is much higher than
the variation lengthscales in the vertical direction, therefore we can
consider $f_{t}$ and $f_{nt}$ constant when considering fluid motions at much
smaller scale than the radius. In order to emphasize the nontraditional
effect, we consider in our work the particular case $f_{t}=0$ (equatorial
region), which allows one to discard the meridional dependence of the problem
(i.e. equation (18) is compatible with $v$, $\partial_{y}=0$). The problem
thus conveniently reduces to a four-by-four matrical equation with parameters
depending on $z$ only. Thus follows the Schrödinger-like form ${\rm
i}\partial_{t}\psi=\hat{\mathcal{H}}[z,\partial_{x},\partial_{z}]\psi$,
introduced in equation (1).

Wavebands and Chern numbers – We introduce the symbol $H[z,k_{x},k_{z}]$ of
the operator $\hat{\mathcal{H}}[z,\partial_{x},\partial_{z}]$, defined by the
following relation, called the Weyl correspondence Moyal (1949); Onuki (2020):

$\begin{split}&\hat{\mathcal{H}}[z,\partial_{x},\partial_{z}]\Psi(x,z)=\frac{1}{(2\pi)^{2}}\int
H\left[\frac{z+z^{\prime}}{2},k_{x},k_{z}\right]\\\
&\times\Psi(x^{\prime},z^{\prime}){\rm e}^{{\rm
i}[k_{x}(x-x^{\prime})+k_{z}(z-z^{\prime})]}\>{\rm d}x^{\prime}{\rm
d}z^{\prime}{\rm d}k_{x}{\rm d}k_{z}\ .\end{split}$ (20)

To clarify formula (20), the symbol is a 4-by-4 matrix-valued function and
$\hat{\mathcal{H}}$ is an operator applying on test functions $\Psi$ defined
over the $(x,z)$ space. For a given $(x,z)$, $\Psi(x,z)\in\mathbb{C}^{4}$.
Moreover, in this case, both $H$ and $\hat{\mathcal{H}}$ are Hermitian for
their respective canonical inner products.

Since the scalar operator $-{\rm i}c_{s}(z)\partial_{z}$ is associated to the
scalar symbol $c_{s}k_{z}+\frac{\rm i}{2}\frac{{\rm d}c_{s}}{{\rm d}z}$, we
get

$\begin{split}&H[z,k_{x},k_{z}]=\\\ &\begin{pmatrix}0&-{\rm
i}f_{nt}&0&c_{s}(z)k_{x}\\\ {\rm i}f_{nt}&0&-{\rm i}N(z)&c_{s}(z)k_{z}+{\rm
i}S(z)\\\ 0&{\rm i}N(z)&0&0\\\ c_{s}(z)k_{x}&c_{s}(z)k_{z}-{\rm
i}S(z)&0&0\end{pmatrix}\end{split}$ (21)

When all coefficients $N$, $S$ and $c_{s}$ are constant, one may naturally
look at the plane-wave solutions of the Schrödinger-like equation, i.e.
$\psi(x,z,t)\propto{\rm e}^{{\rm i}(k_{x}x+k_{z}z-\omega t)}$. In this case
the dispersion relation is given by the eigenvalues of the symbol $H$, and the
polarization relations between the different fields for an eigenmode are given
by the corresponding eigenvector of $H$.

Roughly speaking, the index theorems tell us that a singularity in the
spectrum of the operator $\hat{\mathcal{H}}$ shall be generated by a
corresponding topological singularity of the symbol $H$. While the first may
be quite difficult to identify, the second is more straightforward: for a
given point $(z,k_{x},k_{z})$ the symbol $H$ has four eigenvalues (except at
degeneracy points), thus defining four wavebands (indexed by $n=-1,-2,1,2$) in
this space. We can also define four bundles with the (normalized) eigenvectors
$\Psi_{n}[z,k_{x},k_{z}]$ defined on the parameter space. A topological
singularity point in parameter space may then arise for a waveband whenever
the corresponding family of eigenvectors or bundle can’t be continuously
defined on any closed surface surrounding this point. These singularities
correspond to the degeneracy or crossing points between the wavebands, which
are much easier to identify. As explained in the main text, a pair of
degeneracy points appear in the bulk spectrum when $z$ is such that $S(z)=\pm
f_{nt}$: one between bands $1$ and $2$ at
$(k_{x},k_{z})=(\mp\frac{\sqrt{N^{2}+f_{nt}^{2}}}{c_{s}},0)$ and its symmetric
$-\omega$ between bands $-2$ and $-1$ at
$(\pm\frac{\sqrt{N^{2}+f_{nt}^{2}}}{c_{s}},0)$. Since the degeneracy condition
involves the variable $z$ only through $S$, we will work for convenience in
the parameter space $(S,k_{x},k_{z})$ instead of $(z,k_{x},k_{z})$ (this does
not change the values of the Chern numbers).

In order to quantify the topological charge associated to a degeneracy point
between two bands, we first define the Berry curvature Berry (1984) of the
waveband $n$ as

$\mathbf{F}^{(n)}[S,k_{x},k_{z}]=\nabla\times\left({\rm
i}\Psi_{n}^{\dagger}\nabla\Psi_{n}\right)$ (22)

where $\nabla$ is the gradient in parameter space $(S,k_{x},k_{z})$. One can
check that this quantity is real. In our case the parameter space is three-
dimensional, therefore the Berry curvature (22) can be expressed as a three-
dimensional vector. As the external derivative of the $n^{th}$ bundle’s
connection ${\rm i}\Psi_{n}^{\dagger}\nabla\Psi_{n}$, the Berry curvature is a
relevant quantity to identify topological singularities. More precisely, one
can show that the flux $(1/2\pi)\oint\mathbf{F}^{(n)}\cdot{\rm
d}\mathbf{\Sigma}$ through any closed surface in parameter space is an integer
(${\rm d}\mathbf{\Sigma}$ pointing in the outward direction). Much like the
Gauss formula, this integer counts the degree of a singularity (or
multivaluedness) enclosed by the integration surface. Therefore we can
associate a topological charge, called a Chern number Berry (1984), to an
isolated degeneracy point of the waveband $n$, as

$\mathcal{C}_{n}=\frac{1}{2\pi}\oint\mathbf{F}^{(n)}\cdot{\rm
d}\mathbf{\Sigma}\ ,$ (23)

where the integral is performed on a closed surface enclosing only the
crossing point considered. As a degeneracy point involves two wavebands (say
$n$ and $m$), topological charges appear by pair at this point, with
$\mathcal{C}_{n}=-\mathcal{C}_{m}$.

In flow models, as explained in the main text, the index theorem manifests
through the so-called bulk-boundary correspondence Delplace  _et al._ (2017);
Perrot  _et al._ (2019); Venaille and Delplace (2021), relating the
topological charge $\mathcal{C}_{n}$ to the number of modes transiting between
the corresponding wavebands of operator $\hat{\mathcal{H}}$ when crossing a
degeneracy point.

## References

  * Büttiker (1988) M. Büttiker, Physical Review B 38, 9375 (1988).
  * Parker  _et al._ (2020) J. B. Parker, J. Marston, S. M. Tobias, and Z. Zhu, Physical Review Letters 124, 195001 (2020).
  * Fu and Qin (2021) Y. Fu and H. Qin, Nature Communications 12, 1 (2021).
  * Wang  _et al._ (2020) H. Wang, S. Gupta, and B. e. a. Xie, Front. Optoelectron. 13, 50 (2020).
  * Shankar  _et al._ (2020) S. Shankar, A. Souslov, M. J. Bowick, M. C. Marchetti, and V. Vitelli, arXiv preprint arXiv:2010.00364 (2020).
  * Nash  _et al._ (2015) L. M. Nash, D. Kleckner, A. Read, V. Vitelli, A. M. Turner, and W. T. Irvine, Proceedings of the National Academy of Sciences 112, 14495 (2015).
  * Ozawa  _et al._ (2019) T. Ozawa, H. M. Price, A. Amo, N. Goldman, M. Hafezi, L. Lu, M. C. Rechtsman, D. Schuster, J. Simon, O. Zilberberg, _et al._ , Reviews of Modern Physics 91, 015006 (2019).
  * Delplace  _et al._ (2017) P. Delplace, J. Marston, and A. Venaille, Science 358, 1075 (2017).
  * Wang  _et al._ (2015) Y.-T. Wang, P.-G. Luan, and S. Zhang, New Journal of Physics 17, 073031 (2015).
  * Iga (1995) K. Iga, Journal of Fluid Mechanics 294, 367 (1995).
  * Iga (2001) K. Iga, Fluid Dynamics Research 28, 465 (2001).
  * Godin (2012) O. A. Godin, Physical review letters 108, 194501 (2012).
  * Vallis (2017) G. K. Vallis,  _Atmospheric and oceanic fluid dynamics_ (Cambridge University Press, 2017).
  * Gerkema  _et al._ (2008) T. Gerkema, J. Zimmerman, L. Maas, and H. Van Haren, Reviews of Geophysics 46 (2008).
  * Tort  _et al._ (2014) M. Tort, T. Dubos, F. Bouchut, and V. Zeitlin, Journal of fluid mechanics 748, 789 (2014).
  * Matsuno (1966) T. Matsuno, Journal of the Meteorological Society of Japan. Ser. II 44, 25 (1966).
  * Venaille and Delplace (2021) A. Venaille and P. Delplace, Physical Review Research, in press, arXiv:2011.03440 (2021).
  * Lamb (1911) H. Lamb, Proceedings of the Royal Society of London. Series A, Containing Papers of a Mathematical and Physical Character 84, 551 (1911).
  * Perrot  _et al._ (2019) M. Perrot, P. Delplace, and A. Venaille, Nature Physics 15, 781 (2019).
  * Moyal (1949) J. E. Moyal, in  _Mathematical Proceedings of the Cambridge Philosophical Society_ , Vol. 45 (Cambridge University Press, 1949) pp. 99–124.
  * Nakahara (2003) M. Nakahara,  _Geometry, Topology and Physics_ (CRC Press, 2003).
  * Berry (1984) M. V. Berry, Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences 392, 45 (1984).
  * Burns  _et al._ (2020) K. J. Burns, G. M. Vasil, J. S. Oishi, D. Lecoanet, and B. P. Brown, Physical Review Research 2, 023068 (2020).
  * (24) A. Leclerc, G. Laibe, P. Delplace, A. Venaille, and N. Perez, in preparation .
  * Takeda  _et al._ (2008) Y. Takeda, S. Kawanomoto, and N. Ohishi, The Astrophysical Journal 678, 446 (2008).
  * Onuki (2020) Y. Onuki, Journal of Fluid Mechanics 883 (2020).