3.1. Forced simulations

For a significant duration, the vorticity fields display only random fluctuation at the size of the forcing scale. This is seen in figure 1, (a, run 1) and (c, run 2), both from $ t = 10$ , where time is measured in units $ \epsilon ^{-1/3} a^{2/3}$ . At $ t \sim 100$ , four coherent vortices emerge from the random field in both runs and keep growing in magnitude, but not in width. The vortices arrange themselves on a great circle with $ 90^\circ$ separation and like signed vortices at opposite poles. In figure 1, (b, run 1) and (d, run 2), we see the four vortices in final states, where the vortex quartet is shown rotated to the equatorial plane. The great circle defining their configuration undergoes a random solid body rotation around the origin, with a characteristic velocity that is far slower than the typical fluid velocities. It may seem paradoxical that the circle moves in this way, since the mean angular momentum of the flow is zero, which is ascertained by the numerical method. Noting that the fluid velocity at the points defining the centres of the vortices is always zero, makes it less paradoxical. The circle is not a material line moving with the fluid, but moves as a solid body through the fluid which has no motion as a solid body. Animations of the vorticity fields in stereographic projections are available at Run1, Run2.

Figure 1. Vorticity fields. Run 1: (a) t = 10, (b) t = 1256; run 2: (c) t = 10, (d) t = 1067. In (b) and (d) the fields have been rotated so that the vortices lie on the equator.

Figure 2. Mean enstrophy, $\Omega$ , normalised by $l_f(1+l_f)$ versus time; (a) run 1, (b) run 2. Evolution of compensated energy spectra; (c) run 1, (d) run 2. The colour of each spectrum is associated with the time where it is computed, as indicated in (a) and (b). The dashed lines in (c–d) indicate the Kolmogorov constant $ C = 6$ .

Figure 3. Azimuthal enstrophy spectra normalised by mean enstrophy, $\Omega$ , for different $ l$ from run 2 at t = 1067. Panels show (a) $ l=2,4,6,8,10$ , (b) $ l = 20,30,40$ .

In figure 2(a,b), we see the time evolution of mean enstrophy. The cascades are set up at $ t \sim l_{f}^{-2/3}$ , after which enstrophy is quasistationary and then enters into a stage of slow growth. Also in this stage, the enstrophy growth is less than 1 % of the enstrophy injection, implying that virtually all injected enstrophy goes into a forward cascade. After the set up of the cascades, mean energy is growing at a constant rate which is close to $ 0.5$ in both simulations, implying that approximately half of the injected energy is transferred to large wavenumbers where it is dissipated. In figure 2(c,d), we see the evolution of the compensated energy spectrum, $ E(l) |\varPi _E |^{-2/3} l^{5/3}$ , where $ |\varPi _E |$ is taken as the energy growth. In the early stage, the spectra are similar to Fourier spectra from planar simulations (Smith & Yakhot Reference Smith and Yakhot1994; Boffetta, Celani & Vergassola Reference Boffetta, Celani and Vergassola2000; Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007), especially in run 1, where the Kolmogorov constant can be estimated as $ C \approx 6$ , indicated by a dashed line. In a later stage, energy accumulates in low wavenumbers of even order, leading to a zig-zag shaped spectrum with high/low energy in even/odd modes. This is first visible at $ l \lt 10$ and subsequently appears at $ l \gt 10$ . The zig-zag spectrum is a manifestation of the symmetric arrangement of the vortices. To see this, we expand the vorticity field in a spherical harmonic system whose polar axis coincides with the polar axis of the great circle of the configuration

(3.1) \begin{equation} \omega (\theta , \phi ) = \sum _{l=1}^{l_{\textit{max}}} \sum _{m=-l}^{l} \omega _l^m Y_l^{m}(\theta , \phi ) , \end{equation}

where $ \theta$ is the colatitude and $ \phi$ is the longitude. A symmetric configuration of four vortices of equal strength on the equator satisfies the symmetries $ \omega (90^\circ + \theta , \phi ) = \omega (90^\circ - \theta , \phi )$ and $ \omega (\theta ,\phi + 90^\circ ) = - \omega (\theta , \phi )$ . The spherical harmonics that satisfy these symmetries are those with $ \{ l = \mathrm{even}, m= 2+4i \}$ , where $ i$ is an integer. In figure 3, we see the enstrophy $ m$ -spectra from run 2 for (a) $ l = 2, 4, 6, 8,10$ and (b) $ l = 20,30,40$ . In both plots, essentially all enstrophy is found in $ m = 2+ 4i$ . A close inspection of figure 2 reveals that there is a secondary period in $ E(l)$ , with extra strong peaks in $ l = 2+4i$ .

Clearly, the spherical harmonic energy spectrum of the condensate is very different from the continuous Fourier spectrum of the two coherent vortices in the condensate simulated by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007). However, also in this case it can be argued that energy accumulates in certain modes, due to the overall symmetry, although this is not revealed by the Fourier spectrum, $ E(k)$ , which only contains information of the averaged energy as a function of the magnitude, $ k$ , of the wavenumber vector $ \boldsymbol k$ . A vorticity field of two coherent vortices of equal strength and opposite sign in a $ 2\pi$ -periodic domain, whose centres are located along a diagonal, obeys a symmetry of the form $ \omega (x+\pi , y-\pi ) = - \omega (x, y)$ (see figure 1c in Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007). All energy of such a field is contained in Fourier modes, $ {\boldsymbol k} = (k_x, k_y)$ , satisfying $ k_x + k_y = \mathrm{odd}$ , while there is no energy in modes satisfying $ k_x+k_y = \mathrm{even}$ (see Appendix A).

Figure 4. Energy and enstrophy fluxes from run 1. Panels show (a,b) t = 5; (c,d) t = 1256. Black: total fluxes. Blue: fluxes generated by the first sums in (1.1) and (1.2). Red: fluxes generated by the second sums. Green: fluxes generated by local triads.

Figure 5. Energy and enstrophy fluxes from run 2. Panels show (a,b) t = 5; (c,d) t = 752. Black: total fluxes. Blue: fluxes generated by the first sums in (1.1) and (1.2). Red: fluxes generated by the second sums. Green: fluxes generated by local triads.

A striking feature of the spectra seen in figure 2(c,d) is the attenuation of energy close to $ l_{\!f}$ , which is most clearly seen in run 1 but is also seen in run 2. In the end state of run 1, the zig-zag spectrum is visible all the way up to $ l = l_{\!f}$ , where the energy content is one order of magnitude smaller than in the early stage. A similar attenuation of energy close to the forcing wavenumber is also seen in the Fourier energy spectrum of the end state of the simulation by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007) (see their figure 3a), although it is not as strong as in the end state of our run 1. Most likely, had the simulation been continued, an equally strong attenuation would have been observed in the Fourier spectrum as in the spherical harmonic spectrum. The observation that energy is attenuated close to the forcing wavenumber suggests that the presence of the vortices qualitatively alters the way energy is transferred in wavenumber space. More specifically, it suggests that condensation is a process in which interactions between widely separated scales of motion play a crucial role. To investigate this, we divide the set of all triads into local, defined as $ n \lt 5p$ , where $ n$ is the middle and $ p$ the smallest wavenumber, and non-local, defined as the complementary set. This is the same division of interactions into local and non-local as was introduced by Kraichnan (Reference Kraichnan1971) in Fourier space. Then we use (1.1) and (1.2) to calculate the fluxes generated by each subset. Since each triad conserves both energy and enstrophy, the same is true for each subset. In figures 4 and 5(a,b) we see the energy and enstrophy fluxes at an early time, $t = 5$ , from both simulations. The contributions from the first sums in (1.1) and (1.2), both denoted by $\varPi _1$ , are positive, while the contributions from the second sums, both denoted by $\varPi _2$ , are negative. This is in agreement with the prediction by Fjørtoft (Reference Fjørtoft1953) that the middle wavenumber in each triad interaction acts as a source of energy and enstrophy for the other two wavenumbers. At $t = 5$ , the local energy and enstrophy fluxes are approximately half of the total fluxes close to $ l_{\!f}$ , while the local energy flux is less than half of the total at $ l \sim 100$ in run 2 (figure 5 a,b). This is in general agreement with predictions by Kraichnan (Reference Kraichnan1971) and Chen et al. (Reference Chen, Ecke, Eyink, Rivera, Wan and Xiao2006) that non-local interactions are relatively strong in 2-D turbulence. During the formation of the condensate non-local interactions become increasingly dominant. In figures 4 and 5(c,d) we see the fluxes in the end states of the simulations. In run 1, the fluxes generated by the local triads are essentially zero in the end state, and in run 2 they are considerably smaller in the end state than in the early state. In particular, it is striking that there is such a weak contribution to the enstrophy flux. Since the enstrophy flux through $ l$ is generated by triads where $ p \lt n \lt l \leqslant s$ , the flux through $ l$ is totally dominated by triads such that $ p \lt l/5$ . Thus, the condensate not only absorbs small-scale energy generated by the forcing, it also drains the flow from small-scale enstrophy.

Figure 6. (a) Normalised vorticity as function of normalised distance from the vortex centre, for the four vortices at $ t= 500$ (black) and $ t = 1000$ (magenta) from run 1. (b,c) Vorticity field from run 1 at $ t= 500$ and $ t = 1000$ . (d,e,f) Corresponding plots for run 2.

In figure 6, we see the normalised mean vorticity profiles for the four vortices from both runs at $ t = 500$ and $ t = 1000$ and vorticity fields at these times. Kolokolov & Lebedev (Reference Kolokolov and Lebedev2016) developed a theory for the structures of vortices that are generated by an inverse energy cascade that is balanced by large-scale friction. The theory is not applicable to our fields, since we have no large-scale friction. All of the four vortices are more or less equal in shape in both runs, and they are not wider at $ t = 1000$ compared with $ t = 500$ . The width of the vortices in run 1 is approximately four times that of the vortices in run 2, reflecting the fact that $ l_{\!f}$ is four times larger in run 2 than in run 1. In the limit of large $ l_{\!f}$ , we can thus expect that the condensate would consist of four point vortices of equal strength, configured in the way we observe the vortices in our simulations. Indeed, this is a quite mysterious result. However, it is fully consistent with the simulations on a plane, exhibiting a Fourier energy spectrum $ E(k) \sim k^{-3}$ , which is the form of spectrum generated by two point vortices on a periodic plane (see Appendix A). We will not pretend to be able to explain why the condensate on the sphere consists of four point vortices of equal strength that are configured in this way. However, a comparison with some results derived from the theory of interacting point vortices on a sphere may provide some clues to an explanation.

3.2. Four point vortices on a sphere

The classical theory of a dynamical system of point vortices on a plane (see Eyink & Sreenivasan Reference Eyink and Sreenivasan2006 for a review) was extended to a system of point vortices on a sphere by Polvani & Dritschel (Reference Polvani and Dritschel1993). The theory has been applied to investigate the stability and equilibria of a system of four point vortices on a sphere (Dritschel Reference Dritschel2020) and further extended to general curved surfaces (Dritschel & Boatto Reference Dritschel and Boatto2015). The equations of motion of a system of point vortices on the unit sphere can be written as

(3.2) \begin{equation} \boldsymbol{u}_j = \frac {\mathrm{d} \boldsymbol{n}_j}{\mathrm{d} t} = \sum _{i \ne j} \kappa _i \frac {\boldsymbol{n}_i \times \boldsymbol{n}_j}{1-\boldsymbol{n}_i \boldsymbol{\cdot }\boldsymbol{n}_j} \, , \end{equation}

where $ \boldsymbol{n}_i$ is the normal unit vector of the sphere at the position of the $ i$ th vortex and $ \kappa _i =$ is the strength of the $ i$ th vortex, which is related to its circulation, $ \varGamma _i$ , as $ \kappa _i = \varGamma _i/4\pi$ . The system (3.2) conserves angular momentum and energy defined as

(3.3) \begin{align}&\qquad\qquad {\boldsymbol L} = \sum _{i} \kappa _i \boldsymbol{n}_i , \end{align}

(3.4) \begin{align}& E = -\sum _{i \ne j} \kappa _i \kappa _j \ln (1- \boldsymbol{n}_i \boldsymbol{\cdot }\boldsymbol{n}_j ) , \end{align}

and describes the motions of the point vortices in an inertial frame where the mean angular momentum of the flow is $\boldsymbol L$ and the energy is $E$ . The energy can either can be interpreted as a potential energy of interaction between the point vortices or a kinetic energy of the flow. The expression (3.4) can be deduced from the derivation found in Dritschel & Boatto (Reference Dritschel and Boatto2015).

Since the position of each vortex is specified by two coordinates and the conservation laws provide four holonomic constraints, it appears that a system with $ k$ vortices has $ 2k-4$ degrees of freedom. In the following, we develop a symmetry argument showing that the number is reduced to $ 2k-7$ in the case $ {\boldsymbol L} = {0}$ . In this case, the motion can only depend on how the vortices are configured relative to each other. To count the number of coordinates that are needed to specify such a configuration, fix the first vortex on the sphere. That leaves us with $ 2k-2$ coordinates. However, one of these can be removed, because the whole configuration can be rotated around a main axis going through the point where we have fixed the first vortex, without changing the relative configuration. That leaves us with $ 2k-3$ coordinates. The conservation laws provide four holonomic constraints. Thus, there are $ 2k-7$ degrees of freedom. The four vortex system is thus the smallest possible, and this system has a single degree of freedom. In Appendix B, we formalise this argument. In the case $ {\boldsymbol L} \ne 0$ , which corresponds to a rotating sphere where the flow has zero mean angular momentum, a similar symmetry argument shows that the numbers of degrees of freedom is $ 2k-5$ and the smallest number of vortices is three.

The motion of the four vortex system with zero angular momentum can be represented by constant energy trajectories in a 2-D parameter space. Since the system has a single degree of freedom, there can be little doubt that such trajectories correspond to non-chaotic periodic motions. Extreme points, if they exist, correspond to stable and stationary solutions. To find the constant energy trajectories and possible extreme points of such a system we let $ (\theta _i, \phi _i), \, i = 1,2,3,4$ be the coordinates of the vortices, where $ \theta$ is the colatitude and $ \phi$ is the longitude. As argued above, without loss of generality we can set $ \theta _1 = \pi /2, \phi _1= 0$ and $ \theta _2 = \pi /2$ . In Appendix B, we show that these coordinates specify the positions of the vortices in a rotating frame in which the first vortex is at rest and the second vortex has a single degree of freedom. The remaining five coordinates, $ \phi _2, \theta _3, \phi _3, \theta _4$ and $ \phi _4$ , are constrained by angular momentum and energy conservation, expressed as

(3.5) \begin{align} \kappa _3 \cos \theta _3 + \kappa _4 \cos \theta _4 & = 0 , \end{align}

(3.6) \begin{align} \kappa _1 + \kappa _2 \cos \phi _2 + \kappa _3 \sin \theta _3 \cos \phi _3 + \kappa _4 \sin \theta _4 \cos \phi _4 & = 0 , \end{align}

(3.7) \begin{align} \kappa _2 \sin \phi _2 + \kappa _3 \sin \theta _3 \sin \phi _3 + \kappa _4 \sin \theta _4 \sin \phi _4 & = 0 , \end{align}

(3.8) \begin{align} \nonumber E & = - \kappa _1 \kappa _2\ln (1- \cos \phi _2) -\kappa _1 \kappa _3 \ln (1-\sin \theta _3\cos \phi _3) - \kappa _1\kappa _4 \ln (1-\sin \theta _4 \cos \phi _4) \\ \nonumber &\quad - \kappa _2 \kappa _3 \ln (1-\sin \theta _3 \cos (\phi _3-\phi _2)) - \kappa _2 \kappa _4 \ln (1- \sin \theta _4 \cos (\phi _4 - \phi _2)) \\ &\quad - \kappa _3 \kappa _4 \ln (1-\cos \theta _3 \cos \theta _4 - \sin \theta _3 \sin \theta _4 \cos (\phi _4 - \phi _3)) . \end{align}

In figure 7(a), we plot the constant energy trajectories in the space of $ \phi _2$ and $ \theta _3$ for the case $ \kappa _1 = \kappa _3 = 1, \, \kappa _2=\kappa _4 = -1$ . In figure 7(b) we plot the corresponding constant kinetic energy trajectories for a system of four Rankine vortices, for which the vorticity is constant and equal to $ \pm 1/4r_v^2$ out to a distance, $ r_v= 0.1 a$ , from the vortex centre and equal to zero for larger distance. The strengths of the vortices are then $ \kappa _i = \pm 1$ . The energy levels were calculated numerically (see Appendix C). In both plots, the energy levels have been shifted so that the maximum energy is zero. Apart from this, there is no renormalisation. It is reassuring to see that we obtain very similar results in the two cases. There is a single maximum at $ \phi _2 = \theta _3 = \pi /2$ , corresponding to the configuration we observe. This is consistent with the demonstration given by Dritschel (Reference Dritschel2020) that this is the only steady and stable configuration of two pair of point vortices of equal but opposite strength, in the case $ {\boldsymbol L } = \boldsymbol{0}$ . That the configuration is stationary can intuitively be grasped by noting that the net motion of one vortex induced by the other three is zero. That the configuration is an energy extreme, may be intuitively grasped by regarding $ E$ as a potential energy of interaction.

Figure 7. Constant energy trajectories for the four vortex system. (a) Interaction energy for four point vortices, calculated using the point vortex model. (b) Kinetic energy of four localised vortices, calculated numerically. The energy level has been shifted in both plots, so that the maximum energy is zero. Apart from this, the energy levels have not been renormalised.

Is there an interesting interpretation of this? First, it provides a clue to why the vortices are fixed in their configuration. They cannot move out from it, because that would decrease the energy of the system. Second, it may also provide a clue with deeper significance. In some sense, the flow seeks a maximum energy solution. This observation makes it very tempting to seek an explanation based on a statistical mechanics variational principle applied to an inviscid fluid, such as the maximum entropy principle (Robert Reference Robert1991; Robert & Sommeria Reference Robert and Sommeria1991). However, noting that the formation takes place under the action of a stochastic forcing, leading to a state with constant energy growth and almost constant enstrophy, it is difficult to see how such a principle can be applied. Enstrophy dissipation is a crucial element of the formation process. The configuration seems to be the one that maximises the energy for a given enstrophy. To prove or disprove this, we would need to investigate whether the observed field is the one of maximum energy out of all possible fields with a prescribed enstrophy, which seems to be an overwhelming task. To be more modest, we investigate if it is the maximum energy configuration out of all four point vortex configurations. To do this, we note that the sum of the squares of the strengths is an enstrophy measure. Strictly speaking, the mean enstrophy is infinite for a point vortex system if the strengths are finite. However, if we regard the vortices as Rankine vortices with constant vorticity out to a radius $ d$ , the mean enstrophy is proportional to $ d^{-2} \sum _i \kappa _i^2$ . For every finite $ d$ , the sum of the squares of the strengths will thus be an enstrophy measure. Thus, we numerically determine the maximum energy configuration of all configurations satisfying (3.5)–(3.7) and all combinations of $ \kappa _1, \kappa _2, \kappa _3, \kappa _4$ satisfying

(3.9) \begin{align}& \kappa _1 + \kappa _2 + \kappa _3 + \kappa _4 = 0 , \end{align}

(3.10) \begin{align}& \kappa _1^2 + \kappa _2^2 + \kappa _3^2 + \kappa _4^2 = 4 . \end{align}

For a given set of parameters $\kappa _1$ , $\kappa _2$ , $\theta _3$ and $\phi _2$ , the associated solutions of the vortex positions and strengths satisfying equations (3.5)–(3.7) and (3.9)–(3.10) were identified and the energy was calculated by means of equation (3.8). The search in the 4-D space was first made by a random search algorithm to seek the highest local maximum, identified in the neighbourhood of the great circle configuration $\kappa _1=1$ , $\kappa _2={-}1$ , $\theta _3=\pi /2$ and $\phi _2=\pi /2$ . After that, a numerical optimisation based on the Broyden–Fletcher–Goldfarb–Shanno algorithm was performed converging to the same maximum. Indeed, there is a single maximum configuration, which is exactly the configuration of four vortices of equal strength that we observe.

The fact that the configuration is an energy maximum, suggests that it will remain steady if the forcing is turned off. We have carried out two experiments to investigate this. In the first one, we took the full end state field of run 1 as the initial field. In the second one, we made an effort to initialise the field in such a way that it had a maximum degree of symmetry, by setting all modes for which $ l \gt 30$ to zero as well as all remaining modes for which $ l = {\textrm {odd}}$ or $ l = {\textrm {even} }$ and $ m \ne 2+4i$ . In both cases, the configuration remained steady, however, with a spectacular difference that we cannot explain. In the first case, the great circle continued to rotate in a stationary solid body rotation around an axis in its own plane. In the second case, the flow remained completely steady.

3.3. Freely decaying simulations from random initial conditions

Dritschel, Qi & Marston (Reference Dritschel, Qi and Marston2015) undertook simulations of freely decaying 2-D turbulence on a sphere, leading to late time solutions with four coherent vortices, quite similar to the solutions we observe. The simulations were carried out to test statistical mechanics theories based on the 2-D Euler equation (Miller Reference Miller1990; Robert Reference Robert1991; Robert & Sommeria Reference Robert and Sommeria1991), predicting that the solutions will relax to steady equilibrium states. The four vortices in the late time solutions reported by Dritschel et al. (Reference Dritschel, Qi and Marston2015) were of unequal strength and underwent unsteady oscillatory motions. They also reported evidence of small-scale vortices emerging in the flow fields at high resolution, and concluded that the flow remains persistently unsteady with no indications of relaxation to a steady state. Evidently, similar conclusions are not applicable to the end states of our forced simulations, which display a high degree of symmetry and steadiness. To explore the similarities between the solutions of our forced simulations and the late time solutions of decaying simulations, we have carried out three decaying simulations from random initial conditions with a prescribed energy spectrum. In runs 3 and 4, the initial spectrum is centred around $ l_c = 8$ , while the resolution is four times higher in run 4 compared with run 3. The initial conditions of run 3 and run 4 are not identical although the initial spectrum is the same, since the random phases are different in the initial fields. In run 5, the initial spectrum is centred around $ l_c = 100$ .

Figure 8. Evolution of (a) mean energy and (b) enstrophy from runs 3, 4 and 5.

Figure 9. Energy spectra from (a) run 4, (b) run 5. Enstrophy fluxes from (c) run 4, (d) run 5.

In figure 8, we see the evolution of mean energy and mean enstrophy from the three decaying runs. Time is measured in units $ a/\sqrt {E_0}$ , where $ E_0$ is the initial mean energy which is set to unity, as is $ a$ . As seen in the figure, the energy is more or less constant in runs 3 and 4, implying that there is essentially no energy dissipation at all in these runs, while there is some energy dissipation early in run 5, reflecting the fact that it is initialised with a spectrum peaking at $ l_c = 100$ . In all three runs, the flow relaxes to a state in which mean enstrophy, as well as mean energy, can be regarded as constant. In figure 9(a,b), we see the evolution of the energy spectra from runs 4 and 5. In the early stage of run 4, the initial spectrum centred at $ l_c = 8$ , broadens to a wide spectrum of the form $ E(l) \sim l^{-3}$ , in accordance with the prediction of Kraichnan (Reference Kraichnan1967) for the forward enstrophy cascade. This is clearly seen at $ t = 15$ . In a corresponding way, the initial spectrum in run 5, centred at $ l_c = 100$ , broadens to a wide spectrum peaking at $ l \lt 10$ . This broadening is caused by the inverse energy cascade. For comparison, we have added a line with slope $-5/3$ , although it is unclear if the spectrum conforms to this shape. In figure 9(c,d), we see the spectral enstrophy flux from the same runs. Using the value $ \varPi _{\varOmega } = 0.5$ and fitting the spectrum in run 4 at $ t=15$ to $ E(l) = K \varPi _{\varOmega }^{2/3} l^{-3}$ we obtain $ K \approx 2$ , not too far from $ K = 1.3$ , as found by Lindborg & Alvelius (Reference Lindborg and Alvelius2000) in Fourier space. At late times, the spectra approach a steady shape and the enstrophy fluxes become very small. The late time solutions are similar to the late time solutions of Dritschel et al. (Reference Dritschel, Qi and Marston2015), exhibiting four coherent vortices.

Figure 10. Normalised mean vorticity as function of normalised distance from the vortex centre and vorticity fields at $ t= 500$ and $ t = 900$ from (a–c) run 3 and (d–f) run 4. In (d) the profiles at $ t = 500$ (black) and $ t= 900$ (magenta) are indistinguishable.

In figure 10, we see the vorticity profiles for the four vortices and the vorticity fields at $ t = 500$ and $ t= 900$ from runs 3 and 4. The profiles are computed as the average vorticity along a circle with distance $ d$ from the vortex centre. In both cases, there are no signs of other vortices in the field and the profile of each vortex has not changed in the time span $ t \in [500, 900]$ . In the higher resolved run, the profiles at $ t = 500$ and $ t= 900$ are completely indistinguishable. The field consists of four persistent vortices that are slowly advected over the sphere, without any changes in their internal structure. Unlike the forced simulations, the four vortices are not of equal strength. They are wider and less symmetric. Some of the profiles are step-like, as also found by Dritschel et al. (Reference Dritschel, Qi and Marston2015).

In figure 11, we see time series of the geodesic distances between like signed, $ d^{++}$ and $ d^{{-}{-}}$ , and opposite signed vortices, $ d^{+-}$ , from runs 3 and 4. In both cases, the distances undergo periodic motions that are in phase. In run 3, the period of $d^{+-}$ is half the period of $d^{++}$ and $d^{{-}{-}}$ . In run 4, it is the other way around. The motions bear every sign of a dynamical system with a few degrees of freedom, possibly a single one, as in the four point vortex model. It is noteworthy that the configuration is not very far from lying on a great circle in run 3, where the amplitudes of the motions are quite small, while the configuration is further away from a great circle in run 4, where the amplitudes are larger. Interpreting these results in the light of the four point vortex model, the configuration in run 3 is closer to the energy maximum of the four vortex system than the configuration in run 4. We can conclude that increased resolution does not move us closer to the maximum.

Figure 11. Normalised geodesic distances between vortex centres versus time. Distances between opposite signed vortices denoted as $ d^{+-}/a$ and distances between like signed vortices denoted as $ d^{++}/a$ and $ d^{{-}{-}}/a$ . (a,b) Run 3, (c,d) run4.

The theory of Robert (Reference Robert1991) and Robert & Sommeria (Reference Robert and Sommeria1991) predicts that solutions to the Euler equations will relax to states in which there is a functional relation between the vorticity and the streamfunction. Such states are stationary solutions to the Euler equations, since $ {\boldsymbol{u}} \boldsymbol{\cdot }\boldsymbol{\nabla }\omega \equiv 0$ , when $ \omega = F(\psi )$ . The prediction can be tested by making a scatter plot of $ (\psi , \omega )$ . If all points fall on the same line, there is a functional dependence, otherwise not. In figure 12, we present such plots from the forced simulations, run 1 and run 2, and the decaying simulations, run 3 and run 4.

The plot where the points are closest aggregated in the vicinity of a single line is the one from run 1. There are two branches of high and low vorticity, where we see two distinctive lines representing the four vortex cores. Apparently, the vortex cores are very stable. At lower values of vorticity, the points do not fall onto a single line but onto a broader band. The plot from run 2 is similar to the plot from run 1, with some differences. The two main branches of high and low vorticity are steeper and contain higher magnitude values. This difference is a reflection of the fact that the vortices in run 2 are more narrow than in run 1. A close inspection of the two main branches in run 2 reveals that each of them is divided into two sub-branches, corresponding to two vortices of positive/negative sign. This is particularly clear if we look at the negative branch. At lower values of vorticity, the band in which the points aggregate is somewhat broader in run 2 than in run 1. These differences can be attributed to the fact that the effective time span of condensation was longer in run 1 than in run 2. A continuation of run 2 would probably result in a scatter plot coming closer to a stationary solution. We conclude that the condensate is formed in a close vicinity of a stationary solution to the unforced Euler equation.

In the decaying simulations, the points fall onto four branches of bands, corresponding to the four vortices. That fact that the bands have limited width can be interpreted as a sign of a relatively high degree of local stationarity, while the fact that the bands are separated can be interpreted as a reflection of global unsteadiness. The difference between the structure of the two plots from run 3 and run 4 can be attributed to the difference in the initial conditions. Our plots are similar to the corresponding plots of Dritschel et al. (Reference Dritschel, Qi and Marston2015) (see their figure 12), with an important difference. Dritschel et al. (Reference Dritschel, Qi and Marston2015) used two different numerical methods, denoted as Geodesic Grid Method and Combined Lagrangian Advection Method. Using the second method, they found that the scatter plot (their figure 12b) exhibited a thin mist of points structured in vertical lines, falling outside the main branches, while the scatter plot of the first method (their figure 12a) did not show any such feature. Putting a high degree of confidence in the second method, Dritschel et al. (Reference Dritschel, Qi and Marston2015) interpreted this feature as evidence of a multitude of small-scale vortices emerging with increased resolution. Since none of our plots show any such feature, we are sceptical to this interpretation. We see no signs of small-scale vortices emerging with increased resolution in our scatter plots, in our plots of the vorticity fields in physical space or in spectral space.

Figure 12. Scatter plots of $ (\omega , \psi )$ from the end states of (a) run 1, (b) run 2, (c) run 3, (d) run 4.

Figure 13. Evolution of the vorticity field over time in run 5.

We now turn our attention to our last simulation, run 5, that was carried out mainly to investigate if a field that is initialised with a spectrum centred at a large wavenumber ( $ l_c = 100$ ) will display a late time solution consisting of four localised vortices. The answer is yes, although we had to run the simulation for a very long time to confirm this. In figure 13, we see how the vorticity field develops over time as we run the simulation. A multitude of vortices are formed at a quite early time in the simulation. Their number is subsequently reduced as they interact in a random way and undergo vortex pairings. In the end, there are four vortices of unequal strength undergoing periodic motions. In figure 14, we see the motions of the geodesic distances between the vortices, plotted in a similar way as in figure 11. Again, we see regular periodic motions that are in phase, but the pattern is a little bit more complex than in runs 3 and 4. The main period of the modulations in $ d^{+-}$ is the same as the main period of the modulations in $ d^{++}$ and $ d^{{-}{-}}$ . There is also a period of smaller amplitude modulations in $ d^{++}$ and $ d^{{-}{-}}$ which is just as long. On top of this, there is a longer period that is modulating the magnitude of the peaks. A future study will reveal if these motions can be described by the four point vortex model with a single degree of freedom.

Figure 14. Geodesic distances between two vortex centres versus time, from run 5, plotted in the same way as in figure 11.