# Equilibrium and nonequilibrium properties of Euler turbulence

Mahendra K Verma  Department of Physics, Indian Institute of Technology
Kanpur, Kanpur 208016, India<EMAIL_ADDRESS>Soumyadeep Chatterjee  Department
of Physics, Indian Institute of Technology Kanpur, Kanpur 208016, India
<EMAIL_ADDRESS>Shashwat Bhattacharya  Institute for Thermodynamics
and Fluid Dynamics, Technische Universität Ilmenau, Postfach 100565, D-98684
Ilmenau, Germany

###### Abstract

In this article, we report the equilibrium and nonequilibrium features of two-
dimensional (2D) and three-dimensional (3D) Euler turbulence. To obtain a full
range of equilibrium spectra, we perform pseudo-spectral simulations of Euler
turbulence using $\delta$-correlated velocity field as an initial condition.
These simulations provide zero energy flux and Maxwell-Boltzmann distribution
for the velocity field, thus providing direct verification of the absolute
equilibrium theory of turbulence. However, for ordered initial condition, 2D
Euler turbulence remains out of equilibrium, with flow getting more ordered
with time. We show that the hydrodynamic entropy of 2D Euler turbulence
decreases with time, even though the system is isolated.

Arrow of time, Nonequilibrium phenomena, Energy transfers, Multiscale physics,
Turbulence.

##  I INTRODUCTION

The present article is based on an invited talk delivered in APPC15. It
contains review of past works on Euler turbulence, as well as several new
results, e.g., nonequilibrium behaviour of two-dimensional (2D) Euler
turbulence, and hydrodynamic entropy.

Physical processes are either in equilibrium or out of equilibrium [1, 2].
Thermodynamics provides many examples of equilibrium processes, e.g., thermal
gas, Bose gas, magnetic systems in heat bath. In a gas or liquid under
equilibrium, there is no net flow of energy or matter from one region to
another statistically. This property is called detailed balance. Besides, the
average energy and entropy of an equilibrium system remain invariant in time
[3, 4].

On the other hand, nonequilibrium systems are time-dependent with detailed
balance broken [1, 2]. Earth’s atmosphere, turbulent convection, hydrodynamic
turbulence, and earthquakes are examples of such systems. Kolmogorov’s theory
of turbulence [5, 6, 7, 8] provides valuable insights into dissipative
hydrodyanmic turbulence. In this theory, a viscous fluid is forced at large
scales. The energy injected at the large scale is transferred to intermediate
scale (called inertial range) and then to small scales, where the injected
energy is dissipated. Under a steady state, the inertial-range energy spectrum
is $E(k)=K_{\mathrm{Ko}}\epsilon^{2/3}k^{-5/3}$, where $K_{\mathrm{Ko}}$ is
Kolmogorov’s constant, $\epsilon$ is the energy flux in the inertial range,
and $k$ is the wavenumber.

In this paper, we focus on turbulence in incompressible Euler equation, which
is the hydrodynamic equation with zero external force and zero viscosity. As
we describe below, turbulence in Euler equation, referred to as Euler
Turbulence, is very different from Kolmogorov’s model of turbulence, which
applies to viscous flows. Kraichnan [9] and Lee [10] argued that Euler
turbulence has similarities with equilibrium thermodynamics, and constructed
absolute equilibrium theory of Euler turbulence. By invoking Liouville’s
theorem for the Fourier modes of Euler equation, Kraichnan and Lee derived
equilibrium solution of the three-dimensional (3D) Euler equation with a
finite number of Fourier modes, also called truncated Euler equations. For
this solution, the kinetic energy flux vanishes, and the kinetic energy
spectra are

$\displaystyle E(k)\propto\frac{k^{2}}{\beta-\gamma
k^{2}}~{}\mathrm{for~{}3D},$ (1) $\displaystyle
E(k)\propto\frac{k}{\beta+\gamma k^{2}}~{}\mathrm{for~{}2D},$ (2)

where $\beta$ and $\gamma$ are constants, which are associated with the
conserved quantities: the kinetic energy ($\int d{\bf r}u^{2}/2$) and kinetic
helicity ($\int d{\bf r}({\bf u}\cdot\bm{\omega})$) in 3D, and kinetic energy
and enstrophy ($\int d{\bf r}\omega^{2}/2$) in 2D, where $\bm{\omega}$ is the
vorticity field. Onsager [11] modelled 2D Euler flow using a collection of
point vortices that interact with logarithmic potential. Further, Onsagar
analyzed the above system and predicted negative temperature and a large
cluster of same-circulation vortices for large energies.

There have been numerous efforts to verify the aforementioned predictions of
Kraichnan [9] and Lee [10]. Most numerical works start with ordered initial
condition that asymptotically approach the above equilibrium state. For
example, Cichowlas  _et al._ [12] simulated 3D Euler turbulence using a large-
scale Taylor-Green vortex as an initial condition. They observed Kolmogorov’s
spectrum for the intermediate wavenumbers and $k^{2}$ spectrum at large
wavenumbers, which is a combination of nonequilibrium and equilibrium states.
In later works, Krstulovic et al. [13] simulated truncated Euler equation with
a large-scale helical flow as an initial condition, and obtained Kraichnan’s
helical absolute equilibrium state at small scales. Dallas et al. [14] and
Alexakis and Brachet [15, 16] studied Kolmogorov flow where the forcing is
employed at intermediate scales. They observed that the flow at scales larger
than the forcing scale reaches a thermal equilibrium and exhibits $k^{2}$
energy spectrum. Similar behaviour has been observed for time-dependent
projected Gross-Pitaevskii equation [17], as well as for truncated
dissipation-less Burgers equation [18, 19]. Also refer to earlier reviews by
Orszag [20], Kraichnan and Montgomery [21].

Recently, Verma [22] employed $\delta$-correlated velocity field as initial
condition and performed preliminary simulation of 3D Euler turbulence. They
observed $k^{2}$ energy spectrum from the beginning itself. In the present
paper, we perform detailed simulations with and without kinetic helicty and
report equilibrium solutions for both 2D and 3D Euler equation, i.e., Eqs. (2,
1). In these solutions, the velocity field is as random as in thermodynamic
gas. In addition, we show that the kinetic energy flux is zero and that
detailed balance in energy transfers is respected for the equilibrium
solutions. Interestingly, similar solutions have been observed for
dissipation-less Burgers equation and Korteweg–de Vries (KdV) equation when
they were started with $\delta$-correlated fields [23].

Euler turbulence exhibits interesting nonequilibrium behaviour that are
different for 2D and 3D. The nonequilibrium states are obtained for ordered
initial condition, as in Cichowlas  _et al._ [12]. Customarily, 3D
conservative systems thermalize when they are started with a nonequilibrium
configuration. As shown by Cichowlas  _et al._ [12], 3D Euler turbulence
thermalizes in this spirit. Regarding 2D Euler turbulence, for coherent
velocity field as an initial condition, $E(k)$ differs significantly from Eq.
(2). For example, Fox and Orszag [24] reported deviations from Eq. (2) at
small wavenumbers for enstrophy-dominated 2D Euler turbulence. Seyler  _et
al._ [25] observed large vortex structures, similar to those in a discrete
vortex system [26]. Verma and Chatterjee [27] showed that for several ordered
initial condition, 2D Euler remains out of equilibrium throughout its
evolution. These works indicate 2D Euler turbulence is out of equilibrium.
Even though, 2D Euler turbulence is in nonequilibrium state, Robert and
Sommeria [28], and Bouchet and Venaille [29] analyzed such structures in the
framework of equilibrium statistical mechanics.

Euler equation is time reversible due to an absence of viscous dissipation.
This is another reason why thermodynamic entropy of Euler turbulence is
constant. Note, however, that the Euler equation exhibits irreversibility due
to its chaotic nature. The hydrodynamic entropy, to be describe in this paper,
captures the irreversibility and disorder of Euler turbulence quite well.

The structure of the paper is as follows. In Sections II we review the
analytical works on equilibrium states of Euler turbulence. Section III
summarizes the hydrodynamic entropy formalism, whereas Section IV covers the
numerical procedure for solving Euler equations. Sections V covers the
numerical results on equilibrium properties of Euler turbulence, whereas
Sections VI and VII cover the nonequilibrium properties of 3D and 2D Euler
turbulence respectively. We conclude in Section VIII.

##  II GOVERNING EQUATIONS AND EQUILIBRIUM STATES

The incompressible Euler equation, which is the hydrodynamic equation with
zero external force and zero viscosity, is

$\displaystyle\partial_{t}\mathbf{u}+(\mathbf{u}\cdot\nabla)\mathbf{u}$
$\displaystyle=$ $\displaystyle-\nabla p,$ (3)
$\displaystyle\nabla\cdot\mathbf{u}$ $\displaystyle=$ $\displaystyle 0,$ (4)

where ${\bf u}$ and $p$ are the velocity and pressure fields respectively. In
Fourier space, Eq. (3) gets transformed to [30, 7]

$\frac{d}{dt}\hat{\bf u}({\bf k})+i\sum_{\bf r}{\bf k}\cdot\hat{\bf u}({\bf
q})\hat{\bf u}({\bf r})=-i{\bf k}p({\bf k})$ (5)

where ${\bf q=k-r}$, and $\hat{\bf u}({\bf k}),p({\bf k})$ are the Fourier
transforms of $\mathbf{u}({\bf x})$ and $p({\bf x})$ respectively. The modal
energy corresponding to ${\bf\hat{u}(k)}$ is defined as $E({\bf
k})=|{\bf\hat{u}(k)}|^{2}/2$, whose evolution equation is [31, 32]

$\frac{d}{dt}E({\bf k})=\sum_{\bf p}\Im[\\{{\bf
k\cdot\hat{u}(q)\\}\\{\hat{u}(r)\cdot\hat{u}^{*}(k)}\\}]=T({\bf k}),$ (6)

where $\Im[.]$ represents the imaginary part of the argument, and $T({\bf k})$
represents the nonlinear energy transfer to Fourier mode ${\bf\hat{u}(k)}$
from all other modes. Note that the pressure term of the Euler equation does
not contribute to the evolution of $E({\bf k})$.

Under a statistical steady state, the average rate of change of kinetic energy
is zero, i.e., $d\left\langle E({\bf k})\right\rangle/dt=0$. Hence, Eq. (6)
yields the following relation for all ${\bf k}$’s,

$\left\langle T({\bf k})\right\rangle=0.$ (7)

The energy flux $\Pi(k)$ is defined as the net nonlinear energy transfer from
all the modes residing inside the sphere of radius $k$ to the modes outside
the sphere. In terms of $T({\bf k})$, the energy flux $\Pi(k)$ is defined as
[31, 32]

$\Pi(k)=-\sum_{|{\bf}k^{\prime}|=0}^{k}T({\bf k^{\prime}}).$ (8)

Therefore, under statistical steady state, [23]

$\left\langle\Pi(k)\right\rangle=0.$ (9)

Thus, for a steady state, the energy flux for Euler equation is zero. The
condition $\left\langle T({\bf k})\right\rangle=0$ for all ${\bf k}$’s
indicates that the system respects detailed balance in energy transfers. Thus,
we can claim that the system is under equilibrium, as in thermodynamic
equilibrium. This general condition is expected to work for conservative
systems. Later in this paper we show that Eq. (9) is satisfied for 3D Euler
turbulence, but not for 2D Euler turbulence.

Liouville’s theorem is often invoked to characterize the statistical
properties of a system in equilibrium. For example, Gibbs measure is an
invariant measure of a system satisfying Liouville’s theorem [33, 18]. On this
basis, Lee [10] and Kraichnan [9] derived that $E(k)\sim k^{2}$ for nonhelical
Euler turbulence, and

$E(k)\sim\frac{k^{2}}{\beta^{2}-\gamma^{2}k^{2}}$ (10)

for helical turbulence, where $\gamma$ and $\beta$ are constants. Note,
however, that Eqs. (6, 7, 8), which are related to the flux formalism, provide
alternative framework for classifying equilibrium systems.

In the following discussion, we derive equilibrium energy distribution for 3D
Euler turbulence using helical basis [34, 32]. The helical Fourier modes
associated with wavenumber ${\bf k}$ are [34, 35, 32]

$u_{+}(\mathbf{k}),\quad u_{-}(\mathbf{k}).$ (11)

We denote the corresponding energies as $E_{+}(\bf k)$ and $E_{-}(\bf k)$
respectively. Under equilibrium, $E_{+}$ is distributed among all the positive
helicity modes, $E_{+}(\bf k)$. Similar distribution takes place for
$E_{-}({\bf k})$. The average values of $E_{+}(\bf k)$ and $E_{-}(\bf k)$ are
$\left\langle E_{+}(\bf k)\right\rangle$ and $\left\langle E_{-}(\bf
k)\right\rangle$ respectively. Note that

$\displaystyle E({\bf k})$ $\displaystyle=$ $\displaystyle E_{+}({\bf
k})+E_{-}({\bf k}),$ (12) $\displaystyle H({\bf k})$ $\displaystyle=$
$\displaystyle k(E_{+}({\bf k})-E_{-}({\bf k})),$ (13)

where $E({\bf k})$ and $H({\bf k})$ are modal energy and modal kinetic
helicity respectively. Hence, $\left\langle E_{+}(\bf
k)\right\rangle\neq\left\langle E_{-}({\bf k})\right\rangle$ for helical
flows.

For a given k, $E_{+}(\bf k)$ and $E_{-}(\bf k)$ fluctuate around their
respective mean values. We denote the fluctuations in $E_{+}(\bf k)$ and
$E_{-}(\bf k)$ using $e_{+}$ and $e_{-}$ respectively. The probability for
observance of $(e_{+},e_{-})$ is denoted by $P(e_{+},e_{-})$. Using
conservation laws, we arrive at the following constraints:

$\displaystyle\int_{0}^{\infty}P(e_{+},e_{-})de_{+}de_{-}$ $\displaystyle=$
$\displaystyle 1$ (14)
$\displaystyle\int_{0}^{\infty}e_{\pm}P(e_{+},e_{-})de_{+}de_{-}$
$\displaystyle=$ $\displaystyle\left\langle E_{\pm}(\bf k)\right\rangle,$ (15)

To obtain the distribution for the fluctuations $e_{\pm}$, we extremize the
following function:

$\displaystyle S$ $\displaystyle=$ $\displaystyle-\sum
P(e_{+},e_{-})\log(P(e_{+},e_{-}))$ (16) $\displaystyle-\alpha(\sum
P(e_{+},e_{-})-1)$ $\displaystyle-\beta_{+}(\sum
P(e_{+},e_{-})e_{+}-\left\langle E_{+}(\bf k)\right\rangle)$
$\displaystyle-\beta_{-}(\sum P(e_{+},e_{-})e_{-}-\left\langle E_{-}(\bf
k)\right\rangle),$

where $\sum$ is a shorthand for $\int_{0}^{\infty}de_{+}de_{-}$. Taking
derivative of $S$ with relative to $P(e_{+},e_{-})$ yields

$-\log(P(e_{+},e_{-}))-1-\alpha-\beta_{+}e_{+}-\beta_{-}e_{-}=0.$ (17)

Hence,

$P(e_{+},e_{-})=\exp[-(\alpha+1)]\exp(-\beta_{+}e_{+}-\beta_{-}e_{-}).$ (18)

Substitution of Eq. (18) in Eq. (14) yields

$\exp[-(\alpha+1)]=\beta_{+}\beta_{-}.$ (19)

Hence,

$P(e_{+},e_{-})=\beta_{+}\beta_{-}\exp(-\beta_{+}e_{+}-\beta_{-}e_{-}).$ (20)

Substitution of the above in Eq. (15) yields

$\left\langle E_{\pm}(\bf k)\right\rangle=1/\beta_{\pm}.$ (21)

Note that $\beta_{\pm}$ can be a function of $k$, and that

$e_{\pm}=\frac{1}{2}(e\pm h/k),$ (22)

where $e$ and $h$ are the kinetic energy and kinetic helicity respectively.
Substitution of the above expressions in Eq. (20) yields

$P(e,h)=(\beta^{2}-\gamma^{2}k^{2})\exp(-\beta e-\gamma h),$ (23)

where

$\beta_{\pm}=\beta\pm\gamma k.$ (24)

Using Eqs. (12, 13, 21, 24) we derive

$\displaystyle\left\langle E(\bf k)\right\rangle$ $\displaystyle=$
$\displaystyle\frac{2\beta}{(\beta^{2}-\gamma^{2}k^{2})},$ (25)
$\displaystyle\left\langle H(\bf k)\right\rangle$ $\displaystyle=$
$\displaystyle-\frac{2\gamma k^{2}}{(\beta^{2}-\gamma^{2}k^{2})}.$ (26)

A nonhelical flow is a special case with $\gamma=0$ or $\beta_{+}=\beta_{-}$.
For this case, $\left\langle E(\bf k)\right\rangle=2/\beta$ and $\left\langle
H(\bf k)\right\rangle=0$, and the velocity components in Craya-Herring basis
are independently random. However, for helical case, correlations develop
among the Fourier modes. As is evident from Eq. (21, 24, 25), the total
kinetic energy $E$, as well as $E_{\pm}$, are not equipartitioned among all
the Fourier modes. Also note that the shell spectra for the energy and kinetic
helicity are

$\displaystyle\left\langle E(k)\right\rangle=4\pi k^{2}\left\langle E({\bf
k})\right\rangle,$ (27) $\displaystyle\left\langle H(k)\right\rangle=4\pi
k^{2}\left\langle H({\bf k})\right\rangle.$ (28)

In our future discussion, we will drop the angular brackets for $E(k)$ and
$E_{\pm}(k)$. Also note that for nonhelical 3D Euler turbulence,

$E(k)=\frac{8\pi k^{2}}{\beta}.$ (29)

We will compare the above expression with our numerical data.

The equilibrium energy spectrum for 2D Euler turbulence has been derived
following similar lines as above [9, 8]. Using the fact that the energy and
enstropy are conserved for 2D Euler turbulence, one can derive that

$P(e)=(\beta+\gamma k^{2})\exp[-(\beta+\gamma k^{2})e],$ (30)

and

$E({\bf k})=\frac{1}{\beta+\gamma k^{2}}.$ (31)

Hence, the shell spectrum for the energy for 2D Euler turbulence is $E(k)=2\pi
kE({\bf k})=2\pi k/(\beta+\gamma k^{2})$.

The connection mentioned above between the $\delta$-correlated velocity field
(white noise) and $k^{2}$ spectrum provides a hint that we should choose white
noise as the initial condition for the equilibrium solution of Euler
turbulence. We follow this strategy in the present paper.

##  III Hydrodynamic Entropy

The thermodynamic entropy of Euler turbulence is constant due to an absence of
viscosity [36]. However, the disorder in Euler turbulence varies with time for
nonequilibrium scenario. Hence, the thermodynamic entropy is not suitable for
quantifying order in Euler turbulence. Therefore, Verma and Chatterjee [27]
defined “hydrodynamic entropy” using Shannon’s formula [37]. They postulated
that the probability of occurrence of a Fourier mode with wavenumber ${\bf k}$
is $p_{\bf k}=E({\bf k})/E$, where $E({\bf k})$ is the modal energy, and $E$
is the total energy. In terms of these quantities, the hydrodynamic entropy of
the flow was defined as [27]

$S=-\sum_{\bf k}p_{\bf k}\log_{2}(p_{\bf k}).$ (32)

The above entropy can be used to quantify hydrodynamic order of a snapshot of
a fluid flow. Note that hydrodynamic entropy differs from the thermodynamic
entropy, which depends on temperature and volume of the system.

##  IV Numerical Details

We perform pseudo-spectral simulation [38, 39] of Euler flow using our code
TARANG [40, 41]. We simulate both equilibrium and nonequilibrium
configurations of Euler turbulence. A key issue in Euler turbulence simulation
is energy conservation, which is not conserved while using standard time
stepping schemes, e.g., Runge-Kutta method. In this paper, we employ position-
extended Forest-Ruth-like (PEFRL) scheme [42, 43] (see Appendix A) that
conserves energy to high precision. We also dealise the code using two-third
rule [38, 39].

For studying nonequililbrium behaviour of Euler turbulence, we employ large-
scale flow structures as initial condition. For example, Taylor-Green vortices
are employed for 3D Euler equation [12]. However, we employ
$\delta$-correlated velocity field as initial condition to study equilibrium
properties of Euler turbulence.

We implement the above random initial condition using Craya-Herring basis [44,
45, 35, 32] whose unit vectors for a wavenumber ${\bf k}$ are

$\displaystyle\hat{e}_{1}(\textbf{k})$ $\displaystyle=$
$\displaystyle(\hat{k}\times\hat{n})/|\hat{k}\times\hat{n}|,$ (33)
$\displaystyle\hat{e}_{2}(\textbf{k})$ $\displaystyle=$
$\displaystyle\hat{k}\times\hat{e}_{1}(\textbf{k}),$ (34)

where $\hat{n}$ is chosen as any direction, and $\hat{k}$ is the unit vector
along k. In this basis, the 3D incompressible velocity field is

${\bf
u}(\textbf{k})=u_{1}(\textbf{k})\hat{e}_{1}(\textbf{k})+u_{2}(\textbf{k})\hat{e}_{2}(\textbf{k}),$
(35)

while 2D incompressible velocity field is

${\bf u}(\textbf{k})=u_{1}(\textbf{k})\hat{e}_{1}(\textbf{k}).$ (36)

For simulating 3D nonhelical flows (zero kinetic helicity), we start with
$u_{1}({\bf k})=0$ and $u_{2}({\bf k})=\sqrt{2E/M}\exp(i\phi_{2}({\bf k}))$,
where $E$ is the total kinetic energy, $M$ is the total number of dealiased
modes, and the phase $\phi_{2}({\bf k})$ is chosen to be a random number from
uniform distribution in a band of $[0,2\pi]$. For the 2D simulation, we take
$u_{1}({\bf k})=\sqrt{2E/M}\exp(i\phi_{1}({\bf k}))$ with random phase for
$\phi_{1}({\bf k})$. For simulating helical Euler turbulence with kinetic
helicity $H$, we choose $|u_{1}({\bf k})|=|u_{2}({\bf k})|=\sqrt{E/M}$ and
random $\phi_{2}({\bf k})$ from a uniform distribution in $[0,2\pi]$. The
other phase is $\phi_{1}({\bf k})=\phi_{2}({\bf k})-\sin^{-1}\sigma_{c}({\bf
k})$, where $\sigma_{c}({\bf k})=\Re[{\bf u^{*}({\bf k})}\cdot\bm{\omega}({\bf
k)}]/(k|\textbf{u}(\textbf{k})|^{2})$. We use $\sigma_{c}({\bf k})$ to control
$H$ [46]. Note that the total number of active dealiased modes $M=(2N/3)^{3}$
for 3D and $M=(2N/3)^{2}$ for 2D, where $N$ is the number of grid points in
each direction.

##  V EQUILIBRIUM Behaviour of Euler turbulence

In this section, we report the equilibrium nature of 2D and 3D Euler
turbulence. We discuss the randomness of the flows using real space density
plots of vorticity, the probability distribution functions, as well as energy
spectra and fluxes.

###  V.1 For 3D Euler equation

We perform spectral simulation using $\delta$-correlated velocity field as an
initial condition. See Sec. IV for details. We employ $128^{3}$ grid and PEFRL
scheme for time stepping with $dt=10^{-4}$. We perform two simulations, one
without kinetic helicity, and the other with kinetic helicity. The total
kinetic energy is 0.019 for the nonhelical (without kinetic helicity) flow.
For the helical flow, the total kinetic energy and kinetic helicity are 0.038
and 2.8 respectively. We perform these simulations up to 25 and 30 eddy
turnover time units respectively.

For the nonhelical Euler turbulence, the total kinetic energy is conserved to
13 significant digits ($E=0.019\pm 9\times 10^{-13}$), while the kinetic
helicity remains negligible throughout (rms value of $6\times 10^{-12}$). For
the helical turbulence, we choose $\sigma_{c}(\textbf{k})=0.9$ for all k’s to
inject significant kinetic helicity. Here, the total kinetic energy and total
kinetic helicity for the run are $0.038\pm 7\times 10^{-12}$ and $2.8\pm
7\times 10^{-9}$ respectively; thus, they remain conserved throughout the run.

Now we describe the equilibrium properties of Euler turbulence using the
results of the above runs. We observe that the nonhelical flow remains random,
as in white noise, at all times. Note, however, that the amplitudes and phases
of all the modes vary randomly with time. In Fig. 1(a,b), we exhibit the
density plots of the perpendicular vorticity components of the flows in the
horizontal and vertical mid planes for a snapshot. The plots clearly
demonstrate the random nature of flow.

Figure 1: (Color online) For 3D nonhelical flow with $\delta$-correlated
velocity as initial condition: Density plots of the perpendicular component of
vorticity on (a) horizontal midplane ($z=\pi$), and (b) vertical midplane
($y=\pi$). These plots indicate random velocity configuration.

We substantiate the randomness of the flows by computing the probability
distribution function (PDFs) of the magnitude of the real-space velocity field
($u$) of a snapshot, and test whether it obeys Maxwell-Boltzmann distribution,
which for a 3D flow is

$P(u)=\sqrt{2/\pi}~{}a^{-3}u^{2}\exp(-u^{2}/2a^{2}),$ (37)

where $a$ is the scale parameter. The numerical $P(u)$ for nonhelical 3D run,
which is exhibited in Fig. 2, matches quite accurately with Eq. (37) with
$a=0.11$. Hence, we claim that the velocity field of Euler turbulence is as
random as the velocity distribution of gas molecules in thermodynamics.

Figure 2: (Color online) For the nonhelical Euler simulations: Probability
distribution functions ($P(u)$) of the velocity magnitude for 2D (red squares)
and 3D (blue circles) real-space flows. The numerical PDFs match closely with
Maxwell-Boltzmann distribution for 3D and for 2D (black curves) thermodynamic
systems.

Next, we compute the energy spectrum and flux for the nonhelical simulation.
As shown in Fig. 3, the energy flux is zero (apart from fluctuations).
Besides, the normalized energy spectrum, $E(k)/k^{2}$, is flat around
$8\pi/\beta$ with $\beta=0.64\times 10^{8}$ [see Eq. (29)]. Hence, we claim
that $E(k)$ varies as $k^{2}$ for the whole range of wavenumbers. These
results are consistent with the $\delta$-correlated (white noise) nature of
the real-space velocity field, thus validating the predictions of absolute
equilibrium theory [9, 10]. The vanishing energy flux follows from the fact
that $\delta$-correlated velocity field has zero triple correlation. In
addition, using field-theoretic arguments, Verma [47, 32, 48] has shown that
the equipartitioned Fourier modes yield zero kinetic energy flux.

Figure 3: (Color online) Plot of normalized energy spectra $E(k)/(k^{2})$ for
the nonhelical simulation [blue curve] and the helical run [magenta curve].
The numerical plots match reasonably well with Kraichnan’s predictions [Eq.
(25)], which are shown as dashed lines. The inset shows the energy fluxes for
the these runs using the same colour.

The probability $P(e,h)$ of 3D Euler turbulence is given by Eq. (23). Since
$H=0$ for the nonhelical run, we set $\gamma=0$ in the equation and compute
$P(e)$. Unfortunately, the numerical data is not sufficient for the
computation of $P(e)$ for a given $E({\bf k})$, hence we sample all $E({\bf
k})$’s of a snapshot to compute $P(e)$. Figure 4 exhibits thus computed
$P(e)$. We observe that the numerical data is described by Eq. (23) quite well
with $\beta=0.64\times 10^{8}$, which is consistent with the energy spectrum
computation.

Figure 4: (Color online) For nonhelical 3D Euler and 2D Euler runs, plots of
$P(e)$ at $t=20$. Here, $\tilde{e}=\beta e$ with $\beta=0.64\times 10^{8}$ for
3D and $0.8\times 10^{7}$ for 2D simulation.

Now, we report the energy spectra, $E_{\pm}(k),E(k)$, and energy flux for the
helical run. We observe that $E_{\pm}(k)$ are consistent with Eqs. (21,24)
with $\beta=0.7\times 10^{8}$ and $\gamma=1.1\times 10^{6}$. We illustrate
$E_{\pm}(k)$ in Fig. 5. We also compute $E(k)$ and plot it in Fig. 3. We find
that $E(k)$ matches with Eq. (25) quite well with $\beta,\gamma$ computed
above using $E_{\pm}(k)$. The deviation of $E(k)$ from $k^{2}$ at large $k$ is
consistent with Eq. (25). Also, as shown in Fig. 3, the energy flux vanishes
for the helical 3D Euler turbulence as well. Based on the above results, it is
evident that the velocity field of the helical run too is under equilibrium.
For small wavenumbers, a small deviation of $E(k)$ from Eq. (25) may be due to
an asymmetry in the energy transfers. For example,
$\textbf{u}(\textbf{k}_{0}=1)$ has no Fourier mode with wavenumber less than
$\textbf{k}_{0}$.

Figure 5: (Color online) For 3D helical simulation, plots of $k^{2}/E_{+}(k)$
(blue curve) and $k^{2}/E_{-}(k)$ (red curve) for the helical simulation. The
best fit curves, shown as black dashed and solid lines, are consistent with
Eq. (21).

###  V.2 For 2D Euler equation

We simulate 2D Euler turbulence on a $1024^{2}$ grid using the method
described in Sec. IV. We carry out the 2D simulation up to 30 time units with
$dt=10^{-3}$. The total energy and total enstrophy of the flow are conserved
with good accuracy, for example, $E=(0.058\pm 9\times 10^{-9})$ and
$E_{\Omega}=4534\pm 10^{-3}$.

For 2D Euler run too, the flow is random as a thermodynamic gas. In Figs. 2
and 4, we illustrate the PDF of $u$ and $e$, both of which show equilibrium
properties similar to 3D Euler turbulence. Note that the Maxwell-Boltzmann
distribution for a 2D flow is

$P(u)=a^{-2}u\exp(-u^{2}/2a^{2}),$ (38)

where $a=0.24$. For our simulation with $\delta$-correlated u, $\gamma\approx
0$ because $k_{\Omega}$, the centroid for enstrophy, exceeds the grid size.
Hence, $E(k)=2\pi k/\beta$, as is evident from Fig. 6. We obtain
$\beta=0.8\times 10^{7}$, which is consistent with $P(e)$ plot of Fig. 2. In
addition, the energy flux is zero for this case, indicating detailed balance
of energy transfer.

Figure 6: (Color online) For 2D Euler simulation, plot of the normalized
energy spectrum $\beta E(k)/(2\pi k)$ with $\beta=0.8\times 10^{7}$. The
energy flux is in the inset.

We also study the phase space projection of the phase space trajectory on the
$\\{\Re(u(\mathbf{q})),\Im(u(\mathbf{q}))\\}$ plane for wavenumber
$\mathbf{q}=(128,128)$. We observe a random scatter of the trajectory whose
extent is around two times
$\sqrt{\left\langle|\hat{u}(k)|^{2}\right\rangle}=0.5\times 10^{-3}$. We
believe that the trajectory will wander off to a larger distance from the
origin if we wait for a longer time. These observations clearly indicate that
2D Euler turbulence exhibits equilibrium behaviour.

Figure 7: (Color online) For 2D Euler simulation, phase space projection of
phase space trajectory on the $\\{\Re(u(\mathbf{q})),\Im(u(\mathbf{q}))\\}$
plane for $\mathbf{q}=(128,128)$. The radius of the black circle in the figure
is $\sqrt{\left\langle|\hat{u}(k)|^{2}\right\rangle}=0.5\times 10^{-3}$.

In the next two sections, we describe nonequilibrium behaviour of 3D and 2D
Euler turbulence.

##  VI Nonequilibrium behaviour and thermalization of 3D Euler turbulence

Cichowlas et al. [12] reported thermalization for 3D Euler turbulence, which
is an important topic of research in nonequilibrium statistical mechanics,
both classical and quantum. They presented a model of thermalization that
relates the energy of the thermalized modes to the transition wavenumber
between the nonequilibrium and thermal modes. Krstulovic and Brachet [49]
constructed a two-fluid model of truncated Euler equation and determined the
effective viscosity and thermal diffusion. They employed EDQNM closure and
Monte-Carlo scheme for their derivation. Note that the intermediate stage of
3D Euler turbulence is in a mix state of equilibrium (large $k$) and
nonequilibrium (intermediate $k$). The flow thermalizes after tens of eddy
turnover times.

In the following discussion, we estimate the time required for thermalization
in 3D Euler turbulence [12]. We denote the wavenumber shells in 3D Euler
turbulence as $k_{0},k_{1},...,k_{N-1},k_{N}$, and assume that the initial
condition is a large-scale vortex with wavenumber $k_{0}$, as in [12].
Nonlinear interactions transfer energy from $k_{0}$ to $k_{1}$, from $k_{1}$
to $k_{2}$, …, $k_{N-1}$ to $k_{N}$. The cascade however stops at $k=k_{N}$
where the energy piles up. After sufficiently large accumulation of energy at
$k_{N}$, the energy starts to grow at wavenumbers shell $k_{N-1}$, and then at
$k_{N-2}$, and so on. This is how the large wavenumber shells acquire $k^{2}$
spectrum, as reported by Cichowlas  _et al._ [12] and Krstulovic and Brachet
[49]. See Fig. 8 for an illustration.

Figure 8: A schematic diagram exhibiting the evolution of energy spectrum
$E(k)$ of 3D Euler turbulence during thermalization. The thin red curve
represents $E(k)$ during the early phase, while the thick red curve at an
intermediate state. At the transition wavenumber $k_{\mathrm{tr}}$, $E(k)$
changes from $k^{-5/3}$ to $k^{2}$. The two regimes, nonequilibrium and
equilibrium, are represented by yellow and blue colors respectively.

Following Kolmogorov’s theory of turbulence, the energy cascade rate to the
large-wavenumber modes can be estimated as $\epsilon_{u}=U^{3}/L\sim
U^{3}k_{0}$, where $L,U$ are the large-scale length and velocity respectively
[5, 6, 7, 8]. This energy flux accumulates at large wavenumbers and builds up
$A(t)k^{2}$ spectrum from the transition wavenumber $k_{\mathrm{tr}}$ to
$k_{\mathrm{max}}$ (see Fig. 8). Therefore, in time $t$,

$\epsilon_{u}t\sim\int_{k_{\mathrm{tr}}}^{k_{\mathrm{max}}}A(t)k^{2}dk,$ (39)

or

$U^{3}k_{0}t\sim A(t)[k_{\mathrm{max}}^{3}-(k_{\mathrm{tr}}(t))^{3}].$ (40)

Over time, $A(t)$ increases and $k_{\mathrm{tr}}(t)$ decreases. Using Eq. (40)
we can deduce the total time taken for thermalization ($T$) as follows. During
the final stage, $k_{\mathrm{tr}}\rightarrow k_{0}\ll k_{\mathrm{max}}$ and
$A(T)\sim E/N^{3}$. Hence,

$T\sim\frac{E}{N^{3}U^{3}k_{0}}k_{\mathrm{max}}^{3}\sim\frac{L}{U}$ (41)

because $k_{\mathrm{max}}\approx N/2$. Thus, a 3D Euler flow with large-scale
vortex as an initial condition is expected to thermalize in order of one eddy
turnover time [12, 50]. In practice, this process takes tens of eddy turnover
time [12].

In Euler turbulence with coherent large-scale structures, the energy in the
inertial and dissipation range is converted incoherent (random) energy at
small scales [49, 51]. In the language of statistical mechanics, the yellow
and blue regions of Fig. 8 could represent system and heat bath respectively.
In this process, coherent energy is being converted to incoherent energy, thus
giving a semblence of frictional effect [49, 13]. Verma [51] argued that
dissipation can emerge in conservative systems in a similar manner. Thus, the
nonequilibrium and equilibrium states of Euler turbulence yield valuable
insights into the thermalization process in conservative systems.

Before we end this section, we describe the hydrodynamic entropy of 3D Euler
turbulence, as reported by Verma and Chatterjee [27]. They simulated 3D Euler
turbulence on a $(2\pi)^{3}$ box with a $N^{3}$ grid, where $N=128$, using
Taylor-Green vortex ($k_{0}=1$) as an initial condition. The simulation was
run up to 180 nondimensional time units using $dt=10^{-4}$ (see Sec. IV). The
total energy for the run is $0.125$ and it is conserved up to 12 decimal
places. The time evolution of hydrodynamic entropy computed using Eq. (32) is
shown in Fig. 9. We observe that the hydrodynamic entropy increases
monotonically and asymptotes to 18, which is near the maximum possible entropy
of 18.3.

Figure 9: (Color online) For 3D Euler turbulence, hydrodynamic entropy
exhibits monotonic growth, with exponential increase in the beginning, and
saturation in the end.

##  VII Nonequilibrium behaviour of 2D Euler turbulence

In a recent work, Verma and Chatterjee [27] simulated 2D Euler turbulence with
ordered initial condition and observed nonequilibrium behaviour. In this
section, we briefly review their results.

Verma and Chatterjee [27] simulated 2D Euler turbulence on a $(2\pi)^{2}$ box
using $512^{2}$ grid. See Sec. IV for details. Two of Verma and Chatterjee
[27]’s runs are given below:

  1. 1.

Run A: The initial velocity profile is taken as $(\sin 11x\cos
11y+\eta_{x},-\cos 11x\sin 11y+\eta_{y})$, where $(\eta_{x},\eta_{y})$ is
random noise. We take $|\eta_{x}|\ll 1$ and $|\eta_{y}|\ll 1$.

  2. 2.

Run B: The initial nonzero velocity Fourier modes are ${\bf u}(1,0)=(0,1)$,
${\bf u}(0,1)=(1,0)$, and ${\bf u}(1,1)=(-i,i)$.

These runs were time advanced up to 170 and 30 turnover times
($2\pi/U_{\mathrm{rms}}$) respectively. The Runs A and B reach steady states
after 100 to 10 eddy turnover times respectively. For Runs A and B, the total
energy and enstrophy are (0.2500954, 4) and (62.17, 6) respectively, and they
are conserved to many significant digits. The initial and final states of the
two runs are shown on the top and bottom panels of Fig. 10. Here, the velocity
field is superposed over the density plots of the vorticity field. The
aysmptotic states of Runs A and B are a vortex-antivortex pair [11] and a
unidirectional flow (shear layer) respectively, which are embedded in small-
scale noisy flow.

Figure 10: (Color online) For Runs A and B of 2D Euler turbulence: (a,c) the
initial states, (b,d) the final states respectively. Here we plot the velocity
field over the density plots of the vorticity field.

In Fig. 11(a,b), we plot the averaged energy spectra and fluxes of the steady
states of the two runs. The two runs have the following energy spectra:

$\displaystyle E(k)$ $\displaystyle=$
$\displaystyle\frac{k}{-2518+237k^{2}}~{}\mathrm{RunA:for}~{}k>10,$ (42)
$\displaystyle E(k)$ $\displaystyle=$
$\displaystyle\frac{k}{-6357840+9361k^{2}}~{}\mathrm{RunB:for}~{}k>40,$

For Runs A and B, $\sqrt{\Omega/E}=15.8,1.2$. Since $\sqrt{\Omega/E}\leq
k_{\Omega}$ (centroid of enstrophy), the wavenumbers far beyond $k_{\Omega}$
are dominated by enstrophy leading to $E(k)\propto k^{-1}$, which corresponds
to an equipartition of enstrophy [52, 27].

The average energy flux $\left\langle\Pi(k)\right\rangle\approx 0$ for
intermediate and large $k$’s. But, $\left\langle\Pi(k)\right\rangle<0$ for
small $k$’s. In particular,
$\min[\left\langle\Pi(k)\right\rangle]\approx-3\times 10^{-4},-10^{-3}$ for
Runs A and B respectively. Thus, the nonzero $\Pi(k)$ breaks the detailed
balance of energy transfers. Hence, 2D Euler turbulence is out of equilibrium.
See Verma and Chatterjee [27] for more details.

Figure 11: (Color online) For the Runs A, B of 2D Euler turbulence: plots of
the averaged energy spectra, $E(k)$, and fluxes, $\Pi(k)$, of the asymptotic
states. Figure (b) exhibits $\Pi(k)$ for small $k$’s, while inset shows
$\Pi(k)$ for the whole range.

Verma and Chatterjee [27] computed the hydrodynamic entropies of the 2D Euler
flows of Runs A and B. The plots of the entropy time series are shown in Fig.
12. For Runs A and B, after initial fluctuations, the entropies decrease
exponentially to asymptotic values of 4.9 and 1.2 respectively. These values
are smaller than the maximum possible value, which is $\log_{2}(M)\approx
16.5$, where $M\approx\pi(512/3)^{2}$. Thus, Verma and Chatterjee [27] showed
that the hydrodynamic entropy of 2D Euler turbulence, an isolated system,
decreases with time for a significant duration. Hence, 2D Euler turbulence is
a rare isolated system that exhibits evolution from disorder to order.

Figure 12: (Color online) Plots of the temporal evolution of hydrodynamic
entropies for Runs A, B (a, b) of 2D Euler turbulence respectively. In each
case, after initial transients, the entropy decreases with time and asymptotes
to an approximate constant value.

##  VIII DISCUSSIONS and CONCLUSIONS

In this paper, we review the equilibrium and nonequilibrium properties of 2D
and 3D Euler turbulence. For $\delta$-correlated velocity field as an initial
condition, both 2D and 3D Euler turbulence exhibit equilibrium behaviour,
predicted by Lee [10] and Kraichnan [9]. However, for ordered initial
condition, 3D Euler turbulence evolves from order to disorder. The above
thermalization of 3D Euler turbulence follows a generic path for energy
conserving system. Here, the nonlinear energy transfer $T(k,t)\rightarrow 0$
as $t\rightarrow\infty$. The energy flux too vanishes asymptotically.

The thermodynamic entropy of Euler turbulence remains constant. Hence, Verma
and Chatterjee [27] constructed hydrodynamic entropy to quantify the
variations in order of Euler turbulence. They showed that the hydrodynamic
entropy of 3D Euler turbulence increases monotonically with time.

Interestingly, we can extrapolate the above thermalization process to Navier-
Stokes equation that includes viscous dissipation, as well as external forcing
at large scales. A fluid is composed of molecules whose total energy is
conserved. However, we can separate the system into two parts: (a) coherent
flow, which is described by the flow equation under continuum approximation,
and (b) random or thermal motion of the molecules in the co-moving frame of
the flow. In Kolmogorov’s picture of turbulence, hydrodynamic range of scales
includes forcing, inertial, and dissipation range, while the random motion of
molecules is described by thermodynamics. These scales are exhibited by
different colors in Fig. 13. The transition wavenumber between the
hydrodynamic and thermodynamic ranges may be approximated by the Kolmogorov
wavenumber, $k_{d}=(\epsilon/\nu^{3})^{1/4}$, where $\nu$ is the kinematic
viscosity. The above picture is similar to findings in recent works by Shukla
et al. [53] and Bandak  _et al._ [54].

Figure 13: A schematic diagram exhibiting the energy spectrum $E(k)$ of 3D
forced hydrodynamic turbulence (includes viscous effects). The yellow and
green regions represent the forcing, inertial, and dissipation ranges of
hydrodynamic description, while the blue region represents the thermal motion
of molecules. The former regions are out of equilibrium, while the latter one
is in quasi-equilibrium.

The above arguments can be extended to quantum systems, at least to
superfluids and Bose-Einstein gas. Many experiments and numerical simulations
of such systems yield Kolmogorov-like $k^{-5/3}$ spectrum (Krstulovic and
Brachet [50], Madeira  _et al._ [55], Fonda  _et al._ [56], Shukla  _et al._
[57], Skrbek and Sreenivasan [58], and references therein) that requires
dissipation at small scales. Small-scale dissipation in such systems are
attributed to interactions of condensate with thermal clouds, or to decay of
vortical motion into phonon excitations (Barenghi  _et al._ [59] and
references therein). This feature may appear odd because quantum systems are
energy conserving. But, the multiscale energy transfer in Euler turbulence
provides an interesting framework to introduce quantum dissipation and
thermalization [63, 60, 61]. This framework could be an alternative to other
approaches that are typically based on modeling the interactions between the
system and the heat bath [60, 61, 62].

In contrast, 2D Euler turbulence remains out of equilibrium, at least for
several ordered initial conditions. This is contrary to what we expect for a
conservative system with a large degree of freedom. Here, the nonlinear energy
transfer and energy flux do not vanish asymptotically. Thus, for such a
scenario, the final states of 2D Euler turbulence are not stationary.
Interestingly, the hydrodynamic entropy of 2D Euler turbulence decreases in
the asymptotic regime. Thus, 2D Euler turbulence is a unique isolated system
that exhibits evolution form disorder to order.

###### Acknowledgements.

The authors thank Arul Lakshminarayan, Stephan Fauve, Marc Brachet, Alex
Alexakis, Hal Takasi, Anurag Gupta, Saikat Ghosh, Franck Plunian, Rodion
Stepanov, and Giorgio Krstulovic for useful discussions. We also thank the
organizers of APPC15 for hosting the meeting. This work is supported by the
project 6104-1 from the Indo-French Centre for the Promotion of Advanced
Research (IFCPAR/CEFIPRA). Soumyadeep Chatterjee is supported by INSPIRE
fellowship (IF180094) from Department of Science & Technology, India.

##  Appendix A Modified PEFRL algorithm

Omelyan  _et al._ [42] extended Forest-Ruth (FR) algorithm [43] to solve
differential equations associated with conservative systems. Their scheme is
called position-extended Forest-Ruth-like (PEFRL). The PEFRL algorithm was
devised for molecular dynamics simulations. We modify it to simulate Euler
flow. The evolution of velocity of an inviscid fluid is an explicit function
of ${\bf u}$. In the following, we write equations for $u_{i}$, the $i$th
component of velocity:

$\frac{du_{i}}{dt}=R(u_{i}).$ (44)

Various steps involved in time advance from $t$ to $t+h$ are as follows (see
Fig. 14):

$\displaystyle(1)$ $\displaystyle u_{i}=u_{i}(t)+\xi hR(u_{i}(t)),$ (45)
$\displaystyle(2)$
$\displaystyle\widetilde{u_{i}}=u_{i}(t)+(1-2\lambda)\frac{h}{2}R(u_{i}),$
$\displaystyle(3)$ $\displaystyle u_{i}=u_{i}+\chi hR(\widetilde{u_{i}}),$
$\displaystyle(4)$ $\displaystyle\widetilde{u_{i}}=\widetilde{u_{i}}+\lambda
hR(u_{i}),$ $\displaystyle(5)$ $\displaystyle
u_{i}=u_{i}+(1-2(\chi+\xi))hR(\widetilde{u_{i}}),$ $\displaystyle(6)$
$\displaystyle\widetilde{u_{i}}=\widetilde{u_{i}}+\lambda hR(u_{i}),$
$\displaystyle(7)$ $\displaystyle u_{i}=u_{i}+\chi hR(\widetilde{u_{i}}),$
$\displaystyle(8)$
$\displaystyle\widetilde{u_{i}}=\widetilde{u_{i}}+(1-2\lambda)\frac{h}{2}R(u_{i}),$
$\displaystyle(9)$ $\displaystyle u_{i}(t+h)=u_{i}+\xi hR(\widetilde{u_{i}}).$

Here, $\widetilde{u}$ is an intermediate velocity placeholder, whereas $\xi$,
$\lambda$, and $\chi$ are constants. Omelyan  _et al._ [42] computed the
optimized values of these constants as

$\displaystyle\xi$ $\displaystyle=$ $\displaystyle 0.1786178958448091$
$\displaystyle\lambda$ $\displaystyle=$ $\displaystyle-0.2123418310626054$
$\displaystyle\chi$ $\displaystyle=$ $\displaystyle-0.06626458266981849.$ (46)

Figure 14: (Color online) Schematic diagram of modified PEFRL algorithm.

## References

  * Zwanzig [2001] R. Zwanzig,  _Nonequilibrium Statistical Mechanics_ (Oxford University Press, 2001).
  * Livi and Politi [2017] R. Livi and P. Politi,  _Nonequilibrium Statistical Physics: A Modern Perspective_ (Cambridge University Press, 2017).
  * Landau and Lifshitz [1980] L. D. Landau and E. M. Lifshitz,  _Statistical Physics_ , 3rd ed., Course of Theoretical Physics (Elsevier, Oxford, 1980).
  * Reif [1965] F. Reif,  _Fundamentals of Statistical and Thermal Physics_ (McGraw-Hill, New York, 1965).
  * Kolmogorov [1941a] A. N. Kolmogorov, Dissipation of Energy in Locally Isotropic Turbulence, Dokl Acad Nauk SSSR 32, 16 (1941a).
  * Kolmogorov [1941b] A. N. Kolmogorov, The local structure of turbulence in incompressible viscous fluid for very large Reynolds numbers, Dokl Acad Nauk SSSR 30, 301 (1941b).
  * Frisch [1995] U. Frisch,  _Turbulence: The Legacy of A. N. Kolmogorov_ (Cambridge University Press, Cambridge, 1995).
  * Lesieur [2008] M. Lesieur,  _Turbulence in Fluids_ (Springer-Verlag, Dordrecht, 2008).
  * Kraichnan [1973] R. H. Kraichnan, Helical turbulence and absolute equilibrium, J. Fluid Mech. 59, 745 (1973).
  * Lee [1952] T. D. Lee, On some statistical properties of hydrodynamical and magneto-hydrodynamical fields, Quart. Appl. Math. 10, 69 (1952).
  * Onsager [1949] L. Onsager, Statistical hydrodynamics, Il Nuovo Cimento 6, 279 (1949).
  * Cichowlas  _et al._ [2005] C. Cichowlas, P. Bonaïti, F. Debbasch, and M. E. Brachet, Effective Dissipation and Turbulence in Spectrally Truncated Euler Flows, Phys. Rev. Lett. 95, 264502 (2005).
  * Krstulovic  _et al._ [2009] G. Krstulovic, P. D. Mininni, M. E. Brachet, and A. G. Pouquet, Cascades, thermalization, and eddy viscosity in helical Galerkin truncated Euler flows, Phys. Rev. E 79, 889 (2009).
  * Dallas  _et al._ [2015] V. Dallas, S. Fauve, and A. Alexakis, Statistical Equilibria of Large Scales in Dissipative Hydrodynamic Turbulence, Phys. Rev. Lett. 115, 204501 (2015).
  * Alexakis and Brachet [2019] A. Alexakis and M. E. Brachet, On the thermal equilibrium state of large-scale flows, J. Fluid Mech. 872, 594 (2019).
  * Alexakis and Brachet [2020] A. Alexakis and M. E. Brachet, Energy fluxes in quasi-equilibrium flows, J. Fluid Mech. 884, 87 (2020).
  * Davis  _et al._ [2001] M. J. Davis, S. A. Morgan, and K. Burnett, Simulations of Bose Fields at Finite Temperature, Phys. Rev. Lett. 87, 160402 (2001).
  * Majda and Timofeyev [2000] A. J. Majda and I. Timofeyev, Remarkable statistical behavior for truncated Burgers-Hopf dynamics, PNAS 97, 12413 (2000).
  * Ray  _et al._ [2011] S. S. Ray, U. Frisch, S. V. Nazarenko, and T. Matsumoto, Resonance phenomenon for the Galerkin-truncated Burgers and Euler equations, Phys. Rev. E 84, 016301 (2011).
  * Orszag [1973] S. A. Orszag, Lectures on the statistical theory of turbulence in fluid dynamics, in  _Les Houches Summer School of Theoretical Physics_ , edited by R. Balian and J. L. Peube (1973) p. 235.
  * Kraichnan and Montgomery [1980] R. H. Kraichnan and D. C. Montgomery, Two-dimensional turbulence, Rep. Prog. Phys. 43, 547 (1980).
  * Verma [2020] M. K. Verma, Boltzmann equation and hydrodynamic equations: their equilibrium and non-equilibrium behaviour, Phil. Trans. R. Soc. A. 378, 20190470 (2020).
  * Verma  _et al._ [2022] M. K. Verma, S. Chatterjee, A. Sharma, and A. Mohapatra, Equilibrium states of Burgers and Korteweg–de Vries equations, Phys. Rev. E 105, 034121 (2022).
  * Fox and Orszag [1973] D. G. Fox and S. A. Orszag, Inviscid dynamics of two-dimensional turbulence, Phys. Fluids 16, 169 (1973).
  * Seyler  _et al._ [1975] C. E. Seyler, Y. Salu, D. Montgomery, and G. Knorr, Two-dimensional turbulence in inviscid fluids or guiding center plasmas, Physics of Fluids 18, 803 (1975).
  * Joyce and Montgomery [1973] G. Joyce and D. Montgomery, Negative temperature states for the two-dimensional guiding-centre plasma, J. Plasma Phys. 10, 107 (1973).
  * Verma and Chatterjee [2022] M. K. Verma and S. Chatterjee, Hydrodynamic Entropy and Emergence of Order in Two-dimensional Euler Turbulence [under review in Phys. Rev. Fluids], arXiv (2022), 2210.06445v1 .
  * Robert and Sommeria [1991] R. Robert and J. Sommeria, Statistical equilibrium states for two-dimensional flows, J. Fluid Mech. 229, 291 (1991).
  * Bouchet and Venaille [2012] F. Bouchet and A. Venaille, Statistical mechanics of two-dimensional and geophysical flows, Phys. Rep. 515, 227 (2012).
  * Leslie [1973] D. C. Leslie,  _Developments in the theory of turbulence_ (Clarendon Press, Oxford, 1973).
  * Kraichnan [1959] R. H. Kraichnan, The structure of isotropic turbulence at very high Reynolds numbers, J. Fluid Mech. 5, 497 (1959).
  * Verma [2019a] M. K. Verma,  _Energy transfers in Fluid Flows: Multiscale and Spectral Perspectives_ (Cambridge University Press, Cambridge, 2019).
  * Huang [1987] K. Huang,  _Statistical Mechanics_ , 2nd ed. (Wiley, New York, 1987).
  * Waleffe [1992] F. Waleffe, The nature of triad interactions in homogeneous turbulence, Phys. Fluids A 4, 350 (1992).
  * Sagaut and Cambon [2018] P. Sagaut and C. Cambon,  _Homogeneous turbulence dynamics_ , 2nd ed. (Cambridge University Press, Cambridge, 2018).
  * Landau and Lifshitz [1987] L. D. Landau and E. M. Lifshitz,  _Fluid Mechanics_ , 2nd ed., Course of Theoretical Physics (Elsevier, Oxford, 1987).
  * Shannon [1948] C. E. Shannon, A mathematical theory of communication, Bell Labs Tech. J. 27, 379 (1948).
  * Boyd [2003] J. P. Boyd,  _Chebyshev and Fourier Spectral Methods_ , 2nd ed. (Dover Publications, New York, 2003).
  * Canuto  _et al._ [1988] C. Canuto, M. Y. Hussaini, A. Quarteroni, and T. A. Zang,  _Spectral Methods in Fluid Dynamics_ (Springer-Verlag, Berlin Heidelberg, 1988).
  * Verma  _et al._ [2013] M. K. Verma, A. G. Chatterjee, R. K. Yadav, S. Paul, M. Chandra, and R. Samtaney, Benchmarking and scaling studies of pseudospectral code Tarang for turbulence simulations, Pramana-J. Phys. 81, 617 (2013).
  * Chatterjee  _et al._ [2018] A. G. Chatterjee, M. K. Verma, A. Kumar, R. Samtaney, B. Hadri, and R. Khurram, Scaling of a Fast Fourier Transform and a pseudo-spectral fluid solver up to 196608 cores, J. Parallel Distrib. Comput. 113, 77 (2018).
  * Omelyan  _et al._ [2002] I. P. Omelyan, I. M. Mryglod, and R. Folk, Optimized Forest-Ruth- and Suzuki-like algorithms for integration of motion in many-body systems, Computer Physics Communications 146, 188 (2002).
  * Forest and Ruth [1990] E. Forest and R. D. Ruth, Fourth-order symplectic integration, Physica D 43, 105 (1990).
  * Craya [1958] A. Craya,  _Contribution à l’analyse de la turbulence associée à des vitesses moyennes_ , Ph.D. thesis, Université de Granoble (1958).
  * Herring [1974] J. R. Herring, Approach of axisymmetric turbulence to isotropy, Phys. Fluids 17, 859 (1974).
  * Sadhukhan  _et al._ [2019] S. Sadhukhan, M. K. Verma, R. Stepanov, F. Plunian, and R. Samtaney, Kinetic helicity and enstrophy transfers in helical hydrodynamic turbulence, Phys. Rev. Fluids 4, 84607 (2019).
  * Verma [2004] M. K. Verma, Statistical theory of magnetohydrodynamic turbulence: recent results, Phys. Rep. 401, 229 (2004).
  * Verma [2022] M. K. Verma, Variable energy flux in turbulence, Journal of Physics A: Mathematical and Theoretical 55, 013002 (2022).
  * Krstulovic and Brachet [2008] G. Krstulovic and M. E. Brachet, Two-fluid model of the truncated Euler equations, Physica D 237, 2015 (2008).
  * Krstulovic and Brachet [2011] G. Krstulovic and M. E. Brachet, Dispersive Bottleneck Delaying Thermalization of Turbulent Bose-Einstein Condensates, Phys. Rev. Lett. 106, 115303 (2011).
  * Verma [2019b] M. K. Verma, Asymmetric energy transfers in driven nonequilibrium systems and arrow of time, Eur. Phys. J. B 92, 190 (2019b).
  * Nazarenko [2011] S. V. Nazarenko,  _Wave Turbulence_ (Springer-Varlog, Berlin, 2011).
  * Shukla  _et al._ [2019a] V. Shukla, B. Dubrulle, S. V. Nazarenko, G. Krstulovic, and S. Thalabard, Phase transition in time-reversible Navier-Stokes equations, Phys. Rev. E 100, 043104 (2019a).
  * Bandak  _et al._ [2022] D. Bandak, N. Goldenfeld, A. A. Mailybaev, and G. Eyink, Dissipation-range fluid turbulence and thermal noise, Phys. Rev. E 105, 065113 (2022).
  * Madeira  _et al._ [2020] L. Madeira, M. A. Caracanhas, F. E. A. dos Santos, and V. S. Bagnato, Quantum turbulence in quantum gases, Annu. Rev. Condens. Matter Phys. 11, 37 (2020).
  * Fonda  _et al._ [2019] E. Fonda, K. R. Sreenivasan, and D. P. Lathrop, Reconnection scaling in quantum fluids., PNAS 116, 1924 (2019).
  * Shukla  _et al._ [2019b] V. Shukla, P. D. Mininni, G. Krstulovic, P. C. Di Leoni, and M. E. Brachet, Quantitative estimation of effective viscosity in quantum turbulence, Phys. Rev. A 99, 043605 (2019b).
  * Skrbek and Sreenivasan [2012] L. Skrbek and K. R. Sreenivasan, Developed quantum turbulence and its decay, Phys. Fluids 24, 011301 (2012).
  * Barenghi  _et al._ [2014] C. F. Barenghi, L. Skrbek, and K. R. Sreenivasan, Introduction to quantum turbulence, PNAS 111 Suppl 1, 4647 (2014).
  * Mohsen [2017] R. Mohsen,  _Classical And Quantum Dissipative Systems_ , 2nd ed. (World Scientific, 2017).
  * Weiss [1999] U. Weiss,  _Quantum Dissipative Systems_ (World Scientific, 1999).
  * D’Alessio  _et al._ [2016] L. D’Alessio, Y. Kafri, A. Polkovnikov, and M. Rigol, From quantum chaos and eigenstate thermalization to statistical mechanics and thermodynamics, Adv. Phys. , 1 (2016).
  * Caldeira and Leggett [1983] A. O. Caldeira and A. J. Leggett, Quantum tunnelling in a dissipative system, Annals of Physics 149, 374 (1983).


# SurGNN: Explainable visual scene understanding and assessment of surgical
skill using graph neural networks

Shuja Khalid
Department of Computer Science
University of Toronto
Toronto, ON M5S
<EMAIL_ADDRESS>
&Frank Rudzicz
Faculty of Computer Science
Dalhousie University
Halifax, NS B3H 4R2
<EMAIL_ADDRESS>
Corresponding author

###### Abstract

This paper explores how graph neural networks (GNNs) can be used to enhance
visual scene understanding and surgical skill assessment. By using GNNs to
analyze the complex visual data of surgical procedures represented as graph
structures, relevant features can be extracted and surgical skill can be
predicted. Additionally, GNNs provide interpretable results, revealing the
specific actions, instruments, or anatomical structures that contribute to the
predicted skill metrics. This can be highly beneficial for surgical educators
and trainees, as it provides valuable insights into the factors that
contribute to successful surgical performance and outcomes. SurGNN proposes
two concurrent approaches – one supervised and the other self-supervised. The
paper also briefly discusses other automated surgical skill evaluation
techniques and highlights the limitations of hand-crafted features in
capturing the intricacies of surgical expertise. We use the proposed methods
to achieve state-of-the-art results on EndoVis19, and custom datasets. The
working implementation of the code can be found at
https://github.com/<redacted>.

##  1 Introduction

Explainable visual scene understanding for surgical skill assessment using
graph neural networks (GNNs) involves developing models that can interpret and
reason about the complex visual data generated during surgical procedures. By
representing the procedure as a graph structure and training a GNN on this
data, we can model the dependencies and relationships between different
surgical actions, instruments, and anatomical structures [34]. This allows us
to extract relevant features from the visual data and predict various aspects
of surgical skill, such as the overall quality of the procedure or the
proficiency of individual actions [32]. This can help surgical educators and
trainees better understand the underlying factors that contribute to
successful surgical outcomes and improve training and performance [25, 5, 28].

Traditional assessment methods often rely on subjective expert ratings or
global performance metrics [10, 26, 13] which can be biased, inconsistent, or
difficult to interpret [13, 36, 31, 14]. By contrast, GNN-based models can
provide a more objective and fine-grained analysis of surgical performance by
interpreting the rich visual data available from surgical videos [34, 32].
GNN-based models can also provide explanations for their predictions, which
can help build trust in the models and facilitate learning and improvement
[38, 15, 20]. Additionally, we use 3D representations in our analysis. To the
best of our knowledge, we are the first paper to consider explainable 3D
representations for surgical skill analysis directly from visual cues.

##  2 Related work

#### Automated surgical skill evaluation

Automated surgical skill evaluation is an active research area, with various
approaches proposed for analyzing surgical data using computer vision and
machine learning [40, 7, 21, 16, 3]. A popular approach is to use hand-crafted
features, such as motion analysis or tool usage, to train classifiers for
surgical skill assessment. For example, Gao et al [9] developed a feature-
based model that used motion analysis to classify suturing and knot-tying
skill in surgical videos. However, these hand-crafted features are often
limited by their ability to capture the complex and subtle movements that
characterize surgical expertise.

To address these limitations, researchers have also explored deep learning
techniques, such as convolutional neural networks (CNNs) and recurrent neural
networks (RNNs), for automated surgical skill evaluation. For example, Cheng
et al [4] developed a CNN-based model that used visual cues to assess the
proficiency of laparoscopic suturing and knot-tying in a simulator. Similarly,
Zia et al [41] proposed an RNN-based model that used motion analysis to
classify surgical skill in a simulation environment. However, these deep
learning models often lack interpretability, making it difficult to understand
the underlying factors that contribute to their predictions. Anastasiou et al
[1] described an action-aware Transformer with multi-head attention producing
inter-video contrastive features to regress the manually-assigned skill
scores. This approach, although novel, requires consensus on what is
considered good performance and has only been tested in the JIGSAWS [9]
dataset, which consists of very short, simulated procedures.

To address this interpretability issue, recent research focuses on explainable
AI, including attention mechanisms and GNNs, for automated surgical skill
evaluation. For example, Hira et al [12] developed an attention-based model
that used a hierarchical CNN to identify relevant visual features for surgical
skill assessment. Meanwhile, Ban et al [2] proposed a GNN-based model that
represented surgical actions as nodes in a graph structure, allowing for a
more fine-grained analysis of surgical skill. These approaches represent
promising directions towards automated surgical skill evaluation systems that
are both accurate and interpretable.

#### Scene understanding using graph neural networks

Scene understanding using GNNs is an active research area, with numerous
approaches for leveraging the rich visual data available in complex scenes
[24, 6]. One popular application of GNNs in scene understanding is semantic
segmentation, where the goal is to assign a semantic label to each pixel in an
image [35]. For example, Wang et al [35] proposed a GNN-based model that used
graph convolutions to incorporate contextual information from neighboring
pixels, improving the accuracy of semantic segmentation in complex scenes.

Another popular application of GNNs in scene understanding is object detection
and recognition, where the goal is to identify and localize objects in an
image. Li et al [22] developed a GNN-based model that used spatial and
semantic information to perform object detection and semantic segmentation
simultaneously. Qi et al [30] proposed a GNN-based model that used message-
passing to refine object proposals and improve object detection performance.

GNNs have also been used for more complex tasks in scene understanding, such
as 3D object detection and reconstruction. For example, Zhang et al [39]
developed a GNN-based model that used graph convolutions to reason about the
3D geometry of objects and their relationships in a scene. Similarly, Yin et
al [37] proposed a GNN-based model that used message passing to refine 3D
object proposals and improve 3D object detection performance.

Overall, these works demonstrate the versatility and effectiveness of GNNs in
scene understanding, and highlight their potential for improving the
performance of computer vision systems in a wide range of applications,
including surgical skill evaluation.

##  3 SurGNN: The model

Our architecture consists of two self-supervised graph attention layers
followed by a global mean pooling layer and, finally, a linear classification
layer for training the model in a supervised manner with the available labels.
To generate the graph embeddings, we consider the following approaches: 2D We
follow the feature extraction approach, presented in [16] where each
instrument is segmented and tracked over the length of a clip. A set of
features are calculated from this temporal data to generate the node
embeddings. 3D Traditional approaches, discussed in section 2 of this paper
use either 2D representations, or robotic kinematic data for 3D analysis.
Since kinematic ground truth data isn’t available for laparoscopic videos, we
use neural radiance fields [8, 29, 23, 27, 17, 18] to generate dynamic scene
renderings. The neural scene renderings are advantageous as they allow for
fixed camera renderings. By fixing the position of the camera, we aim to
remove the relative motion of the camera with respect to the scene. Using the
same approach presented in [16], we generate 3D representations without the
additional noise associated with extraneous camera motion.

Figure 1: Architecture of the SurGNN model. Feature nodes are extracted in the
pre-processing stage and the adjacency matrix generated. The resulting model
is trained with a supervised objective using the provided labels. The self-
supervised approach uses a spectral objective.

###  3.1 Self-supervised Architecture

The self-supervised architecture is inspired by Tsitsulin et al. [33], who
used a spectral loss, which is commonly used in GNNs to optimize the spectral
properties of the learned representation. The spectral loss is defined as the
Frobenius norm of the difference between the spectral embeddings of the
original graph and the reconstructed graph, and is typically used in
conjunction with other loss functions, such as cross-entropy loss. This is
represented mathematically using the following equations:

The normalized Laplacian matrix L is calculated as follows:

$L=I-D^{-1/2}AD^{-1/2}$ (1)

We compute the eigenvectors $\\{v_{1},v_{2},\ldots,v_{n}\\}$ and eigenvalues
$\\{\lambda_{1},\lambda_{2},\ldots,\lambda_{n}\\}$ of the Laplacian matrix
$L$.

$Lv_{i}=\lambda_{i}v_{i},\quad i=1,2,\ldots,n$ (2)

For the self-supervised task, we randomly mask some nodes or edges in the
graph. The goal is to predict the missing parts of the graph from its spectral
representation:

$\text{Original Laplacian Matrix:}\quad L_{\text{original}}=L$ (3)
$\text{Reconstructed Laplacian Matrix:}\quad
L_{\text{reconstructed}}=\text{GNN\\_Decoder}(v_{1},v_{2},\ldots,v_{n})$ (4)

The spectral loss can be defined using the mean squared error (MSE) between
the original and reconstructed Laplacian matrices.

$\text{Spectral
Loss}=\frac{1}{n}\sum_{i=1}^{n}\|L_{\text{original}}^{(i)}-L_{\text{reconstructed}}^{(i)}\|_{F}^{2}$
(5)

The first part of the model architecture in Figure 1 serves as the encoder
which takes the graph representation of the input and creates an embedding of
length 32. Masking out each embedding ans using this loss function for
training encourages the embeddings to capture the spectral properties of the
graph, such as smoothness and locality, while also preserving the modularity
structure of the network. By using this loss, we are able to systematically
train the model without labels.

Figure 2: Our content extraction pipeline for defining the nodes of the
proposed GNN architecture.

##  4 Experiments

###  4.1 Datasets

#### EndoVis19

The EndoVis19 dataset was created for endoscopic video analysis in surgical
scenarios. The dataset consists of over 60 hours of endoscopic video
recordings of 10 different surgical procedures, including cholecystectomy,
appendectomy, and hernia repair. The videos were captured using various
endoscopic cameras, and include both raw and annotated videos, as well as
surgical tool and instrument segmentation masks. It is notable for its
relatively large size and diverse set of surgical procedures, which enables
the development and evaluation of algorithms for a wide range of endoscopic
surgical scenarios. The dataset also includes expert annotations of surgical
phases and actions, which can be used for developing and evaluating algorithms
for automated surgical phase detection and action recognition.

#### Custom

Our original dataset contains 309 sample clips, evenly distributed among the
three major categories of novice, intermediate, and expert. This dataset was
procured as part of data-sharing agreements that [ANON] has with its partner
institutions and is limited to laparoscopic cholecystectomy procedures. The
videos are double-rated and annotated by skilled clinical analysts.

###  4.2 Pre-processing

We direct the reader to Figure 4 and describe the steps to create a surgical
graph from a surgical video:

#### Video segmentation

: Each surgical video needs to be segmented into different phases or steps,
such as the initial incision, tissue dissection, and wound closure. This can
be done manually or using automated techniques such as motion detection or
deep learning algorithms. For the purposes of this paper, this was done
manually and each procedure was broken down into two major components – the
calot and the dissection phases.

#### Graph construction

: Once the different surgical phases have been identified, they are
represented as nodes in a graph, and the relationships between them are
captured as edges. For example, an edge created between calot and dissection
phases indicates that the tissue has been cut and is about to be dissected.

#### Feature extraction

: Various features can be extracted from the surgical graph, such as the
frequency and duration of each action, the number of edges between different
actions, and the overall structure of the graph. We extract this information
directly by extracting instrument motion statistics across frames such as the
position and speed of the instrument over time, illustrated in figure 2. For
the purposes of this paper, we compare and contrast 2 approaches, namely
segmentation assisted 2D [16] and radiance fields assisted 3D feature
extraction [8, 29, 23, 27, 17, 18] as described in Figure 3. For the 3D case,
video clips were limited to 1 minute segments. This is because of the
extremely computationally intensive nature of the reconstructions. Please
refer to the Appendix for further details.

(a) GT

(b) recon.

(c) depth

(d) mask

(e) recon.

(f) depth

(g) mask

Figure 3: 3D representations: We include the reconstructed images at times t=0
left and t=1 right for 3 cases, along with the corresponding instrument
segmentation masks and depth masks. Note: Since this is a fixed camera view
rendering, the ground truth image corresponds to the image at t=0.

#### Data balancing

: The extracted features and their associated features are balanced by
creating synthetic data using the ADASYN [11] algorithm with seven
$k$-neighbours.

#### Graph-based analysis

: The features extracted from the surgical graph can be used to assess various
aspects of surgical skill, such as the speed and accuracy of the surgeon’s
movements, the level of coordination between different actions, and the
complexity of the procedure.

###  4.3 Training

#### Supervised learning

The architecture is shown in Figure 1. We use the Adam optimizer [19] with an
initial learning rate of 0.0025 which decays at 200 and 400 epochs with a
batch size of 32. We train each model for a total of 1000 epochs on one Nvidia
RTX 3090 card.

#### Self-supervised learning

One approach to training graph neural networks (GNNs) is to use spectral self-
supervised learning. This technique involves using the graph Laplacian matrix
to generate embeddings for each node in the graph, which can then be used to
train the GNN.

To implement spectral self-supervised learning for GNNs, we can use the
following steps:

  * •

Construct the adjacency matrix $A$ for the graph of interest.

  * •

Compute the graph Laplacian matrix $L$ from the adjacency matrix using the
formula $L=D-A$, where $D$ is the diagonal matrix of node degrees.

  * •

Compute the eigenvectors and eigenvalues of the Laplacian matrix $L$.

  * •

Select the $k$ smallest eigenvalues and corresponding eigenvectors to form the
spectral embedding matrix $H$.

  * •

Train the GNN using the spectral embedding matrix $H$ as illustrated in
equation 5.

where $\mathbf{X}_{i}$ and $\mathbf{X}j$ are the node features of nodes $i$
and $j$ in the graph, $w_{ij}$ is the weight of the edge connecting nodes $i$
and $j$, and $\mathbf{Z}_{i}$ and $\mathbf{Z}_{j}$ are the corresponding
embeddings generated by the GNN. The spectral loss function aims to minimize
the distance between the node embeddings in the learned feature space while
taking into account the graph structure.

An important advantage of spectral self-supervised learning is that it can be
used to learn embeddings that capture the global structure of the graph,
rather than just local node features. This can be particularly useful in
domains where the relationships between nodes are complex and non-local, as in
the case of surgical videos.

##  5 Explainability

Explainability is the ability of a model to provide understandable and
transparent reasoning for its predictions. In the case of GNNs, explainability
is particularly important because the structure of the input data is often
essential to the prediction. Therefore, understanding how the model arrives at
its prediction can provide valuable insights into the underlying data and the
decision-making process. We visualize the learned node embeddings or feature
representations, which are the intermediate outputs of SurGNN. By visualizing
these embeddings, we gain insights into the relationships between nodes in the
graph and the importance of different features for the prediction. We
visualize the structure of the embeddings and project them in 2D space in
Figure 4. The assigned classifications are thus accompanied by supporting
visualizations which can be invaluable for surgeons interested in improving
their skills by comparing their performance to that of experienced surgeons.

##  6 Results and discussion

We present our findings on two datasets – the custom dataset and the EndoVis19
dataset – in Tables LABEL:table:sst and LABEL:table:endovis respectively. In
both, each method is evaluated for different categories, and we report
performance metrics such as Pearson correlation coefficients, Spearman rank
correlation coefficients, Kendall’s tau, precision, recall, and F1 scores. For
each dataset, the baseline results are calculated by sampling from a Gaussian
distribution and we average over 10 unique runs. To the best of our knowledge,
there aren’t any published baselines in literature that we can compare our
approach against. We thus present the results of our supervised and self-
supervised approaches.

#### [ANON] dataset

The baseline method performs poorly in terms of Spearman and Kendall
correlation coefficients, with a mean value of $0.045$ for the overall score
category. We include another approach, "SurGNN (SS)", which is a self-
supervised approach. In contrast, our proposed method, represented by "SurGNN
(S)" in the table, shows a significant improvement across all performance
metrics and outperforms the baseline method in terms of precision, recall, and
F1.

#### EndoVis19 dataset

The results show that our proposed method, represented by "SurGNN (S)," again
outperforms the baseline method across all performance metrics. For instance,
for the Overall category, our proposed method achieves Pearson, Spearman, and
Kendall’s tau coefficients of 0.709, 0.711, and 0.678, respectively, while the
baseline method has Spearman correlation of only $0.088$. Similarly, for the
Time and Motion category, our proposed method outperforms the baseline method
across all metrics.

Overall, our results suggest that our proposed method, "SurGNN (S)," is
effective in improving performance on both datasets for most categories and
performance metrics. However, the method’s performance may vary depending on
the category and the dataset. For example, the 2D analysis shows signifcantly
better results than the 3D analysis. This could be a results of a few factors:
The 3D analysis is limited to a length of 30s or approximately 60 frames in
total (sampling 2 frames/s). The 30s segment chosen is also a potential source
of error as it was chosen using clinical judgement. Longer samples were shown
to lead to image degeneration, decreasing the quantity of available samples to
211 and 50 samples for the ANON and EndoVis19 datasets respectively.

Table 1: Model performance on [ANON] dataset. SurGNN (S) refers to our results
using the supervised approach whereas SurGNN (SS) refers to the proposed Self-
supervised clustering approach detailed in the previous section.

Method | Category | Mode | N | Pearson | Spearman | Kendall | Precision | Recall | F1
---|---|---|---|---|---|---|---|---|---
- | Overall (baseline) | - | - | - | 0.045 | - | 0.140 | 0.140 | 0.140
SurGNN (S) | Overall | 2D | 309 | 0.724 | 0.717 | 0.671 | 0.690 | 0.690 | 0.680
SurGNN (S) | Precision of Operating technique | 2D | 309 | 0.731 | 0.724 | 0.675 | 0.670 | 0.680 | 0.675
SurGNN (S) | Economy of movements | 2D | 309 | 0.682 | 0.669 | 0.642 | 0.681 | 0.659 | 0.670
SurGNN (SS) | Overall | 2D | 309 | 0.424 | 0.419 | 0.411 | 0.530 | 0.580 | 0.550
SurGNN (SS) | Precision of Operating technique | 2D | 309 | 0.444 | 0.413 | 0.397 | 0.580 | 0.520 | 0.550
SurGNN (SS) | Economy of movements | 2D | 309 | 0.415 | 0.412 | 0.406 | 0.560 | 0.540 | 0.550
SurGNN (S) | Overall | 3D | 211 | 0.355 | 0.325 | 0.316 | 0.430 | 0.440 | 0.435
SurGNN (S) | Precision of Operating technique | 3D | 211 | 0.328 | 0.315 | 0.321 | 0.430 | 0.440 | 0.435
SurGNN (S) | Economy of movements | 3D | 211 | 0.409 | 0.401 | 0.389 | 0.410 | 0.410 | 0.410
SurGNN (SS) | Overall | 3D | 211 | 0.288 | 0.273 | 0.279 | 0.320 | 0.320 | 0.320
SurGNN (SS) | Precision of Operating technique | 3D | 211 | 0.261 | 0.254 | 0.255 | 0.300 | 0.310 | 0.300
SurGNN (SS) | Economy of movements | 3D | 211 | 0.214 | 0.197 | 0.196 | 0.310 | 0.310 | 0.310

Table 2: Model performance on EndoVis19 dataset - SurGNN (S) refers to our
results using the supervised approach whereas SurGNN (SS) refers to the
proposed self-supervised clustering approach detailed in the previous section.

Method | Category | Mode | N | Pearson | Spearman | Kendall | Precision | Recall | F1
---|---|---|---|---|---|---|---|---|---
- | Overall (baseline) | - | - | - | 0.088 | - | 0.180 | 0.180 | 0.180
SurGNN (S) | Overall | 2D | 69 | 0.709 | 0.711 | 0.678 | 0.700 | 0.700 | 0.690
SurGNN (S) | Time and Motion | 2D | 69 | 0.712 | 0.705 | 0.686 | 0.800 | 0.800 | 0.790
SurGNN (S) | Instrument Handling | 2D | 69 | 0.663 | 0.677 | 0.628 | 0.660 | 0.670 | 0.650
SurGNN (SS) | Overall | 2D | 69 | 0.409 | 0.394 | 0.381 | 0.500 | 0.520 | 0.510
SurGNN (SS) | Time and Motion | 2D | 69 | 0.387 | 0.368 | 0.365 | 0.470 | 0.550 | 0.510
SurGNN (SS) | Instrument Handling | 2D | 69 | 0.430 | 0.381 | 0.675 | 0.510 | 0.520 | 0.520
SurGNN (S) | Overall | 3D | 50 | 0.311 | 0.297 | 0.281 | 0.360 | 0.350 | 0.360
SurGNN (S) | Time and Motion | 3D | 50 | 0.324 | 0.312 | 0.309 | 0.360 | 0.360 | 0.360
SurGNN (S) | Instrument Handling | 3D | 50 | 0.317 | 0.316 | 0.297 | 0.320 | 0.320 | 0.320
SurGNN (SS) | Overall | 3D | 50 | 0.156 | 0.144 | 0.121 | 0.250 | 0.250 | 0.250
SurGNN (SS) | Time and Motion | 3D | 50 | 0.169 | 0.154 | 0.142 | 0.275 | 0.275 | 0.275
SurGNN (SS) | Instrument Handling | 3D | 50 | 0.165 | 0.155 | 0.140 | 0.270 | 0.270 | 0.270

Figure 4: Explainable clusters in of surgical procedures allows for surgeries
to be assessed and ranked.

###  6.1 Limitations and future directions

We aim to provide a structured basis for research in visual scene
understanding for surgical applications. Below are some items that should be
considered before widespread application of these systems.

#### Inter-observer variability

: Different expert surgeons may have varying opinions on what constitutes good
or bad surgical skill, leading to inconsistent ratings across different
evaluators. Each dataset used in this analysis has been rated by multiple
analysts and has thus gone through a rigorous quality assurance process.

#### Lack of granularity

: Subjective labels may provide only a high-level assessment of surgical
performance, such as overall performance or proficiency in a particular task,
without providing detailed feedback on specific actions or movements.

#### Bias

: Human evaluators may unconsciously or consciously introduce bias into their
ratings, for example, by favoring or discriminating against certain
individuals or groups based on factors such as gender, race, or ethnicity.

#### Limited feedback

: Subjective labels may not provide specific feedback on how to improve
performance, which can be crucial for surgical trainees to learn and improve.

#### Time-consuming

: Collecting subjective ratings from expert evaluators can be time-consuming
and expensive, especially when evaluating large datasets.

#### Concurrent validity

: The datasets presented in this paper consist primarily of cholecystectomy
procedures. Testing the efficacy of our proposed method on datasets with
different types of demographics is an important line of future work.

To address these issues, objective and quantitative measures of surgical
performance, such as motion analysis or tool usage, have been proposed. These
measures can provide more granular and unbiased assessments of surgical skill,
and can be automated using computer vision and machine learning techniques.
Additionally, such objective measures can provide specific feedback on how to
improve performance and can facilitate more efficient evaluation of large
datasets.

##  7 Conclusion

We present two distinct approaches using graph neural networks (GNNs) for
explainable visual scene understanding in surgical skill assessment. We
demonstrate that, by representing surgical procedures as graph structures and
training them using supervised and self-supervised techniques, relevant
features can be extracted from the complex visual data and used to predict
various aspects of surgical skill. GNNs also allow for interpretable results
by identifying the specific actions, instruments, or anatomical structures
that contribute to the predicted skill metrics. We highlight the limitations
of traditional assessment methods and how GNN-based models can provide a more
objective and fine-grained analysis of surgical performance. Explainable
visual scene understanding models for surgical skill assessment represents an
exciting opportunity to leverage cutting-edge AI techniques to improve
surgical education, training, and outcomes.

## References

  * [1] Dimitrios Anastasiou, Yueming Jin, Danail Stoyanov, and Evangelos Mazomenos.  Keep your eye on the best: Contrastive regression transformer for skill assessment in robotic surgery.  IEEE Robotics and Automation Letters, 8(3):1755–1762, 2023.
  * [2] Yutong Ban, Jennifer A Eckhoff, Thomas M Ward, Daniel A Hashimoto, Ozanan R Meireles, Daniela Rus, and Guy Rosman.  Concept graph neural networks for surgical video understanding.  arXiv preprint arXiv:2202.13402, 2022.
  * [3] François Chadebecq, Laurence B Lovat, and Danail Stoyanov.  Artificial intelligence and automation in endoscopy and surgery.  Nature Reviews Gastroenterology & Hepatology, 20(3):171–182, 2023\.
  * [4] Ke Cheng, Jiaying You, Shangdi Wu, Zixin Chen, Zijian Zhou, Jingye Guan, Bing Peng, and Xin Wang.  Artificial intelligence-based automated laparoscopic cholecystectomy surgical phase recognition and analysis.  Surgical endoscopy, 36(5):3160–3168, 2022.
  * [5] Charity H Evans and Kimberly D Schenarts.  Evolving educational techniques in surgical training.  Surgical Clinics, 96(1):71–88, 2016.
  * [6] Kanglong Fan, Wei Liu, Xiaowen Chen, Bharath Ramesh, and Cheng Xiang.  An interpretable scene understanding framework via graph learning.  Available at SSRN 4238333, 2018.
  * [7] Mahtab J Fard, Sattar Ameri, R Darin Ellis, Ratna B Chinnam, Abhilash K Pandya, and Michael D Klein.  Automated robot-assisted surgical skill evaluation: Predictive analytics approach.  The International Journal of Medical Robotics and Computer Assisted Surgery, 14(1):e1850, 2018.
  * [8] Chen Gao, Ayush Saraf, Johannes Kopf, and Jia-Bin Huang.  Dynamic view synthesis from dynamic monocular video.  In Proceedings of the IEEE/CVF International Conference on Computer Vision, pages 5712–5721, 2021.
  * [9] Yixin Gao, S Swaroop Vedula, Carol E Reiley, Narges Ahmidi, Balakrishnan Varadarajan, Henry C Lin, Lingling Tao, Luca Zappella, Benjamın Béjar, David D Yuh, et al.  Jhu-isi gesture and skill assessment working set (jigsaws): A surgical activity dataset for human motion modeling.  In MICCAI workshop: M2cai, volume 3, 2014.
  * [10] Jean D Gray.  Global rating scales in residency education.  Academic Medicine, 71(1):S55–63, 1996.
  * [11] Haibo He, Yang Bai, Edwardo A Garcia, and Shutao Li.  Adasyn: Adaptive synthetic sampling approach for imbalanced learning.  In 2008 IEEE international joint conference on neural networks (IEEE world congress on computational intelligence), pages 1322–1328. IEEE, 2008\.
  * [12] Sanchit Hira, Digvijay Singh, Tae Soo Kim, Shobhit Gupta, Gregory Hager, Shameema Sikder, and S Swaroop Vedula.  Video-based assessment of intraoperative surgical skill.  International Journal of Computer Assisted Radiology and Surgery, 17(10):1801–1811, 2022.
  * [13] Nancy J Hogle, Yi Liu, R Todd Ogden, and Dennis L Fowler.  Evaluation of surgical fellows’ laparoscopic performance using global operative assessment of laparoscopic skills (goals).  Surgical endoscopy, 28:1284–1290, 2014.
  * [14] Shiva Jayaraman et al.  The objective structured assessment of technical skills (osats) and its role in surgical education.  Annals of the Royal College of Surgeons of England, 91(2):100–105, 2009.
  * [15] Shuiwang Ji et al.  Graph convolutional networks with explainable node importance scores.  In Advances in Neural Information Processing Systems (NeurIPS), 2019\.
  * [16] Shuja Khalid, Vanessa Palter, Teodor Grantcharov, and Frank Rudzicz.  Or vision: Objective, explainable assessment of surgical skill with deep learning.  https://doi.org/10.21203/rs.3.rs-1978829/v1, 2022.
  * [17] Shuja Khalid and Frank Rudzicz.  wildnerf: Complete view synthesis of in-the-wild dynamic scenes captured using sparse monocular data.  arXiv preprint arXiv:2209.10399, 2022.
  * [18] Shuja Khalid and Frank Rudzicz.  Refinerf: Modelling dynamic neural radiance fields with inconsistent or missing camera parameters.  arXiv preprint arXiv:2303.08695, 2023.
  * [19] Diederik P Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [20] Thomas N. Kipf and Max Welling.  Semi-supervised classification with graph convolutional networks.  In International Conference on Learning Representations (ICLR), 2017\.
  * [21] Marc Levin, Tyler McKechnie, Shuja Khalid, Teodor P Grantcharov, and Mitchell Goldenberg.  Automated methods of technical skill assessment in surgery: a systematic review.  Journal of surgical education, 76(6):1629–1639, 2019.
  * [22] Fuxian Li, Jie Feng, Huan Yan, Depeng Jin, and Yong Li.  Crowd flow prediction for irregular regions with semantic graph attention network.  ACM Transactions on Intelligent Systems and Technology (TIST), 13(5):1–14, 2022.
  * [23] Zhengqi Li, Simon Niklaus, Noah Snavely, and Oliver Wang.  Neural scene flow fields for space-time view synthesis of dynamic scenes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6498–6508, 2021.
  * [24] Zhijun Liang, Junfa Liu, Yisheng Guan, and Juan Rojas.  Visual-semantic graph attention networks for human-object interaction detection.  In 2021 IEEE international conference on robotics and biomimetics (ROBIO), pages 1441–1447. IEEE, 2021.
  * [25] Heidi Maertens, A Madani, T Landry, Frank Vermassen, Isabelle Van Herzeele, and R Aggarwal.  Systematic review of e-learning for surgical training.  Journal of British Surgery, 103(11):1428–1437, 2016.
  * [26] JA Martin, Glenn Regehr, Richard Reznick, Helen Macrae, John Murnaghan, Carol Hutchison, and M Brown.  Objective structured assessment of technical skill (osats) for surgical residents.  British journal of surgery, 84(2):273–278, 1997.
  * [27] Keunhong Park, Utkarsh Sinha, Peter Hedman, Jonathan T Barron, Sofien Bouaziz, Dan B Goldman, Ricardo Martin-Brualla, and Steven M Seitz.  Hypernerf: A higher-dimensional representation for topologically varying neural radiance fields.  arXiv preprint arXiv:2106.13228, 2021.
  * [28] Brian Andrew Parsons, Natalie S Blencowe, Andrew D Hollowood, and Janet R Grant.  Surgical training: the impact of changes in curriculum and experience.  Journal of surgical education, 68(1):44–51, 2011.
  * [29] Albert Pumarola, Enric Corona, Gerard Pons-Moll, and Francesc Moreno-Noguer.  D-nerf: Neural radiance fields for dynamic scenes.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 10318–10327, 2021.
  * [30] Charles R Qi, Yin Zhou, Mahyar Najibi, Pei Sun, Khoa Vo, Boyang Deng, and Dragomir Anguelov.  Offboard 3d object detection from point cloud sequences.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 6134–6144, 2021.
  * [31] Osama Samargandi et al.  The objective structured assessment of technical skills (osats) for surgical trainees: A systematic review and meta-analysis.  Journal of Surgical Education, 74(2):388–399, 2017.
  * [32] Lalithkumar Seenivasan, Mobarakol Islam, Chi-Fai Ng, Chwee Ming Lim, and Hongliang Ren.  Biomimetic incremental domain generalization with a graph network for surgical scene understanding.  Biomimetics, 7(2):68, 2022.
  * [33] Anton Tsitsulin, John Palowitch, Bryan Perozzi, and Emmanuel Müller.  Graph clustering with graph neural networks.  arXiv preprint arXiv:2006.16904, 2020.
  * [34] Natalia Valderrama, Paola Ruiz Puentes, Isabela Hernández, Nicolás Ayobi, Mathilde Verlyck, Jessica Santander, Juan Caicedo, Nicolás Fernández, and Pablo Arbeláez.  Towards holistic surgical scene understanding.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2022: 25th International Conference, Singapore, September 18–22, 2022, Proceedings, Part VII, pages 442–452. Springer, 2022\.
  * [35] Hongzhen Wang, Ying Wang, Qian Zhang, Shiming Xiang, and Chunhong Pan.  Gated convolutional neural network for semantic segmentation in high-resolution images.  Remote Sensing, 9(5):446, 2017.
  * [36] Christy Yan et al.  Evaluation of the objective structured assessment of technical skills (osats) tool for assessment of technical skills in laparoscopic surgery.  Journal of Surgical Education, 71(4):495–501, 2014.
  * [37] Junbo Yin, Jianbing Shen, Chenye Guan, Dingfu Zhou, and Ruigang Yang.  Lidar-based online 3d video object detection with graph-based message passing and spatiotemporal transformer attention.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 11495–11504, 2020.
  * [38] Rex Ying et al.  Understanding graph neural networks with graph logits.  In Proceedings of the 37th International Conference on Machine Learning (ICML), 2020.
  * [39] Yanan Zhang, Di Huang, and Yunhong Wang.  Pc-rgnn: Point cloud completion and graph neural network for 3d object detection.  In Proceedings of the AAAI Conference on Artificial Intelligence, volume 35, pages 3430–3437, 2021.
  * [40] Aneeq Zia and Irfan Essa.  Automated surgical skill assessment in rmis training.  International journal of computer assisted radiology and surgery, 13:731–739, 2018.
  * [41] Aneeq Zia, Andrew Hung, Irfan Essa, and Anthony Jarc.  Surgical activity recognition in robot-assisted radical prostatectomy using deep learning.  In Medical Image Computing and Computer Assisted Intervention–MICCAI 2018: 21st International Conference, Granada, Spain, September 16-20, 2018, Proceedings, Part IV 11, pages 273–280. Springer, 2018\.