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. (Reference Bordemann, Meinrenken and Schlichenmaier1994); 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 (Reference Hoppe1982). 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. (Reference Bordemann, Hoppe, Schaller and Schlichenmaier1991, Reference Bordemann, Meinrenken and Schlichenmaier1994)).

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 (Reference Hoppe and Yau1998)

(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 Reference Rios and Straume2014, 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 Reference McLachlan, Modin and Verdier2014, Reference McLachlan, Modin and Verdier2016). 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. Reference Hairer, Lubich and Wanner2006)). Since (2.4) is a Hamiltonian isospectral flow we can apply the Lie–Poisson schemes developed by Modin & Viviani (Reference Modin and Viviani2019). 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 (Reference Viviani2019).

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.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. (Reference Dritschel, Qi and Marston2015). 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. (Reference Bordemann, Meinrenken and Schlichenmaier1994, 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. (Reference Dritschel, Qi and Marston2015) 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 (Reference Lang and Schwab2015).

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.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. (Reference Dritschel, Qi and Marston2015). For a detailed treatment of point vortex dynamics, we refer to the monograph of Newton (Reference Newton2016) or the survey paper by Aref (Reference Aref2007b).

Already Helmholtz (Reference Helmholtz1858) 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 Reference Bogomolov1977). 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.

We run the point vortex dynamics simulation with data from (4.5) using the symplectic Lie–Poisson integrator by McLachlan et al. (Reference McLachlan, Modin and Verdier2016). 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. (Reference Dritschel, Qi and Marston2015).

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 (Reference Laurent-Polz, Montaldi and Roberts2011). 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. Reference Laurent-Polz, Montaldi and Roberts2011, § 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 (Reference Sakajo2007) 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. Reference Marsden, Misiołek, Ortega, Perlmutter and Ratiu2007). 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. (Reference Dritschel, Qi and Marston2015) 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 Reference Sakajo2007). 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. (Reference Dritschel, Qi and Marston2015), 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 Reference Sevryuk1995), 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).