1. Introduction
Two-dimensional turbulence is the study of incompressible hydrodynamics at large (including infinite) Reynolds numbers. It is a vibrant field of both mathematics and physics that began with Euler (Reference Euler1757), who derived the basic equations of motion. Turbulent flows in two space dimensions do not exist as classical fluids in nature. Rather, they constitute basic models of intermediate-scale flows in ‘almost’ two-dimensional (thin) domains (Majda & Bertozzi Reference Majda and Bertozzi2002; Cullen Reference Cullen2006; Majda & Wang Reference Majda and Wang2006; Dolzhansky Reference Dolzhansky2013; Pedlosky Reference Pedlosky2013; Zeitlin Reference Zeitlin2018).
The conditions of two-dimensional (2-D) turbulence can be emulated in experiments. One set-up is a soap film flowing rapidly through a fine comb (Couder Reference Couder1984). Another is a conducting fluid confined to a thin layer and driven into turbulence by a temporally varying magnetic field (Sommeria Reference Sommeria1988). When such a ‘quasi-2-D’ flow is released it self-organises into blob-like condensates. The progression is depicted in figure 1 for a spherical domain. Heuristically, the mechanism is driven by the merging of equally signed vorticity regions. This large-scale fusion is balanced by fine-scale dissipation. In many ways, 2-D turbulence is propelled by the quest to understand the resulting scale separation.
Figure 1. Evolution of vorticity for Euler's equations on the sphere. Vorticity regions of equal sign undergo merging to form stable, interacting vortex condensates.
To make theoretical progress, Onsager (Reference Onsager1949) applied statistical mechanics to a large but finite number 
$N$ of point vortices. They are weak (distributional) solutions of Euler's equations where vorticity is a sum of weighted Dirac pulses. Onsager realised that a fixed number of positive and negative point vortices, confined to a bounded domain, can have energies from 
$-\infty$ to 
$\infty$. The phase volume function, therefore, has an inflexion point at some finite energy. At this energy, the thermodynamical temperature is zero. If the energy exceeds this point, so the temperature is negative, then equally signed vortices should cluster according to thermodynamics. This theory of statistical hydrodynamics is a prominent, although lesser-known, part of Onsager's legacy (Eyink & Sreenivasan Reference Eyink and Sreenivasan2006). Its validation is a long-standing open problem (Marchioro & Pulvirenti Reference Marchioro and Pulvirenti2012, discussion in chap. 7). On the mathematical side, Caglioti et al. (Reference Caglioti, Lions, Marchioro and Pulvirenti1992, Reference Caglioti, Lions, Marchioro and Pulvirenti1995) and Kiessling (Reference Kiessling1993) gave rigorous results on clustering of point vortices in the negative temperature regime as 
$N\to \infty$ under the assumption of ergodicity. Their work has fostered much theoretical progress (Marchioro & Pulvirenti Reference Marchioro and Pulvirenti2012, and references therein). On the experimental side, seventy years after Onsager presented his theory, the conditions of negative temperature point vortex dynamics were experimentally realised in planar Bose–Einstein condensates (Gauthier et al. Reference Gauthier, Reeves, Yu, Bradley, Baker, Bell, Rubinsztein- Dunlop, Davis and Neely2019; Johnstone et al. Reference Johnstone, Groszek, Starkey, Billington, Simula and Helmerson2019). As predicted, persisting vortex clusters emerge.
Onsager's theory cannot be applied to continuous vorticity fields (corresponding to smooth solutions). Consequently, it is natural to look for a statistical theory of continua. One approach is to expand the vorticity field in a Fourier series and then truncate it (Kraichnan Reference Kraichnan1975). The truncated, finite-dimensional system preserves phase volume and quadratic invariants, but not higher-order invariants (Casimir functions). To account also for those invariants, another approach is to maximise the entropy of a probability distribution of coarse-grained vorticity fields (Lynden-Bell Reference Lynden-Bell1967; Miller Reference Miller1990; Robert & Sommeria Reference Robert and Sommeria1991).
All theories based on statistical mechanics assume an ergodic dynamics. Rigorous results on ergodicity are available for 2-D Navier–Stokes on a doubly periodic domain (flat torus) with added regular-in-space noise proportional to the square root of the viscosity 
$\nu$. In this setting there exists a unique stationary measure 
$\mu _\nu$ (Kuksin & Shirikyan Reference Kuksin and Shirikyan2012, Reference Kuksin and Shirikyan2017). As 
$\nu \to 0$, one obtains a stationary measure 
$\mu _0$ for the 2-D Euler equations, but it is not expected to be unique (Kuksin & Shirikyan Reference Kuksin and Shirikyan2017).
Statistical mechanics is not the only approach. Another possibility is to study energy and enstrophy spectra. The inspiration comes from Kolmogorov's (Reference Kolmogorov1941) theory of 3-D turbulence. Notably, Kraichnan (Reference Kraichnan1967) argued that viscous 2-D fluids in forced equilibrium, where energy at an intermediate scale is fed into the system, exhibit a forward cascade of enstrophy into fine scales and a simultaneous backward cascade of energy into large scales. Direct numerical simulations typically support Kraichnan's theory (Xiao et al. Reference Xiao, Wan, Chen and Eyink2009, and references therein).
For 2-D systems with no energy dissipation at the large scale (so that vortex condensation occurs), numerical simulations develop a ‘broken line’ energy spectrum with a steep slope at the large scale, typically steeper than 
$k^{-3}$ where 
$k$ is the wavenumber, and then a swift switch at an intermediate scale to a less steep slope, typically between 
$k^{-5/3}$ and 
$k^{-1}$ (Boffetta & Ecke Reference Boffetta and Ecke2012, and references therein). An approximate broken line energy spectrum is also observed in zonal and meridional wind measurements on Earth over the scales 3–10 000 km (Nastrom, Gage & Jasperson Reference Nastrom, Gage and Jasperson1984).
To better understand characteristic energy spectra it is natural to impose a splitting of the vorticity field 
$\omega = \omega _s + \omega _r$ into a large-scale component 
$\omega _s$ and a small-scale component 
$\omega _r$. A wavelet-based vorticity splitting is proposed by Farge, Schneider & Kevlahan (Reference Farge, Schneider and Kevlahan1999) and applied to numerical simulation on the doubly periodic square (where condensation occurs) by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007). The results (Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007, figure 1f) show an energy spectrum slope of 
$k^{-3}$ for the large-scale component and of 
$k^{-1}$ for the intermediate-to-small-scale component. It is a powerful technique to analyse energy spectra, but it depends on a choice of wavelet basis and parameters identifying the different scales. Therefore, the method cannot give insights into the mechanisms behind vortex condensation or broken line energy spectra.
In this paper, we give a new, canonical decomposition of vorticity. By ‘canonical’ we mean parameter free, determined solely in terms of the data for the 2-D Euler equations: the Poisson bracket and the Hamiltonian function. The decomposition has the following properties:
(i) The vorticity field
$\omega = \omega _s + \omega _r$ is a steady state if and only if $\omega _r = 0$.(ii) Under numerical simulation of Euler's equations,
$\omega _s$ and $\omega _r$ evolve into a separation of scales. The component $\omega _s$ traps large-scale condensates whereas the component $\omega _r$ captures small-scale fluctuations.(iii) After a short transient time, the energy spectrum slope of
$\omega _s$ is approximately $k^{-3}$ and that of $\omega _r$ is between $k^{-5/3}$ and $k^{-1}$.(iv) Over time, the component
$\omega _r$ displays an average enstrophy increase (quantifying the forward enstrophy cascade) and an average energy decrease (quantifying the backward energy cascade).
The coupled equations governing the dynamics of 
$\omega _s$ and 
$\omega _r$ embody a new line of attack for 2-D turbulence. Our standpoint is that a detailed study of these equations may explain the mechanisms behind vortex condensation and broken line energy spectra, or at least give deep insights. The numerical simulations we present suggest so.
1.1. Two-dimensional Euler equations
Our starting point is Euler's equations for an inviscid, incompressible fluid on a 2-D closed surface. Throughout the paper, we take the surface to be the unit sphere 
$\mathbb {S}^2\subset \mathbb {R}^3$. It makes our arguments more explicit and enables numerics. Most concepts are readily transferable to arbitrary closed surfaces (in particular to the flat torus, which is the most studied example in the literature albeit less relevant than the sphere in applications).
In vorticity formulation, Euler's equations on 
$\mathbb {S}^2$ are
\begin{equation} \dot{\omega} =\lbrace\psi,\omega\rbrace ,\quad \varDelta\psi = \omega, \end{equation}
where 
$\{{\cdot },{\cdot }\}$ is the Poisson bracket on 
$\mathbb {S}^2$, 
$\omega$ is the vorticity function of the fluid, related to the fluid velocity 
${\boldsymbol{v}}$ via 
$\omega = \operatorname {curl}\boldsymbol {v}$, and 
$\psi$ is the streamfunction, related to the vorticity function via the Laplace–Beltrami operator 
$\varDelta$. Geometrically, (1.1a,b) constitute an infinite-dimensional Lie–Poisson system (Arnold & Khesin Reference Arnold and Khesin1998, and references therein). The phase space consists of vorticity fields and is equipped with the following infinite-dimensional Poisson bracket:
\begin{equation} \prec \!F, G \!\succ\! (\omega) = \int_{\mathbb{S}^2} \left\{ \frac{\delta F}{\delta \omega}, \frac{\delta G}{\delta \omega} \right\} \omega . \end{equation}The Hamiltonian (total energy) for the vorticity equation (1.1a,b) is
\begin{equation} H ={-}\tfrac{1}{2}\int_{\mathbb{S}^2}\psi\omega . \end{equation}In addition to total energy, there is an infinite number of conservation laws: total angular momentum
\begin{equation} \boldsymbol{L} = \int_{\mathbb{S}^2} \omega\boldsymbol{n} , \quad \boldsymbol{n}\ \text{unit normal on}\ \mathbb{S}^2, \end{equation}and Casimir functions
\begin{equation} C_f(\omega) = \int_{\mathbb{S}^2}f(\omega),\quad \text{for any smooth}\ f\colon \mathbb{R}\to \mathbb{R}. \end{equation}These conservation laws are fundamental for long-time behaviour. In particular, the presence of infinitely many Casimir functions sets apart 2-D from 3-D fluids.
1.2. Overview of the paper
Zeitlin's truncated model of Euler's equations originates from the vorticity formulation (1.1a,b). In the spirit of quantisation theory, the space of vorticity functions is replaced by the space 
$\mathfrak {u}(N)$ of skew-Hermitian complex matrices, while the Poisson bracket is replaced by the matrix commutator. The size 
$N$ of the matrices is the spatial discretisation parameter, related to ‘Planck's constant’ in quantum theory via 
$\hbar = 1/N$. We present an overview of how functions and matrices are related in § 2.
There are at least two advantages of Zeitlin's model. First, it yields a spatial discretisation that preserves all the underlying geometry of the Euler equations: the Hamiltonian structure and the conservation laws. Combined with a symplectic time-integration scheme one can obtain fully structure-preserving numerical methods (Modin & Viviani Reference Modin and Viviani2020a). Second, complicated topological or geometrical properties of the Euler equations can be described in terms of standard tools from linear algebra and matrix Lie groups. Indeed, our splitting naturally arises from the standard spectral decomposition of Hermitian matrices applied to the stream matrix. The splitting has a precise geometric meaning in terms of Lie algebras, but also a dynamical interpretation as the steady and unsteady vorticity components. The details are given in § 3.
Numerical experiments for the canonical splitting are presented and discussed in § 4. We demonstrate that the components in the splitting converge into a separation of scales. They also capture broken line energy spectra. To rigorously answer to what degree simulations based on Zeitlin's model capture the behaviour of Euler's equations is an open problem.
In § 5 we translate our results about Zeitlin's model to the continuous Euler equations. All matrix concepts used for the canonical splitting have classical, fluid dynamical counterparts. But to define them rigorously requires prudence. It is essential to use a weak formulation. The section sets the foundation for an 
$L^\infty$-based theory of canonical vorticity splitting, independent of Zeitlin's model. It is the most mathematical part of the paper. Even so, a heuristic explanation is straightforward. Let us give it here. For a smooth vorticity function 
$\omega$, the aim is to construct a splitting 
$\omega = \omega _s + \omega _r$. Let 
$\boldsymbol{\gamma} (\tau )$ be a closed level curve of 
$\psi$. We define 
$\omega _s$, restricted to 
$\boldsymbol{\gamma}$, to take the constant value
\begin{equation} \left.\omega_s\right|_{\boldsymbol{\gamma}} = \frac{1}{\operatorname{length}(\boldsymbol{\gamma})}\int_{\boldsymbol{\gamma}} \omega\, {\rm d}\tau . \end{equation}
In other words, 
$\omega _s$ is obtained via averaging of 
$\omega$ along the level curves, or streamlines, of 
$\psi$. The mathematical difficulties arise where the level curves contain bifurcations.
2. Background to Euler–Zeitlin equations
In this section we give background to the ‘consistent truncation’ of Euler's equations introduced by Zeitlin (Reference Zeitlin1991, Reference Zeitlin2004). The model relies on quantisation of the Poisson algebra of smooth functions (Hoppe Reference Hoppe1982; Fairlie & Zachos Reference Fairlie and Zachos1989). Recall that the Poisson bracket between two smooth functions 
$f,g\in C^\infty (\mathbb {S}^2)$ is
\begin{equation} \lbrace f,g\rbrace(\boldsymbol{x}) = \boldsymbol{x} \boldsymbol{\cdot} (\boldsymbol{\nabla} f\times\boldsymbol{\nabla} g), \quad \boldsymbol{x}\in\mathbb{S}^2\subset\mathbb{R}^3. \end{equation}
Quantisation, in our context, means to find a projection from smooth functions 
$C^\infty (\mathbb {S}^2)$ to complex skew-Hermitian matrices 
$\mathfrak {u}(N)$ such that the Poisson bracket (2.1) under this map is approximated by the matrix commutator.
On the sphere, quantisation is explicit, obtained as follows. Hoppe & Yau (Reference Hoppe and Yau1998) gave an operator 
$\varDelta _N\colon \mathfrak {u}(N)\to \mathfrak {u}(N)$ with the same spectrum as the Laplace–Beltrami operator (up to truncation 
$N$). The eigenvectors 
$\boldsymbol{\mathsf{T}}^N_{lm}$ of 
$\varDelta _N$ are the quantised analogues of the spherical harmonics 
$Y_{lm}$. It is therefore natural to define the projection 
$\varPi _N\colon C^\infty (\mathbb {S}^2)\to \mathfrak {u}(N)$ mapping functions to matrices via
\begin{equation} C^\infty(\mathbb{S}^2) \ni \omega = \sum_{l=0}^\infty\sum_{m={-}l}^l\omega^{lm} Y_{lm} \mapsto \sum_{l=0}^{N}\sum_{m={-}l}^l \omega^{lm}\boldsymbol{\mathsf{T}}^N_{lm} = W \in \mathfrak{u}(N). \end{equation}For more details and explicit formulae we refer to Hoppe & Yau (Reference Hoppe and Yau1998), Zeitlin (Reference Zeitlin2004) and Modin & Viviani (Reference Modin and Viviani2020a).
The vorticity formulation (1.1a,b) uses only the Laplacian 
$\varDelta$ and the Poisson bracket 
$\{{\cdot },{\cdot }\}$. Their quantised analogues 
$\varDelta _N$ and 
$[{\cdot },{\cdot }]$ give rise to the Euler–Zeitlin equations
\begin{equation} \dot{\boldsymbol{\mathsf{W}}} =[\boldsymbol{\mathsf{P}},\boldsymbol{\mathsf{W}}],\quad \varDelta_N \boldsymbol{\mathsf{P}} = \boldsymbol{\mathsf{W}}, \end{equation}
where 
$\boldsymbol{\mathsf{W}}\in \mathfrak {su}(N)$ is the vorticity matrix and 
$\boldsymbol{\mathsf{P}}\in \mathfrak {su}(N)$ is the stream matrix. The condition
\begin{equation} \boldsymbol{\mathsf{W}} \in \mathfrak{su}(N) = \{\boldsymbol{\mathsf{A}}\in \mathfrak{u}(N)\mid \operatorname{tr} \boldsymbol{\mathsf{A}}=0 \} \end{equation}
corresponds to vanishing total circulation 
$\int _{\mathbb {S}^2}\omega = 0$.
The Euler–Zeitlin equations (2.3a,b) have been studied in various contexts, primarily on the flat torus (Zeitlin Reference Zeitlin1991; McLachlan Reference McLachlan1993; Abramov & Majda Reference Abramov and Majda2003), but also on the sphere (Zeitlin Reference Zeitlin2004; Modin & Viviani Reference Modin and Viviani2020a). Their main feature is that they preserve the rich phase space geometry of the original equations (1.1a,b), namely the Lie–Poisson structure (Modin & Viviani Reference Modin and Viviani2020a,Reference Modin and Vivianic, and references therein). In turn, this structure implies conservation of total energy 
$H(\boldsymbol{\mathsf{W}}) = \mbox{Tr}(\boldsymbol{\mathsf{PW}})/2$, (quantised) Casimirs 
$C_k(\boldsymbol{\mathsf{W}}) = \mbox {{Tr}}(\boldsymbol{\mathsf{W}}^k)$ and angular momentum 
$\boldsymbol {L}=(L_x,L_y,L_z)$. Conventional discretisations do not maintain these conservation laws.
Remark 2.1 Quantisation on the sphere is more complicated to work with than quantisation on the torus. Even so, the Euler–Zeitlin equations are more accurate on the sphere than on the torus. This has a deep geometric reason: quantisation exactly preserves the rotational symmetry of the sphere, but the translational symmetry of the torus is only approximately captured. The quantised Poisson equation on the torus, therefore, suffers from the Gibbs phenomenon, which is not present on the sphere.
In previous work (Modin & Viviani Reference Modin and Viviani2020a,Reference Modin and Vivianic) we develop a Lie–Poisson-preserving numerical method for the Euler–Zeitlin equations on the sphere and we use it to study the long-time behaviour. Our numerical results, along with earlier ones by Dritschel, Qi & Marston (Reference Dritschel, Qi and Marston2015), give strong evidence against the predictions of statistical mechanics theories, derived for the sphere by Herbert (Reference Herbert2013); see also Bouchet & Venaille (Reference Bouchet and Venaille2012). Rather, the results suggest the existence of near-integrable parts of phase space that act as barriers for the statistical predictions to be reached. Those near-integrable solutions take the form of interacting vortex blobs (3 or 4 depending on the total angular momentum), perfectly reflecting integrability conditions for the Hamiltonian dynamics on the sphere (Modin & Viviani Reference Modin and Viviani2021).
During our work with Zeitlin's model, a new point of view emerged. More than a spatial discretisation, the Euler–Zeitlin equations themselves provide new mathematical tools 2-D hydrodynamics. Those tools include, in particular, Lie theory for 
$\mathfrak {u}(N)$, which is exceptionally well understood from quantum theory, representation theory and linear algebra. In particular, through the lens of Lie algebra theory, it is natural to split the vorticity matrix 
$\boldsymbol{\mathsf{W}}$ by orthogonal projection onto the stabiliser of the stream matrix 
$\boldsymbol{\mathsf{P}}$ – the underpinning point of this paper. By simulating the Euler–Zeitlin equations (2.3a,b) using our Lie–Poisson integrator, and then for each output computing the canonical splitting, we see that it captures the dynamics of vortex condensation and scale separation, directly related to the theory of Kraichnan (Reference Kraichnan1967) for an inverse energy cascade.
3. Canonical splitting of the vorticity matrix
In this section, we present and discuss canonical vorticity splitting for Zeitlin's model. Here, ‘canonical’ means that the splitting only depends on the Lie algebra structure, the vorticity matrix 
$\boldsymbol{\mathsf{W}}$ and the stream matrix 
$P$. It does not require any ad hoc choice of scale as previous methods do. The scale separation is a result of the dynamics itself.
Consider again the Euler–Zeitlin equations (2.3a,b). Equip 
$\mathfrak {su}(N)$ with the Frobenius inner product. The canonical splitting of the vorticity matrix
\begin{equation} \boldsymbol{\mathsf{W}}=\boldsymbol{\mathsf{W}}_s + \boldsymbol{\mathsf{W}}_r \end{equation}
is obtained by taking 
$\boldsymbol{\mathsf{W}}_s$ to be the orthogonal projection of 
$\boldsymbol{\mathsf{W}}$ onto the stabiliser of the stream matrix 
$\boldsymbol{\mathsf{P}}$
\begin{equation} \mbox{stab}_\boldsymbol{\mathsf{P}} = \lbrace \boldsymbol{\mathsf{A}}\in \mathfrak{su}(N)\mid [\boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{P}}]=0\rbrace . \end{equation}
If 
$\boldsymbol{\mathsf{P}}$ is generic, then all its eigenvalues are distinct, and 
$\mbox {stab}_\boldsymbol{\mathsf{P}}$ is equivalently given by
\begin{equation} \mbox{stab}_\boldsymbol{\mathsf{P}} = \lbrace \boldsymbol{\mathsf{A}}\in \mathfrak{su}(N)\mid \boldsymbol{\mathsf{A}},\boldsymbol{\mathsf{P}}\mbox{ simultaneously diagonalisable}\rbrace. \end{equation}
In this case 
$\boldsymbol{\mathsf{W}}_s$ is obtained via the spectral decomposition: first find 
$\boldsymbol{\mathsf{E}}\in \mathrm {SU}(N)$ that diagonalises 
$\boldsymbol{\mathsf{P}}$, then define 
$\varPi _\boldsymbol{\mathsf{P}}\colon \mathfrak {su}(N)\rightarrow \mbox {stab}_\boldsymbol{\mathsf{P}}$ as
\begin{equation} \boldsymbol{\mathsf{W}}_s := \varPi_\boldsymbol{\mathsf{P}}(\boldsymbol{\mathsf{W}})=\boldsymbol{\mathsf{E}}\,\mbox{diag}(\boldsymbol{\mathsf{E}}^{\dagger} \boldsymbol{\mathsf{WE}})\boldsymbol{\mathsf{E}}^{\dagger}. \end{equation}Remark 3.1 Relative to the splitting (3.1), the Euler–Zeitlin equations (2.3a,b) can be written
\begin{equation} \dot{\boldsymbol{\mathsf{W}}} =[\boldsymbol{\mathsf{P}},\boldsymbol{\mathsf{W}}_r]. \end{equation}
In fact, 
$\boldsymbol{\mathsf{W}} = \boldsymbol{\mathsf{W}}_s + \boldsymbol{\mathsf{W}}_r$ is a steady solution (equilibrium) if and only if 
$\boldsymbol{\mathsf{W}}_r = 0$. So, in a way, the dynamics emerges from the residual part 
$\boldsymbol{\mathsf{W}}_r$.
3.1. Dynamics of $\boldsymbol{\mathsf{W}}_s$ and $\boldsymbol{\mathsf{W}}_r$
For insight into the splitting (3.1) we derive the Euler–Zeitlin equations in the variables 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$. Consider first a general flow on 
$\mathfrak {su}(N)$ of the form
\begin{equation} \dot{\boldsymbol{\mathsf{P}}} = F(\boldsymbol{\mathsf{P}}), \end{equation}
for some smooth vector field 
$F\colon \mathfrak {su}(N)\to \mathfrak {su}(N)$. Let 
$\boldsymbol{\mathsf{E}} \in \mathrm {SU}(N)$ and 
${\boldsymbol{\mathsf{\Lambda}}} \in \mathrm {diag}_N$ denote an eigenbasis and corresponding eigenvalues for 
$\boldsymbol{\mathsf{P}}$. Given (3.6) we first determine the evolution of 
$\boldsymbol{\mathsf{E}}$ and 
$\boldsymbol{\mathsf{\Lambda}}$. The Lie algebra 
$\mathfrak {su}(N)$ is foliated in orbits (cf. Kirillov Reference Kirillov2004) given by
\begin{equation} {O}_{\boldsymbol{\mathsf{P}}} = \lbrace \boldsymbol{\mathsf{UPU}}^{\dagger}{\mid} \boldsymbol{\mathsf{U}}\in \mathrm{SU}(N)\rbrace. \end{equation}
In the generic case, when all eigenvalues are distinct, the tangent space 
$T_{\boldsymbol{\mathsf{P}}}{O}_{\boldsymbol{\mathsf{P}}}$ is spanned by 
$\{ \mathrm {i}e_k e_l^{\dagger} \}_{k\neq l}$, where 
$\boldsymbol{\mathsf{E}}=[e_1,\ldots,e_N]$ is an orthonormal eigenbasis of 
$\boldsymbol{\mathsf{P}}$. The orthogonal directions 
$T_{\boldsymbol{\mathsf{P}}}{O}_{\boldsymbol{\mathsf{P}}}^\bot = \mathrm {span}\{\mathrm {i}e_ke_k^{\dagger} \} = \mbox {stab}_{\boldsymbol{\mathsf{P}}}$ are the linear subspace of matrices in 
$\mathfrak {su}(N)$ sharing the same eigenbasis (they are simultaneously diagonalisable). Thus, the two projections
\begin{equation} \varPi_{\boldsymbol{\mathsf{P}}}\colon \mathfrak{su}(N)\to \mathrm{stab}_{\boldsymbol{\mathsf{P}}} \quad\text{and} \quad \varPi_{\boldsymbol{\mathsf{P}}}^\bot := \mathrm{Id} - \varPi_{\boldsymbol{\mathsf{P}}} \colon \mathfrak{su}(N)\to \mbox{stab}_{\boldsymbol{\mathsf{P}}}^\bot \end{equation}
correspond to decomposition in the basis 
$\{\mathrm {i}e_ke_k^{\dagger} \}_k$ and 
$\{\mathrm {i}e_ke_l^{\dagger} \}_{k\neq l}$. Notice, as expected, that neither 
$\varPi _\boldsymbol{\mathsf{P}}$ nor 
$\varPi _{\boldsymbol{\mathsf{P}}}^\bot$ depend on the eigenvalues of 
$\boldsymbol{\mathsf{P}}$, only on the eigenbasis.
We can now write (3.6) as
\begin{equation} \dot{\boldsymbol{\mathsf{P}}} = \varPi_{\boldsymbol{\mathsf{P}}} F({\boldsymbol{\mathsf{P}}}) + \varPi_{\boldsymbol{\mathsf{P}}}^\bot F({\boldsymbol{\mathsf{P}}}) . \end{equation}
The first part of the flow changes the eigenbasis but not the eigenvalues and vice versa. The question is: What is the generator of 
${\boldsymbol{\mathsf{P}}}\mapsto \varPi _{\boldsymbol{\mathsf{P}}}^\bot F({\boldsymbol{\mathsf{P}}})$? Since it is isospectral it should be of the form 
${\boldsymbol{\mathsf{P}}}\mapsto [B({\boldsymbol{\mathsf{P}}}),{\boldsymbol{\mathsf{P}}}]$ for some 
$B({\boldsymbol{\mathsf{P}}})\in \mathfrak {su}(N)$. Let us denote 
$X=F({\boldsymbol{\mathsf{P}}})$. If the eigenvalues 
$p_1,\ldots,p_N$ of 
${\boldsymbol{\mathsf{P}}}$ are distinct, then
\begin{align} \varPi_{\boldsymbol{\mathsf{P}}}^\bot \boldsymbol{\mathsf{X}} &= \sum_{k\neq l} x_{kl}\mathrm{i} e_ke_l^{\dagger}{=} \sum_{k\neq l} (p_k-p_l)b_{kl}e_ke_l^{\dagger}\nonumber\\ &= \left[\sum_{k\neq l} b_{kl}e_ke_l^{\dagger}, \boldsymbol{\mathsf{P}}\right] = [\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{P}}], \end{align}
where 
$x_{kl}$ are the components of 
$\boldsymbol{\mathsf{X}}$ in the basis 
$\boldsymbol{\mathsf{E}}$, and 
$b_{kl} := x_{kl}/(p_k-p_l)$. Thus, in the generic case the generator 
$B(\boldsymbol{\mathsf{P}})$ is constructed from the eigenvalues 
$p_1,\ldots,p_N$ and the eigenbasis 
$e_1,\ldots,e_N$ of 
$\boldsymbol{\mathsf{P}}$. This observation allows us to write (3.6) in terms of the eigenvalues and eigenbasis of 
$\boldsymbol{\mathsf{P}}$
\begin{equation} \left.\begin{gathered} \dot{p}_k = e_k^{\dagger} \boldsymbol{\mathsf{X}} e_k , \quad \boldsymbol{\mathsf{X}} = F\left(\sum_{k=1}^N p_k e_k e_k^{\dagger} \right)\\ \dot{e}_k = \boldsymbol{\mathsf{B}} e_k, \quad \boldsymbol{\mathsf{B}} = \sum_{k\neq l} \frac{e_k^{\dagger} \boldsymbol{\mathsf{X}} e_l}{p_k-p_l}e_ke_l^{\dagger} . \end{gathered}\right\} \end{equation}
The matrices 
$\boldsymbol{\mathsf{E}}_{kk}=\mathrm {i}e_ke_k^{\dagger} \in \mathfrak {su}(N)$, forming the eigenbasis of 
$\boldsymbol{\mathsf{P}}$, are the quantised analogues of the level curves of the streamfunction 
$\psi$.
We now apply the Lie theory machinery to obtain the dynamics of 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$. From the definition of 
$\boldsymbol{\mathsf{W}}_s$ we have
\begin{align} \dot{\boldsymbol{\mathsf{W}}}_{s} &= \dfrac{{\rm d}}{{\rm d}t}(\boldsymbol{\mathsf{E}}\,\mathrm{diag}(\boldsymbol{\mathsf{E}}^{\dagger} \boldsymbol{\mathsf{WE}})\boldsymbol{\mathsf{E}}^{\dagger})\nonumber\\ &=[\dot{\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger},\boldsymbol{\mathsf{W}}_s] - \varPi_P([\dot{\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger},\boldsymbol{\mathsf{W}}_r]), \end{align}
where we used 
$\varPi _P(\dot {\boldsymbol{\mathsf{W}}})=0$ and 
$\dot {\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger} =-\boldsymbol{\mathsf{E}}\dot {\boldsymbol{\mathsf{E}}^{\dagger} }$. The dynamics for 
$\boldsymbol{\mathsf{P}}$ is similar
\begin{equation} \dot{\boldsymbol{\mathsf{P}}} = [\dot{\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger},\boldsymbol{\mathsf{P}}] + \boldsymbol{\mathsf{E}}\dot{\boldsymbol{\mathsf{\Lambda}}} \boldsymbol{\mathsf{E}}^{\dagger}. \end{equation}
Hence, a formula for 
$\dot {\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger}$ is needed. But we know that the dynamics of 
$\boldsymbol{\mathsf{P}}$ can be orthogonally decomposed as
\begin{equation} \dot{\boldsymbol{\mathsf{P}}} = \varPi_{\boldsymbol{\mathsf{P}}}^\perp(\varDelta_N^{{-}1}[\boldsymbol{\mathsf{P}},\varDelta \boldsymbol{\mathsf{P}}]) + \varPi_{\boldsymbol{\mathsf{P}}}(\varDelta_N^{{-}1}[\boldsymbol{\mathsf{P}},\varDelta \boldsymbol{\mathsf{P}}]), \end{equation}so
\begin{equation} {}[\dot{\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger},\boldsymbol{\mathsf{P}}]= \varPi_{\boldsymbol{\mathsf{P}}}^\perp(\varDelta_N^{{-}1}[\boldsymbol{\mathsf{P}},\varDelta \boldsymbol{\mathsf{P}}]). \end{equation}
Notice that 
$\dot {\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger}$ can be taken in 
$\mbox {stab}_P^\perp$. In fact, the dynamics of 
$\boldsymbol{\mathsf{W}}_s$ remains the same for any 
$\dot {\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger} +\boldsymbol{\mathsf{S}}$ with 
$\boldsymbol{\mathsf{S}}\in \mbox {stab}_{\boldsymbol{\mathsf{P}}}$. The map
\begin{equation} {}[{\cdot},\boldsymbol{\mathsf{P}}]:\mbox{stab}_{\boldsymbol{\mathsf{P}}}^\perp\rightarrow\mbox{stab}_{\boldsymbol{\mathsf{P}}}^{\perp} \end{equation}
is invertible so 
$\dot {\boldsymbol{\mathsf{E}}}\boldsymbol{\mathsf{E}}^{\dagger}$ is uniquely determined in 
$\mbox {stab}_{\boldsymbol{\mathsf{P}}}^\perp$. In conclusion, we have derived the following result for the dynamics of 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$.
Theorem 3.2 Let 
$\boldsymbol{\mathsf{W}} = \boldsymbol{\mathsf{W}}(t)$ and 
$\boldsymbol{\mathsf{P}}=\boldsymbol{\mathsf{P}}(t)$ be the vorticity and stream matrix for a solution to the Euler–Zeitlin equations (2.3a,b). Let 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ respectively be the orthogonal projections of 
$\boldsymbol{\mathsf{W}}$ onto 
$\mbox {stab}_{\boldsymbol{\mathsf{P}}}$ and its orthogonal complement as in (3.1). Then 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ satisfy the following system of equations:
\begin{equation} \left.\begin{gathered} \dot{\boldsymbol{\mathsf{W}}}_{s} = [\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{W}}_s] - \varPi_{\boldsymbol{\mathsf{P}}}[\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{W}}_r]\\ \dot{\boldsymbol{\mathsf{W}}}_r ={-}[\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{W}}_s] + \varPi_{\boldsymbol{\mathsf{P}}}[\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{W}}_r] + [\boldsymbol{\mathsf{P}},\boldsymbol{\mathsf{W}}_r], \end{gathered}\right\} \end{equation}
where 
$\boldsymbol{\mathsf{P}}=\varDelta _N^{-1}(\boldsymbol{\mathsf{W}}_s+\boldsymbol{\mathsf{W}}_r)$ and 
$\boldsymbol{\mathsf{B}}$ is the unique solution to
\begin{equation} \left[\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{P}}\right] = \varPi_{\boldsymbol{\mathsf{P}}}^\bot\varDelta_N^{{-}1}[\boldsymbol{\mathsf{P}},\boldsymbol{\mathsf{W}}_r], \quad \boldsymbol{\mathsf{B}}\in\mbox{stab}_{\boldsymbol{\mathsf{P}}}^\perp. \end{equation} From Theorem 3.2 we can deduce properties of 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$. First, if 
$\boldsymbol{\mathsf{W}}_r=0$ then 
$\boldsymbol{\mathsf{B}}\in \mbox {stab}_{\boldsymbol{\mathsf{P}}}\cap \mbox {stab}_{\boldsymbol{\mathsf{P}}}^\perp$ so 
$\boldsymbol{\mathsf{B}}=0$. Conversely, if 
$\boldsymbol{\mathsf{B}}=0$ then 
$\dot {\boldsymbol{\mathsf{W}}}_s=0$ and
\begin{equation} \dot{\boldsymbol{\mathsf{W}}}_r =[\varDelta_N^{{-}1}(\boldsymbol{\mathsf{W}}_s+\boldsymbol{\mathsf{W}}_r),\boldsymbol{\mathsf{W}}_r]. \end{equation}
Hence, in that case 
$\boldsymbol{\mathsf{W}}_s$ plays the role of a fixed topography, and 
$\boldsymbol{\mathsf{W}}_r$ satisfies the Euler–Zeitlin equation with forcing (3.19). From (3.18a,b) we deduce that 
$\boldsymbol{\mathsf{B}}=0$ also implies
\begin{equation} \mbox{{Tr}}(\varDelta_N^{{-}1}\boldsymbol{\mathsf{W}}_s[\varDelta_N^{{-}1}\boldsymbol{\mathsf{W}}_r,\boldsymbol{\mathsf{W}}_r])=0. \end{equation} Another observation is that if 
$[\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{W}}_s]=0$ then 
$\boldsymbol{\mathsf{B}}=0$ so 
$\dot {\boldsymbol{\mathsf{W}}}_s=0$. Vice versa, 
$\varPi _{\boldsymbol{\mathsf{P}}}[\boldsymbol{\mathsf{B}},\boldsymbol{\mathsf{W}}_r]=0$ implies (3.20). This means that it is possible to have an evolution of the eigenvectors of 
$P$ without any change of the eigenvalues, but not the other way around.
Remark 3.3 The projection 
$\varPi _{\boldsymbol{\mathsf{P}}}\colon \mathfrak {su}(N)\rightarrow \mbox {stab}_P$ has rank 
$N$ in the generic case. Hence, the dynamics of 
$\boldsymbol{\mathsf{W}}_s$ in the moving frame 
${\boldsymbol{\mathsf{E}}}_{kk}$ can be described by only 
$N$ components, namely its eigenvalues. Therefore, the vorticity splitting can be interpreted as a reduced dynamics. As we shall see in § 5 below, the projection 
$\varPi _{\boldsymbol{\mathsf{P}}}$ is a quantised version of a mixing operator. Such operators were used by Shnirelman (Reference Shnirelman1993) to characterise stationary flows.
3.2. Energy and enstrophy splitting
We now study how the canonical splitting (3.1) couples with energy and enstrophy. Since 
$\mbox {{Tr}}(\boldsymbol{\mathsf{P}}\boldsymbol{\mathsf{W}}_r)=0$, the energy, given by the square of the energy norm, fulfils
\begin{align} H(\boldsymbol{\mathsf{W}}) &= \tfrac{1}{2}\mbox{{Tr}}(\boldsymbol{\mathsf{W}}_s\varDelta_N^{{-}1}(\boldsymbol{\mathsf{W}}_s+\boldsymbol{\mathsf{W}}_r))\nonumber\\ &= \tfrac{1}{2}\mbox{{Tr}}(\boldsymbol{\mathsf{W}}_s\varDelta_N^{{-}1}\boldsymbol{\mathsf{W}}_s)-\tfrac{1}{2}\mbox{{Tr}}(\boldsymbol{\mathsf{W}}_r\varDelta_N^{{-}1}\boldsymbol{\mathsf{W}}_r). \end{align}Yet, the enstrophy, corresponding to the enstrophy norm, fulfils
\begin{equation} E(\boldsymbol{\mathsf{W}}) ={-}\mbox{{Tr}}(\boldsymbol{\mathsf{W}}_s^2) - \mbox{{Tr}}(\boldsymbol{\mathsf{W}}_r^2). \end{equation}These equalities furnish the interesting relations
\begin{equation} \left.\begin{gathered} H(\boldsymbol{\mathsf{W}}) = H(\boldsymbol{\mathsf{W}}_s) - H(\boldsymbol{\mathsf{W}}_r)\\ E(\boldsymbol{\mathsf{W}}) = E(\boldsymbol{\mathsf{W}}_s) + E(\boldsymbol{\mathsf{W}}_r). \end{gathered}\right\} \end{equation}
Notice that 
$H(\boldsymbol{\mathsf{W}}_s)$ is always larger than 
$H(\boldsymbol{\mathsf{W}})$ and that the energies of 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ have to increase or decrease at the same rate. On the other hand, if the enstrophy of 
$\boldsymbol{\mathsf{W}}_s$ decrease at some rate then the enstrophy of 
$\boldsymbol{\mathsf{W}}_r$ must increase at the same rate. The canonical splitting is thus coherent with Kraichnan's (Reference Kraichnan1967) description of an inverse energy cascade and a forward enstrophy cascade.
The energy–enstrophy splitting (3.23) has a geometric interpretation. It says that 
$\boldsymbol{\mathsf{W}}$ and 
$\boldsymbol{\mathsf{W}}_r$ are orthogonal in the energy norm, whereas 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ are orthogonal in the enstrophy norm. Consider the plane spanned by 
$\boldsymbol{\mathsf{W}}$ and 
$\boldsymbol{\mathsf{W}}_r$, and let 
$H_r=H(\boldsymbol{\mathsf{W}}_r)$ and 
$H_0=H(\boldsymbol{\mathsf{W}})$. Then, since 
$\boldsymbol{\mathsf{W}}_r = (0,\sqrt {H_r})$ and 
$\boldsymbol{\mathsf{W}}=(\sqrt {H_0},0)$ in this plane, energy corresponds to the Euclidean norm on 
$\mathbb {R}^2$. We want to express the enstrophy norm relative to the energy norm. First, observe that 
$\boldsymbol{\mathsf{W}}_s = (\sqrt {H_0},-\sqrt {H_r})$. The positive definite matrix 
$\boldsymbol{\mathsf{G}}$ for the enstrophy inner product restricted to the 
$(\boldsymbol{\mathsf{W}},\boldsymbol{\mathsf{W}}_r)$-plane can be written 
$\boldsymbol{\mathsf{G}}=\boldsymbol{\mathsf{C}}^\top \boldsymbol{\mathsf{C}}$, where the matrix 
$\boldsymbol{\mathsf{C}}\in \mathbb {R}^{2\times 2}$ is determined by 
$\boldsymbol{\mathsf{C}}\cdot (0,\sqrt {H_r}) = (0,\sqrt {E_0}\sin \alpha )$ and 
$\boldsymbol{\mathsf{C}}\cdot (\sqrt {H_0},0) = \sqrt {E_0}(\cos \alpha,\sin \alpha )$. Here, 
$\alpha$ is the angle between 
$\boldsymbol{\mathsf{W}}$ and 
$\boldsymbol{\mathsf{W}}_s$ in the enstrophy norm, and 
$E_0=E(\boldsymbol{\mathsf{W}})$. Then we have
\begin{equation} \boldsymbol{\mathsf{G}} = \begin{bmatrix}
\dfrac{E_0}{H_0} &
\dfrac{E_0}{\sqrt{H_r}\sqrt{H_0}}\sin^2\alpha \\
\dfrac{E_0}{\sqrt{H_r}\sqrt{H_0}}\sin^2\alpha &
\dfrac{E_0}{H_r}\sin^2\alpha \end{bmatrix}.
\end{equation}
Proposition 3.4 Let 
$\boldsymbol{\mathsf{W}}\neq 0$. Then
\begin{equation} 0< H(\boldsymbol{\mathsf{W}})\leqslant H(\boldsymbol{\mathsf{W}}_s) < E(\boldsymbol{\mathsf{W}}_s)\leqslant E(\boldsymbol{\mathsf{W}}). \end{equation}Moreover, with notation as above
\begin{equation} \left.\begin{gathered} H_r\leqslant E_r=E_0 \sin^2\alpha \\ N^2 H_r \geqslant E_r = E_0 \sin^2\alpha. \end{gathered}\right\} \end{equation}
Hence, if 
$\sin \alpha \neq 0$ then
\begin{equation} \dfrac{E_0}{N^2}\leqslant \dfrac{H_r}{\sin^2\alpha}\leqslant E_0. \end{equation}Proof. First notice that 
$\boldsymbol{\mathsf{W}}_s=0$ if and only if 
$\boldsymbol{\mathsf{W}}=0$ since 
$\sqrt {H({\cdot })}$ is a norm. The inequalities (3.25) then follow from (3.23) and the fact that the enstrophy is always larger than the energy. To get the second inequality of (3.26), we use that the largest eigenvalue of the discrete Laplacian 
$\varDelta _N$ is 
$(N-1)N$.
Remark 3.5 In the limit 
$N\rightarrow \infty$ the ratio 
$\sin ^2\alpha /H_r$ in (3.27) is potentially unbounded. It could, and does, happen that the enstrophy norm of 
$\boldsymbol{\mathsf{W}}_r$ is far from zero while its energy norm tends to zero. This corresponds to 
$\boldsymbol{\mathsf{W}}_r$ being shifted towards small scales while not decreasing in enstrophy. It is a manifestation of Kraichnan's theory of forward enstrophy and inverse energy cascades.
4. Dynamically emerging scale separation
This section contains numerical evidence that the canonical vorticity splitting (3.1) captures the dynamics of scale separation. We also provide theoretical arguments to support these results. To give a complete mathematical proof that fully explains the numerical observations is a challenge not addressed in this paper.
A good reason to study the vorticity splitting 
$\boldsymbol{\mathsf{W}} = \boldsymbol{\mathsf{W}}_s + \boldsymbol{\mathsf{W}}_r$ is that unsteadiness precisely means non-vanishing of 
$\boldsymbol{\mathsf{W}}_r$. From an analytical point of view, 
$\boldsymbol{\mathsf{W}}_s$ represents a projection of 
$\boldsymbol{\mathsf{W}}$ onto a smoother subspace. Indeed, the relation via the Laplace operator between 
$\boldsymbol{\mathsf{P}}$ and 
$\boldsymbol{\mathsf{W}}$ says that 
$\boldsymbol{\mathsf{P}}$ admits two derivatives more than 
$\boldsymbol{\mathsf{W}}$ does. Hence, since 
$\boldsymbol{\mathsf{W}}_s$ is related to 
$\boldsymbol{\mathsf{P}}$ via a polynomial relationship, 
$\boldsymbol{\mathsf{W}}_s$ is, in general, more regular than 
$\boldsymbol{\mathsf{W}}$. Vice versa, 
$\boldsymbol{\mathsf{W}}_r$ contains the rougher part of 
$\boldsymbol{\mathsf{W}}$. The tempting conjecture is that 
$\boldsymbol{\mathsf{W}}_s$ represents the low-dimensional, large-scale dynamics, whereas 
$\boldsymbol{\mathsf{W}}_r$ represents the noisy, small-scale dynamics. To assess this conjecture we present two numerical simulations, both with randomly generated, smooth initial data. These two simulations represent the two generic behaviours described by Modin & Viviani (Reference Modin and Viviani2020a): formation of either 3 or 4 coherent vortex formations, strongly correlated to the momentum–enstrophy ratio. (We have run many more simulation with randomly generated initial conditions. The two simulations presented here capture the universal behaviour in all simulations.)
4.1. Vanishing momentum simulation
This simulation starts from smooth, randomly generated initial data. Each spherical harmonic coefficient 
$\omega ^{lm}$ in the range 
$2\leqslant l \leqslant 10$ was drawn from the standard normal distribution. The remaining coefficients were set to zero. The data were then transformed to a vorticity matrix of size 
$N=512$ (see Modin & Viviani (Reference Modin and Viviani2020a) for details and explicit formulae). For reproducibility, the initial conditions are available in the supplementary material at https://doi.org/10.1017/jfm.2022.457.
For time discretisation of the Euler–Zeitlin equations (2.3a,b) we use the second-order isospectral midpoint method (Modin & Viviani Reference Modin and Viviani2020b). The time step length is 
$1.2239\times 10^{-4}\ \textrm {s}$. This corresponds to 
$0.2$ time units in Zeitlin's model, which scales time by 
$4\sqrt {{\rm \pi} }/N^{3/2}$.
To visualise the fluid motion stages, we sample at the initial time (
$t=0$), at an intermediate time during mixing (
$t=13\ \textrm {s}$) and at a long time well after the large-scale condensates are formed (
$t=318\ \textrm {s}$). The vorticity matrices at these outputs are then transformed to functions in azimuth–elevation coordinates. The result is shown in figure 1.
Due to vanishing momentum (
$\omega ^{1m}=0$ in the initial data), integrability theory for the Hamiltonian dynamics on the sphere suggests that 4 vortex blobs should appear (Modin & Viviani Reference Modin and Viviani2020a,Reference Modin and Vivianib). Indeed they do appear and then pass into quasi-periodic orbits. The first movie of the supplementary material captures the entire process. On top of the large-scale condensates, a noisy, fine-scale structure emerges. In essence, we see a separation of scales.
Figure 2 displays azimuth–elevation fields for the stream matrix 
$\boldsymbol{\mathsf{P}}$ and the vorticity matrix components 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ at the same sampled output times. After some time of initial mixing, the large scales of the vorticity are all contained in 
$\boldsymbol{\mathsf{W}}_s$, whereas 
$\boldsymbol{\mathsf{W}}_r$ collects the small-scale fluctuations. The long-time 
$\boldsymbol{\mathsf{W}}_r$ state resembles noise but, as such, is not completely uniform. In fact, the non-uniform character captures the quasi-periodic dynamics since 
$\dot {\boldsymbol{\mathsf{W}}} = [\boldsymbol{\mathsf{P}},\boldsymbol{\mathsf{W}}_r]$. An open problem is to model the noise by a stochastic process.
Figure 2. Vanishing momentum simulation. Progression of the stream matrix 
$\boldsymbol{\mathsf{P}}$ and the components 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ for the same simulation as in figure 1. Initially, 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ are similar in nature, but they evolve so 
$\boldsymbol{\mathsf{W}}_s$ traps the large-scale condensates whereas 
$\boldsymbol{\mathsf{W}}_r$ captures the small-scale fluctuations.
Figure 3 shows the evolution of the energy and enstrophy of 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$. Over time, 
$E_r$ increases whereas 
$H_r$ decreases. Also notice that 
$E_s$, 
$E_r$, 
$H_s$ and 
$H_r$ fulfil the energy and enstrophy relations (3.23). The residual vorticity 
$\boldsymbol{\mathsf{W}}_r$ decreases in energy norm but increases in enstrophy norm. This is a quantification of Reference KraichnanKraichnan's (Reference Kraichnan1967) qualitative description.
Figure 3. Vanishing momentum simulation. Evolution of the decomposed enstrophies 
$E_s$ and 
$E_r$ (a), and decomposed energies 
$H_s$ and 
$H_r$ (b). The dashed, vertical lines indicate the sample times in figures 1, 2 and 4. The energy 
$H_r$ decays almost to zero so that most of the energy is contained in 
$H_s$ (reflecting the inverse energy cascade). On the other hand, the enstrophy 
$E_r$ increases over time (reflecting the forward enstrophy cascade).
Figure 4. Vanishing momentum simulation. Spectrum in log–log scale for the energies 
$H$, 
$H_s$ and 
$H_r$ at the initial (a), intermediate (b) and long times (c). The dashed lines indicate the slopes 
$l^{-3}$ and 
$l^{-1}$. The slope of 
$H_s$ is almost settled at the intermediate time. The slope of 
$H_r$ takes much longer to settle. At long times, the broken line spectrum of 
$H$ is captured well by the components 
$H_s$ and 
$H_r$, which themselves have almost the same average slope at each scale.
The scale separation of the vorticity is even more clear in the spectral domain. Figure 4 contains energy spectra for 
$\boldsymbol{\mathsf{W}}$, 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ at the sampled output times. The energy level 
$H(l)$, corresponding to the wavenumber 
$l=1,\ldots, N$, contains the energy of the modes for the spherical harmonics 
$Y_{lm}$ with 
$m=-l,\ldots,l$. We notice that the energy spectrum of 
$\boldsymbol{\mathsf{W}}$ is similar in nature to that described by Boffetta & Ecke (Reference Boffetta and Ecke2012). Indeed, the ‘broken line’ slope in the energy spectrum of 
$\boldsymbol{\mathsf{W}}$ originates from an 
$l^{-3}$ slope of 
$\boldsymbol{\mathsf{W}}_s$ and an 
$l^{-1}$ slope of 
$\boldsymbol{\mathsf{W}}_r$. Thus, the vorticity splitting yields a scale separation of the vorticity field that exactly reflects the broken line spectra previously observed in numerical simulations and empirical observations. The broken line spectrum is not reached at the intermediate time before mixing has settled.
4.2. Non-vanishing momentum simulation
In this simulation, the initial data were generated much like before, but now the range of non-zero spherical harmonics coefficients 
$\omega ^{lm}$ is 
$1 \leqslant l \leqslant 10$. Consequently, the total angular momentum is no longer zero. For reproducibility also of this simulation, the generated initial conditions are available in the supplementary material. Time discretisation, step size lengths, etc. are selected as in the previous simulation.
Figure 5 shows azimuth–elevation fields corresponding to the stream matrix 
$\boldsymbol{\mathsf{P}}$, the vorticity matrix 
$\boldsymbol{\mathsf{W}}$ and the components 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$, sampled at initial time (
$t=0$), intermediate time (
$t=13$ s) and long time (
$t=344$ s). The entire motion is captured in the second movie of the supplementary material. Three vortex blobs emerge. The formation of these, from initial data with non-vanishing momentum, is again predicted and demonstrated by Modin & Viviani (Reference Modin and Viviani2020a). It reflects the integrability of the low-dimensional Hamiltonian dynamics on the sphere. As before, the large-scale vorticity patterns are contained in the 
$\boldsymbol{\mathsf{W}}_s$ component. The 
$\boldsymbol{\mathsf{W}}_r$ component swiftly develops noisy fluctuations. At long times, it is less uniform than in the vanishing momentum simulation. The reason is that the 3 blobs here move faster than the 4 blobs in the previous simulation.
Figure 5. Non-vanishing momentum simulation. Progression of the stream matrix 
$\boldsymbol{\mathsf{P}}$ and the components 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ of the vorticity matrix 
$\boldsymbol{\mathsf{W}}$. Initially, 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$ are similar in nature, but they evolve so that 
$\boldsymbol{\mathsf{W}}_s$ contains the large-scale condensates whereas 
$\boldsymbol{\mathsf{W}}_r$ contains the small-scale fluctuations.
The scale separation of the vorticity is again evident from the energy spectrum of 
$\boldsymbol{\mathsf{W}}$. Figure 6 gives energy spectra for the three vorticity fields 
$\boldsymbol{\mathsf{W}}$, 
$\boldsymbol{\mathsf{W}}_s$ and 
$\boldsymbol{\mathsf{W}}_r$. The results are analogous to those in figure 4.
Figure 6. Non-vanishing momentum simulation. Spectrum in log–log scale for the energies 
$H$, 
$H_s$ and 
$H_r$ at the initial (a), intermediate (b) and long times (c). The dashed lines indicate the slopes 
$l^{-3}$ and 
$l^{-1}$. The slope of 
$H_s$ is almost settled at the intermediate time. The slope of 
$H_r$ takes much longer to settle. At long times, the broken line spectrum of 
$H$ is captured well by the components 
$H_s$ and 
$H_r$, which themselves have almost the same average slope at each scale.
4.3. Streamfunction–vorticity branching and blobs
In the literature on 2-D turbulence, steady solutions are often characterised by a functional relation between the streamfunction and vorticity. Branching in such relations has been used as a measure of unsteadiness (Dritschel et al. Reference Dritschel, Qi and Marston2015).
Since 
$\boldsymbol{\mathsf{W}} = \boldsymbol{\mathsf{W}}_s$ if and only if 
$\boldsymbol{\mathsf{W}}$ is a steady solution, it is natural to consider scatter plots between values of 
$\boldsymbol{\mathsf{P}}$ and 
$\boldsymbol{\mathsf{W}}_s$. Such diagrams are given at the initial, intermediate and long times, in figure 7 for the vanishing momentum simulation, and in figure 8 for the non-vanishing momentum simulation. The following interpretation of branches is more fundamental than interpretations related to unsteadiness. Each branch represents and characterises a specific blob in 
$\boldsymbol{\mathsf{W}}_s$. Upward branches represent blobs with positive values. Vice versa for downward branches. This is particularly clear at the long times, where there are fewer blobs. However, the interpretation is valid also at the initial and intermediated times, as revealed by carefully comparing the branching diagrams with the values of 
$\boldsymbol{\mathsf{P}}$ and 
$\boldsymbol{\mathsf{W}}_s$ in figures 2 and 5. The end of each specific branch (sometimes they overlap) is then readily identified with a specific blob. The steepness of the branch corresponds to 
$\textrm {d}\boldsymbol{\mathsf{W}}_s/\textrm {d}\boldsymbol{\mathsf{P}}$. This can be used to determine the shape of the blob, assuming axisymmetry. For example, in figure 8 at long times, the steeper of the two upward branches corresponds to the left-most, sharper of the two positive blobs in figure 5. The two negative blobs at long times in figure 2 are almost indistinguishable, which is reflected as overlapping downward branches in figure 7.
Figure 7. Vanishing momentum simulation. Values of 
$\boldsymbol{\mathsf{P}}$ vs values of 
$\boldsymbol{\mathsf{W}}_s$. The end of each branch corresponds to a blob in the 
$\boldsymbol{\mathsf{W}}_s$ plot in figure 2. Upward if the blob is positive, otherwise downward. For example, at intermediate time, the sharp upward branch close to the 
$y$-axis matches the small positive blob of 
$\boldsymbol{\mathsf{W}}_s$ in the lower-right corner of the corresponding plot in figure 2.
Figure 8. Non-vanishing momentum simulation. Values of 
$\boldsymbol{\mathsf{P}}$ vs values of 
$\boldsymbol{\mathsf{W}}_s$. The end of each branch corresponds to a blob in figure 5. For example, at the initial time, a careful study reveals that the steep, overlapping branches just to the right of the 
$y$-axis match the two positive–negative blob pairs seen in the corresponding plot in figure 5.
5. Canonical splitting of the vorticity function
In this section, we map our results for the Euler–Zeitlin equations to the original Euler equations (1.1a,b). Indeed, all the concepts needed in the canonical splitting for the Euler–Zeitlin equations have classical counterparts; some are listed in table 1. However, one has to be careful to rigorously define these concepts: the transition from the quantised to the classical equations is valid only in a weak sense. Mathematically, the correct framework is 
$L^\infty$ weak-star convergence. Formally, we may nevertheless proceed as follows, keeping in mind that concepts are transferable only in the weak sense. The key ingredient is that the projection 
$\varPi _{\boldsymbol{\mathsf{P}}}$ onto the stabiliser of 
$\boldsymbol{\mathsf{P}}$ corresponds to averaging along the level sets of the streamfunction 
$\psi$. This gives the canonical splitting for the vorticity function via the projector 
$\varPi _\psi$ as
\begin{equation} \omega = \varPi_\psi(\omega) + \varPi_\psi^\bot(\omega) = \omega_s + \omega_r . \end{equation}
The projection 
$\varPi _\psi$ is time dependent, since the level curves of 
$\psi$ change with time.
Table 1. Dictionary between Euler's and Zeitlin's models of hydrodynamics.
Let us now proceed with more mathematical details on this construction. First, recall again Euler's equations in vorticity formulation
\begin{equation} \dot{\omega} =\lbrace\psi,\omega\rbrace, \quad \Delta \psi = \omega. \end{equation}
To define the continuous analogue of the vorticity matrix splitting, we have to understand equations (5.2a,b) in the weak sense. Indeed, in general, the projections 
$\varPi _\psi$ and 
$\varPi _\psi ^\bot$ cannot preserve the smoothness of 
$\omega$. But for any 
$p\in [1,\infty ]$ they are continuous operators from 
$C^0(\mathbb {S}^2)$ to 
$L^p(\mathbb {S}^2)$ with operator norm one. Since continuous functions are dense in 
$L^p(\mathbb {S}^2)$, we can extend 
$\varPi _\psi$ and 
$\varPi _\psi ^\bot$ to continuous operators on 
$L^\infty (\mathbb {S}^2)$. This result fits well with the global well posedness of (5.2a,b), which precisely requires a vorticity function in 
$L^\infty$ (Yudovich Reference Yudovich1963; Majda & Bertozzi Reference Majda and Bertozzi2002; Marchioro & Pulvirenti Reference Marchioro and Pulvirenti2012).
Let us first give (5.2a,b) in the weak sense. For any 
$p \geqslant 2$, if 
$\omega \in L^p(\mathbb {S}^2)$ then 
$\psi \in W^{2,p}(\mathbb {S}^2)\subset H^1(\mathbb {S}^2)$. We define the weak Poisson bracket as
\begin{equation} \int_{\mathbb{S}^2}\lbrace\psi,\omega\rbrace\phi ={-}\int_{\mathbb{S}^2}\omega\lbrace\psi,\phi\rbrace , \end{equation}
for any test function 
$\phi \in C^\infty (\mathbb {S}^2)$. Hence, we define the stabiliser of 
$\psi$ as
\begin{equation} \mbox{stab}_\psi:=\lbrace f\in L^2(\mathbb{S}^2) \mid \lbrace f,\psi\rbrace=0\rbrace. \end{equation}
Next, we define the 
$L^2$ orthogonal projection 
$\varPi _\psi$ onto 
$\mbox {stab}_\psi$.
Proposition 5.1 If 
$p\geqslant 2$ and 
$\psi \in W^{2,p}(\mathbb {S}^2)$ then 
$\mbox {stab}_\psi$ is a closed subspace of 
$L^2(\mathbb {S}^2)$.
Proof. Let 
$\lbrace f_n\rbrace \subset \mbox {stab}_\psi$, such that 
$f_n\rightarrow f$ in 
$L^2$, for 
$n\rightarrow \infty$. We want to show that 
$f\in \mbox {stab}_\psi$. Let 
$\phi$ be a function in 
$C^1(\mathbb {S}^2)$, then we get
\begin{equation} \left|\int \lbrace f,\psi\rbrace\phi \right|= \left|\int \lbrace \phi,\psi\rbrace f \right|=\left|\int \lbrace \phi,\psi\rbrace (f-f_n) \right|\leqslant \|\lbrace \phi,\psi\rbrace\|_0\|f-f_n\|_0\rightarrow 0, \end{equation}
for 
$n\rightarrow \infty$.
The operator 
$\varPi _\psi$ has an explicit form when evaluated on continuous functions. To state it, we first make the following assumption on the critical points of 
$\psi$.
Assumption 5.2 Let 
$\psi \in C^1(\mathbb {S}^2)$ be the streamfunction. Then the critical points of 
$\psi$ define a set of zero Lebesgue measure on 
$\mathbb {S}^2$, such that it is never dense in any neighbourhood of one of its points.
We say that 
$\psi$ is generic whenever it satisfies Assumption 5.2. Consider now some 
$f\in C^1(\mathbb {S}^2)$. Then 
$f\in \mbox {stab}_\psi$ if and only 
$\boldsymbol {\nabla } f$ and 
$\boldsymbol {\nabla }\psi$ are parallel. Since we take 
$\psi$ to be generic, the points where 
$\{ \boldsymbol{x}\mid \boldsymbol {\nabla }\psi (\boldsymbol{x}) = 0\}$ lie on a set of zero measure, nowhere dense. Therefore, since 
$f$ is continuous, 
$f\in \mbox {stab}_\psi$ if it is constant on the connected components of the level curves of 
$\psi$. Then the projection of 
$f$ onto 
$\mbox {stab}_\psi$ can be defined by evaluating 
$f$ on the level curves of 
$\psi$, i.e. the streamlines. Let 
$\boldsymbol{\gamma}$ be a connected component of a streamline, then define the projection 
$\varPi _\psi :C^1(\mathbb {S}^2)\rightarrow \mbox {stab}_\psi$ as
\begin{equation} \left.\varPi_\psi(f)\right|_{\boldsymbol{\gamma}}=\frac{1}{\mbox{length}(\boldsymbol{\gamma})}\int_{\boldsymbol{\gamma}} f \,{\rm d}s. \end{equation}
In the limit when 
$\boldsymbol{\gamma}$ approaches a single point, clearly 
$\varPi _\psi (\omega )_{|\boldsymbol{\gamma} }=f(\boldsymbol{\gamma} )$.
The operator 
$\varPi _\psi$ does not, in general, preserve the continuity of 
$f$. Indeed, consider a bifurcation saddle point 
$\boldsymbol{x}\in \mathbb {S}^2$, namely a saddle point of 
$\psi$ such that the streamline passing through 
$x$ contains a bifurcation point. We then have the following result.
Proposition 5.3 Let 
$\psi$ be generic and 
$\varPi _\psi$ the projection defined in (5.6). If 
$\boldsymbol{x}\in \mathbb {S}^2$ is a bifurcation saddle point for 
$\psi$, then there exists 
$f\in C^1(\mathbb {S}^2)$ such that 
$\varPi _\psi (f)$ is discontinuous at the streamline passing through 
$\boldsymbol{x}$. Vice versa, given 
$f\in C^1(\mathbb {S}^2)$, if 
$\varPi _\psi (f)$ is discontinuous at some point 
$\boldsymbol{x}\in \mathbb {S}^2$, then the streamline passing through 
$\boldsymbol{x}$ contains a bifurcation saddle point.
Proof. (Sketch) The issue about the continuity of 
$f$ can be treated locally. Hence, let us work in Cartesian coordinates. Let 
$\boldsymbol{x}\in \mathbb {S}^2$ be a bifurcation saddle point for 
$\psi$ and 
$\boldsymbol{\gamma}$ the streamline passing through 
$\boldsymbol{x}$. Then, let 
$\beta$ be a curve intersecting 
$\boldsymbol{\gamma}$ only in 
$\boldsymbol{x}$ and let 
$f$ be a smooth function positive at one side of 
$\beta$ and negative at the other one, such that 
$\int _{\boldsymbol{\gamma}} f \,\textrm {d}s=0$. Then, with 
$\boldsymbol{x}$ a bifurcation point, for any neighbourhood 
$U$ of 
$\boldsymbol{x}$, there exist streamlines totally contained in one or another side of 
$\beta$. Then, the average of 
$f$ on those streamlines is either strictly positive or negative, creating a discontinuity at 
$\boldsymbol{\gamma}$.
Vice versa, let 
$f\in C^1(\mathbb {S}^2)$, such that 
$\varPi _\psi (f)$ is discontinuous at some point 
$\boldsymbol{x}\in \mathbb {S}^2$. Then, the streamline passing through 
$x$ cannot be homeomorphic to any of those in some tubular neighbourhood. Hence, the streamline passing through 
$\boldsymbol{x}$ must contain a critical point for 
$\psi$, which also is a bifurcation saddle point.
However, we have the following regularity for 
$\varPi _\psi$.
Proposition 5.4 For any 
$p\in [1,\infty ]$, the operator 
$\varPi _\psi$ can be extended to a bounded operator with norm one on 
$L^p(\mathbb {S}^2)$.
Proof. Let us first notice that 
$\varPi _\psi$ is defined from 
$C^0(\mathbb {S}^2)$ to 
$L^1(\mathbb {S}^2)$, and satisfies 
$\|\varPi _\psi f\|_{L^1}\leqslant \|f\|_{L^1}$, for any 
$f\in C^0(\mathbb {S}^2)$. Since 
$C^0(\mathbb {S}^2)$ is dense in 
$L^1(\mathbb {S}^2)$, it is possible to extend 
$\varPi _\psi$ to a bounded operator on 
$L^1(\mathbb {S}^2)$. Secondly, 
$\varPi _\psi$ is well defined from the space of simple functions to 
$L^\infty (\mathbb {S}^2)$, and satisfies 
$\|\varPi _\psi f\|_{L^\infty }\leqslant \|f\|_{L^\infty }$, for any 
$f$ simple. Since the space of simple functions is dense in 
$L^\infty (\mathbb {S}^2)$, it is possible to extend 
$\varPi _\psi$ to a bounded operator on 
$L^\infty (\mathbb {S}^2)$. Furthermore, since 
$\varPi _\psi$ fixes the constant functions, it must be that 
$\|\varPi _\psi \|_{L^1}=\|\varPi _\psi \|_{L^\infty }=1$. Hence, by the Riesz–Thorin theorem, we conclude that 
$\|\varPi _\psi \|_{L^p}=1$, for any 
$p\in [1,\infty ]$.
Hence, from now on, let us consider (5.2a,b) in the weak form, for 
$\omega \in L^\infty (\mathbb {S}^2)$. It is clear that 
$\varPi _\psi ^2=\varPi _\psi$. Moreover, we can formally define the operator 
$\varPi _\psi$ via the kernel
\begin{equation} K(\boldsymbol{x},\boldsymbol{y}) = \frac{1}{\mbox{length}(\boldsymbol{\gamma}_{\boldsymbol{x}})}\delta_{\boldsymbol{\gamma}_{\boldsymbol{x}}}(\boldsymbol{y}),\quad \text{for any}\ \boldsymbol{x},\boldsymbol{y}\in\mathbb{S}^2 \end{equation}
where 
$\boldsymbol{\gamma} _{\boldsymbol{x}}$ is the connected component of the streamline passing through 
$\boldsymbol{x}$. In this way we get that 
$\varPi _\psi$ is self-adjoint with respect to the 
$L^2$ inner product, i.e. for any 
$f,g\in C^\infty (\mathbb {S}^2)$
\begin{align} \int_{\mathbb{S}^2} f\varPi_\psi g &= \int_{\mathbb{S}^2} f(\boldsymbol{x})\int_{\mathbb{S}^2}K(\boldsymbol{x},\boldsymbol{y}) g(\boldsymbol{y})\nonumber\\ &= \int_{\mathbb{S}^2} g(\boldsymbol{y})\int_{\mathbb{S}^2}K(\boldsymbol{x},\boldsymbol{y}) f(\boldsymbol{x})\nonumber\\ &=\int_{\mathbb{S}^2} g\varPi_\psi f . \end{align}
By extension these equalities are valid for 
$f,g\in L^p(\mathbb {S}^2)$ whenever 
$p\in [1,\infty ]$.
Remark 5.5 We notice that the operator 
$\varPi _\psi$, defined by the kernel 
$K(\boldsymbol{x},\boldsymbol{y})$, is a mixing operator (or polymorphism or bistochastic operator), as introduced by Shnirelman (Reference Shnirelman1993, Reference Shnirelman2013). Such operators give rise to a partial ordering on 
$L^2(\mathbb {S}^2)$: for any 
$f,g\in L^2(\mathbb {S}^2)$, 
$f\dashv g$ if there exists a mixing operator, defined via the kernel 
$K$, such that 
$f=K\ast g$. In his work, Shnirelman shows that stationary flows are characterised as minimal elements of this ordering. In a way, our work shows that it is enough to consider mixing operators of the form 
$\varPi _\psi$. Within this class, 
$\omega$ is minimal if there exists a streamfunction such that 
$\varPi _\psi \omega =\omega$. As we see next, this in turn implies that 
$\omega$ is stationary.
Proposition 5.6 Let 
$\psi \in C^1(\mathbb {S}^2)$ be generic. For 
$\omega \in L^\infty (\mathbb {S}^2)$ we then have
\begin{equation} \omega\in\mbox{stab}_\psi \iff \varPi_\psi \omega=\omega . \end{equation}Proof. We prove the result for 
$\omega \in C^1(\mathbb {S}^2)$, then conclude by extension. Let 
$\omega \in \mbox {stab}_\psi$. Then, 
$\boldsymbol {\nabla }\omega$ is parallel to 
$\boldsymbol {\nabla }\psi$ almost everywhere. Hence, the gradient of 
$\omega$ along any streamline must vanish, and so on each connected component it is constant. By continuity of 
$\omega$ we deduce that 
$f$ must be constant also on the streamlines containing critical points. Therefore, 
$\varPi _\psi \omega =\omega$. Assume now that 
$\varPi _\psi \omega =\omega$. Then 
$\omega$ must be constant on each connected component of a streamline. Hence, 
$\boldsymbol {\nabla } \omega$ is orthogonal to the streamlines. Since 
$\boldsymbol {\nabla }\psi$ is also orthogonal to the streamlines, we conclude that 
$\lbrace \omega,\psi \rbrace =0$, i.e. 
$\omega \in \mbox {stab}_\psi$.
We are now in the position to derive continuous analogues of the results in § 3 (which, remember, are based on Lie theory for matrices). First, the streamfunction 
$\psi$ satisfies the equation
\begin{equation} \dot{\psi} =\varDelta^{{-}1}\lbrace\psi,\Delta \psi\rbrace. \end{equation}
This equation is not Hamiltonian. But we can split the right-hand side into a Hamiltonian and non-Hamiltonian part via the projection 
$\varPi _\psi$
\begin{equation} \dot{\psi} = \varPi_\psi^\bot\varDelta^{{-}1}\lbrace\psi,\Delta \psi\rbrace + \varPi_\psi\varDelta^{{-}1}\lbrace\psi,\Delta \psi\rbrace. \end{equation}
Analogous to the quantised case, we seek a generator for 
$\varPi _\psi ^\bot \varDelta ^{-1}\lbrace \psi,\Delta \psi \rbrace$. That is, a function 
$b\colon \mathbb {S}^2\to \mathbb {R}$ such that
\begin{equation} \lbrace b,\psi\rbrace = \varPi_\psi^\bot\varDelta^{{-}1}\lbrace\psi,\Delta \psi\rbrace. \end{equation}
It is clear that a necessary condition for the equation 
$\lbrace b,\psi \rbrace = f$ to have a solution 
$b$ is that 
$f\in \mbox {stab}_\psi ^\perp$. Indeed, if 
$b\in C^1(\mathbb {S}^2)$ then 
$\lbrace b,\psi \rbrace \in \mbox {stab}_\psi ^\perp$ since
\begin{equation} \int_{\mathbb{S}^2} \lbrace b,\psi\rbrace g ={-}\int_{\mathbb{S}^2} \lbrace g,\psi\rbrace b = 0 \quad \text{for any} \ g\in\mbox{stab}_\psi . \end{equation}
However, in general (5.12) can be solved only where 
$\boldsymbol {\nabla }\psi \neq 0$. Around the critical points of 
$\psi$ the gradient of 
$b$ is potentially unbounded. Moreover, the right-hand side in (5.12) can be discontinuous at the level curves of 
$\psi$ containing saddle points of 
$\psi$. Hence, 
$b$ can only be defined almost everywhere. Furthermore, we have the following:
Proposition 5.7 Let 
$f\in C^0(\mathbb {S}^2)$ and 
$\psi$ be generic. Then 
$f\in \mbox {stab}_\psi ^\perp$ if and only if there exists 
$b$ almost everywhere smooth, such that 
$\lbrace b,\psi \rbrace = f$ on the set 
$\boldsymbol {\nabla }\psi \neq 0$.
Proof. The if part is clear. Let instead take 
$f\in \mbox {stab}_\psi ^\perp$. Then, for any point 
$\boldsymbol{x}\in \mathbb {S}^2$, we have to solve the partial differential equation (PDE) for 
$b$
\begin{equation} \boldsymbol{\nabla}^\perp\psi\boldsymbol{\cdot}\boldsymbol{\nabla} b = f, \end{equation}
where 
$\boldsymbol {\nabla }^\perp \psi =\boldsymbol{x}\times \boldsymbol {\nabla }\psi$. If 
$\boldsymbol {\nabla }\psi$ does not vanish, (5.14) can be solved by integration in the direction of 
$\boldsymbol {\nabla }^\perp \psi$. At the points where 
$\boldsymbol {\nabla }\psi$ does not vanish, 
$\boldsymbol {\nabla } b$ is not defined by (5.14), and it can be unbounded around those points. Hence, the field 
$b$ is almost everywhere smooth and satisfies 
$\lbrace b,\psi \rbrace = f$, where 
$\boldsymbol {\nabla }\psi \neq 0$.
5.1. Dynamics of $\omega _s$ and $\omega _r$
To derive the dynamical equations for 
$\omega _s$, we cannot directly define the field 
$b$ corresponding to the quantised field 
$B$ above. Instead, we consider the volume-preserving vector field
\begin{equation} X[\psi]:=\varPi_\psi^\bot\varDelta^{{-}1}\lbrace\psi,\Delta \psi\rbrace. \end{equation}
We note that 
$X$ corresponds to the infinitesimal action of a map 
$\varphi _t$ which transports 
$\psi$ by deforming its level curves. Hence, 
$\varphi _t$ acts naturally on 
$\mbox {stab}_\psi$. Let us write 
$\varPi _\psi ^t$ for the projection onto 
$\mbox {stab}_\psi$ at time 
$t$. Then, for any point 
$\boldsymbol{x}\in \mathbb {S}^2$, let
\begin{equation} \boldsymbol{\gamma}(t)=\lbrace \boldsymbol{y}\mid \psi(\boldsymbol{y})=\psi(\boldsymbol{x})\mbox{ and } \boldsymbol{y} \hbox{ in connected component of } \boldsymbol{x} \rbrace , \end{equation}
and 
$\textrm {d}\widehat {s_t}$ be the normalised Lebesgue measure on 
$\gamma (t)$. We have then the formal identity
\begin{align} \omega_s(t,\boldsymbol{x}) &= \varPi_\psi^t\omega(t,\boldsymbol{x}) \nonumber\\ &= \int_{\boldsymbol{\gamma}(t)}\omega(t,\boldsymbol{y}) \,{\rm d}\widehat{s_t}(\boldsymbol{y}) \nonumber\\ &= \int_{\boldsymbol{\gamma}(0)}(\varphi_t^*\omega)(t,\boldsymbol{y}) \,{\rm d}\widehat{s_0}(\boldsymbol{y})\nonumber\\ &= \varPi_\psi^0(\varphi_t^*\omega)(t,\varphi_t^{{-}1}(x)). \end{align}
Hence, for any test function 
$\phi \in C^\infty (\mathbb {S}^2)$
\begin{align} \dfrac{{\rm d}}{{\rm d}t}\int_{\mathbb{S}^2}\omega_s(t,\boldsymbol{x})\phi(\boldsymbol{x}) &= \dfrac{{\rm d}}{{\rm d}t}\int_{\mathbb{S}^2}\varPi_\psi^0(\varphi_t^*\omega)(t,\varphi_t^{{-}1}(\boldsymbol{x}))\phi(\boldsymbol{x}) \nonumber\\ &=\int_{\mathbb{S}^2}\left(\varPi_\psi^0(\varphi_t^*\mathcal{L}_X\omega)(t,\varphi_t^{{-}1}\boldsymbol{x}) - \mathcal{L}_X\varPi_\psi^0(\varphi_t^*\omega)(t,\varphi_t^{{-}1}\boldsymbol{x})\right)\phi(\boldsymbol{x})\nonumber\\ &={-}\int_{\mathbb{S}^2}\omega(t,\boldsymbol{x}) \mathcal{L}_X\varPi_\psi^t\phi(\boldsymbol{x}) + \mathcal{L}_X\varPi_\psi^t\omega(t,\boldsymbol{x})\phi(\boldsymbol{x}) , \end{align}
where 
$\mathcal {L}_X$ is the Lie derivative, which simply acts on functions as 
$\mathcal {L}_Xf=X[f]$. We notice from the previous calculations that 
$\mathcal {L}_X$ has to be evaluated only on elements in 
$\mbox {stab}_\psi$. Hence, the time derivative of 
$\omega _s$ is well defined in the weak sense.
Let us now formally denote 
$X:=-\lbrace b,{\cdot }\rbrace$. Then, interpreting the Poisson bracket in the weak sense, we can write the dynamical system for 
$\omega _s$ and 
$\omega _r$ as
\begin{equation} \left.\begin{gathered} \dot{\omega_s} = \lbrace b,\omega_s\rbrace - \varPi_\psi\lbrace b,\omega_r\rbrace\\ \dot{\omega_r} ={-}\lbrace b,\omega_s\rbrace + \varPi_\psi\lbrace b,\omega_r\rbrace + \lbrace\psi,\omega_r\rbrace\\ \lbrace b,\psi\rbrace = \varPi_\psi^\bot\varDelta^{{-}1}\lbrace\psi,\omega_r\rbrace, \end{gathered}\right\} \end{equation}
where 
$b$ is implicitly defined by the third equation and 
$\psi =\varDelta ^{-1}(\omega _s+\omega _r)$. We notice that the equations of motion for 
$\omega _s$ can also be written in a more compact form as
\begin{equation} \dot{\omega_s}=[\varPi_\psi,\mathcal{L}_X]\omega, \end{equation}where the square bracket is the commutator of operators.
Finally, notice that the energy and enstrophy splitting is valid also in the classical setting
\begin{equation} \left.\begin{gathered} H(\omega) = H(\omega_s) - H(\omega_r)\\ E(\omega) = E(\omega_s) + E(\omega_r). \end{gathered}\right\} \end{equation}6. Conclusions and outlook
Based on the Euler–Zeitlin equations we have reported on a new technique for studying 2-D turbulence via canonical splitting of vorticity. In numerical simulations, this splitting dynamically develops into a separation of scales. These numerical results are supported by some theoretical evidence. To develop a full mathematical understanding of these phenomena, even completely within the setting of Zeitlin's model, is an interesting open problem. We have further presented mathematical foundations for a weak 
$L^\infty$ theory in the continuous setting, independent of Zeitlin's model.
As the numerical simulations so strikingly capture the scale separation process, and as the inverse relations of the corresponding energy–enstrophy splitting reflect the stationary theory of Kraichnan, further numerical and theoretical studies of the canonical vorticity splitting shall likely unveil more details on the mechanism behind vortex condensation. Furthermore, the splitting into large scales 
$\omega _s$ and small scales 
$\omega _r$ suggests using these variables as a basis for stochastic model reduction (cf. Jain, Timofeyev & Vanden-Eijnden Reference Jain, Timofeyev and Vanden-Eijnden2015) of the 2-D Euler equations, with 
$\omega _r$ modelled as multiplicative noise.
Supplementary material and movies
Supplementary material and movies are available at https://doi.org/10.1017/jfm.2022.457.
Acknowledgements
The computations were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC) at C3SE partially funded by the Swedish Research Council through grant agreement no. 2018-05973.
Funding
This work was supported by the Swedish Research Council (K.M., grant number 2017-05040); the Knut and Alice Wallenberg Foundation (K.M., grant number WAF2019.0201), (M.V., grant number KAW2020.0287); and the research grant Junior Visiting Fellowship by the Scuola Normale Superiore of Pisa, Italy.
Declaration of interests
The authors report no conflict of interest.
Open access
