Two-dimensional flow on the sphere

Equilibrium statistical mechanics predicts that inviscid, two-dimensional, incompressible flow on the sphere eventually reaches a state in which spherical harmonic modes of degrees $n=1$ and $n=2$ hold all the energy. By a separate theory, such flow is static in a reference frame rotating at angular speed $2\Omega/3$ with respect to the inertial frame. The vorticity field in the static frame is an accident of the initial conditions, but, once established, it lasts forever under the stated assumptions. We investigate the possibility of such behavior with a stereographic-coordinate model that conserves energy and enstrophy when the viscosity vanishes.


Introduction
Jack Herring made fundamental contributions to the theory of two-dimensional turbulence and to quasigeostrophic turbulence. The former is the prototype of the latter, which adds the key dynamical ingredients-rotation, density stratification, and topography-required to make the theory relevant to the Earth's atmosphere and oceans. Much of the theory is based based upon planar geometry and the hope that the beta-plane approximation captures the most significant effects of curvature and coordinate system rotation. A primary purpose of this paper will be to show that this is not the case: Two-dimensional incompressible flow on the sphere differs significantly from two-dimensional flow on the plane even in the case of vanishing rotation, or, as we prefer to say, even in the case of vanishing angular momentum.
A second primary purpose is to illustrate the advantages of stereographic coordinates for both analytic and numerical work on the sphere. In stereographic coordinates the equations governing two-dimensional turbulence on the sphere take a form, (31) or (37)(38)(39), that is very similar to the corresponding formulation in Cartesian coordinates on the plane. Only the appearance of the smoothly varying metric coefficient h distinguishes the two dynamics. From a mathematical point of view, h is responsible for all of the differences between flow in the two geometries.
To quickly appreciate the effects of spherical geometry, first consider inviscid, incompressible, twodimensional flow on the plane. Unlike on the sphere, which has no boundary, the prescribed boundary on the plane is determinative. Let it be the rectangle shown in Figure 1.  Let the flow within this rectangle be truncated to exclude all the Fourier modes with wavenumbers greater than some prescribed cutoff k c . Let the initial energy be concentrated at a wavenumber k 0 that lies between k c and the lowest wavenumber as determined by the boundary of the system.
Equilibrium statistical mechanics (Kraichnan, 1967) predicts the statistically steady state attained by this system after an undetermined period of adjustment. Interest attaches to the sequence of equilibrium states that occur as k c → ∞. In this limit the theory predicts that all of the energy ends up in the lowest-wavenumber mode, and that all of the enstrophy not contained in that mode appears in modes near k c . As k c → ∞ this excess enstrophy is expelled to infinity, disappearing in much the same way as if it had been removed by viscosity. The tendency for all the energy to crowd into a single lowest mode is often referred to as 'Bose-Einstein condensation', after a closely analogous result in quantum mechanics.
Suppose that the initial conditions are such that the conserved circulation around the rectangular boundary vanishes. Then the lowest-wavenumber mode, which has a single extremum within the rectangle, remains unexcited. All of the energy flows into the next-lowest mode, the one depicted in Figure 1. Because the boundary is a rectangle (and not, say, a square or a circle), there is only one next-lowest mode. If all of the energy ends up in this mode, it must appear in one of the two states depicted in Figure 1. Moreover, once either of these states is established, the system cannot flip suddenly to the opposite state, because that would require the energy in one of these states to flow temporarily into higher-wavenumber modes, and that, according to statistical mechanics, is extremely unlikely. Thus the two states depicted in Figure 1 act as a single-register binary memory. The state that actually occurs is an accident of the initial conditions, but once established it lasts for a very long time.
Still considering flow on the plane, suppose that there are two lowest modes available to receive all the energy. For the statistical mechanics we take the microcanonical ensemble as seems appropriate. Let the energy be E = x 2 + y 2 where x and y are the amplitudes of the two modes. If x were the only mode, then, as already noted, its probability distribution would be where, here and below, we ignore normalization constants. However, if there are two modes, then the probability distribution of either mode is if x 2 < E and zero otherwise. Like (1), the distribution (2) is sharply peaked at ± √ E but now there is significant probability of finding x anywhere in between. If there are three lowest modes, then the analogous calculation predicts that x is uniformly distributed between ± √ E. The implication of these calculations is clear: If the lowest available mode is degenerate with two or more members, then each member sees the other members as a reservoir from which energy may be borrowed or loaned, allowing leakage between states and a breakdown of memory. In all of these examples the average flow predicted by statistical mechanics is misleading and beside the point. For example, if the ensemble consists of the two equal-probability states depicted in Figure 1, then the average flow vanishes, but a vanishing velocity field is very unrepresentative of either state. Such considerations demonstrate how even the simplest calculations with probability distributions require careful interpretation. Now consider the corresponding equilibrium state on the sphere (Frederiksen and Sawford, 1980). Spherical harmonics replace Fourier modes, and the harmonic degree n with cutoff n c replaces the wavenumber k with cutoff k c . The lowest mode Y 0 0 is irrelevant because the stream function is arbitrary by a constant, and because the vorticity integrates to zero over the sphere (Gauss's constraint). The next lowest modes are three: Y 0 1 , Y ±1 1 . However, the amplitudes of these three modes are determined by the three components of angular momentum, which is conserved by our dynamics even with nonzero viscosity. Thus the lowest-degree modes into which the energy can pile up are those of degree n = 2, namely, Y 0 2 , Y ±1 2 , Y ±2 2 . By the reasoning of the previous paragraph, the Bose-Einstein condensate would be one in which these 5 modes freely exchange energy. Instead, it turns out, the 5 modes of degree 2 lock together in a pattern that is perfectly static in a uniformly rotating reference frame. This behavior is dramatically different from flow on the plane, and is solely a consequence of spherical geometry.
Let the sphere be centered at x = y = z = 0, and suppose, without loss of generality, that the angular momentum (if nonzero) points toward positive z. Then the amplitude of the Y 0 1 -mode is a nonzero constant, and the amplitudes of Y ±1 1 vanish. The Y 0 1 -mode corresponds to a solid-body flow about the z-axis with constant angular speed Ω, and the corresponding vorticity is called the Coriolis parameter. As shown in Section 3, the flow consisting of the Y 0 1 -mode and the 5 Y 2 -modes with arbitrary amplitudes represents an exact solution to the inviscid vorticity equation on the sphere. If Ω = 0 this flow field is steady in the inertial frame. If Ω = 0 the flow field is steady in a reference frame rotating at speed 2Ω/3 with respect to the inertial frame. The amplitude Ω of the Y 0 1 -mode is determined by the conserved angular momentum. The amplitudes of the Y 2 -modes are accidents of the initial conditions in the same way that initial conditions determine which of the two states in Figure 1 eventually emerge. However, once these amplitudes are determined via Bose-Einstein condensation, they persist indefinitely under the stated assumptions. Thus Bose-Einstein condensation into Y 2 -modes represent a potential memory storage that does not, as in planar geometry, require a 'probability well' for its preservation. For the plane, as for the sphere, triads with the same wavenumber magnitude (the analog of spherical harmonic degree) do not interact, either with each other or with modes of different wavenumber magnitude. However, typical planar geometries forbid the occurrence of many modes with the same wavenumber magnitude. For instance, in the case of an infinitely periodic box, in which the wave vector corresponds to integer pairs (n, m), modes with the same n 2 + m 2 are rare, as is seen from a list of Pythagorean Triples.
More generally, flow consisting only of Y 1 modes and modes of a single spherical harmonic degree n is perfectly static in a reference frame rotating at angular speed (3) Ω − 2Ω n(n + 1) with respect to the inertial frame. According to Kochin et al. (1964) this result was known to Ertel. Section 3 supplies details. There are 2n + 1 spherical harmonics with the same degree n. These facts invite us to think of the stream function for flow on the sphere as the sum, of contributions ψ n of degree n, where each ψ n is itself the sum of 2n + 1 spherical harmonics. The coordinates (ξ, η) are arbitrary, but we strongly prefer them to be stereographic coordinates. Again, if only ψ 1 and one ψ n with n > 1 are present, then ψ is steady in the frame rotating at angular speed (3). In the general case where all ψ n are present, the various ψ n interact to produce a turbulent flow. However, the interaction between ψ n with very different autorotation rates (3)-very different n-is likely to be weak. And since the difference between autorotation rates is greatest and varies most rapidly at small n, we expect energy transfer to large spatial scales to slow dramatically as the energy reaches low n, for sufficiently large Ω. This is essentially Rhines (1975)'s 'beta-arrest' phenomenon for flow on the beta plane, but now taking a form particular to the sphere. See also Vallis and Maltrud (1993). This study was begun with the hope that the special static solutions described above might be a factor in intermediate-range weather prediction. Although flickers of that hope remain, there are significant differences between the Earth's atmosphere and the very idealized dynamics considered in this paper. First and foremost, angular momentum conservation applies only to the entire Earth-atmosphere-ocean system. Mountain torque provides the main coupling between the atmosphere and the solid planet. The latter, with its enormously greater moment of inertia, experiences only small changes in its rotation rate, but fluctuations in the Y 0 1 component of atmospheric flow are very significant and lead to a peak in the atmospheric energy spectrum at n = 1 (Boer and Shepherd, 1983). Second, because the Rhines scale corresponding to the Earth's atmosphere is smaller than the planetary radius, beta-arrest prevents the leftward movement of energy through the spectrum from reaching n = 2. As a result of these two factors, the atmospheric energy spectrum actually shows a minimum at n = 2. See Sawford and Frederiksen (1983) for the statistical mechanics including mountain torque and Egger et al. (2007) for a thorough review of atmospheric angular momentum.
The plan of the paper is as follows. In Section 2 and Appendix A we derive the Navier Stokes equations for flow on the sphere in stereographic coordinates. The form of Navier-Stokes viscosity for flow on curved surfaces has been the subject of disagreement. For instance, the viscosity given by Gill (1982) does not conserve angular momentum. We follow the recommendation of Gilbert et al. (2014) that the viscosity conserve angular momentum, and in Appendix A we show that their recommended viscosity takes an especially simple form in stereographic coordinates.
In Section 3 we prove the special solutions described briefly above. Our proof using stereographic coordinates is substantially simpler than the proofs given by Thompson (1982) or Verkley (1984) using spherical coordinates.
In Section 4 and Appendix B we use the method of Salmon and Talley (1989) to construct an Arakawatype model in stereographic coordinates. The model conserves energy and enstrophy when the viscosity vanishes.
Section 5 demonstrates the accuracy of the model and tests our speculation that Bose-Einstein condensation on the sphere produces a static flow. Our results suggest that the time required to achieve the static state is probably infinite, but that the flow rather quickly becomes quasi-static.
It is a pleasure to dedicate this paper to the memory of Jack Herring. Rick Salmon acknowledges the importance of Jack's mentorship at an early stage of his career and holds the happy memories of a lifelong friendship. Rick cannot think of Jack without also remembering his wife Betty, who shared Jack's kind and generous spirit, and greatly enjoyed the foibles and eccentricities of his scientific friends.

Sterographic coordinates
We solve the vorticity equation for incompressible flow on the unit sphere, x 2 + y 2 + z 2 = 1.
We refer to the point (x, y, z) = (0, 0, 1) as the 'north pole', and (x, y, z) = (0, 0, −1) as the 'south pole', even when the angular momentum vanishes, the case commonly referred to as 'non-rotating.' The 'equator' is the intersection of the plane z = 0 with (5). The conserved angular momentum is a vector in the direction of (0, 0, 1). Besides the Cartesian embedding coordinates (x, y, z), we shall refer to three other coordinate systems. The first of these are the usual spherical coordinates, defined by x = cos θ cos λ y = cos θ sin λ z = sin θ where λ is the longitude and θ the latitude. The (λ, θ) coordinates cover the sphere except for singularities at the two poles. We also use the stereographic coordinates (7) ξ = x 1 − z , η = y 1 − z defined by a line, emanating from the north pole, that intersects the equatorial plane at (ξ, η, 0) and the sphere at (x, y, z). (ξ, η) cover the sphere except for a singularity at infinity, corresponding to the north pole. We also use defined by the intersection with the equatorial plane of a line emanating from the south pole. (ξ,η) cover the sphere except for singularity at infinity corresponding to the south pole. The stereographic systems (7) and (8) prove to be more useful than the spherical coordinates; the primary challenge is to match them together. Our strategy will be to solve the southern hemisphere dynamics in (ξ, η) within the unit circle to solve the northern hemisphere dynamics in (ξ,η) within the unit circle (10)ξ 2 +η 2 < 1, and to match the two solutions together at the equator, The transformation equations between the two stereographic systems are A useful relation is Points within the unit circle on the (ξ,η) plane transform to points outside the unit circle on the (ξ, η) plane, and vice versa. On the equator itself, (ξ, η) = (ξ,η). It is also useful to record the inverse transformations from either set of stereographic coordinates to the Cartesian embedding coordinates, For a thorough introduction to stereographhic coordinates, see Needham et al. (1997); Needham (2021).
Infinitesimal displacements on the surface of the sphere satisfy ds 2 = dx 2 + dy 2 + dz 2 = cos 2 θdλ 2 + dθ 2 = 4 (1 + ξ 2 + η 2 ) 2 dξ 2 + dη 2 = 4 1 +ξ 2 +η 2 2 dξ 2 + dη 2 where ds is infinitesimal distance tangent to the surface of the sphere. Thus all of our systems fit the general form where (ξ 1 , ξ 2 ) are general orthogonal coordinates, and diag(h 2 1 , h 2 2 ) is the metric tensor. In arbitrary orthogonal coordinates, the incompressibility condition takes the form where the overdot denotes the time derivative following a fluid particle. Since the velocity components tangent to the sphere in the directions of ξ 1 and ξ 2 are respectively In the case of spherical coordinates (ξ 1 , ξ 2 ) = (λ, θ), (h 1 , h 2 ) = (cos θ, 1), and (23-25) take the familiar forms To obtain the dynamics in stereographic coordinates, we proceed from the general equations (24-25). Some care must be taken concerning the handedness of the coordinate systems. Viewed from the exterior of the unit sphere, (ξ,η) and (λ, θ) constitute right-handed systems, whereas (ξ, η) constitutes a left-handed system. That is The general covariant point of view cares nothing about handedness, but we want our computed vorticity and stream function to be the same as would be obtained by solving the entire problem in spherical coordinates, and, in particular, we want these quantities to be continuous at the equator; it would be a great nuisance to keep track of differently defined stream functions in the two hemispheres. With this in mind, the northern hemisphere equations take the forms (31)ĥ 2 q t + ∂(ψ, q) ∂(ξ,η) = 0,ĥ 2 q =∇ 2 ψ,ĥ = 2 1 +ξ 2 +η 2 within the unit circle (10); and the southern hemisphere equations take the forms within the unit circle (9). Here, The Laplace Beltrami operator in (26) or (29) is equivalent toĥ −2∇2 and h −2 ∇ 2 . Apart from the factors of h andĥ, (31-33) have the same form as the corresponding equations for planar geometry in Cartesian coordinates. This greatly facilitates numerical solution. The sign change between (31) and (32) reflects the change in handedness. To see that this makes sense, recall that the first equation in (31) and (32) also applies to a passive scalar advected by the velocity field represented by ψ. The sign change between (31) and (32) ensures that the (cross-equatorial) flux of q out of the unit circle (9) matches the flux of q into the unit circle (10), with no need to redefine or re-interpret ψ or q.
Alternatively, by replacing the vorticity q in (31) and (32) by ξ 2 +η 2 + 1 we obtain the equations of motion in a frame rotating at angular speed Ω.
We solve (31) on (10), and (32) on (9), with the matching conditions that ψ and q be continuous at the equator (11). If the initial condition is given by prescribing q, then q must satisfy the Gauss constraint, where the integrals are over the two unit circles (i.e. the whole sphere), but the constraint (35) is automatically maintained thereafter. If the initial condition is given by prescribing ψ, the Gauss constraint is automatically satisfied. Although we work primarily with (31) and (32), it is useful to record the corresponding forms of the momentum equations. Introducing the notation, the southern hemisphere equations take the forms The northern hemisphere equations take the same form as (37-40) but with hats applied to all the variables except q. Again, if h 2 = 1, (37-40) are formally identical to the equations for two-dimensional incompressible motion in planar geometry. However, no choice of variables can make h 2 equal unity on the unit sphere; instead it must obey the requirement that the Gaussian curvature be uniform and equal to unity. The effects of curvature can never be transformed away, but it is interesting that they can be confined to a single factor in (40). This is what motivates the definition (36). However, it must be emphasized that the symbols (u, v) invite misinterpretation: the velocity of fluid particles tangent to the sphere is and not (u, v).
The inviscid dynamics (31-32) or (37-40) conserve energy in the form, total vorticity in the form and enstrophy in the form The integrations in (43-45) are over the whole sphere, with the understanding that they are actually computed as the sum of integrations over the two unit circles. The vanishing of (44) corresponds to the Gauss constraint. More generally, the inviscid dynamics conserves Casimirs of the form where F (q) is any function of the vorticity for which the integral in (46) converges; (44) and (45) are particularly important cases of (46). The inviscid dynamics also conserves angular momentum. Angular momentum is a vector with 3 components, each of which is conserved. However, since the location of the axis that determines the origin of our stereographic system has no particular significance, we may assume, without loss of generality, that the angular momentum vector points in the z-direction. Then the sole nonvanishing component of angular momentum is It remains to discuss the viscosity. Gilbert et al. (2014) examine three proposed definitions of viscosity for incompressible flow on a curved surface. All three are covariant, and none of the three can be rejected on purely fundamental grounds. However, only one of the three forms of viscosity conserves angular momentum. We agree with Gilbert et al. (2014) that this angular-momentum conserving property is the proper choice for applications in which the viscosity is actually an eddy viscosity representing the mixing-but not the destruction-of angular momentum by unresolved motions. In Appendix A we show that the viscosity suggested by Gilbert et al. (2014) takes the surprisingly simple form ∂u ∂t + · · · = −ν ∂q ∂η (48) ∂v ∂t + · · · = +ν ∂q ∂ξ (49) in the dynamics (37-40), and (50)ĥ 2 q t + · · · = ν∇ 2 q, h 2 q t + · · · = ν∇ 2 q in the dynamics (31-32). The ellipses represent the conservative terms already discussed. When the viscosity in (50) is added to the dynamics (31-32) we find that Thus viscosity dissipates energy and enstrophy, but not angular momentum. The importance of angular momentum conservation by viscosity seems not to be widely appreciated.
Many modelers prefer a more scale-selective hyperviscosity to any form of Navier-Stokes viscosity. However, hyperviscosity introduces unphysical effects. For example, hyperviscosity creates artificial local extrema in vorticity. For this reason, we believe that Navier-Stokes viscosity is well worth the extra cost in spatial resolution.
We emphasize that, in our dynamics, the sphere plays no role except to shape the flow. Without drag or mountain torque, only curvature and angular momentum conservation distinguish our dynamics from the flow in a box on the plane.

Exact solutions
Let Y n be any solution of (54) ∇ 2 LB Y n = −n(n + 1)Y n where n is a positive integer. The coordinate system used to represent Y n is arbitrary. In a spherical coordinate system, the general solution of (54) is where P m n is a Legendre function of the first kind, and A m and γ m are arbitrary constants. We recognize (55) as the sum of all spherical harmonics of degree n with arbitrary amplitudes A m and phases γ m . The customary normalization coefficients have been lumped into the A m . The representation (55) applies to any system of spherical coordinates. That is, if we transform (55) to an alternative set of spherical coordinates-obtained, for example, by changing the polar axis with respect to which the spherical harmonics are defined-we obtain an expression identical to (55) but with a different set of A m 's and γ m 's. We may in fact use any system of coordinates, including the more convenient stereographic coordinates, because the following proof relies only on the covariant property (54).
Let ψ Ω (l,m,n) be the streamfunction corresponding to solid-body flow at angular speed Ω about an axis in the direction of the unit vector (l, m, n). Let Y n be a solution of (54), and let Y ω n be the pattern Y n rotating at the angular speed in spherical coordinates. However, the axis (l, m, n) is arbitrary, and, more importantly, the spherical harmonics can be defined with respect to any axis whatsoever; it only matters that they are all of degree n. In fact, the solution (57), properly interpreted, holds in any system of coordinates. Viewed from a reference frame rotating with the solid-body flow-that is, viewed in rotating coordinates-the solution (57) corresponds to a retrograde, 'pseudo-westward', propagation of the Y n -pattern at an angular velocity equal to the last term in (56), which differs from the frequency of Rossby-Haurwitz waves in its lack of a factor m in the numerator. This m-factor is missing because we here define ω as angular frequency rather than wave frequency. That is, in (58) we write m(λ − ωt) rather than mλ − ωt. Kochin et al. (1964) attribute the solution (57) to Hans Ertel in 1945, but they do not give a reference. See also Rochas (1984). It was re-discovered by Herring's colleague Thompson (1982) for the special case in which (l, m, n) coincides with the z-axis, as in (58), and Y n consists of a single spherical harmonic whose axis of symmetry is inclined with respect to the z-axis by a prescibed angle. Thompson proposed this special exact solution as a test of numerical codes based on spherical coordinates. Verkley (1984) matched Thompson's solution in one region of the sphere with an analogous solution, in another region, in which Y n is replaced by the solutions of ∇ 2 LB Y n = dY n with d a positive constant. The matching conditions are that the two solutions precess about the z-axis at the same angular speed and have continuous derivatives at the boundary between regions. The family of solutions discovered by Ertel, Thompson, and Verkley are noteworthy in that they seem to be the only analytic solutions of the vorticity equation on the sphere (besides 'zonal' flows in which ψ = F (θ) with F arbitrary) that do not include point vortices. See also Neven (1992) and references therein.
The last term in (65) is needed to satisfy the Gauss constraint, but it does not contribute to the vorticity equation, in which only derivatives of q appear. By the chain rule, h 2 q t = −n(n + 1)h 2 ∂ ∂t Y n (s,s) = −n(n + 1)iωh 2 −s ∂Y n ∂s +s ∂Y n ∂s where we have used the identities For the second term in the vorticity equation we have The sum of (66) and (68) vanishes for ω given by (56). Thus (62) satisfies our dynamics for any positive integer n. This proof is simpler than the proofs given by Thompson (1982) or Verkley (1984) using spherical coordinates.

Numerical model
We solve (31-32) by the method of Salmon and Talley (1989), a generalization of Arakawa (1997)'s method that conserves discrete analogues of the energy (43), total vorticity (44), and enstrophy (45) when ν = 0. This method has been further generalized to the shallow water equations; see Salmon (2009) and references therein.
We begin by noting that the vorticity equation in (32) is equivalent to the statement that for any function α(ξ, η). The integration is over the entire sphere. By parts integrations, (69) takes the more useful form Our strategy is to discretize all the variables in (70) in such a way that the conservation of energy, total vorticity, and enstrophy correspond to purely algebraic cancellations between the terms. The conservation of enstrophy (45) corresponds to the choice α = q in (70). This requires that the discrete form of (71) dξdη q ∂(ψ, q) ∂(ξ, η) + ψ ∂(q, q) ∂(ξ, η) + q ∂(q, ψ) ∂(ξ, η) vanish, and this occurs automatically if the discrete Jacobian is antisymmetric. Conservation of total vorticity corresponds to α = 1, and requires that the discrete form of (72) dξdη ∂(ψ, q) ∂(ξ, η) vanish. Again, this occurs automatically. Energy conservation corresponds to α = −ψ, and we discuss it further below.
We discretize (70) by covering the interior of both stereographic circles with quadrilateral finite elements. The placement of the elements is the same in both circles. Most of the elements are perfect squares, as illustrated in Figure 2 for a very low resolution case, but near the unit circle itself the quadrilaterals are deformed, with one or more nodes lying on the circle. The nodal values of ψ and q are the dependent variables of the model. Each interior node corresponds to a point on one of the hemispheres, and is shared by 4 elements. Each equatorial node is shared by up to 6 elements, counting elements on both sides of the equator. At every time, the state of the system is determined by the values of ψ or q at N = 2N i + N eq nodes, where N i is the number of interior nodes within each of the two unit circles, and N eq is the number of (shared) equatorial nodes. Appendix B offers further details.
In claiming conservation properties, we regard time derivatives as exact. That is, we ignore errors that result from the finite time step. Experience shows these errors to be small in comparison to those that result from space discretization. This issue arises when we consider the finite-element discretization of the left-hand side of (70), namely where the summation is over all N nodes, A i is the area of the ξη-plane associated with node i (in a sense to be made precise), and α i , h i , q i are the nodal values of α, q and h. Note that the metric factor (74) h i = 2 1 + ξ 2 i + η 2 i is determined by the arrangement of the quadrilateral elements, and is independent of time. Replacing the left-hand side of (70) with (73), and the right-hand side of (70) with the finite-element representation described briefly above and more thoroughly in Appendix B, we obtain the model dynamics by requiring that the discrete form of (70) vanish for arbitrary nodal values α i . Setting α i = q i and treating the time derivative as exact, we conclude that the discrete enstrophy Energy conservation corresponds to the choice α i = −ψ i . Again the right-hand side of (70) vanishes if the discrete Jacobian is antisymmetric, and hence To make (76) the time derivative of a discrete energy, we must consider the Poisson equation in (32) relating ψ and q. It is equivalent to the statement that (77) − dξdη βh 2 q = dξdη ∇β · ∇ψ for any function β(ξ, η). The arbitrary function β(ξ, η) is analogous to α(ξ, η). Let be the discrete approximation to (77), where w ij is the set of weights that arise by discretizing the right side of (77). We require (78) to vanish for any set of β i . The β i are completely arbitrary, and can take different values at different times. By first regarding the β as time independent, we obtain the general result, Then, setting β i = ψ i (t), treating time derivatives as exact, and assuming that the weights w ij are symmetric (w ij = w ji ), we obtain Since by (76) this must vanish, we conclude that the discrete energy In overall summary, enstrophy is conserved if (70) is discretized in such a way that its Jacobians are antisymmetric, and energy is conserved if (77) is discretized in such a way that the dot product is symmetric. Salmon (2009) shows how these two steps may be interpreted more generally in terms of the Poisson bracket and the Hamiltonian of the system. The model does not conserve other Casimirs of the form (46). Appendix B gives further details of the discretizations.

Numerical solutions
For flow with rms velocity U on a sphere of radius a with angular momentum corresponding to the solid body rotation rate Ω, the Rossby number Ro = U/Ωa and the Rhines scale Rh = Ro 1/2 a. For values typical of the Earth's atmosphere-U =10 m/sec, a=6400 km, and Ω = 2π/day-we have Ro = 0.02 and Rh = 0.15a. For Earth's ocean and the giant planets the ratio of the Rhines Rh scale to the planetary radius a is even smaller. This regime, which features zonal jets with widths on the order of Rh, has attracted the greatest interest from modelers, but the non-rotating case Ω = 0 is also interesting, because it reveals the effects of curvature by itself.
In our experiments U = a = 1, where U is the rms velocity at the intial time in the rotating frame, that is, the rms velocity excluding solid body rotation. The Rossby number Ro = 0.02 corresponding to the Earth's atmosphere then corresponds to Ω = 50. We present solutions of the two cases Ω = 0 and Ω = 50. In the solutions with Ω = 50, the displayed stream function and vorticity do not include the n = 1 modes associated with the angular momentum and therefore depict the fields in the rotating frame. In all of our experiments, the viscosity, if nonzero, has the value ν = q rms ∆ 2 , where q rms is the rms vorticity and ∆ 2 is the area of the square gridboxes on the stereographic plane. The number N = 241, 000 of nodes corresponds to the cutoff n c = 490 in spherical harmonic degree n. Figure 3 shows the streamfunction and vorticity in 3 views of the sphere at time t = 2.0 in an experiment with Ω = 0 (the non-rotating case). The initial energy is confined to spherical harmonics of degree n = 16, 17, 18 with the amplitudes randomly assigned. Figure 4 shows the same fields in an experiment beginning from the same initial conditions but with Ω = 50. In Figures 3 and 4, (a) is a polar view of the stream function in the northern hemisphere, as seen by an observer at z = ∞ above the north pole, and (d) is the same field within the stereographic cirlê ξ 2 +η 2 < 1 corresponding to the northern hemisphere. Although both views cover the entire northern hemisphere, the polar view (a) severely distorts features near the equator: The ratio between true distance and apparent distance in (a) varies between unity at the pole and infinity at the equator. In contrast, the stereographic view (d) introduces a scale distortion that varies only between 1 and 2. It offers the superior view of the entire northern hemisphere field. Figure 4. The same as Figure 3, but for the case Ω = 50 of Earth-like rotation. This solution differs from that of Figure 3 in that isolated coherent vortices are less prominent, and that elongated, jet-like structures are beginning to appear, especially near the equator in (c), where the β-effect is greatest. Figure 4 shows the same 5 views as in Figure 3, but for the case Ω = 50 of nonvanishing angular momentum. Figure 4 differs from Figure 3 in that its isolated coherent vortices are less prominent than in Figure 3, and elongated, jet-like structures are beginning to appear, especially near the equator (Figure  4c), where the β-effect is greatest. Figures 3 and 4 illustrate the advantage of stereographic projection. Whereas at least six conventional views of the sphere are needed to give a fair representation of a single field, two stereographic projections offer a complete and relatively undistorted depiction.
Next we consider two experiments in which the initial state is sharply confined to spherical harmonics of degree n = 6. Angular momentum corresponding to a solid body rotation rate Ω = 50 is present, but now it points in the direction of the y-axis in the Cartesian embedding space. This is achieved by replacing the Coriolis parameter f in (34) by (82) 2Ωy = 2Ω 2η ξ 2 + η 2 + 1 = 2Ω 2η ξ 2 +η 2 + 1 .
According to Section 3, the n = 6 initial state should rotate rigidly, in a retrograde manner, about the y-axis at the angular rate ω ≡ 2Ω/(6(6 + 1)). Figure 5 shows this motion as seen by an observer at x = ∞. From its perspective, the pattern rotates upward in the figure, around the sphere, completing a full cycle of revolution in the recurrence time 2π/ω = 2.639 shown on the far right of the figure. The solutions depicted in Figure 5 are both a check on the theory and a strong test of the numerical code, because the rotation carries the pattern from the interior of one stereographic circle to the other and then back again, twice crossing the irregularly shaped elements near the bounding circles (11). The inviscid experiment (bottom row of Figure 5) exhibits small, grid-scale oscillations in the vorticity field, which do not appear in the viscous experiment (top row). Finally, we explore the suggestion in Section 1 that solutions lock into a static state as Bose-Einstein condensation occurs. Theory predicts Bose-Einstein condensation only as n c → ∞ and t → ∞-two limits that are impossible to achieve numerically. The interesting question is whether quasi-static states appear with finite resolution and in finite time. According to equilibium statistical mechanics (Frederiksen and Sawford, 1980), the equilibrium energy in spherical harmonic degree n is (83) E n = 2n + 1 α + βn(n + 1) , n ≥ 2 where the constants α and β are determined by the energy and enstrophy. The equilibrium spectrum (83) is the analogue of Kraichnan (1967)'s spectrum (84) E(k) = k α + βk 2 for flow on the plane. The numerator in (83) is the number of modes of degree n. Bose-Einstein condensation corresponds to the limit α → −2(2 + 1)β as n c → ∞.
Several authors investigate freely decaying turbulence on the sphere and its resemblance to the predictions of equilibrium statistical mechanics, including the more recondite Miller-Sommeria-Robert-Montgomery (MRSM) theory. See Qi and Marston (2014), Dritschel et al. (2015), Modin and Viviani (2020), and Jagad and Samtaney (2021). For a summary of the MRSM theory, see Eyink and Sreenivasan (2006). The requirement that t → ∞ has proved formidable. MRSM theory claims to predict the structure of coherent vortices by the principle of maximum entropy. Salmon (2018) offers a different picture, in which coherent vortices are low entropy sites that exist in order to maximize the rate of entropy production in the area outside the vortices.
Here we consider the case Ω = ν = 0 of nonrotating, inviscid flow with a modest resolution corresponding to n c = 240. The initial condition corresponds to energy distributed equally among the three degrees n = 4, 5, 6, that is, E 4 = E 5 = E 6 with all other E n = 0. The amplitudes of the spherical harmonics within each band are randomly assigned, and the initial state is normalized such that u rms = 1. A calculation based upon (83) predicts that 99.6% of the energy will end up in n = 2, the lowest available mode. Again, E 1 = 0 at all time by the conservation of angular momentum. Figure 6 summarizes the evolution of the flow from t = 0 to t = 60. The time t = 60 is the time required for a fluid particle to circumnavigate a great circle 60/2π = 9.5 times at the rms velocity of the flow. The dashed curve in Figure 6 is the fraction of energy in n = 2. Beginning from zero, it only gradually approaches its predicted value of 0.996.
The solid curves in Figure 6 measure the approach to a static configuration. Let A 2,m be the amplitude of Y m 2 in the spherical harmonic representation of ψ, and define (85) p = (A 2,0 , A 2,1 , A 2,−1 , A 2,2 , A 2,−2 ) [(A 2,0 ) 2 + (A 2,1 ) 2 + (A 2,−1 ) 2 + (A 2,2 ) 2 + (A 2,−2 ) 2 ] 1/2 . Thus p is a unit vector in the 5-dimensional space spanned by the amplitudes of the 5 spherical harmonics comprising n = 2. If the flow were indeed to become static, then the five components of p would become time-independent at values that are accidents of the initial conditions. Our choice of conventionally defined spherical harmonics as amplitudes to be tested is both arbitrary and irrelevant; any other choice of mutually orthogonal functions that span ψ 2 would serve as well. Figure 6. Evolution of a non-rotating, inviscid flow in which the initial energy is confined to spherical harmonic degrees n = 4, 5, 6. The dashed curve (right scale) is the fraction of energy in n = 2, which vanishes initially, but, according to theory, will eventually reach 0.996. The other curves, labeled 2, m, represent the normalized coefficient of Y m 2 in the spherical harmonic representation of the flow, which vary between ±1 (left scale). Figure 6 shows that the approach to the putative equilibrium state is extremely slow. Although E 2 generally increases, its increase is not monotonic, but instead exhibits persistent, relatively rapid oscillations. The components of p, labeled n, m in the figure oscillate over their full range of ±1, but the periods of their oscillations appear to lengthen as time increases. It is conceivable that this flow eventually reaches a static state, but this evidently requires a very long time. On the other hand, since the time unit is the time required for a fluid particle to move a distance equal to the radius of the sphere, the flow at t = 60 is already quasi-static in the sense that it requires at least several such time units for the components of p to change appreciably.
Substituting (B6), (B7) and the corresponding expressions for the other Jacobians back into (B5), we obtain (B8) dξdη ∇β · ∇ψ = 4 A 1234 4 i,j,n,m=1 Γ ijnm [β i η j ψ n η m + ξ i β j ξ n ψ m ] = 4 i,j=1 The C ij may be calculated analytically computed and stored. The final sum in (B8) represents the contribution of a single element. To discretize the integral (B4) over the entire sphere, we sum over all the elements to obtain (78). The weight w ij vanishes unless node i and node j share at least one element. Thus, up to 6 elements contribute to each w ij . However, since C ij is symmetric for each contributing element, w ij is symmetric as well, and therefore the discrete energy (81) is conserved.
Summing over all the elements gives the discrete approximation to (70). The left-hand sides of (69) and (77) take the discrete form where f i = dq i /dt in the case of (69) and f i = q i in the case of (77). In (B15) the summation is over nodes, and A i is the area of the ξη-plane that is 'assigned' to, i.e. closest to, node i. The precise rule of assignment is arbitrary, but it must obey the normalization requirement that the sum of the A i 's equals the area within a unit circle, or that the sum of the h 2 i A i 's on both unit circles equals the area of the unit sphere.
By requiring the coefficient of each β i in (78) to vanish, we obtain N equations representing a logical discrete approximation to h 2 q = ψ ξξ + ψ ηη . By requiring the coefficient of each α i in the above-described discrete analogue of (32) to vanish, we obtain N equations representing a logical discrete approximation to (32).