A Casimir preserving scheme for long-time simulation of spherical ideal hydrodynamics

Abstract

The incompressible two-dimensional Euler equations on a sphere constitute a fundamental model in hydrodynamics. The long-time behaviour of solutions is largely unknown; statistical mechanics predicts a steady vorticity configuration, but detailed numerical results in the literature contradict this theory, yielding instead persistent unsteadiness. Such numerical results were obtained using artificial hyperviscosity to account for the cascade of enstrophy into smaller scales. Hyperviscosity, however, destroys the underlying geometry of the phase flow (such as conservation of Casimir functions), and therefore might affect the qualitative long-time behaviour. Here, we develop an efficient numerical method for long-time simulations that preserve the geometric features of the exact flow, in particular conservation of Casimirs. Long-time simulations on a non-rotating sphere then reveal three possible outcomes for generic initial conditions: the formation of either 2, 3 or 4 coherent vortex structures. These numerical results contradict the statistical mechanics theory and show that previous numerical results, suggesting 4 coherent vortex structures as the generic behaviour, display only a special case. Through integrability theory for point vortex dynamics on the sphere we present a theoretical model which describes the mechanism by which the three observed regimes appear. We show that there is a correlation between a first integral $\unicode[STIX]{x1D6FE}$ (the ratio of total angular momentum and the square root of enstrophy) and the long-time behaviour: $\unicode[STIX]{x1D6FE}$ small, intermediate and large yields most likely 4, 3 or 2 coherent vortex formations. Our findings thus suggest that the likely long-time behaviour can be predicted from the first integral $\unicode[STIX]{x1D6FE}$.

1 Introduction

The motion of an ideal fluid restricted to the surface of a sphere is a fundamental model in oceanography, meteorology and astrophysics (see Majda & Bertozzi (2002); Dolzhansky (2012); Pedlosky (2013); Zeitlin (2018), and references therein). The equations of motion, first studied by Euler in 1757, encode a rich geometry – a Lie–Poisson structure – which results in conservation of energy, momentum and Casimir functions (see Arnold 1966; Marsden & Weinstein 1983; Arnold & Khesin 1998).

The ultimate ‘fate’ of two-dimensional (2-D) fluid motion in a bounded domain is largely unknown (Newton 2016). Statistical mechanics theories, such as developed by Miller (1990) and Robert & Sommeria (1991), are based on maximizing the entropy of a coarse-grain probability distribution of the macroscopic states under conservation of energy and (at least some of) the Casimirs. Such models predict a steady equilibrium of large-scale coherent vortex structures, with a functional relation between vorticity and streamfunction.

To test the statistical model of Miller, Robert & Sommeria (MRS) a natural approach is to use long-time numerical simulations. A serious complication is the ‘inverse energy cascade’ where energy from small scales is eventually fed into large scales whereas enstrophy cascades in the forward direction towards smaller scales. This process was first described by Kraichnan (1967). Of course, in a numerical simulation the spatial resolution is finite, so one can never fully resolve the fine-scale structure. As a remedy, a common approach is to adopt a subgrid model, most often hyperviscosity, to account for the enstrophy cascade to smaller scales (see Qi & Marston (2014) and references therein). The inverse energy cascade is related to the conservation of Casimirs, although the exact relation is unknown. In addition to energy, circulation (linear Casimir) and enstrophy (quadratic Casimir), there are several numerical investigations reporting that cubic and possibly higher-order Casimirs also play a role in the formation of large-scale coherent vortex structures (Abramov & Majda 2003; Dubinkina & Frank 2010). On the non-rotating sphere, Dritschel, Qi & Marston (2015) provided numerical evidence that, for randomly generated initial data, the long-time behaviour results in a non-steady interaction largely between two positive and two negative coherent vortex structures (referred to as vortex blobs in this paper) essentially governed by finite-dimensional point vortex dynamics. Seemingly persistent unsteadiness in numerical solutions of 2-D Euler fluids was also reported by Segre & Kida (1998) but for special initial conditions. Dritschel et al. (DQM) argue that, in fact, the unsteady four vortex blob behaviour is generic. This statement is in stark contrast to the previous notion that a steady equilibrium is the generic behaviour. However, DQM used methods with hyperviscosity and in their simulations the percentage decay in enstrophy is between 30 % and 60 %, so hyperviscosity clearly comes into play, but precisely how and if it affects the long-time result is unclear.

In this paper, based on a new numerical method that exactly conserves all Casimir functions thereby eliminating the need for hyperviscosity, we give strong evidence that neither MRS nor DQM are correct. Or, in a way, they are both correct – it all depends on the regime of the initial conditions. Based on the non-dimensional non-negative number $\unicode[STIX]{x1D6FE}$ given by the quotient between the total angular momentum and the total enstrophy, we identify three different regimes: generically, (‘generic’ here means that the initial vorticity is sampled as a random field in the space of $L^{2}$ functions, as described in § 3.2 below) if $\unicode[STIX]{x1D6FE}\lesssim 0.15$ then most likely 4 vortex blobs form (the behaviour observed by DQM), if $\unicode[STIX]{x1D6FE}\gtrsim 0.4$ then most likely 2 vortex blobs form (the behaviour suggested by MRS) and if $0.15\lesssim \unicode[STIX]{x1D6FE}\lesssim 0.4$ we have found a new, intermediate regime where most likely 3 vortex blobs form. The 2 vortex blob formation is steady (or at least almost steady), whereas the 3 and 4 blob formations are unsteady. Furthermore, through point vortex dynamics, we suggest a theoretical mechanism which qualitatively explains the three regimes. This theory, which also predicts results observed on the torus, is not based on statistical mechanics (i.e. maximizing entropy, like MRS) but rather on integrability theory (results on quasi-periodicity) for point vortex dynamics.

As mentioned already, the central tool in the discovery of the three regimes is a new numerical scheme for ideal fluids on rotating or non-rotating spheres that encapsulate the full Lie–Poisson geometry (in particular conservation of associated Casimirs). (It is clear from the definition of $\unicode[STIX]{x1D6FE}$ that a scheme with hyperviscosity, such as those used by Dritschel et al. (2015) with 30 %–60 % decay in enstrophy but no decay in total angular momentum, could never be used to correctly identify the three regimes.) It is based on geometric quantization theory developed by Hoppe (1982), Hoppe & Yau (1998) in conjunction with the Lie–Poisson preserving numerical time discretization developed by Modin & Viviani (2019). The method can be seen as a spherical analogue of the spatial discretization of the Euler equations on the torus suggested by Zeitlin (1991) and the associated numerical time discretization suggested by McLachlan (1993).

We now continue the introduction with a more detailed exposition of the equations of motion, an overview of the space and time discretization and a summary of our main findings.

Consider a homogeneous, incompressible, inviscid, two-dimensional fluid, constrained to the unit sphere $\mathbb{S}^{2}$ embedded in Euclidean $\mathbb{R}^{3}$ and possibly rotating with constant angular speed about a fixed axis. The equations of motion are given by Euler’s equations of hydrodynamics

(1.1)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\dot{\boldsymbol{v}}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\boldsymbol{v}=-\unicode[STIX]{x1D735}p-2\widetilde{\unicode[STIX]{x1D734}}\times \boldsymbol{v},\\ \unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v}=0,\end{array}\right\}\end{eqnarray}$$

where $v$ is the velocity vector field of the fluid, $p$ is its internal pressure and $\widetilde{\unicode[STIX]{x1D734}}=(\unicode[STIX]{x1D734}\boldsymbol{\cdot }\boldsymbol{n})\boldsymbol{n}$ is the projection of the angular rotation vector $\unicode[STIX]{x1D6FA}\in \mathbb{R}^{3}$ to the normal $n$. The term $-2\widetilde{\unicode[STIX]{x1D734}}\times \boldsymbol{v}$ is due to the Coriolis force. Equivalent to (1.1) is the barotropic vorticity equation (also called the quasi-geostrophic equation in the case $\widetilde{\unicode[STIX]{x1D734}}\neq 0$), formulated in terms of the vorticity variable $\unicode[STIX]{x1D714}=(\unicode[STIX]{x1D735}\times \boldsymbol{v})\boldsymbol{\cdot }\boldsymbol{n}$. By Stokes’ theorem we necessarily have $\int \unicode[STIX]{x1D714}=0$ corresponding to zero circulation. Euler’s equations (1.1) can now be written

(1.2)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\dot{\unicode[STIX]{x1D714}}=\{\unicode[STIX]{x1D713},\unicode[STIX]{x1D714}\},\\ \unicode[STIX]{x0394}\unicode[STIX]{x1D713}=\unicode[STIX]{x1D714}-f,\end{array}\right\}\end{eqnarray}$$

where $f:=2\unicode[STIX]{x1D734}\boldsymbol{\cdot }\boldsymbol{n}$ is the Coriolis parameter, $\unicode[STIX]{x1D6E5}$ is the Laplace–Beltrami operator, $\{\cdot ,\cdot \}$ is the Poisson bracket and the streamfunction $\unicode[STIX]{x1D713}$ is unique by the additional condition

(1.3)$$\begin{eqnarray}\int \unicode[STIX]{x1D713}=0.\end{eqnarray}$$

The vorticity equation (1.2) constitutes an infinite-dimensional Lie–Poisson system (cf. Arnold & Khesin 1998) on the space of smooth zero mean functions

(1.4)$$\begin{eqnarray}C_{0}^{\infty }(\mathbb{S}^{2})=\left\{\unicode[STIX]{x1D714}\in C^{\infty }(\mathbb{S}^{2})\left|\int \unicode[STIX]{x1D714}=0\right.\right\}.\end{eqnarray}$$

The Hamiltonian is

(1.5)$$\begin{eqnarray}H(\unicode[STIX]{x1D714})={\textstyle \frac{1}{2}}\int (\unicode[STIX]{x1D714}-f)\unicode[STIX]{x1D713},\end{eqnarray}$$

and the (infinitely many) Casimir functions are given, for any smooth real function $g\in C^{\infty }(\mathbb{R})$, by

(1.6)$$\begin{eqnarray}{\mathcal{C}}(\unicode[STIX]{x1D714})=\int g(\unicode[STIX]{x1D714}).\end{eqnarray}$$

Often $g$ is chosen as monomials, and the corresponding Casimirs

(1.7)$$\begin{eqnarray}{\mathcal{C}}_{k}(\unicode[STIX]{x1D714})=\int \unicode[STIX]{x1D714}^{k},\quad k=1,2,\ldots\end{eqnarray}$$

are called linear, quadratic, cubic, etc. Each Casimir (1.6) is indeed a first integral

(1.8)$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\int g(\unicode[STIX]{x1D714})=-\int g^{\prime }(\unicode[STIX]{x1D714})\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}=-\int \boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}g(\unicode[STIX]{x1D714})=\int (\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{v})g(\unicode[STIX]{x1D714})=0,\end{eqnarray}$$

where we have used that

(1.9)$$\begin{eqnarray}\{\unicode[STIX]{x1D713},\unicode[STIX]{x1D714}\}_{p}=\boldsymbol{p}\boldsymbol{\cdot }(\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{p}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D714}_{p})=(\boldsymbol{p}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}_{p})\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}_{p}=-\boldsymbol{v}_{p}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D714}_{p}\end{eqnarray}$$

for any $\unicode[STIX]{x1D713},\unicode[STIX]{x1D714}\in C^{\infty }(\mathbb{S}^{2})$ and any $\boldsymbol{p}\in \mathbb{S}^{2}$. Notice, in particular, that the Casimirs are conserved for any choice of Hamiltonian; this reflect the underlying Lie–Poisson geometry which is foliated in co-adjoint orbits preserved by any Hamiltonian flow (cf. Marsden & Ratiu 1999, ch. 13–14).

The traditional approach to numerical discretization of a partial differential equation (PDE) is to construct schemes of high local order of accuracy, using for example finite element or finite volume schemes. Rather than focusing on local accuracy, we take here conservation of the Casimir functions (1.6) and the underlying geometric structure as a guiding principle for spatial discretization: we wish to replace the infinite-dimensional Lie–Poisson structure $(C_{0}^{\infty }(\mathbb{S}^{2}),\{\cdot ,\cdot \})$ by a finite-dimensional analogue. We require the number of conserved quantities to increase with the size of the spatial discretization. This cannot be achieved by a truncated spectral decomposition of the vorticity, essentially because the space spanned by a truncated spectral basis is not closed under the Poisson bracket. Instead, we take the approach proposed by Zeitlin (2004) based on the theory of geometric quantization studied in Hoppe (1982), Bordemann et al. (1991) and Bordemann, Meinrenken & Schlichenmaier (1994). It provides a sequence, indexed by $N=1,2,\ldots \,$, of finite-dimensional Lie algebras, that converges to the infinite-dimensional Lie algebra of smooth functions on the sphere as $N\rightarrow \infty$. The sequence is given explicitly by the Lie algebra $\mathfrak{s}\mathfrak{u}(N)$ (or $\mathfrak{s}\mathfrak{l}(N,\mathbb{C})$) ($\mathfrak{s}\mathfrak{u}(N)$ is the Lie algebra of $N\times N$ skew-Hermitian complex matrices with trace zero, $\mathfrak{s}\mathfrak{l}(N,\mathbb{C})$ is the Lie algebra of $N\times N$ complex matrices with trace zero) for $N=1,2,\ldots \,$. For any choice of $N$ we get an ordinary differential equation (ODE) which is a finite-dimensional analogue of (1.2)

(1.10)$$\begin{eqnarray}\dot{\unicode[STIX]{x1D652}}=[\unicode[STIX]{x1D6E5}_{N}^{-1}(\unicode[STIX]{x1D652}-\unicode[STIX]{x1D641}),\unicode[STIX]{x1D652}]_{N},\end{eqnarray}$$

where $\unicode[STIX]{x1D652}\in \mathfrak{s}\mathfrak{u}(N)$ (corresponding to the vorticity $\unicode[STIX]{x1D714}$), $\unicode[STIX]{x1D641}\in \mathfrak{s}\mathfrak{u}(N)$ (corresponding to the Coriolis parameter $f$), $\unicode[STIX]{x1D6E5}_{N}:\mathfrak{s}\mathfrak{u}(N)\rightarrow \mathfrak{s}\mathfrak{u}(N)$ is the discrete Laplace–Beltrami operator (corresponding to $\unicode[STIX]{x1D6E5}$) and $[\cdot ,\cdot ]_{N}$ is the rescaled matrix commutator (corresponding to $\{\cdot ,\cdot \}$). The matrix differential equation (1.10) is an isospectral flow, meaning that the eigenvalues of $W$ are invariant in time. The conservation of these eigenvalues corresponds to the conservation of the Casimirs. Exactly how $\unicode[STIX]{x1D652}$ in (1.10) approximates $\unicode[STIX]{x1D714}$ in (1.2) is described in a complicated (but explicit) linear change of coordinates between $W$ and a truncated spherical harmonics basis. Details are given in § 2.1. A feature of the spatial discretization is that $\unicode[STIX]{x1D652}\mapsto \unicode[STIX]{x1D6E5}_{N}^{-1}(\unicode[STIX]{x1D652}-\unicode[STIX]{x1D641})$ can be computed in only $O(N^{2})$ operations. Thus, the main computational complexity is due to matrix multiplications in the bracket $[\cdot ,\cdot ]$ (which has complexity $O(N^{3})$). Details on the computational complexity are given in § 2.3.

To discretize (1.10) in time we apply a Lie–Poisson preserving isospectral symplectic Runge–Kutta integrator (Modin & Viviani 2019). These numerical methods exactly conserve (i.e. up to rounding errors) the discrete Casimirs (eigenvalues), they nearly conserve the Hamiltonian (‘nearly’ in the sense of backward error analysis of symplectic integrators, cf. Hairer, Lubich & Wanner (2006)), and they exactly conserve the Lie–Poisson flow structure (in short, this means that the time discretized system corresponds to a continuous Lie–Poisson flow on $\mathfrak{s}\mathfrak{u}(N)$ for a slightly modified Hamiltonian). The IsoSRK integrators are necessarily implicit, thus requiring nonlinear root solving at each time step. As a comparison, we also employ the standard explicit Heun method for time discretization of (1.10).

In § 3 we present numerical simulations on a non-rotating sphere ($\unicode[STIX]{x1D641}=0$). First, in § 3.1, we use the same randomly generated initial data as suggested by Dritschel et al. (2015). Long-time simulations are carried out for both types of time discretizations (IsoSRK and Heun) and various levels of spatial discretization. Our numerical results verify, but now without hyperviscosity, the formation of a quadruple of vortex blobs moving quasi-periodically with no sign of reaching steadiness. However, although the DQM initial conditions were randomly generated, we claim they cannot represent the generic behaviour because the total angular momentum is zero. The motivation by Dritschel et al. to set momentum to zero was ‘to avoid starting with a flow organized at the largest possible scale’. Herein lies the implicit assumption that the value of the momentum does not affect the qualitative behaviour. On the doubly periodic square (i.e. the flat torus) the assumption is correct: momentum does not influence the dynamics and can therefore safely be set to zero. On the sphere, however, the momentum strongly affects the dynamics. In fact, our results suggest that the generic qualitative behaviour on a non-rotating sphere is essentially governed by the value of the total angular momentum. Indeed, in § 3.2 we generate 16 sets of initial vorticity as samples from a Gaussian random field on the space of $L^{2}$-functions. In the corresponding 16 long-time simulations we observe the following qualitative behaviour: 5 of them give 4 vortex blobs, 9 of them give 3 vortex blobs and 2 of them give 2 vortex blobs. We also observe that the non-dimensional number $\unicode[STIX]{x1D6FE}=\Vert \boldsymbol{L}\Vert /(R\sqrt{{\mathcal{C}}_{2}})$ (total angular momentum divided by the radius of the sphere times the square root of enstrophy) gives a probabilistic indication of which ‘qualitative regime’ the fluid configuration develops into: small values (approximately less than 0.15) result in 4 vortex blobs, large values (approximately larger than 0.4) result in 2 vortex blobs and intermediate values result in 3 vortex blobs. The number $\unicode[STIX]{x1D6FE}$, computable from the initial conditions, is thus implicated in predicting the fluid’s long-time qualitative behaviour. Of the three regimes, only the 2 vortex formation is steady (up to a constant speed rotation about the momentum axis).

It is natural to ask for a theoretical model explaining the three observed regimes. Clearly, the statistical mechanics based MRS theory is insufficient; it incorrectly predicts steadiness and does not predict or offer insights into why there should be three regimes. Instead, we have found a different theory which explains the mechanisms by which the regimes appear: it is closely related to integrability theory for point vortex dynamics (PVD). Recall that a Hamiltonian system is called integrable if there is a local change of variables in which the dynamics is described by quasi-periodic linear motion on tori. (Equivalently, integrability of a $2n$-dimensional Hamiltonian system can be characterized by the existence of $n$ first integrals in involution (cf. Arnold 1989).) PVD constitute a class of Hamiltonian $N$-particle systems that describe, at least formally, special solutions to the Euler equations (1.1) in the non-rotating case ($\unicode[STIX]{x1D6FA}=0$). Aref (2007a) refers to PVD as ‘a classical mathematics playground’: although the connection to fluid mechanics has always remained in the background, mathematicians have studied these finite-dimensional Hamiltonian systems in their own right, observing that ‘many strands of classical mathematical physics come together’ (Aref 2007a, § I). A frequently addressed question is whether a particular number of point vortices on some given geometry (for example the sphere) yields integrable dynamics or not. In § 4 of this paper we (re)connect the mathematical theory for integrability of PVD to the long-time behaviour of a continuous, generic incompressible fluid, thereby obtaining an explanation of the three observed regimes. This is briefly how the mechanism works:

1. (i) Smaller vortex formations of the same sign merge to larger formations when their trajectories come close enough (the inverse energy cascade).

2. (ii) The motion of $N$ vortex blobs is accurately described by $N$ point vortices as long as the blobs are well separated (so that no merging occurs). A careful, numerical evaluation of this assumption is given in § 4.1.

3. (iii) If the motion of $N$ vortex blobs is not integrable, then, sooner or later, two vortex blobs of equal sign will reach a point in phase space where they are close enough to merge.

4. (iv) If, however, the motion of the $N$ vortex blobs is integrable (or at least close enough to integrable in the Kolmogorov–Arnold–Moser (KAM) sense, see § 4.2) then the motion remains quasi-periodic with well-separated trajectories and no further merging occurs (integrability acts as a barrier in phase space, preventing further merging of blobs). (From a mathematical viewpoint, the integrability prevents the dynamical system from being ergodic. Ergodicity is assumed in statistical mechanics theories such as MRS.)

To summarize, vortex blobs of equal sign continue to merge until integrability blocks them from doing so. Thus, in order to find the long-time behaviour, one has to find the largest possible number of point vortices for which the dynamics is integrable. Here is the key point: on the non-rotating sphere integrability depends on the total angular momentum. A 4-PVD system on the sphere is integrable if the momentum is zero, but non-integrable if the momentum is non-zero (Sakajo 2007). If momentum is close to zero one still obtains ‘integrable-like’ dynamics since integrable systems are stable in the sense of KAM theory for small perturbations (the small momentum configuration can be viewed as a perturbation of a zero momentum configuration). This explains why 4 vortex blobs is the stable long-time regime for fluid configurations with a small $\unicode[STIX]{x1D6FE}$ parameter. If the momentum in a 4 blob configuration is above the threshold where KAM can be applied, the dynamics is chaotic and sooner or later two of the blobs will merge into a 3 blob configuration. Since 3-PVD systems on the sphere are integrable (regardless of the momentum), this explains the intermediate 3 blob regime. It remains to explain why 2 blobs are sometimes formed. If $\unicode[STIX]{x1D6FE}$ is large enough, there are already two dominant vortex blobs from the start, so the smaller vortex formations are directly merged with these two without passing through the stable 3 vortex blob regime. We thereby have an explanation of the mechanism leading to the three observed regimes.

Conclusions and an outlook to future research are presented in § 5. Although our main focus is on the non-rotating sphere, we have included in appendix A numerical examples of Rossby–Haurwitz waves on a rotating sphere, to illustrate the usability of the new method also in the rotating case (relevant for quasi-geostrophic flows in atmospheric dynamics).

2 Numerical integration algorithm

For spatial discretization we use the system of differential equations developed by Zeitlin (2004), based on the work of Hoppe et al. on the approximation of infinite-dimensional Lie algebras (Bordemann et al. 1991, 1994). The Poisson algebra of smooth functions on the sphere is approximated by the finite-dimensional matrix Lie algebras $\mathfrak{s}\mathfrak{l}(N,\mathbb{C})$, for the Poisson algebra $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$, and $\mathfrak{s}\mathfrak{u}(N)$ for the Poisson algebra $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{R})$. To discretize the equations in time we use the class of isospectral symplectic Runge–Kutta methods developed by Modin & Viviani (2019).

2.1 Spatial discretization via geometric quantization

This section is devoted to the technique used to get a finite dimension analogue of the Euler equations on a sphere. The main theoretical concept behind the approach is the so called $L_{\unicode[STIX]{x1D6FC}}$-approximation.

2.1.1 The $L_{\unicode[STIX]{x1D6FC}}$-approximation

Consider a Lie algebra $(\mathfrak{g},[\cdot ,\cdot ])$ and a family of labelled Lie algebras $(\mathfrak{g}_{\unicode[STIX]{x1D6FC}},[\cdot ,\cdot ]_{\unicode[STIX]{x1D6FC}})_{\unicode[STIX]{x1D6FC}\in \unicode[STIX]{x1D644}}$, where $\unicode[STIX]{x1D6FC}\in I=\mathbb{N}$ or $\mathbb{R}$. Furthermore, assume that, to any element of this family, a distance $d_{\unicode[STIX]{x1D6FC}}$ and a surjective projection map $p_{\unicode[STIX]{x1D6FC}}:\mathfrak{g}\rightarrow \mathfrak{g}_{\unicode[STIX]{x1D6FC}}$ are associated. Then an $L_{\unicode[STIX]{x1D6FC}}$-approximation $(\mathfrak{g}_{\unicode[STIX]{x1D6FC}},[\cdot ,\cdot ]_{\unicode[STIX]{x1D6FC}})_{\unicode[STIX]{x1D6FC}\in \unicode[STIX]{x1D644}}$ of $(\mathfrak{g},[\cdot ,\cdot ])$ should fulfil:

1. (i) if $x,y\in \mathfrak{g}$ and $d_{\unicode[STIX]{x1D6FC}}(p_{\unicode[STIX]{x1D6FC}}(x),p_{\unicode[STIX]{x1D6FC}}(y))\rightarrow 0$ as $\unicode[STIX]{x1D6FC}\rightarrow \infty$, then $x=y$;

2. (ii) for all $x,y\in \mathfrak{g}$ we have $d_{\unicode[STIX]{x1D6FC}}(p_{\unicode[STIX]{x1D6FC}}([x,y]),[p_{\unicode[STIX]{x1D6FC}}(x),p_{\unicode[STIX]{x1D6FC}}(y)]_{\unicode[STIX]{x1D6FC}})\rightarrow 0$ as $\unicode[STIX]{x1D6FC}\rightarrow \infty$;

3. (iii) for $\unicode[STIX]{x1D6FC}\gg 0$ the projections $p_{\unicode[STIX]{x1D6FC}}$ are surjective.

The above definition is given in Bordemann et al. (1994); it is a weak requirement to obtain a limit for a sequence of Lie algebras.

Consider now the smooth complex functions on the sphere with vanishing mean, denoted $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$. This vector space is endowed with a Poisson structure $\{\cdot ,\cdot \}$ given by the skew symmetric bilinear form on $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$

(2.1)$$\begin{eqnarray}\{f,g\}(\boldsymbol{x})=|\boldsymbol{X}_{f}(x),\boldsymbol{X}_{g}(\boldsymbol{x}),\boldsymbol{x}|,\end{eqnarray}$$

where $\boldsymbol{X}_{h}(\boldsymbol{x})=\boldsymbol{x}\times \unicode[STIX]{x1D735}h(\boldsymbol{x})$ is the Hamiltonian vector field associated with the Hamiltonian function $h\in C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$. With this bracket, $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$ becomes an infinite-dimensional Poisson algebra; in particular, it is an infinite-dimensional Lie algebra.

A basis for $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$ is given by the complex spherical harmonics, expressed in the standard azimuthal-inclination coordinates $(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})$ by

(2.2)$$\begin{eqnarray}Y_{lm}(\unicode[STIX]{x1D719},\unicode[STIX]{x1D703})=\sqrt{{\displaystyle \frac{2l+1}{4\unicode[STIX]{x03C0}}}{\displaystyle \frac{(l-m)!}{(l+m)!}}}P_{l}^{m}(\cos \unicode[STIX]{x1D703})\text{e}^{\text{i}m\unicode[STIX]{x1D719}},\quad l\geqslant 1,\,m=-l,\ldots ,l,\end{eqnarray}$$

where $P_{l}^{m}$ are the associated Legendre polynomials (i.e. solutions to the general Legendre equation). Using this basis, an explicit approximating sequence for $C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C})$ was constructed by Hoppe (1982). The sequence is given by the matrix Lie algebras $(\mathfrak{s}\mathfrak{l}(N,\mathbb{C}),[\cdot ,\cdot ]_{N})_{N\in \mathbb{N}}$, where $[\cdot ,\cdot ]_{N}:=N^{3/2}[\cdot ,\cdot ]$ is a rescaling of the matrix commutator $[\cdot ,\cdot ]$. The distances $d_{N}$ are given by suitable matrix norms, and the projections $p_{N}$ are obtained by associating with each spherical harmonic $Y_{lm}$ a matrix $\text{i}\unicode[STIX]{x1D64F}_{lm}^{N}\in \mathfrak{s}\mathfrak{l}(N,\mathbb{C})$ defined by

(2.3)$$\begin{eqnarray}(\unicode[STIX]{x1D64F}_{lm}^{N})_{m_{1}m_{2}}=(-1)^{[(N-1)/2]-m_{1}}\sqrt{2l+1}\left(\begin{array}{@{}ccc@{}}{\displaystyle \frac{N-1}{2}} & l & {\displaystyle \frac{N-1}{2}}\\ -m_{1} & m & m_{2}\end{array}\right),\end{eqnarray}$$

where the bracket denotes the Wigner 3j-symbols. The following $L_{\unicode[STIX]{x1D6FC}}$-convergence result for this approximating sequence have been established:

Theorem 1 (Bordemann et al. (1991, 1994)).

Consider the Poisson algebra $(C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C}),\{\cdot ,\cdot \})$ with Poisson bracket defined by (2.1). Then, for the projections $p_{N}$ and any choice of matrix norms $d_{N}$, $(\mathfrak{s}\mathfrak{l}(N,\mathbb{C}),[\cdot ,\cdot ]_{N})_{N\in \mathbb{N}}$ is an $L_{\unicode[STIX]{x1D6FC}}$-approximation of $(C_{0}^{\infty }(\mathbb{S}^{2},\mathbb{C}),\{\cdot ,\cdot \})$.

2.1.2 The quantized system

We can now derive the spatial discretization of the Euler equations via the $L_{\unicode[STIX]{x1D6FC}}$-approximation in Theorem 1, thereby obtaining a finite-dimensional ‘quantized’ system. We begin without the Coriolis parameter.

For any $N\in \mathbb{N}$ an analogue of the Euler equations (1.2) is the following flow of matrices

(2.4)$$\begin{eqnarray}\dot{\unicode[STIX]{x1D652}}=[\unicode[STIX]{x1D6E5}_{N}^{-1}\unicode[STIX]{x1D652},\unicode[STIX]{x1D652}]_{N},\end{eqnarray}$$

where $\unicode[STIX]{x1D652}\in \mathfrak{s}\mathfrak{l}(N,\mathbb{C})$ and $\unicode[STIX]{x1D6E5}_{N}^{-1}$ is the inverse of the discrete Laplacian, given by the following formula of Hoppe & Yau (1998)

(2.5)$$\begin{eqnarray}\unicode[STIX]{x1D6E5}_{N}=\frac{N^{2}-1}{2}\left([\unicode[STIX]{x1D653}_{3}^{N},[\unicode[STIX]{x1D653}_{3}^{N},\boldsymbol{\cdot }]]-\frac{1}{2}[\unicode[STIX]{x1D653}_{+}^{N},[\unicode[STIX]{x1D653}_{-}^{N},\boldsymbol{\cdot }]]-\frac{1}{2}[\unicode[STIX]{x1D653}_{-}^{N},[\unicode[STIX]{x1D653}_{+}^{N},\boldsymbol{\cdot }]]\right),\end{eqnarray}$$

where $\unicode[STIX]{x1D653}_{\pm }^{N}\propto \unicode[STIX]{x1D64F}_{1\pm 1}^{N}$, $\unicode[STIX]{x1D653}_{3}^{N}\propto \unicode[STIX]{x1D64F}_{10}^{N}$. The crucial property of $\unicode[STIX]{x1D6E5}_{N}^{-1}$ is that $\unicode[STIX]{x1D6E5}_{N}^{-1}\unicode[STIX]{x1D64F}_{lm}^{N}=(-l(l+1))^{-1}\unicode[STIX]{x1D64F}_{lm}^{N}$, for any $l=1,\ldots ,N$, $m=-l,\ldots ,l$. That is, the basis elements $\unicode[STIX]{x1D64F}_{lm}^{N}$ are eigenvectors of the discrete Laplacian $\unicode[STIX]{x1D6E5}_{N}$, which is a direct analogue to the continuous case where the spherical harmonics $Y_{lm}$ are eigenvectors of the Laplace–Beltrami operator $\unicode[STIX]{x1D6E5}$.

Let us again, now explicitly, discuss the connection between the continuous vorticity equation (1.2) and the quantized version (2.4). First, notice that (2.4) is an isospectral flow; it preserves the eigenvalues of $\unicode[STIX]{x1D652}=\unicode[STIX]{x1D652}(t)$. This isospectral property is a direct analogue of preservation of Casimirs in the continuous flow (1.2). Given a continuous vorticity function expanded in the spherical harmonics basis, $\unicode[STIX]{x1D714}=\sum \unicode[STIX]{x1D714}^{lm}Y_{lm}$, the projection operator $p_{N}$ is given by

(2.6)$$\begin{eqnarray}p_{N}(\unicode[STIX]{x1D714})=\mathop{\sum }_{l=1}^{N-1}\mathop{\sum }_{m=-l}^{l}\text{i}\unicode[STIX]{x1D714}^{lm}\unicode[STIX]{x1D64F}_{lm}^{N}.\end{eqnarray}$$

If the continuous vorticity $\unicode[STIX]{x1D714}$ is real valued, then the spherical harmonics coefficients fulfil $\unicode[STIX]{x1D714}^{lm}=(-1)^{m}\unicode[STIX]{x1D714}^{l(-m)}$. The corresponding condition on the matrix $\unicode[STIX]{x1D652}\in \mathfrak{s}\mathfrak{l}(N)$ is $\unicode[STIX]{x1D652}+\unicode[STIX]{x1D652}^{\dagger }=0$, i.e. it belongs to the subalgebra $\mathfrak{s}\mathfrak{u}(N)$ of trace-free skew Hermitian matrices. Thus, we need to restrict the quantized flow (2.4) to $\mathfrak{s}\mathfrak{u}(N)$, which is possible since $\mathfrak{s}\mathfrak{u}(N)$ is a matrix Lie algebra (so it is closed under the matrix commutator $[\cdot ,\cdot ]$) and since the discrete Laplacian $\unicode[STIX]{x1D6E5}_{N}$ restricts to an operator $\mathfrak{s}\mathfrak{u}(N)\rightarrow \mathfrak{s}\mathfrak{u}(N)$ (corresponding to the fact that the continuous Laplace–Beltrami operator $\unicode[STIX]{x1D6E5}$ on $C^{\infty }(\mathbb{S}^{2},\mathbb{C})$ restricts to real functions $C^{\infty }(\mathbb{S}^{2},\mathbb{R})$).

Recall from the introduction that the continuous vorticity equation (1.2) is a Lie–Poisson system with Hamiltonian given by (1.5). Likewise, the quantized equation (2.4) is a Lie–Poisson system on the dual of $\mathfrak{s}\mathfrak{u}(N)$ with Hamiltonian given by

(2.7)$$\begin{eqnarray}H(\unicode[STIX]{x1D652})={\textstyle \frac{1}{2}}\operatorname{Tr}(\unicode[STIX]{x1D6E5}_{N}^{-1}\unicode[STIX]{x1D652}\unicode[STIX]{x1D652}^{\dagger }).\end{eqnarray}$$

The continuous Casimir functions ${\mathcal{C}}_{k}(\unicode[STIX]{x1D714})$ for (1.2) correspond, up to a normalization constant depending on $N$, to the following Casimir functions for (2.4)

(2.8)$$\begin{eqnarray}C_{k}(\unicode[STIX]{x1D652})=\operatorname{Tr}(\unicode[STIX]{x1D652}^{k})\quad \text{for}~k=2,\ldots ,N.\end{eqnarray}$$

As $N\rightarrow \infty$ we have convergence to the corresponding moments ${\mathcal{C}}_{k}(\unicode[STIX]{x1D714})$ of the continuous vorticity (see Rios & Straume 2014, Corollary 8.1.2). We remark that the matrices $\unicode[STIX]{x1D64F}_{lm}^{N}$, with the Frobenius inner product, share the orthogonality properties of $Y_{lm}$, with the $L^{2}(\mathbb{S}^{2},\mathbb{C})$ inner product. Therefore, if the initial vorticity $\unicode[STIX]{x1D714}$ is represented by a finite linear combination of spherical harmonics, then choosing $N$ sufficiently large, the continuous Hamiltonian $H(\unicode[STIX]{x1D714})$ and enstrophy (quadratic Casimir) ${\mathcal{C}}_{2}(\unicode[STIX]{x1D714})$ exactly coincide with the quantized analogues $H(\unicode[STIX]{x1D652})$ and $C_{2}(\unicode[STIX]{x1D652})$.

In the rotating case the quantized system is

(2.9)$$\begin{eqnarray}\dot{\unicode[STIX]{x1D652}}=[\unicode[STIX]{x1D6E5}_{N}^{-1}(\unicode[STIX]{x1D652}-\unicode[STIX]{x1D641}),\unicode[STIX]{x1D652}]_{N},\end{eqnarray}$$

where $\unicode[STIX]{x1D641}=2\unicode[STIX]{x1D6FA}\text{i}\unicode[STIX]{x1D64F}_{10}^{N}$ represents the discrete Coriolis parameter. The Hamiltonian in this case is given by

(2.10)$$\begin{eqnarray}H(\unicode[STIX]{x1D652})={\textstyle \frac{1}{2}}\operatorname{Tr}(\unicode[STIX]{x1D6E5}_{N}^{-1}(\unicode[STIX]{x1D652}-\unicode[STIX]{x1D641})(\unicode[STIX]{x1D652}-\unicode[STIX]{x1D641})^{\dagger }).\end{eqnarray}$$

2.2 Time discretization

To obtain a complete algorithm we also have to discretize time. For this, we use two different schemes. The first is implicit and preserves the Lie–Poisson structure. The second is explicit but does not preserve the Lie–Poisson structure.

2.2.1 Isospectral midpoint method

To take advantage of the quantization of the original equations, it is preferable to solve the quantized system (2.4) in time using a Lie–Poisson integrator, i.e. a time stepping scheme that preserves the Lie–Poisson structure (cf. McLachlan, Modin & Verdier 2014, 2016). This way we obtain exact conservation of the Casimir functions and near conservation of the Hamiltonian (in the sense of backward error analysis of symplectic integrators (cf. Hairer et al. 2006)). Since (2.4) is a Hamiltonian isospectral flow we can apply the Lie–Poisson schemes developed by Modin & Viviani (2019). We use here the second-order isospectral midpoint rule (IsoMP). Given a time step parameter $h>0$ it is given by

(2.11)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D652}_{n}=\left(\unicode[STIX]{x1D644}-\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}\widetilde{\unicode[STIX]{x1D652}}\right)\widetilde{\unicode[STIX]{x1D652}}\left(\unicode[STIX]{x1D644}+\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}\widetilde{\unicode[STIX]{x1D652}}\right),\\ \displaystyle \unicode[STIX]{x1D652}_{n+1}=\left(\unicode[STIX]{x1D644}+\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}\widetilde{\unicode[STIX]{x1D652}}\right)\widetilde{\unicode[STIX]{x1D652}}\left(\unicode[STIX]{x1D644}-\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}\widetilde{\unicode[STIX]{x1D652}}\right),\end{array}\right\}\end{eqnarray}$$

where $\unicode[STIX]{x1D644}$ is the identity matrix. The matrix $\widetilde{\unicode[STIX]{x1D652}}$ is an auxiliary variable implicitly defined (together with $\unicode[STIX]{x1D652}_{n+1}$) by the two equations in (2.11). For further details on the method (2.11) we refer to Viviani (2019).

In presence of the Coriolis parameter $F$ the IsoMP scheme is

(2.12)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \unicode[STIX]{x1D652}_{n}=\left(\unicode[STIX]{x1D644}-\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}(\widetilde{\unicode[STIX]{x1D652}}-\unicode[STIX]{x1D641})\right)\widetilde{\unicode[STIX]{x1D652}}\left(\unicode[STIX]{x1D644}+\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}(\widetilde{\unicode[STIX]{x1D652}}-\unicode[STIX]{x1D641})\right),\\ \displaystyle \unicode[STIX]{x1D652}_{n+1}=\left(\unicode[STIX]{x1D644}+\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}(\widetilde{\unicode[STIX]{x1D652}}-\unicode[STIX]{x1D641})\right)\widetilde{\unicode[STIX]{x1D652}}\left(\unicode[STIX]{x1D644}-\frac{h}{2}\unicode[STIX]{x1D6E5}_{N}^{-1}(\widetilde{\unicode[STIX]{x1D652}}-\unicode[STIX]{x1D641})\right).\end{array}\right\}\end{eqnarray}$$

The IsoMP method (2.11) (and (2.12)) exactly conserves angular momentum and the Casimirs $C_{k}(\unicode[STIX]{x1D652})$, and nearly conserves the Hamiltonian $H(\unicode[STIX]{x1D652})$ (its value oscillates in time without drift).

2.2.2 Heun’s method

As an alternative to the Lie–Poisson preserving time discretization just described, we also consider the explicit Heun method. Explicit methods, such as Heun’s, exhibit linear drift in the first integrals. However, if the linear drift is slow in comparison with the total simulation time, an explicit method might be the most competitive choice since it avoids nonlinear root solving. An efficient implementation of Heun’s method for the quantized system (2.4) is the following:

(2.13)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D646}_{1}=\unicode[STIX]{x1D6E5}_{N}^{-1}\unicode[STIX]{x1D652}_{n}\unicode[STIX]{x1D652}_{n},\\ \displaystyle \widetilde{\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D652}_{n}+h\left(\unicode[STIX]{x1D646}_{1}-\unicode[STIX]{x1D646}_{1}^{\dagger }-\frac{1}{N}\operatorname{Tr}(\unicode[STIX]{x1D646}_{1}-\unicode[STIX]{x1D646}_{1}^{\dagger })\unicode[STIX]{x1D644}\right),\\ \unicode[STIX]{x1D646}_{2}=\unicode[STIX]{x1D646}_{1}+\unicode[STIX]{x1D6E5}_{N}^{-1}\widetilde{\unicode[STIX]{x1D652}}\widetilde{\unicode[STIX]{x1D652}},\\ \displaystyle \unicode[STIX]{x1D652}_{n+1}=\unicode[STIX]{x1D652}_{n}+\frac{h}{2}\left(\unicode[STIX]{x1D646}_{2}-\unicode[STIX]{x1D646}_{2}^{\dagger }-\frac{1}{N}\operatorname{Tr}(\unicode[STIX]{x1D646}_{2}-\unicode[STIX]{x1D646}_{2}^{\dagger })\unicode[STIX]{x1D644}\right).\end{array}\right\}\end{eqnarray}$$

In the presence of the Coriolis parameter $F$ the scheme becomes

(2.14)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D646}_{1}=\unicode[STIX]{x1D6E5}_{N}^{-1}(\unicode[STIX]{x1D652}_{n}-\unicode[STIX]{x1D641})\unicode[STIX]{x1D652}_{n},\\ \displaystyle \widetilde{\unicode[STIX]{x1D652}}=\unicode[STIX]{x1D652}_{n}+h\left(\unicode[STIX]{x1D646}_{1}-\unicode[STIX]{x1D646}_{1}^{\dagger }-\frac{1}{N}\operatorname{Tr}(\unicode[STIX]{x1D646}_{1}-\unicode[STIX]{x1D646}_{1}^{\dagger })\unicode[STIX]{x1D644}\right),\\ \unicode[STIX]{x1D646}_{2}=\unicode[STIX]{x1D646}_{1}+\unicode[STIX]{x1D6E5}_{N}^{-1}(\widetilde{\unicode[STIX]{x1D652}}-\unicode[STIX]{x1D641})\widetilde{\unicode[STIX]{x1D652}},\\ \displaystyle \unicode[STIX]{x1D652}_{n+1}=\unicode[STIX]{x1D652}_{n}+\frac{h}{2}\left(\unicode[STIX]{x1D646}_{2}-\unicode[STIX]{x1D646}_{2}^{\dagger }-\frac{1}{N}\operatorname{Tr}(\unicode[STIX]{x1D646}_{2}-\unicode[STIX]{x1D646}_{2}^{\dagger })\unicode[STIX]{x1D644}\right).\end{array}\right\}\end{eqnarray}$$

2.3 Complexity

At first sight, it looks like the most demanding computational operation is the inversion of the discrete Laplacian $\unicode[STIX]{x1D6E5}_{N}$: it is a linear operator on $\mathfrak{s}\mathfrak{l}(N,\mathbb{C})$ and thus a fourth-order tensor, so dense linear algebra would require $O(N^{4})$ operations. This is clearly not possible, even for moderate values of $N$. However, from the formula (2.5) of Hoppe and Yau one can deduce

(2.15)$$\begin{eqnarray}\displaystyle (\unicode[STIX]{x1D6E5}_{N})_{M_{1}M_{2}}^{M_{1}^{\prime }M_{2}^{\prime }} & = & \displaystyle 2\unicode[STIX]{x1D6FF}_{M_{1}}^{M_{1}^{\prime }}\unicode[STIX]{x1D6FF}_{M_{2}}^{M_{2}^{\prime }}(s(s+1)-M_{1}M_{2})\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FF}_{M_{1}+1}^{M_{1}^{\prime }}\unicode[STIX]{x1D6FF}_{M_{2}+1}^{M_{2}^{\prime }}\sqrt{s(s+1)-M_{1}(M_{1}+1)}\sqrt{s(s+1)-M_{2}(M_{2}+1)}\nonumber\\ \displaystyle & & \displaystyle -\,\unicode[STIX]{x1D6FF}_{M_{1}-1}^{M_{1}^{\prime }}\unicode[STIX]{x1D6FF}_{M_{2}-1}^{M_{2}^{\prime }}\sqrt{s(s+1)-M_{1}(M_{1}-1)}\sqrt{s(s+1)-M_{2}(M_{2}-1)},\qquad\end{eqnarray}$$

for $M_{1},M_{1}^{\prime },M_{2},M_{2}^{\prime }=1,\ldots ,N$ and $s=(N-1)/2$. Notice that the tensor $\unicode[STIX]{x1D6E5}_{N}$ is tridiagonal over the diagonal $M_{1}=M_{1}^{\prime }$ and $M_{2}=M_{2}^{\prime }$, i.e. it is sparse and contains only $O(N^{2})$ non-zero entries; we store $\unicode[STIX]{x1D6E5}_{N}$ as an $N^{2}\times N^{2}$ sparse matrix. Remarkably, this sparse matrix also admits a sparse $\unicode[STIX]{x1D647}\unicode[STIX]{x1D650}$-factorization, i.e. a factorization of upper and lower diagonal matrices $\unicode[STIX]{x1D647}$ and $\unicode[STIX]{x1D650}$ which are also sparse with $O(N^{2})$ non-zero entries. Thus, to compute the inverse $\unicode[STIX]{x1D6E5}_{N}^{-1}\unicode[STIX]{x1D652}$ requires just a single $\unicode[STIX]{x1D647}\unicode[STIX]{x1D650}$-factorization (which is $O(N^{3})$ operations) and thereafter only $O(N^{2})$ operations every time $\unicode[STIX]{x1D6E5}_{N}$ is applied. In essence, since the number of time steps for long-time simulations typically are of the order $O(10^{6})$, this means that inversion of the discrete Laplacian only counts as $O(N^{2})$ operations.

We solve the nonlinear equation (2.11) with Newton iterations. Thus, under the assumption that the average number of iterations per step is independent of $N$, the global complexity of the algorithm per time step is first $O(N^{2})$ (for applying $\unicode[STIX]{x1D6E5}_{N}^{-1}$) and then $O(N^{3})$ (for the two matrix multiplications corresponding to computing the commutator $[\cdot ,\cdot ]$). In summary, this means that the full complexity of the algorithm, per time step, is $O(N^{3})$.

2.4 Time scaling

Recall that the correspondence between the matrix commutator on $\mathfrak{s}\mathfrak{u}(N)$ and the Poisson bracket on $C^{\infty }(\mathbb{S}^{2},\mathbb{R})$ is $N^{3/2}[\cdot ,\cdot ]\approx c\{\cdot ,\cdot \}$ for some constant $c$. The requirement that 1 time unit of the vorticity equation (1.2) corresponds to 1 time unit of the quantized system (2.4) as $N\rightarrow \infty$ implies $c=\sqrt{16\unicode[STIX]{x03C0}}$. In our simulations below we normalize the time scaling of the quantized equations by rescaling the initial conditions by $\Vert \unicode[STIX]{x1D652}_{0}\Vert$ and setting $[\cdot ,\cdot ]_{N}=[\cdot ,\cdot ]$. This way, the non-dimensional time step $h$ corresponds to

(2.16)$$\begin{eqnarray}\unicode[STIX]{x1D6FF}t={\displaystyle \frac{h\sqrt{16\unicode[STIX]{x03C0}}}{N^{3/2}\Vert W_{0}\Vert }}\end{eqnarray}$$

seconds of real time. In all our simulations below we use the non-dimensional time step $h=0.1$. A summary of the complete algorithm is given in Algorithm 1; it is implemented using MATLAB and available online. (The code is available at bitbucket.org/kmodin/euler-sphere-quantization.)

3 Simulation results

3.1 Initial data with zero momentum

We run our method with the same (randomly generated but zero momentum) initial data suggested by DQM, i.e. Dritschel et al. (2015). We use $N=501$, $[\cdot ,\cdot ]_{N}=[\cdot ,\cdot ]$, and a dimensionless time step of $h=0.1$. With these parameters, the simulation time $t_{k}$ at step $k$ in the original units of time is computed by the formula $t_{k}=k/13\,643$, (derived from the formula in § 2.4). We simulate with both the IsoMP and the Heun time integration. For IsoMP, we use Newton-type iterations with a tolerance of $10^{-13}$.

As already discussed in the introduction, the numerical results by DQM show that steady state is not reached, but rather four main vortex formations that move around the sphere, surrounded by smaller-scale vortices. Let us now compare with our results. The vorticity at various output times is displayed, using spherical coordinates, in figure 1 for the two different time integration methods (IsoMP and Heun).

Figure 1. Simulation with two different time integration methods: IsoMP (2.11) and Heun (2.13). Vorticity $\unicode[STIX]{x1D714}(x,t)$ clockwise from the top left at $t=0~\text{s}$, 4 s, 40 s, 400 s, for the initial data in DQM. The horizontal axis is the azimuth $\unicode[STIX]{x1D711}\in [0,2\unicode[STIX]{x03C0}]$ and the vertical axis is minus the inclination $\unicode[STIX]{x1D703}\in [0,\unicode[STIX]{x03C0}]$. The results are visually indistinguishable up to $t=40$. At $t=400$ there are some differences in the positions of the vortex blobs. See also Movies 1, 2 and 3 of the supplementary material.

At time $t=4$ our simulations and those in DQM give visually indistinguishable results. At the early–intermediate vorticity, at time $t=40$, there is already a clear visible difference to DQM. However, there is no visible difference between our two numerical time-integration schemes. This indicates that, for time step lengths in the selected range, the choice of discretization in space, rather than time, dominates the numerical errors.

At $t=400~\text{s}$ all simulations show the same qualitative feature: four large vortices moving about in the domain. The exact positions of the vortices are different between all the simulations (also between IsoMP and Heun). There are two positive and two negative vortex blobs. The exact strengths vary slightly between the blobs (see § 4 for further discussion about the vortex strengths).

When we run the simulation, either IsoMP or Heun, for long times a clear pattern emerges: the 4 vortex formations are moving quasi-periodically. The initial vortex mixing phase, up until the four vortex blobs have been formed at approximately $t=200$, is captured in Movie 1 of the supplementary material, available at https://doi.org/10.1017/jfm.2019.944. (All the movies are also available at bitbucket.org/kmodin/euler-sphere-quantization.) The fast-forward Movie 2 of the whole simulation shows a short emerging phase of vortex mixing followed by a stable but unsteady large-scale quasi-periodic interaction of the four vortices. In § 4 we discuss in detail the relation to stability of quasi-periodic point vortex solutions. Movie 3 shows a simulation with the same initial conditions, but at the higher spatial resolution $N=1001$. The qualitative behaviour is the same, with four vortex blobs forming and then circulating about each other in a quasi-periodic fashion. However, the distribution of vortex formation is different in the high resolution simulation, with the positive instead of the negative blobs closer to the poles.

Let us continue the discussion here with the conservation properties of our method. Figure 2 shows the variation of the energy and enstrophy during the simulation. For IsoMP, the energy is nearly conserved by a factor $10^{-6}$ with no sign of drift, whereas the enstrophy has the same variation as the Newton tolerance we have used, $10^{-13}$. For Heun, we see that, albeit energy and enstrophy have linear drifts from their original values, the variation is quite small and, in particular, the energy changes less than with IsoMP. The negligible drift of energy and enstrophy is likely the reason why Heun perform so well. We stress, however, that there is a drift, so at some point the numerics will break down, whereas with IsoMP such a breakdown will not occur since symplecticity is preserved.

Figure 2. Hamiltonian variation $|H-H_{0}|$ (a,c) and enstrophy variation $|E-E_{0}|$ (b,d) with the IsoMP (a,b) and Heun (c,d) time integrators.

The difference between IsoMP and Heun is more pronounced for the higher-order Casimir functions of (2.4). In fact, computing the maximal absolute variation of the eigenvalues of $W$, after $5\times 10^{6}$ time steps, we get with IsoMP a value of the order $10^{-12}$, whereas with Heun a value of the order 1. Even considering only the third and fourth momenta of the vorticity, the Heun scheme has an absolute variation, after $5\times 10^{6}$ time steps, of the order $10^{-3}$.

In addition to integral invariants, such as energy and enstrophy, the continuous vorticity equation (1.2) also conserves pointwise measures, such as the maximum vorticity

(3.1)$$\begin{eqnarray}\Vert \unicode[STIX]{x1D714}\Vert _{\infty }:=\sup _{x\in \mathbb{S}^{2}}|\unicode[STIX]{x1D714}(x)|.\end{eqnarray}$$

Formally, the conservation of $\Vert \unicode[STIX]{x1D714}\Vert _{\infty }$ follows from conservation of the Casimir functions ${\mathcal{C}}_{k}(\unicode[STIX]{x1D714})$ as $k\rightarrow \infty$. Indeed, since the corresponding Casimir functions $C_{k}(\unicode[STIX]{x1D652})$ of the quantized system approximate ${\mathcal{C}}_{k}(\unicode[STIX]{x1D714})$ one can deduce (formally) that $\Vert \unicode[STIX]{x1D714}\Vert _{\infty }$ is nearly conserved without any drift (just like the energy). In fact, this result follows rigorously from a theorem by Bordemann et al. (1994, Theorem 4.1), who proved that there is a constant $c\geqslant 0$, independent of $N$, such that

(3.2)$$\begin{eqnarray}\Vert \unicode[STIX]{x1D652}\Vert \leqslant \Vert \unicode[STIX]{x1D714}\Vert _{\infty }\leqslant \Vert \unicode[STIX]{x1D652}\Vert +\frac{c}{N},\end{eqnarray}$$

where $\Vert \unicode[STIX]{x1D652}\Vert$ is the matrix (operator) norm of $\unicode[STIX]{x1D652}\in \mathfrak{s}\mathfrak{u}(N)$ and $\unicode[STIX]{x1D714}$ is the vorticity function corresponding to $W$. Since $\Vert \unicode[STIX]{x1D652}\Vert$ is the largest eigenvalue (in magnitude) of $\unicode[STIX]{x1D652}$, and since all the eigenvalues are conserved by the quantized flow (the isospectral property), we get that $\Vert \unicode[STIX]{x1D714}\Vert _{\infty }$ is nearly conserved in the quantized system (i.e. it is an adiabatic invariant for the quantized flow).

To measure the unsteadiness in the simulated flow we look at the relation between the vorticity $\unicode[STIX]{x1D714}$ and the streamfunction $\unicode[STIX]{x1D713}$ at $t=400$. The MRS theory predicts a steady flow determined by a functional relation $\unicode[STIX]{x1D714}=F(\unicode[STIX]{x1D713})$ between the vorticity and streamfunction. Figure 3 contains a scatter plot of $\unicode[STIX]{x1D713}$ and $\unicode[STIX]{x1D714}$ for both IsoMP (2.11) and Heun (2.13). We notice that the shape of the resulting diagrams has branches, similar to those in DQM, indicating unsteadiness. Our branches are more diffuse than those in DQM since no artificial dissipation is added in our model. We also see a slight difference between IsoMP and Heun: the one obtained with IsoMP has more defined branches.

Figure 3. Scatter plots of vorticity $\unicode[STIX]{x1D714}$ versus streamfunction $\unicode[STIX]{x1D713}$ at $t=400$. (a) Using the IsoMP time integrator (2.11), and (b) using the Heun time integrator (2.13).

3.2 Generic initial data

In this section we present the results obtained with our numerical scheme on randomly generated initial conditions. We show that the generic behaviour for long times described by Dritschel et al. (2015) it is not attained for non-zero angular momentum of the fluid. In our simulations we use $N=501$, $[\cdot ,\cdot ]_{N}=[\cdot ,\cdot ]$, and a dimensionless time step of $h=0.1$. Again, with these parameters the simulation time at step $k$ in the original units of time is computed by the formula $t=k/13\,643$. We simulate with the Heun time integration as it is faster; for time evolutions as long as 400 real seconds the decay in enstrophy is negligible (see figure 2).

The generic random initial vorticity is obtained as follows. Consider the expansion of the vorticity function in terms of spherical harmonics

(3.3)$$\begin{eqnarray}\unicode[STIX]{x1D714}(x)=\mathop{\sum }_{l=1}^{\infty }\mathop{\sum }_{m=-l}^{l}\unicode[STIX]{x1D714}^{lm}Y_{lm}(x).\end{eqnarray}$$

Then, $\unicode[STIX]{x1D714}\in L^{2}(\mathbb{S}^{2})$ means that $\sum _{l=1}^{\infty }\sum _{m=-l}^{l}|\unicode[STIX]{x1D714}^{lm}|^{2}<\infty$. We set the level of truncation $l_{max}=N-1=500$ and we generate the coefficients as random variables such that $\unicode[STIX]{x1D714}^{lm}l^{1+\unicode[STIX]{x1D700}}\sim {\mathcal{N}}(0,1)$, where ${\mathcal{N}}(0,1)$ is the normal Gaussian distribution and $\unicode[STIX]{x1D700}=10^{-3}$. We stress that $L^{2}(\mathbb{S}^{2})$ as the space for initial conditions is a natural choice in terms of Fourier analysis. Generating random initial conditions as just described corresponds mathematically to drawing samples from the isotropic Gaussian random field on $L^{2}(\mathbb{S}^{2})$ as described by Lang & Schwab (2015).

Figure 4. Pairs of initial (upper) and final (lower) vorticities for the 16 generic simulations with $L^{2}(\mathbb{S}^{2})$ random initial data. The numbers labelling the simulations correspond to those in figure 9.

In this setting, we run 16 simulations on a cluster for long times. The vorticity for the simulations are given in figure 4 at time $t=0$ and $t=400$ real seconds. From the initially chaotic vorticity we see at time $t=400$ three qualitatively scenarios, with either 2, 3 or 4 persistent coherent vortices. In § 4 we explain this phenomenon in terms of integrability of point vortex dynamics and KAM theory. Movies 5, 6 and 7 of the supplementary material show the complete evolution of simulations 1 (giving 2 blobs), 4 (giving 3 blobs) and 7 (giving 4 blobs).

4 Relation to point vortex dynamics

We now explain the connection between the long-time behaviour of the Euler equations (1.1) on a non-rotating sphere and integrability theory of point vortex dynamics. Recall from the introduction that our theory is based on the following two assumptions:

1. (i) The inverse energy cascade operates in such a way that smaller vortex formations of the same sign merge into larger ones when they get close enough.

2. (ii) PVD describes the motion of vortex blobs well, as long as the blobs are well separated so that no merging occurs.

Based on the simplest, zero momentum case, we first give a detailed numerical verification of the second of these assumptions. We then give the connection to integrability. After that, we address the generic case of non-zero momentum and we show how our simulation results, with the three observed regimes, is a consequence of our theory.

4.1 Zero momentum case

In this section we give a detailed study of the relation of our simulation results to the dynamics of four point vortices on the sphere, following up the brief study in Dritschel et al. (2015). For a detailed treatment of point vortex dynamics, we refer to the monograph of Newton (2016) or the survey paper by Aref (2007b).

Already Helmholtz (1858) knew that the incompressible Euler equations admit solutions with vorticity supported on single points. Such solutions also appear for the vorticity equations (1.2) on a sphere in the non-rotational case (Bogomolov 1977). That is, vorticity is a finite sum of $n$ Dirac delta distributions

(4.1)$$\begin{eqnarray}\unicode[STIX]{x1D714}=\mathop{\sum }_{i=1}^{n}\unicode[STIX]{x1D6E4}_{i}\unicode[STIX]{x1D6FF}_{\boldsymbol{x}_{i}},\end{eqnarray}$$

where $\unicode[STIX]{x1D6E4}_{i}\in \mathbb{R}$ are the strengths and $\boldsymbol{x}_{i}\in \mathbb{S}^{2}$ are the positions of the point vortices. The solutions evolve according to an ordinary differential equation known as the point vortex equation

(4.2)$$\begin{eqnarray}\dot{\boldsymbol{x}_{i}}=\frac{1}{4\unicode[STIX]{x03C0}}\mathop{\sum }_{i\neq j}\unicode[STIX]{x1D6E4}_{j}{\displaystyle \frac{\boldsymbol{x}_{j}\times \boldsymbol{x}_{i}}{1-\boldsymbol{x}_{i}\boldsymbol{\cdot }\boldsymbol{x}_{j}}},\end{eqnarray}$$

for $i=1,\ldots ,n$. Notice that multiplying all $\unicode[STIX]{x1D6E4}_{i}$ by a factor does not change the phase space trajectories (only their speed), so only relative strengths are of importance to us. Our aim is to extract the positions and relative strengths of the vortex blobs in the DQM simulation, to compare their motion with the corresponding system of $n=4$ point vortices.

Figure 5. Behaviour of the vortex blobs for the same initial data as in figure 1. The quasi-periodic motion of the blobs is tracked in the scatter plots.

To extract the trajectories on the sphere of the 4 vortex blobs in the simulation from the previous § 3.1, we use a tracking algorithm based on Python/SciPy. The result is given in figure 5.

Now, for the Euler equation on a non-rotating sphere, the total angular momentum $\boldsymbol{L}$ is conserved

(4.3)$$\begin{eqnarray}\frac{\text{d}}{\text{d}t}\boldsymbol{L}=\frac{\text{d}}{\text{d}t}\int \unicode[STIX]{x1D714}\boldsymbol{n}=0.\end{eqnarray}$$

The DQM simulation is set up with vanishing total angular momentum, $\boldsymbol{L}=0$. Thus, the point vortex solutions should fulfil

(4.4)$$\begin{eqnarray}\int \mathop{\sum }_{i=1}^{n}\unicode[STIX]{x1D6E4}_{i}\unicode[STIX]{x1D6FF}_{\boldsymbol{x}_{i}}\boldsymbol{n}=\mathop{\sum }_{i=1}^{n}\unicode[STIX]{x1D6E4}_{i}\boldsymbol{x}_{i}=0.\end{eqnarray}$$

If we set $\unicode[STIX]{x1D6E4}_{1}=1$ (since we are only looking for relative strengths), then, for generic positions $\boldsymbol{x}_{i}$, this yields a linear set of equations from which $\unicode[STIX]{x1D6E4}_{2},\unicode[STIX]{x1D6E4}_{3},\unicode[STIX]{x1D6E4}_{4}$ can be determined from the positions alone. The computed relative strengths thereby obtained correspond well with those obtained by numerical integration over circular domains covering the blobs. In summary, we have the following extracted positions (expressed in inclination $\unicode[STIX]{x1D711}$ and azimuth $\unicode[STIX]{x1D703}$) and corresponding computed relative strengths

(4.5)$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\unicode[STIX]{x1D711}=[2.3218,-0.9638,-2.5283,0.8511],\\ \unicode[STIX]{x1D703}=[1.3017,1.8837,1.577,1.5896],\\ \unicode[STIX]{x1D6E4}=[1,0.9002,-0.5436,-0.4178].\end{array}\right\}\end{eqnarray}$$

To obtain the absolute strengths, i.e. to determine the scaling, one might use the total energy integral, noting that the point vortex Hamiltonian is quadratic in the strengths.

Remark 1. The fact that the relative strengths of the point vortices are uniquely determined (in the generic case) by the positions given vanishing total angular momentum is interesting. It shows that there is a connection between the strengths and the positions. One may ask to what extent the strengths determine the positions for zero momentum configurations. That is, what is the dimension of the manifold of four point vortices with vanishing total angular momentum. Heuristically, just counting constraints, one gets the dimension of all possible four point vortex configurations, 8, minus the dimension of the 3 angular momentum constraints, which gives 5 dimensions. To investigate this question in detail, one can use symplectic reduction theory (cf. Marsden & Ratiu 1998).

We run the point vortex dynamics simulation with data from (4.5) using the symplectic Lie–Poisson integrator by McLachlan et al. (2016). Let us now compare this point vortex simulation with the tracked blob motion in figure 5.

In the chosen spherical coordinates, the motion of the tracked blobs in figure 5 looks complicated, but, in fact, when plotted on the sphere, one can see that it is almost a steady rigid rotation about a fixed axis. By least squares we find the best approximating rotation axis, and we use new spherical coordinates with the new rotation axis as the north pole. The resulting trajectories, of both the tracked blobs and the computed four point vortex dynamics with data (4.5), are given in figure 6. We see that the motion between the point vortices and the tracked blobs are in good agreement, as also reported by Dritschel et al. (2015).

Figure 6. Motion of point vortices (a,c) versus tracking of vortex blobs in the simulation in figure 5 (b) and in the simulation at low resolution ($N=51$), as described in (4.6) (d).

Looking at the almost pure rotational trajectories in figure 6, it is natural to ask if there is an underlying relative equilibrium, i.e. a close-by solution given by a simultaneous, steady rotation of all the four point vortices. The answer is no: that relative equilibrium is in fact a static equilibrium, because the total angular momentum is zero. Thus, the ‘wobbling’ in figure 6 is necessary for unsteadiness to occur. Continuing this train of thought, we may look for static, non-stable equilibria for zero momentum four point vortex dynamics with arbitrary strengths. Based on symmetry considerations, a general study of equilibria for point vortex dynamics on the sphere is carried out by Laurent-Polz, Montaldi & Roberts (2011). For general strengths there does not seem to exist equilibria, but for pairs of two equal positive and two equal negative strengths there are, given by staggered rings (see Laurent-Polz et al. 2011, § 8). Since the computed strengths (4.5) almost come in such pairs, and since at any instance in time the configuration of the vortex blobs in figure 5 is almost given by such staggered rings, we are, in this sense, always close to equilibria, but they are unstable. That the strength of the vortex blobs almost comes in pairs is likely not a coincidence.

The simulation in § 3.1 generating the blob formation is long enough to cover about 3 revolutions of the blobs about each other. Although this is considered a ‘long-time’ simulation of the Euler equations, it is not very long if one wants to study the stability of the quasi-periodic trajectories. If we run the point vortex simulation for approximately 30 revolutions, we see in figure 6(c) that the trajectories appear to keep on wobbling about the zonal lines. But even 30 revolutions is not much. With a much longer simulation of about 2800 revolutions, we see, as plotted in figure 7, a different pattern emerge: the positions of the vortices are spreading out by a very slow precession. These numerical experiments indicate that the dynamics of the four point vortices restricted to the submanifold of vanishing total angular momentum is integrable, or at least quasi-integrable in the KAM sense. The frequencies would then be the oscillations within each revolution (highest), the rotation (intermediate) and one or two much lower frequencies for the precession. In general, the dynamics of four point vortices is not integrable. However, the submanifold of vanishing total angular momentum is special, as it is the only submanifold of fixed angular momentum that is invariant under arbitrary rotations (if you rotate a configuration of zero momentum, it still has zero momentum). Indeed, Sakajo (2007) showed that the dynamics of four point vortices with zero angular momentum is integrable. As a theoretical approach aiming to prove integrability, one could also proceed by zero momentum Hamiltonian reduction (cf. Marsden et al. 2007). Roughly, it goes as follows. The Lie group $\text{SO}(3)$ of rotations acts on the configuration space $(\mathbb{S}^{2})^{4}$ of point vortices. The corresponding Nöther integrals are the total angular momentum. Since the area form on $\mathbb{S}^{2}$ is preserved by rotations, the action is symplectic. By Poisson reduction we thereby obtain a new Hamiltonian system on the Poisson manifold $(\mathbb{S}^{2})^{4}/\text{SO}(3)$ of dimension 5. Now, every Poisson manifold is foliated in symplectic leaves. In particular, we have the special zero momentum leaf, given by restriction to the zero set of the total angular momentum. We thereby obtain a Hamiltonian system on the symplectic manifold given by the zero momentum leaf of dimension 2, which is always integrable.

Figure 7. Long-time behaviour (about 2800 revolutions) of the four point vortex simulation. The axis of rotation is slowly drifting, yielding a precession of the trajectories. The indicated lines show the motion in the last simulated period. These numerical results verify that the motion is integrable, i.e. quasi-integrable, with three or four frequencies: the wobbling during each rotation (highest frequency), the rotation about the axis (intermediate frequency) and the precession of the rotation axis (lowest frequency/frequencies).

Our findings in this section show that the initial conditions used by DQM, although random in the higher-order spherical harmonics, is special since it has zero angular momentum. That is, one cannot expect the long-time behaviour obtained with the DQM initial conditions to be generic for initial conditions with non-zero angular momentum. Indeed, if four vortex blobs in a non-zero momentum configuration are formed, there might be further mixing, since their motion most likely will be chaotic. For zero momentum, however, the quasi-periodic behaviour acts as a barrier, preventing further mixing. We thus predict that for vanishing angular momentum, quasi-periodic asymptotics is the generic behaviour. To investigate this question in detail is yet another future topic.

We now want to illustrate how the quasi-periodic motion of the blobs can be obtained even for a very coarse spatial discretization $N=51$. The initial vorticity here is:

(4.6)$$\begin{eqnarray}\unicode[STIX]{x1D714}_{0}(\boldsymbol{x})=\mathop{\sum }_{i=1}^{4}\unicode[STIX]{x1D6E4}_{i}\exp (-20|\boldsymbol{x}-\boldsymbol{x}_{i}(\unicode[STIX]{x1D711}_{i},\unicode[STIX]{x1D703}_{i})|^{2})+C(\boldsymbol{x}),\end{eqnarray}$$

for $\unicode[STIX]{x1D6E4},\unicode[STIX]{x1D711},\unicode[STIX]{x1D703}$ as in (4.5), and $C(\boldsymbol{x})$ such that $\unicode[STIX]{x1D714}_{0}$ integrate to zero and has momentum $\boldsymbol{L}=0$. The result is given in figure 8 and Movie 4, and the resulting vortex blob motion tracking (in adapted spherical coordinates) is given in figure 6(d). We obtain good agreement with both the point vortex simulation and the full high resolution simulation with $N=501$. That Gaussian vortex blob simulations can be carried with small discretization parameters $N$ is important, because it opens up for much longer simulations studying the stability of quadruple vortex blob formations. We anticipate that the slow precession seen in point vortex dynamics also happens in vortex blob dynamics.

Figure 8. Motion of initially Gaussian blobs discretized with low spatial discretization ($N=51$). The shape, strengths, initial positions of the blobs are given by (4.6). The tracked motion of the blobs is in good agreement with the $N=501$ simulation and the point vortex dynamics (compare with figure 5).

4.2 Generic case

The ideas presented in the previous paragraph can be extended to the non-zero momentum vorticity. Our simulations in § 3.2 suggest that the four blob formation in Dritschel et al. (2015) is specific for initial conditions with small angular momentum. In fact, as already mentioned in the introduction, there is a correlation between the first integral $\unicode[STIX]{x1D6FE}:=\Vert \boldsymbol{L}\Vert /(R\sqrt{{\mathcal{C}}_{2}})$ and the number of coherent vortices that persist in the final state (see figure 9). As can be seen in the simulations, there exists a finite range for $\unicode[STIX]{x1D6FE}$ (approximately $0.15\lesssim \unicode[STIX]{x1D6FE}\lesssim 0.4$) for which the mixing of vortex blobs is continuous until the dynamics reaches a quasi-periodic motion of three blobs and no more mixing occurs after that. Above this range, the momentum prevails on the other modes, allowing only the persistence of two large vortices. To the best of our knowledge this phenomenon has not been previously described.

Figure 9. Values of $\unicode[STIX]{x1D6FE}=\Vert \boldsymbol{L}\Vert /(R\sqrt{{\mathcal{C}}_{2}})$ for the simulations of § 3.2. The grey scale correspond to the number of vortex blobs observed in the final state: 2, 3 or 4. Notice that the value of $\unicode[STIX]{x1D6FE}$ largely determines the number of vortex blobs in the final state: 4 when $\unicode[STIX]{x1D6FE}\lesssim 0.15$, 3 when $0.15\lesssim \unicode[STIX]{x1D6FE}\lesssim 0.4$ and 2 when $\unicode[STIX]{x1D6FE}\gtrsim 0.4$.

Based on the two assumptions presented at the beginning of this section, an explanation for the observed phenomenon is offered through integrability properties of point vortex dynamics as already laid out in the introduction. Indeed, it is known that for non-zero momentum the three point vortex dynamics is integrable, whereas it is not integrable in general for four point vortices (Sakajo 2007). In § 3.2 our numerical simulations show that when the angular momentum $\boldsymbol{L}$ is non-zero there may occur further mixing from the four vortices found by Dritschel et al. (2015), leading to a final state of three or two vortices. This can be understood in terms of perturbation of an integrable configuration of point vortices. In fact, starting from the zero momentum four blobs, one can understand the modification of momentum $\boldsymbol{L}$ to a non-zero value as a perturbation of the zero-level set. As noticed in § 3.2, up to the critical value of $\unicode[STIX]{x1D6FE}\sim 0.15$, the four point vortex dynamics persists and is quasi-periodic. The reason for such a situation is that the momentum $\boldsymbol{L}$ is a small perturbation, in the sense of the KAM theory (cf. Sevryuk 1995), of an integrable system of point vortices, and small perturbations do not destroy the invariant tori, so the quasi-periodicity is still intact, acting as a barrier for further mixing. However, when $\unicode[STIX]{x1D6FE}\gtrsim 0.15$, the perturbation from zero momentum is large enough to destroy the four point vortex integrable state, leading to chaotic trajectories of the blobs and therefore further mixing up to the next integrable configuration of three point vortices. Eventually, increasing the magnitude of the momentum $\boldsymbol{L}$ over a certain threshold ($\unicode[STIX]{x1D6FE}\gtrsim 0.4$), the final state of the vorticity can be described by two antipodal point vortices only, aligned in the direction of the momentum $\boldsymbol{L}$. These two vortices are now so large from the start that they tend to directly swallow the smaller vortices without passing through the quasi-periodic three blob formation (although, if one looks carefully in simulations 1 and 8 corresponding to the higher values of $\unicode[STIX]{x1D6FE}$, one can trace a small third vortex blob which does not affect the dynamics).

Remark 2. We point out that the relation of the final state of the total vorticity to the point vortex dynamics strongly depends on the manifold where the equations take place. On a torus, in fact, the total angular momentum does not play any role for the integrability of the point vortex dynamics. Instead, the total circulation of the point vortices (i.e. the sum of the vortices’ strengths) is determinant. With zero circulation, three point vortex dynamics on the torus is integrable, whereas for non-zero circulation only two point vortex dynamics is integrable. This explains why the latter configuration of two large blobs has been extensively observed (see for example Qi & Marston 2014), whereas no three large blobs on a torus appear in the simulations. Indeed, prescribing a final state of zero circulation of three blobs (notice, not only a zero circulation of the total vorticity since one has to subtract the constant background circulation) is not possible, unlike prescribing zero momentum. Hence, on a torus, our theory predicts that the generic behaviour is the formation of two steady point vortices, as also observed numerically. We hypothesize that our theory, connecting long-time behaviour with integrability of point vortex dynamics, is valid for the Euler equations on any 2-D surface. To investigate, and possibly verify, this claim in full detail is a future research topic.

5 Conclusions and outlook

We have developed a new numerical algorithm for the Euler equations on a sphere that preserves, up to machine precision, the Casimir functions of (2.4)–(2.9) and nearly conserves the Hamiltonian (see § 2). The spatial discretization is based the work of Hoppe (1982) on geometric quantization of the infinite-dimensional Poisson algebra of smooth functions on the sphere by matrix Lie algebras $\mathfrak{s}\mathfrak{u}(N)$ for an increasing $N$ (corresponding to the spatial discretization parameter), and the suggestion of Zeitlin (2004) to use this quantization for Euler fluids. The resulting finite-dimensional dynamical system on $\mathfrak{s}\mathfrak{u}(N)$ is isospectral, corresponding to conservation of Casimir functions, and preserves a Lie–Poisson structure, corresponding to the Hamiltonian structure of ideal fluids.

On the one hand, long-time simulations on a non-rotating sphere for zero momentum initial vorticity confirm the results in Dritschel et al. (2015) of a quadruple of coherent vortex formations, but now without introducing artificial hyperviscosity into the equations (see § 3.1). On the other hand, for non-zero momentum vorticity, our results show that the generic behaviour suggested in Dritschel et al. (2015) is incorrect: the situation is more complicated yielding either 2, 3 or 4 coherent vortex formations.

In § 4, comparing the motion of the obtained vortex blob formations with point vortex dynamics, we presented a theoretical explanation describing the mechanism for the asymptotic behaviour of the solutions to the Euler equations: the inverse energy cascade continues until two, three or four vortex blobs have been formed, with the number of vortices correlated to the ratio between the magnitude of the momentum and the square root of the enstrophy. After that, the vortex blob formation is blocked from further mixing by the quasi-periodic motion imposed by the integrability of the point vortex dynamics. This way, we establish integrability theory of point vortex dynamics together with KAM perturbation theory as the fundamental theory underlying the formation of unsteady but quasi-periodic coherent vortex formations.

As an outlook, the connections to integrability theory could be studied in much more detail, for example, the relation between the regularity of the generic initial data and the qualitative properties of the final state of the system, e.g. the size of the vortex blobs. Perhaps more pressing is to get better statistics for the correlation between $\unicode[STIX]{x1D6FE}$ and the long-time behaviour: instead of just 16 simulations, we aim to run 512 or more simulations to collect statistics from. One could also try with initial vorticity from the more regular Gaussian random field on the Sobolev space $H^{s}(\mathbb{S}^{2})$. Another aspect to investigate is the long-time behaviour of the vorticity on a rotating sphere. In appendix A, we present some numerical results indicating that quasi-periodic behaviour can also be reached, but is now more complicated than what can be achieved by point vortex dynamics. One could also look deeper at the mechanism behind the inverse energy cascade in the quantized equations. For example, the standard Poisson bracket between two spherical harmonics feeds into harmonics with larger wavenumbers $l$. In the quantized bracket, however, high wavenumber harmonics are fed to lower wave numbers. This might explain why the inverse energy Cascade works well despite spatial truncation. Another benefit of the quantized fluid model is that it is possible to introduce viscosity while still preserving all of the Casimirs. Indeed, one can add a gradient term of (some approximation of) the entropy functional in such a way that isospectrality is preserved; an example is the Brockett flow (Brockett 1991) which is known to correspond to a gradient flow of entropy on the space of multivariate Gaussian probability distributions (Modin 2017). If viscosity is added to the quantized vorticity equations in such a way, the resulting model can be seen as a mix of computational ideal fluid dynamics (corresponding to the conservative part of the dynamics) and the MRS statistical mechanics model (corresponding to the pure spectral preserving entropy maximizing gradient flow).