1. Introduction
There are two outstanding papers in the vast literature on two-dimensional turbulence – the papers by Fjørtoft (Reference Fjørtoft1953) and Kraichnan (Reference Kraichnan1967). Fjørtoft (Reference Fjørtoft1953) showed that energy cannot consistently flow in a single direction in wavenumber space, due to the fact that there is a second conserved quantity in a two-dimensional non-divergent flow, namely squared vorticity or enstrophy. Correctly, he predicted that energy on average will flow from large to small wavenumbers in spherical harmonic space for a two-dimensional flow on a sphere, meaning that there is an energy transfer from small to large physical scales. Kraichnan (Reference Kraichnan1967) developed the theory of the double cascade. If a source is injecting energy in a two-dimensional flow on a plane at spatial scales corresponding to a particular wavenumber 
$k_f$, energy will consistently flow in the direction of smaller wavenumbers and enstrophy in the opposite direction. A constant energy flux range will develop at wavenumbers smaller than 
$k_f$, with a Fourier energy spectrum of the form
\begin{equation} E_F(k) = C|\varPi_{E}|^{2/3} k^{{-}5/3} , \end{equation}
where C is a non-dimensional constant and 
$| \varPi _E |$ is the upscale energy flux, while a constant enstrophy flux range will develop at larger wavenumbers, with an energy spectrum of the form
\begin{equation} E_F(k) = \kappa \varPi_{{\mathcal{E}}} ^{2/3} k^{{-}3} , \end{equation}
where 
$\kappa$ is a non-dimensional constant and 
$\varPi _{{\mathcal {E}}}$ is the downscale enstrophy flux. In a subsequent paper Kraichnan (Reference Kraichnan1971) suggested that (1.2) also should contain a logarithmic correction. Essential parts of the theory have been confirmed in experiments (e.g. Paret & Tabeling Reference Paret and Tabeling1997; Kellay & Goldburg Reference Kellay and Goldburg2002) and numerical simulations (e.g. Boffetta Reference Boffetta2007; Boffetta & Musacchio Reference Boffetta and Musacchio2010). Energy flows upscale and enstrophy flows downscale relative to a source and the energy spectrum often develops close to the predictions, although the 
$k^{-3}$ spectrum in the enstrophy cascade range only can be reproduced in simulations with extremely high resolution (Lindborg & Vallgren Reference Lindborg and Vallgren2010). For an excellent review on experiments and simulations, see Boffetta & Ecke (Reference Boffetta and Ecke2012). One of the most fascinating predictions made by Kraichnan (Reference Kraichnan1967) is that a ‘condensate’ will be formed when the upscale energy cascade reaches the lowest-order wavenumber mode 
$k_0$, corresponding to the physical length scale of the domain confinement. When this happens, the flow undergoes a transition from a chaotic to a more ordered state – a prediction that has been verified by numerical simulations (Smith & Yakhot Reference Smith and Yakhot1994; Chertkov et al. Reference Chertkov, Connaughton, Kolokolov and Lebedev2007).
There are several unresolved issues connected to the double cascade theory. In addition to the studies that seem to confirm the theory, there are a number of numerical studies (Boure Reference Boure1994; Danilov & Gurarie Reference Danilov and Gurarie2001; Tran & Bowman Reference Tran and Bowman2004; Scott Reference Scott2007; Vallgren Reference Vallgren2011) that display conflicting results. In all these simulations, the spectrum is considerably steeper than 
$k^{-5/3}$ and in all simulations a multitude of coherent vortices are formed in the upscale cascade range. A more complete review of these results is given by Burgess & Scott (Reference Burgess and Scott2017). There is also an issue regarding locality of triad interactions. Kraichnan (Reference Kraichnan1967) pointed out that his theory requires that the cascades are sufficiently local in Fourier space and raised concerns as to whether this is true for the enstrophy cascade. Numerical studies (e.g. Ohkitani Reference Ohkitani1990; Maltrud & Vallis Reference Maltrud and Vallis1991) have confirmed that the enstrophy cascade is, indeed, dominated by highly non-local interactions. There is also numerical evidence (Vallgren Reference Vallgren2011) that the upscale energy cascade is driven by highly non-local interactions. It is not clear if such a high degree of non-locality as observed in the simulations is consistent with the double cascade theory. Finally, there are also open questions related to the condensate hypothesis. Given a forcing wavenumber 
$k_f$ and an energy injection rate 
$P$, at what time does the condensate start to form? Will the condensate encompass all wavenumbers from 
$k_0$ to 
$k_f$ as suggested by the simulation by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007), or just a limited wavenumber band close to 
$k_0$, as suggested by the simulation of Smith & Yakhot (Reference Smith and Yakhot1994)?
With few exceptions (e.g. Falkovich & Gaweȩdzki Reference Falkovich and Gaweȩdzki2014) most theoretical and numerical studies of two-dimensional turbulence have investigated the standard case of a flow confined to a plane with double periodic boundary conditions. There are two features of such a flow that makes it a very attractive object of investigation. First, the Fourier description is very transparent from an analytical point of view. Second, using the fast Fourier transform it is relatively easy to design efficient pseudo-spectral codes simulating such a flow. However, it is no longer extremely demanding to carry out highly resolved simulations of a flow on a sphere. Quite recently, it has been shown (Wedi, Hamrud & Mozdzynski Reference Wedi, Hamrud and Mozdzynski2013) that a fast spherical harmonic transform can be effectively implemented in a global circulation model. There is also a GPU accelerated version of the slow transform which is highly competitive (private communication with N. Wedi). This development has not been matched by a corresponding analytical development. Since the study of Boer (Reference Boer1983) there have been few advancements of the spectral analysis of a flow on a sphere. The aim of the present study is to bring the spherical harmonic description to the same degree of transparency as the Fourier description and thereby pave the way for studies in which the unresolved issues of two-dimensional turbulence can be investigated by means of simulations of a flow on a sphere.
2. Kinetic energy spectrum on a sphere
In turbulence theory, the kinetic energy wavenumber spectrum is often defined as a Fourier transform of the covariance function of velocity measured at two different points and an assumption of isotropy is often made to simplify the description. Boer (Reference Boer1983) has advanced such an approach to the analysis of two-dimensional incompressible turbulence on a sphere, with the Fourier transform replaced by the spherical harmonic transform. We will adopt a different approach and argue that the constraint of isotropy does not lead to any simplifications but only to restrictions. Fortunately, the spherical geometry offers another attractive mode of description where the spatial mean value over the entire sphere is the only statistical measure that is used and no assumption of isotropy is invoked.
We let 
$\psi$ be a real scalar function defined on a sphere with radius 
$a$, 
$(\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z)$ the unit vectors associated with a Cartesian coordinate system with origin at the centre of the sphere, and 
$(\theta, \phi )$ the corresponding polar coordinates, where the 
$z$-axis is the polar axis, 
$\theta$ is the colatitude and 
$\phi$ is the longitude. We expand 
$\psi$ in spherical harmonics as
$$\begin{gather} \psi(\theta, \phi) = \sum_{l=0}^{\infty} \psi_{l}(\theta, \phi) , \end{gather}$$
$$\begin{gather}\psi_{l}(\theta, \phi) = \sum_{m={-}l}^{l} a_{lm} Y_{l}^{m}(\theta, \phi) , \end{gather}$$
where we have separated the summation into two different equations to emphasise that the 
$l$- and 
$m$-numbers have different significance. The spherical harmonics are defined as
\begin{equation} Y_l^m(\theta, \phi) = \sqrt{(2l+1) \frac{(l-m)!}{(l+m)!}} \, P_l^m(\cos \theta) {\rm e}^{{\rm i}\phi}, \end{equation}
where 
$P_l^m$ are the associated Legendre functions. They are eigenfunctions of the Laplace operator and orthogonal with respect to integration over the sphere,
$$\begin{gather} \nabla^2 Y_l^m ={-}\frac{l(l+1)}{a^2} Y_l^m , \end{gather}$$
$$\begin{gather}\frac{1}{4 {\rm \pi}} \int Y_l^m {Y_s^p}^* \,{\text{d}} \varOmega = \delta_{ls} \delta_{mp} , \end{gather}$$and have the symmetry property
\begin{equation} Y_l^{{-}m} = ({-}1)^m {Y_l^{m}}^{*} , \end{equation}
where 
$\varOmega$ is solid angle and superscript 
$*$ denotes the complex conjugate. The fact that 
$\psi$ is real imposes the constraint
\begin{equation} a_{l({-}m)} = ({-}1)^m {a_{lm}}^{*} , \end{equation}
on the expansion coefficients for 
$\psi$, implying that 
$\psi _l (\theta, \phi )$ also is a real function.
As pointed out by Tang & Orszag (Reference Tang and Orszag1978), the 
$l$-number is more significant than the 
$m$-number since it is only 
$l$ that enters into the eigenvalue of the Laplace operator. The significance of the 
$l$-number is more clearly revealed by using a result from group representation theory to show that 
$\psi _l (\theta, \phi )$ is invariant under a coordinate transformation between systems with different polar axes. A new Cartesian coordinate system is defined as 
$R [(\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z)] = (\boldsymbol {e}_{x^\prime }, \boldsymbol {e}_{y^\prime }, \boldsymbol {e}_{z^\prime })$, where 
$R(\alpha, \beta, \gamma )$ is a rotation defined by the Euler angles 
$(\alpha, \beta, \gamma )$. We let 
$(\theta ^{\prime }, \phi ^{\prime })$ be the polar coordinates corresponding to the primed system, with the 
$z^\prime$-axis being the new polar axis. For a given 
$R(\alpha, \beta, \gamma )$ there is a coordinate transformation 
$(\theta, \phi ) \rightarrow (\theta ^\prime, \phi ^\prime )$, meaning that 
$(\theta ^\prime, \phi ^\prime )$ and 
$(\theta, \phi )$ are two pairs of coordinates referring to the same points on the sphere. In the primed system we can expand 
$\psi$ as
$$\begin{gather} \psi(\theta^\prime, \phi^\prime) = \sum_{l=0}^{\infty} \psi^\prime_{l}(\theta^\prime, \phi^\prime), \end{gather}$$
$$\begin{gather}\psi^\prime_{l}(\theta^\prime, \phi^\prime) = \sum_{p={-}l}^{l} b_{lp} Y_{l}^{p}(\theta^\prime, \phi^\prime) . \end{gather}$$A powerful result of group representation theory (see chapter 15 in Wigner Reference Wigner1959) states that
$$\begin{gather} Y_{l}^p (\theta^\prime, \phi^\prime) = \sum_{m={-}l}^{l} D_{pm}^{(l)} Y_{l}^{m} (\theta, \phi) , \end{gather}$$
$$\begin{gather}Y_{l}^m (\theta, \phi ) = \sum_{p={-}l}^{l} {D_{pm}^{(l)}}^{*} Y_{l}^{p} (\theta^\prime, \phi^\prime) , \end{gather}$$
where 
$D_{pm}^{(l)}(\alpha, \beta, \gamma )$ is a unitary matrix of order 
$2l+1$, the so called ‘Wigner D-matrix’, which is an irreducible representation of the 
$SO(3)$ group. Using (2.10) and orthogonality we find
\begin{align} b_{lp} &= \frac{1}{4{\rm \pi}} \int Y_{l}^{p}(\theta^\prime, \phi^\prime)^{*} \psi\, \text{d}\varOmega\nonumber\\ &= \frac{1}{4{\rm \pi}} \int \sum_{m={-}l}^{l} {D_{pm}^{(l)}}^* Y_{l}^{m} (\theta, \phi)^{*} \psi \,\text{d} \varOmega = \sum_{m={-}l}^{l} a_{lm}{D_{pm}^{(l)}}^{*} . \end{align}Substituting (2.12) into (2.9) and using (2.11) we obtain
\begin{equation} \psi_{l}^\prime(\theta^\prime, \phi^\prime) = \sum_{p={-}l}^{l} \sum_{m={-}l}^{l} a_{lm}{D_{pm}^{(l)}}^{*} Y_{l}^{p}(\theta^\prime, \phi^\prime) = \sum_{m={-}l}^{l} a_{lm} Y_{l}^{m} (\theta, \phi) = \psi_l(\theta, \phi) . \end{equation}
Since 
$(\theta, \phi )$ and 
$(\theta ^\prime, \phi ^\prime )$ are two pairs of coordinates referring to the same points on the sphere, (2.13) shows that 
$\psi _l^\prime$ is exactly the same function as 
$\psi _l$. In other words, for each 
$l$, 
$\psi _{l}$ is a coordinate invariant scalar field.
We will now let 
$\psi$ play the role of a streamfunction. A standard way of doing this would be to define the velocity field as 
$\boldsymbol {u} = - \boldsymbol {e}_r \times \boldsymbol {\nabla } \psi$, where 
${\boldsymbol {e}}_r$ is the radial unit vector. However, since 
$\boldsymbol {e}_r$ is not an intrinsic vector of the geometry such an approach would lead us into complications when differentiation is carried out on 
${\boldsymbol {u}}$. We will adopt another approach, also used by Eriksson & Nordmark (Reference Eriksson and Nordmark2020), which greatly simplifies differential operations, but comes with a price – the introduction of a notation which is rarely used in fluid mechanics. We let 
$\boldsymbol {\nabla }$ denote the surface gradient operator acting along the tangent of the sphere and 
$\star$ the so-called ‘Hodge star operator’ (see e.g. Bishop & Goldberg Reference Bishop and Goldberg1980; Schutz Reference Schutz1980) which produces the Hodge dual of the element that it is operating on. In our case the star is simply rotating a tangent vector on the sphere counter-clockwise a right angle about 
$\boldsymbol {e}_r$. Equipped with these tools we define the velocity field as
\begin{equation} \boldsymbol{u} \equiv{\star} \boldsymbol{\nabla} \psi , \end{equation}
ensuring that it is divergence free, 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = \boldsymbol {\nabla } \boldsymbol {\cdot } \star \boldsymbol {\nabla } \psi = 0$. With (2.14) we have the opposite sign convention for 
$\psi$ as compared with the standard given by 
$\boldsymbol {u} = - \boldsymbol {e}_r \times \boldsymbol {\nabla } \psi$. The operational rules including 
$\{ \boldsymbol {\nabla }, \star, \boldsymbol {\cdot } \}$ are given in Appendix A. The advantage of the notation becomes clear when the vorticity field is to be defined. Instead of acting on 
$\boldsymbol {u}$ with the curl operator and then project on the radial unit vector – a procedure that would bring us into trouble (see Krishnamurthy Reference Krishnamurthy2019) – we define the vorticity as the divergence of the rotated velocity, that is
\begin{equation} \omega \equiv \boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u} = \boldsymbol{\nabla} \boldsymbol{\cdot} \star{\star} \boldsymbol{\nabla} \psi ={-} \nabla^2 \psi , \end{equation}
where we have used 
$\star \star = -1$. Just as in the case of 
$\psi$, the sign convention of 
$\omega$, as defined in (2.15), will be the opposite compared with the standard. Readers who feel uncomfortable with this can restore the standard sign conventions of 
$\psi$ and 
$\omega$ by introducing a non-standard sign convention of the Hodge star by letting it perform a clockwise instead of a counter-clockwise rotation. That (2.15) is an appropriate definition of vorticity on the sphere is clear from the fact that it permits us to calculate the circulation along any closed curve on the sphere as the corresponding surface integral of 
$\omega$, in accordance with Stokes’ theorem.
The velocity and vorticity fields can be expanded as
\begin{equation} \boldsymbol{u} = \sum_{l=1}^{\infty}\boldsymbol{u}_l , \quad \omega = \sum_{l=1}^{\infty} \omega_l , \end{equation}where
\begin{equation} \boldsymbol{u}_l ={\star} \boldsymbol{\nabla} \psi_l , \quad \omega _l = \frac{l(l+1)}{a^2} \psi_l . \end{equation}
For each 
$l$, 
$\boldsymbol {u}_l$ is a coordinate invariant vector field and 
$\omega _l$ is a coordinate invariant scalar field. That a non-divergent velocity field on a sphere can be completely decomposed as in (2.16a,b) follows from the more general result that vector spherical harmonics form a complete set (Freeden & Gutting Reference Freeden and Gutting2008).
The mean kinetic energy and enstrophy over the sphere can now be expressed as
$$\begin{gather} \bar{E} = \frac{1}{4 {\rm \pi}} \int \frac{1}{2} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \,{\text{d}} \varOmega = \sum_{l=1}^{\infty} E(l) , \end{gather}$$
$$\begin{gather}\bar{{\mathcal{E}}} = \frac{1}{4 {\rm \pi}} \int \frac{1}{2} \omega^2 \,{\text{d}} \varOmega = \sum_{l=1}^{\infty} \frac{l(l+1)}{a^2} E(l) , \end{gather}$$where
\begin{equation} E(l) = \frac{l(l+1)}{2a^2} \sum_{m={-}l}^{l} a_{lm} a_{lm}^* = \frac{l(l+1)}{2a^2} \sum_{m={-}l}^{l} b_{lm} b_{lm}^* , \end{equation}
is the kinetic energy spectrum. From (2.20) it is clear that 
$E(l)$ is invariant under coordinate transformations between systems with different polar axes. It should be emphasised that the result is derived without making any assumption regarding the structure of the flow field. For example, the flow may consist of a single localised structure on the sphere, such as a vortex or an equatorial jet. Energy spectra calculated using spherical harmonic systems whose polar axes are oriented differently relative to the structure will all be equal.
2.1. Isotropy
Since the invariance of 
$E(l)$ with respect to rotations of the polar axis is a property that is satisfied for any flow field, isotropy must constrain the energy distribution even further. The statistics of an isotropic field should contain no information on the direction of the polar axis. Strictly speaking, isotropy therefore implies that the energy distribution among modes with different 
$m$ and the same 
$l$ also must be invariant under rotations. It is easy to see that there is no such distribution for a finite 
$l$ and that isotropy, in a strict sense, therefore can be fulfilled only in the limit 
$l \rightarrow \infty$. This is analogous to a flow field on a plane with periodic boundary conditions, for which isotropy can be only approximately fulfilled for finite wavenumbers. Likewise, for large but finite 
$l$, say 
$l > 100$, the distribution of energy among modes with the same 
$l$ and different 
$m$ may become approximately invariant under rotations, with increasing accuracy for larger and larger 
$l$. According to the analysis of the two-point covariance function by Boer (Reference Boer1983) isotropy implies that the variance of the expansion coefficient should be independent of 
$m$. Expressed in a different language kinetic energy should be equipartitioned among 
$m$-modes with the same 
$l$. Intuitively, it is quite clear that this must be the case. In a neighbourhood of a pole, we can make a flat Earth approximation and pass over to plane polar coordinates, as long as curvature effects are small. The 
$\theta$-coordinate of the spherical description will transform to a radial coordinate 
$r = |\theta |a$, while 
$\phi$ retains the role of an azimuth. Looking at the flow field from the pole, isotropy means statistical independence of the azimuthal coordinate for a fixed radial coordinate, which in spectral space translates to independence of 
$m$ for a fixed 
$l$. Since the polar axis can be chosen freely, this reasoning can be applied generally.
It is evident that isotropy cannot be even approximately fulfilled for low-order modes. Consider the 
$l = 1$ mode. With reference to a specific coordinate system, 
$(\boldsymbol {e}_x, \boldsymbol {e}_y, \boldsymbol {e}_z)$, we have
\begin{equation} \psi_1(\theta, \phi) = a_{1({-}1)} Y_1^{{-}1}(\theta, \phi) + a_{10} Y_1^0 (\theta, \phi) + a_{11} Y_1^1 (\theta, \phi) . \end{equation}
By (2.10) there exists another coordinate system, 
$(\boldsymbol {e}_{x^\prime }, \boldsymbol {e}_{y^\prime }, \boldsymbol {e}_{z^\prime })$, where
\begin{equation} \psi_1^\prime (\theta^\prime, \phi^\prime) = b_{10} Y_1^0 (\theta^\prime, \phi^\prime) = b_{10} \sqrt{3} \cos \theta^\prime , \end{equation}
with 
$b_{1(-1)} = b_{11} = 0$. The corresponding velocity field can be written as
\begin{equation} \boldsymbol{u}_1 = b_{10} \frac{\sqrt{3}}{a} \sin \theta^\prime \boldsymbol{e}_{\phi^\prime} , \end{equation}
which is a solid body rotation around the 
$z^\prime$-axis. It is quite clear that 
$\psi _1$ must be put to zero under the constraint of isotropy, and so must 
$\psi _2$ and 
$\psi _3$ and so forth up to a certain 
$l$ for which it may be assumed that isotropy is approximately fulfilled. The constraint of isotropy would thus force us to refrain from studying large scale dynamics and for this reason we deem it as too restrictive.
3. The spectral energy equation
The incompressible Navier–Stokes equations in a rotating frame of reference on a sphere with radius 
$a$ can be written as
\begin{equation} \frac{\partial \boldsymbol{u}}{\partial t} ={-} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} - \boldsymbol{\nabla} (p/\rho) - f \star \boldsymbol{u} + \nu \left (\Delta \boldsymbol{u} + \frac{\boldsymbol{u}}{a^2} \right) , \end{equation}
where 
$\boldsymbol {\nabla }$ is the gradient operator acting along the tangent of the sphere, 
$f = 2 {\boldsymbol \varOmega } \boldsymbol {\cdot } \boldsymbol {e}_r$ is the Coriolis parameter, 
$\boldsymbol{\varOmega}$ is the rotation vector and
\begin{equation} \Delta {\boldsymbol{u}} ={-} \star \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) + \frac{\boldsymbol{u}}{a^2} , \end{equation}
is the vector Laplace of 
$\boldsymbol {u}$ (see Appendix B). The centrifugal force is included in the pressure term. Compared with the standard form of the Navier–Stokes, there is an extra metric term, 
$\nu {\boldsymbol {u}}/a^2$, on the right-hand side of (3.1). In Appendix B we show how this term arises from the spherical geometry by using tensor analysis on a two-dimensional manifold. At the request of one of the reviewers we also show in Appendix C how the viscous term in (3.1) can be derived from the expression of the viscous term of the three-dimensional equation by considering a fluid layer within a thin spherical shell with stress free boundaries.
The equation for 
$E(l)$ can be derived by taking the scalar product between each term in (3.1) and 
$\boldsymbol {u}_l$ and averaging over the sphere. In this way, we obtain the time derivative of 
$E(l)$ from the left-hand side of (3.1),
\begin{equation} \frac{1}{4{\rm \pi}} \int \boldsymbol{u}_l \boldsymbol{\cdot} \sum_{n= 1}^{\infty} \frac{\partial \boldsymbol{u}_n}{\partial t} \,{\text{d}} \varOmega = \frac{1}{4{\rm \pi}} \int \boldsymbol{u}_l \boldsymbol{\cdot} \frac{\partial \boldsymbol{u}_l}{\partial t} \,{\text{d}} \varOmega = \frac{\partial E(l)}{\partial t} . \end{equation}The spectral term emanating from the nonlinear term in (3.1) we name the transfer term
\begin{equation} T(l) ={-}\frac{1}{4{\rm \pi}} \int (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}) \boldsymbol{\cdot} \boldsymbol{u}_l \,{\text{d}} \varOmega , \end{equation}in accordance with standard vocabulary in isotropic turbulence theory. Energy conservation is expressed as
\begin{equation} \sum_{l=1}^{\infty} T(l) ={-}\frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{u} \frac{1}{2} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u} \right ) \,{\text{d}} \varOmega = 0 . \end{equation}Enstrophy is conserved, because vorticity is conserved at a fluid particle. Using the rules in Appendix A we have
\begin{align} \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} &={-} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) \star \boldsymbol{u} + \tfrac{1}{2} \boldsymbol{\nabla} ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{u}) ,\end{align}
\begin{align} \boldsymbol{\nabla} \boldsymbol{\cdot} \star ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ) &={-} \boldsymbol{\nabla} \boldsymbol{\cdot} \star [ (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) \star \boldsymbol{u} ] \nonumber\\ &= \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) + (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega . \end{align}
Enstrophy conservation imposes a constraint on 
$T(l)$ which can be derived in the following way:
$$\begin{gather} T(l) ={-}\frac{1}{4{\rm \pi}} \int (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}) \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi_l \,{\text{d}} \varOmega = \frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} [{\star} (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}}) \psi_l ] \,{\text{d}} \varOmega \nonumber\\ - \frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \star ( \boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u} ) \psi_l \,{\text{d}} \varOmega ={-} \frac{a^2}{l(l+1)} \frac{1}{4{\rm \pi}} \int (\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \omega ) \omega_l \,{\text{d}} \varOmega , \end{gather}$$
where we have used (3.7), (2.17a,b) and the rules in Appendix A. Multiplying both sides by 
$l(l+1)$ and summing over 
$l$ we find
\begin{equation} \sum_{l=1}^{\infty} l(l+1) T(l) ={-} \frac{a^2}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{u} \frac{1}{2} \omega^2 \right ) \,{\text{d}} \varOmega = 0 , \end{equation}which is the expression for enstrophy conservation.
The contribution from the pressure term to the spectral energy equation is trivially zero,
\begin{equation} -\frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} (p/\rho) \boldsymbol{\cdot} \boldsymbol{u}_l \,{\text{d}} \varOmega ={-}\frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} (p \boldsymbol{u}_l/\rho ) \,{\text{d}} \varOmega = 0 . \end{equation}As found in previous investigations based on the vorticity equation (e.g. Tang & Orszag Reference Tang and Orszag1978) the Coriolis term does not make any contribution to the spectral energy equation. We confirm this by direct integration in the system whose polar axis coincides with the rotation axis,
\begin{align} &-\frac{1}{4{\rm \pi}} \int f \star \boldsymbol{u} \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi_l \,{\text{d}} \varOmega ={-}\frac{1}{4{\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{u} f \psi_l) \,{\text{d}} \varOmega \nonumber\\ &\qquad+\frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} f \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi \psi_l \,{\text{d}} \varOmega ={-}\frac{1}{4 {\rm \pi}} \int 2 | {\boldsymbol \varOmega} | \sin \theta \boldsymbol{e}_{\theta} \boldsymbol{\cdot} \left( -\frac{\partial_{\phi} \psi \boldsymbol{e}_\theta}{a \sin\theta} +\frac{\partial_{\theta} \psi \boldsymbol{e}_\phi }{a} \right ) \psi_l \,{\text{d}} \varOmega \nonumber\\ &\quad= \frac{1}{4{\rm \pi} a} \int 2 | {\boldsymbol \varOmega} | \sum_{n= 1}^\infty \frac{\partial \psi_n}{\partial \phi} \psi_l \,{\text{d}} \varOmega = 0 . \end{align}
That (3.11) is zero can be seen as follows. For 
$n = l$, use 
$\psi _l \partial _\phi \psi _l = 0.5 \partial _\phi (\psi _l^2)$ and integrate in 
$\phi$. For 
$n \ne l$, use 
$\partial _{\phi } Y_{n}^{m} = {\rm i}m Y_{n}^{m}$ and orthogonality of the spherical harmonics. That the Coriolis term does not contribute to the spectral energy equation can also be understood from the fact that it may be incorporated into the 
$l= 1$ mode, by transforming the equation back to a non-rotating frame of reference. The mean angular momentum with respect to the polar axis of a chosen spherical harmonic system can be calculated as
\begin{equation} \overline{L_z} = \frac{1}{4{\rm \pi}} \int_{0}^{2{\rm \pi}} \int_{0}^{\rm \pi} \sqrt{3} \sin^3 \theta a_{10} \,\mathrm{d} \theta \,\mathrm{d} \phi = \frac{2}{\sqrt{3}} a_{10} , \end{equation}
where the contributions from all other modes vanish by Rodrigues’ formula. Clearly, in the absence of drag, 
$\overline {L_z}$ is a conserved quantity and 
$a_{10}$ is therefore conserved and so is 
$\boldsymbol {u}_1$, since the mean angular momentum with respect to any axis is conserved.
The calculation of the contribution from the viscous term is straightforward. First,
\begin{align} \Delta \boldsymbol{u}_l &={-}\star \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star{\star} \boldsymbol{\nabla} \psi_l ) + \frac{\boldsymbol{u}_l}{a^2} ={\star} \boldsymbol{\nabla} (\nabla^2 \psi_l ) + \frac{\boldsymbol{u}_l}{a^2} \nonumber\\ &={-}\frac{l(l+1)}{a^2} \star \boldsymbol{\nabla} \psi_l + \frac{\boldsymbol{u}_l}{a^2} ={-} \frac{[l(l+1)-1]}{a^2} \boldsymbol{u}_l . \end{align}Then, using this expression we find that the contribution from the viscous term in (3.1) to the spectral energy equation is
\begin{equation} -\frac{2\nu [l(l+1)-2]}{a^2} E(l) . \end{equation}
Previous studies (e.g. Fjørtoft Reference Fjørtoft1953; Tang & Orszag Reference Tang and Orszag1978; Boer Reference Boer1983) gave an expression of the form 
$-2\nu l(l+1) E(l)/a^2$. That this cannot be completely accurate can be understood from the fact that there can be no viscous dissipation in the 
$l=1$ mode.
To summarise, we have derived the spectral energy equation on a sphere
\begin{equation} \frac{\partial E(l)}{\partial t} = T(l) -\frac{2\nu [l(l+1)-2]}{a^2}E(l) , \end{equation}where conservation of angular momentum, energy and enstrophy can be expressed as
$$\begin{gather} T(1) = 0 , \end{gather}$$
$$\begin{gather}\sum_{l=2}^{\infty} T(l) = 0 , \end{gather}$$
$$\begin{gather}\sum_{l= 2}^{\infty} l(l+1) T(l) = 0 . \end{gather}$$
The 
$l = 1$ mode represents a stationary solid body rotation. A linear Ekman drag force, 
$-\alpha \boldsymbol {u}$, added to the right-hand side of (3.1), modelling the interaction between the fluid and a solid ground, would give rise to a term 
$-2\alpha E(l)$ on the right-hand side of (3.15). Such a term would act on all modes, including the 
$l=1$ mode, attenuating the mean angular momentum of the fluid. The (3.15) is similar to the corresponding equation for the energy spectrum of two-dimensional homogeneous isotropic turbulence on a plane. An important difference is that (3.15) has been derived without invoking assumptions of homogeneity and isotropy.
4. Triad interactions
The transfer term can be written as a sum of interactions involving triad of modes
\begin{equation} T(l) ={-}\sum_{s =2}^{\infty} \sum_{n=2}^{\infty} \frac{1}{4{\rm \pi}} \int (\boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_n ) \boldsymbol{\cdot} \boldsymbol{u}_l \,{\text{d}} \varOmega . \end{equation}
Using a formula derived by Gaunt (Reference Gaunt1929), involving the integral of the product between three associated Legendre polynomials, Silberman (Reference Silberman1954) brilliantly proved that the nonlinear interaction within a triad of modes, 
$\{l,n,s\}$, is zero, unless
$$\begin{gather} l + n + s = {\text{odd integer}} , \end{gather}$$
$$\begin{gather}| l -n | < s < l+n , \end{gather}$$
where it should be noted that the triangle inequality (4.3) is strict. Below, it will be shown that the interaction is also zero if any two of 
$\{ l, n, s \}$ are identical. Fjørtoft (Reference Fjørtoft1953) derived relations involving energy exchanges within a triad under the assumption that energy exchanges take place only in a single triad, that is at initial time for initial conditions consisting of only three modes. It will now be shown that his relations hold generally, for energy exchanges within any triad of modes independently of all other exchanges.
We define the energy transfer from 
$n$ to 
$l$ in the interaction involving 
$s$ as the third mode, as
\begin{equation} T_{E}(s, n, l) ={-} \frac{1}{4{\rm \pi}} \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{n} ) \boldsymbol{\cdot} {\boldsymbol{u}}_{l } \, \mathrm{d} \varOmega . \end{equation}
The transfer is antisymmetric with respect to interchange of 
$n$ and 
$l$, that is
\begin{equation} T_{E}(s, n, l) ={-} T_{E}(s, l, n) , \end{equation}
meaning that the gain of 
$l$ is the loss of 
$n$. This follows from
\begin{align} 0&= \int \boldsymbol{\nabla} \boldsymbol{\cdot} ({\boldsymbol{u}}_{s} \boldsymbol{u}_{n} \boldsymbol{\cdot} {\boldsymbol{u}}_{l} ) \, \mathrm{d} \varOmega \nonumber\\ &= \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{n} ) \boldsymbol{\cdot} {\boldsymbol{u}}_{l} \, \mathrm{d} \varOmega + \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{l} ) \boldsymbol{\cdot} {\boldsymbol{u}}_{n} \, \mathrm{d} \varOmega , \end{align}
where we have used 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_{s} = 0$. In a similar way we define the transfer of enstrophy from 
$n$ to 
$l$, in the triad interaction involving 
$s$ as the third mode, as
\begin{align} &T_{{\mathcal{E}}}(s,n,l) ={-}\frac{1}{4{\rm \pi}} \int ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega} _{n} ) {\omega}_{l} \, \mathrm{d} \varOmega \nonumber\\ &\quad ={-}\frac{l(l+1)n(n+1)}{a^4} \frac{1}{4{\rm \pi}} \int \star \boldsymbol{\nabla} \psi_s \boldsymbol{\cdot} \boldsymbol{\nabla} \psi _{n} \psi_{l} \, \mathrm{d} \varOmega . \end{align}It is clear that
$$\begin{gather} T_{\mathcal{E}} (s,n,l) ={-} T_{\mathcal{E}}(s,l,n) , \end{gather}$$
$$\begin{gather}T_{\mathcal{E}} (s,s,l) = 0 , \end{gather}$$implying that the transfer is zero if two modes in a triad are identical.
We define the total gain of energy and enstrophy in unit time of mode 
$l$ at the expense of modes 
$n$ and 
$s$ as
$$\begin{gather} G(l, n, s) = T_E (s,n,l) + T_E (n,s,l) , \end{gather}$$
$$\begin{gather}H(l, n, s) = T_{\mathcal{E}} (s,n,l) + T_{\mathcal{E}} (n,s,l) . \end{gather}$$In the next step we show that
\begin{equation} H(l, n, s) = \frac{l(l+1)}{a^2} G(l,n,s) . \end{equation}First, a straightforward calculation gives
\begin{equation} H(l,n,s) ={-} \frac{1}{4{\rm \pi}} \frac{ l(l+1)}{a^2} \int ( \boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_n + \boldsymbol{u}_n \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_s ) \psi_l\, \mathrm{d} \varOmega , \end{equation}where we have used (2.17b). Then, a corresponding calculation gives
\begin{align} G(l,n,s) &={-}\frac{1}{4 {\rm \pi}} \int ( {\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} _{n} + \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{s} ) \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi_l \,\mathrm{d} \varOmega \nonumber\\ & = \frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} [{\star} ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} _{n} + \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{s} ) \psi_l ] -\boldsymbol{\nabla} \boldsymbol{\cdot} \star ({\boldsymbol{u}}_{s} \boldsymbol{\cdot} \boldsymbol{\nabla} {\boldsymbol{u}} _{n} + \boldsymbol{u}_{n} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{s} ) \psi_l \,\mathrm{d} \varOmega \nonumber\\ &= \frac{1}{4 {\rm \pi}} \int \boldsymbol{\nabla} \boldsymbol{\cdot} \star [ (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}_s) \star \boldsymbol{u}_n +(\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}_n) \star \boldsymbol{u}_s ] \psi_l \,\mathrm{d} \varOmega \nonumber\\ &={-} \frac{1}{4{\rm \pi}} \int ( \boldsymbol{u}_n \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_s + \boldsymbol{u}_s \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega}_n ) \psi_l \,\mathrm{d} \varOmega , \end{align}
where we have used the rules in Appendix A and 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_n = \boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u}_s = 0$. Now, (4.12) follows from (4.13) and (4.14). Energy and enstrophy conservation follow from (4.5) and (4.8),
$$\begin{gather} G(l,n,s) + G(n,s,l) + G(s,l,n) = 0 , \end{gather}$$
$$\begin{gather}H(l, n,s) + H(n,s,l) + H(s,l,n) = 0 . \end{gather}$$From (4.12) and (4.16) we also have
\begin{equation} l(l+1) G(l,n,s) + n(n+1) G(n,s,l) + s(s+1) G(s,l,n) = 0 . \end{equation}
For a triad of modes, 
$\{l, n, s \}$, such that 
$l < n < s$, the Fjørtoft (Reference Fjørtoft1953) relations now follow from (4.15) and (4.17),
\begin{equation} \left. \begin{aligned} \frac{G(l,n,s)}{G(s,l,n)} & = \frac{s(s+1) - n(n+1)}{n(n+1)-l(l+1)} > 0 , \\ \frac{G(n,s,l)}{G(l,n,s)} & = \frac{l(l+l) - s(s+1)}{s(s+1)-n(n+1)} < 0 , \\ \frac{G(s,l,n)}{G(n,s,l)} & = \frac{n(n-1) - l(l+1)}{l(l+1)-s(s+1)} < 0 . \end{aligned} \right\} \end{equation}
The relations (4.18) are analogous to the relations for energy exchanges within a triad, 
$\{k,p,q\}$, of Fourier modes, derived by Kraichnan (Reference Kraichnan1967, equation (2.7)). Apart from the different in notation – our 
$G(l,n,s)$ corresponds to 
$T(k,p,q)$ of Kraichnan – the difference is in the eigenvalue of the Laplace operator, which is 
$k^2$ in Fourier space and 
$l(l+1)$ in spherical harmonics space.
Nowhere in the derivation of (4.18) have we used that the elementary mode for each 
$l$ has been defined as a sum over the contributions from all 
$m$, as in (2.16a,b). The relations (4.18) would therefore hold also for a triad of single spherical harmonic modes. Such a triad can be specified as 
$\{ (l,m), (n,p), (s,q) \}$, corresponding to spherical harmonics 
$Y_l^m$, 
$Y_n^p$ and 
$Y_s^q$. However, the distribution of energy between modes with the same 
$l$ and different 
$m$ depends on the chosen coordinate system and for this reason it appears to be of limited interest that (4.18) actually hold even though no summation over 
$m$ is carried out. A related but more interesting result, is that for a chosen coordinate system there is no net enstrophy or energy exchange within a triad where all three modes have the same 
$l$. With a slight generalisation of the notation in (4.7) we can define the enstrophy gain of mode 
$(l,m)$ at the expense of mode 
$(n,p)$ in the triad interaction involving 
$(s, q)$ as the third mode as
$$\begin{gather} T_{{\mathcal{E}}}[(s,q), (n,p), (l,m) ]={-}\frac{1}{4{\rm \pi}} \int ({\boldsymbol{u}}_{sq} \boldsymbol{\cdot} \boldsymbol{\nabla} {\omega} _{np} ) {\omega}_{lm} \, \mathrm{d} \varOmega \nonumber\\ ={-}\frac{l(l+1)n(n+1)}{a^4} \frac{1}{4{\rm \pi}} \int \star \boldsymbol{\nabla} \psi_{sq} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi _{np} \psi_{lm} \, \mathrm{d} \varOmega . \end{gather}$$Clearly, we have
$$\begin{gather} T_{{\mathcal{E}}}[(s,q), (l,m), (n,p) ] ={-} T_{{\mathcal{E}}}[(s,q), (n,p), (l,m) ] , \end{gather}$$
$$\begin{gather}T_{{\mathcal{E}}}[(s,q), (s,p), (l,m) ] ={-} T_{{\mathcal{E}}}[(s,p), (s,q), (l,m) ] . \end{gather}$$
The enstrophy gain of mode 
$(l,m)$ at the expense of the two other modes is
\begin{equation} H[ (l,m), (s,q), (n,p) ] = T_{{\mathcal{E}}}[(s,q), (n,p), (l,m) ]+ T_{{\mathcal{E}}}[(n,p), (s,q), (l,m) ] . \end{equation}It now follows that
\begin{equation} H[(l,m),(l,q), (l,p)] = 0 , \end{equation}
which shows that there is no net enstrophy or energy exchange within a triad where all three modes have the same 
$l$.
5. Vanishing kinetic energy dissipation and cascade directions
Apart from the simultaneous conservation of energy and enstrophy by the nonlinear term there is a rigorous result which is encapsulating the difference between two- and three-dimensional turbulence, and that is the vanishing of kinetic energy dissipation in the limit 
$\nu \rightarrow 0$. The mean kinetic energy and enstrophy equations can be written as
$$\begin{gather} \frac{\mathrm{d} \bar{E}}{\mathrm{d} t } ={-}\epsilon , \end{gather}$$
$$\begin{gather}\frac{\mathrm{d} \bar{\mathcal{E}}}{\mathrm{d} t} ={-}\eta , \end{gather}$$where
$$\begin{gather} \epsilon \equiv \sum_{l = 2}^{\infty} \frac{2\nu ( l(l+1)-2)}{a^2} E(l) , \end{gather}$$
$$\begin{gather}\eta \equiv \sum_{l=2}\frac{2\nu ( l(l+1)-2)l(l+1)}{a^4} E(l) , \end{gather}$$are the mean kinetic dissipation and mean enstrophy dissipation, respectively. We now have
\begin{equation} \epsilon < 2 \nu \bar{\mathcal{E}} , \end{equation}implying that there can be no Richardson energy cascade in two-dimensional turbulence.
In order to say something more on energy and enstrophy transfer in wavenumber space some additional assumption is needed. A deep insight of Fjørtoft (Reference Fjørtoft1953) was that the inequalities (4.18) imply that there are two types of triad interactions. The first type is where the middle wavenumber 
$n$ acts as a source for each of the other two modes 
$l$ and 
$s$, and the second type is where it acts as a sink. He assumed that first type interactions are dominant and argued that as a consequence more energy will flow upscale than downscale. As pointed out by Merilees & Warn (Reference Merilees and Warn1975) his arguments were partly flawed, since they are based on the assumption that more energy is transferred to the smallest wavenumber mode 
$l$ than to the largest wavenumber mode 
$s$ in a triad interaction, which is not generally true. Kraichnan (Reference Kraichnan1967) developed the notion of two separated infinitely broad constant flux ranges – a downscale enstrophy flux range at large wavenumbers where the energy flux is zero and an upscale energy flux range at small wavenumbers where the enstrophy flux is zero. In both ranges, first type interactions are dominant, in accordance with the assumption of Fjørtoft (Reference Fjørtoft1953). Most subsequent studies (e.g. Eyink Reference Eyink1996; Tran & Shepherd Reference Tran and Shepherd2002; Gkioulekas & Tung Reference Gkioulekas and Tung2007) have used an assumption of stationarity in the analysis of cascade directions, which can be invoked only if some extra terms are added to (3.1). Usually, the added terms are forcing acting in a narrow wave number band and dissipation acting at small and large wave numbers. Of particular interest is an inequality derived by Gkioulekas & Tung (Reference Gkioulekas and Tung2007) under the assumption of stationarity and forcing in a wavenumber band 
$[k_1, \, k_2]$. In Fourier space, the ‘Danilov inequality’, named by Gkioulekas & Tung (Reference Gkioulekas and Tung2007) after the researcher who communicated it to them, reads
\begin{equation} \varPi_{{\mathcal{E}}}(k) - k^2 \varPi_{E}(k) > 0 , \end{equation}
where 
$\varPi _{{\mathcal {E}}}(k)$ and 
$\varPi _{E}(k)$ are the spectral enstrophy and energy fluxes respectively, and 
$k$ is the magnitude of the wavenumber vector. As shown by Gkioulekas & Tung (Reference Gkioulekas and Tung2007) (5.6) holds for a stationary flow at all wavenumbers outside 
$[k_1, \, k_2]$. It will now be shown that a corresponding inequality can be derived in spherical harmonic space under the assumption of dominance of first type of interactions for which it holds at all wavenumbers.
In spherical harmonic space, we define spectral energy and enstrophy fluxes as
$$\begin{gather} \varPi_{E}(k) \equiv \sum_{l = k}^{\infty} T(l) ={-} \sum_{l=2}^{k-1} T(l) , \end{gather}$$
$$\begin{gather}a^2\varPi_{\mathcal{E}}(k) \equiv \sum_{l = k}^{\infty} l(l+1)T(l) ={-} \sum_{l=2}^{k-1} l(l+1)T(l) . \end{gather}$$
There are some basic mathematical facts that deserve to be pointed out. First, it is evident that the fluxes cannot be completely constant in a range unless 
$T(l)$ is identically zero in the same range. The notion of a constant flux range is therefore an idealisation which should be adopted with caution. Second, the maxima and minima of 
$\varPi _{E}(k)$ coincide with the maxima and minima of 
$\varPi _{\mathcal {E}}(k)$. If we permit ourselves to think about 
$T(l)$ as a continuous function, the maxima and minima of the fluxes are the zeros of 
$T(l)$.
Following a similar derivation as carried out by Kraichnan (Reference Kraichnan1967) we can use the conservation and symmetry properties of 
$G(l,n,s)$ to rewrite the fluxes as
$$\begin{gather} \varPi_{E} (k) = \sum_{A_1(k,l,n,s)} G(s, n, l)- \sum_{A_2(k,l,n,s)} G(l,n,s) , \end{gather}$$
$$\begin{gather}a^2 \varPi_{\mathcal{E}}(k) = \sum_{A_1(k,l,n,s)} s(s+1)G(s, n, l) - \sum_{A_2(k,l,n,s)} l(l+1)G(l,n,s) , \end{gather}$$where
$$\begin{gather} A_1(k,l,n,s) = \{ \{l,n,s\} ; 1 < l < n < k \le s < l+n ; \; l+n+s = {\text{odd}} \} , \end{gather}$$
$$\begin{gather}A_2(k,l,n,s) = \{ \{l,n,s\} ; 1 < l < k \le n < s < l+n ; \; l+n+s = {\text{odd}} \} . \end{gather}$$
A detailed derivation of (5.9) is carried out in Appendix D and a derivation of (5.10) can be carried out in the same way. We note that in the summation over 
$A_1$, 
$G(s,n,l)$ is the gain of the largest mode 
$s$ and in the summation over 
$A_2$, 
$G(l, n, s)$ is the gain of the smallest wavenumber mode 
$l$. If the first type of interaction is dominant for all wavenumbers, the sums over 
$A_1$ in (5.9) and (5.10) give positive contributions to the fluxes while the sums over 
$A_2$ give negative contributions. For such a flow we will have
\begin{equation} a^2 \varPi_{\mathcal{E}}(k) - k(k+1)\varPi_E(k) \ge 0 , \end{equation}
for all 
$k$. The inequality (5.13) implies that 
$T(l)$ must be positive both in the limit of small and large 
$l$. Therefore, 
$T(l)$ must have at least two maxima, one at small and one at large 
$l$. Assume that there are only two maxima. Then, there is a middle range where 
$T(l)$ is negative and has a minimum. If we permit ourselves to think about the fluxes as continuous functions, they will both have a single zero. The enstrophy flux will be positive at the zero of the energy flux and the energy flux will be negative at the zero of the enstrophy flux. With narrow band forcing and a sufficiently large separation of length scales it can also be shown that the amount of energy which is flowing to the small wavenumber region is larger than the amount of energy flowing to the large wavenumber region and that the opposite is true for enstrophy (Eyink Reference Eyink1996; Gkioulekas & Tung Reference Gkioulekas and Tung2007).
If the second type of interaction were dominant the opposite inequality
\begin{equation} a^2 \varPi_{\mathcal{E}}(k) - k(k+1)\varPi_E(k) \le 0 , \end{equation}
would hold for all 
$k$. In this case, 
$T(l)$ must be negative in the limit of small as well as large wavenumbers and have at least two minima. Assume that there are just two minima. Then there is a middle wavenumber region where 
$T(l)$ is positive and has a maximum. Energy and enstrophy will flow from large and small wavenumbers into the middle region. Even though it is possible that such a transfer function can be generated for a short period of time from some initial conditions – for example an initial 
$E(l)$ having a broad minimum in a middle wavenumber range – it must be highly improbable. In particular, it seems extremely unlikely that enstrophy would flow from the region of very large wavenumbers where viscosity is acting into the middle wavenumber region. Two-dimensional turbulence cannot relax towards a stationary or quasi-stationary state as long as the transfer is dominated by the second type of interaction.
Thus, it seems that the Fjørtoft (Reference Fjørtoft1953) assumption of dominance of the first type of interaction may be the key to a deeper understanding of the problem of cascade directions. Whether the assumption can be rooted in some deeper principle is a very interesting question. Kraichnan (Reference Kraichnan1967) suggested that first type interactions represent ‘a statistical plausible spreading of the excitation in wave number: out of the middle wave number into the extremes’. Waleffe (Reference Waleffe1992) developed a stability argument suggesting that energy is generally flowing from the mode in which the gain/loss is largest. In two dimensions this is always the middle mode. In three dimensions, on the other hand, the dominant triad interaction is the one where the smallest wavenumber mode is loosing energy to the two larger ones. The argument is appealing, because the upscale energy transfer in two dimensions and the downscale transfer in three dimensions are derived from the same principle.
6. The double cascade theory expressed in spherical harmonic space
Kraichnan (Reference Kraichnan1967) formulated the equations of motions in an infinitely large quadratic domain with period boundary conditions in which case the discrete Fourier description passes over to an integral description. Formally, the description is valid only in the limit 
$Lk \rightarrow \infty$, where 
$L$ is the side length of the domain and 
$k$ is the magnitude of the wave vector. Under the assumption of isotropy he carried out a similarity analysis of the expressions for the energy fluxes, in which he made ingenious use of the conservation laws for triadic interactions. Based on this analysis he predicted that there are two possible similarity ranges – an upscale energy flux range in which the enstrophy flux is zero and a downscale enstrophy flux range where the energy flux is zero. It should be realised that such ranges are mathematical idealisations which only can be approached in the limit of large separation of length scales. In the constant energy flux range, the Fourier energy spectrum scales as (1.1) and in the constant enstrophy flux range it scales as (1.2), with a possible logarithmic correction.
From the similarity between the spectral energy equation (3.15), the Fjørtoft (Reference Fjørtoft1953) relations (4.18) and the flux expressions (5.9) and (5.10), on the one hand, and the corresponding expression in Fourier space, on the other hand, it is clear that there is a very close correspondence between the two descriptions in the limit of large wavenumbers. A heuristic demonstration of how they are related in this limit can be developed in the following way. Assume that the spherical harmonic spectrum peaks at 
$l \gg 1$ and the Fourier spectrum peaks at 
$k \gg 1/L$, and both spectra fall off exponentially at large wavenumbers. Then, for 
$n \ge 0$ we have
\begin{equation} \sum_{l=1}^{\infty}a^{{-}2n} [l(l+1)]^n E(l) \approx \sum_{l=1}^{\infty} (l/a)^{2n} E(l) \approx \int_{0}^{\infty} k^{2n} E_F(k) \,{\text{d}} k , \end{equation}
with increasing accuracy for increasing 
$n$. The reason why (6.1) must hold is that the two sides are expressions of the same physical quantity. For 
$n = 0$, (6.1) is the mean energy, for 
$n =1$ the mean enstrophy and for 
$n= 2$ the mean palenstrophy. With a broad band spectrum peaking at 
$l \gg 1$, we can approximate the sum in the middle term of (6.1) with an integral. There is only one straightforward way to do this, and that is
\begin{equation} \sum_{l=1}^{\infty} (l/a)^{2n} E(l) \approx \int_{0}^{\infty} k^{2n} E(ak) a\,\text{d} k \approx \int_{0}^{\infty} k^{2n} E_F(k) \,{\text{d}} k . \end{equation}
Since this must hold for all finite 
$n$, it can be argued that the spherical harmonic spectrum must converge to the Fourier spectrum in the limit of large 
$l$, where we must have
\begin{equation} E(l) = \frac{E_F(l/a)}{a} . \end{equation}
For the same reason it can be argued that 
$G(l,n,s)$ must converge to the corresponding triadic transfer function (see Kraichnan Reference Kraichnan1967), 
$T(k,p,q)$, in Fourier space in the limit of large wavenumbers, where we must have
\begin{equation} G(l,n,s) = \frac{T(l/a, n/a, s/a)}{a^3} . \end{equation}
The triangle inequality (4.3) of Silberman (Reference Silberman1954), which is a necessary condition for non-zero triad interactions in spherical harmonic space, corresponds exactly to the triangle inequality 
$| k-p| < q < k+p$, which is the condition in Fourier space. Thus, it is quite clear that the double cascade theory can equally well be formulated in spherical harmonic space as in Fourier space. If we use the expressions (4.18) in (5.9) and (5.10), make the approximation 
$l(l+1) \approx l^2$ and let the sums pass over to integrals, the similarity analysis of Kraichnan (Reference Kraichnan1967) can be repeated in every detail in spherical harmonic space. It may be objected that the approximation 
$l(l+1) \approx l^2$, which must be made in order for the similarity assumption to be consistent, indicates that the theory cannot be as accurately expressed in spherical harmonic space as in Fourier space. However, the assumptions of isotropy and similarity are only justified in the limit of large wavenumbers and in this limit the approximation 
$l(l+1) \approx l^2$ is also justified. The theory predicts that the spherical harmonic spectrum should scale as
\begin{equation} E(l) = \frac{C |\varPi_E|^{2/3} (l/a)^{{-}5/3}}{a} , \end{equation}in the constant energy flux range and
\begin{equation} E(l) = \frac{\kappa \varPi_{{\mathcal{E}}}^{2/3} (l/a)^{{-}3}}{a} , \end{equation}in the constant enstrophy flux range.
6.1. The two cascades are connected
In a flow with limited scale separation the two cascades are inevitable interconnected to some degree. Numerical simulations of the double cascade in Fourier space indicate that there is a stronger connection than may be anticipated from theory. For a given forcing wavenumber 
$k_f$, a necessary condition for an efficient upscale energy cascade is that there is a relatively broad enstrophy cascade range at 
$k > k_f$. If this condition is not met, as in the simulations of Xia et al. (Reference Xia, Wan, Chen and Eyink2009), more energy is actually flowing in the downscale direction than in the upscale direction. This is an indication of the importance of non-local triad interactions. Numerical investigations (Ohkitani Reference Ohkitani1990; Maltrud & Vallis Reference Maltrud and Vallis1991) show that the enstrophy flux emanates from highly non-local interactions where the smallest wavenumber typically lies outside the constant enstrophy flux range and the two other wavenumbers lie within the same range. This observation is consistent with the hypothesis of Polyakov (Reference Polyakov1993) that the enstrophy cascade is dominated by non-local triad interactions in which the smallest wavenumber mode is infrared. In spherical harmonic space this would mean that it is dominated by triads, 
$\{l, n, s\}$, such that 
$l \ll n < k < s < l+n$, appearing in the sum over 
$A_1$ in (5.10). Now, for each 
$k$ such that 
$l < k < n$ the same triad which is appearing with a positive sign in the expression for the enstrophy flux, is also appearing with a negative sign in the sum over 
$A_2$ in the expression (5.9) for the energy flux. Thus, the negative energy flux and the positive enstrophy flux may, to a large extent, emanate from the same set of triads, which would mean that the cascades are produced by the same physical mechanism. According to Boffetta & Ecke (Reference Boffetta and Ecke2012) the mechanism of the upscale energy cascade is ‘interaction of strain and vortices of different sizes’, while the downscale enstrophy cascade is the result of ‘large-scale straining of vortices through vortex-gradient stretching’. The two descriptions may refer to the same phenomenon.
6.2. The condensate may be reached in two different ways
As pointed out in the introduction, a number of numerical investigations show that the energy spectrum in the constant energy flux range can be considerably steeper than 
$k^{-5/3}$, sometimes as steep as 
$k^{-3}$, as in the simulation of Boure (Reference Boure1994). Common features of these simulations are that viscous damping is weak, both in the enstrophy and the energy cascade range, allowing for built up of non-local interactions, that a multitude of vortices is formed in the upscale cascade range and that the peak of the energy spectrum is moving very slowly towards small wavenumbers. Although the characteristic radius of the coherent vortices is much smaller than the side, 
$L$, of the computational domain, the state is in several ways similar to the end state in the simulation by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007) of a condensate, as also pointed out by Boffetta & Ecke (Reference Boffetta and Ecke2012). The condensate consists of two counter-rotating vortices of size 
${\sim }L$ with energy spectrum of the form 
$k^{-3}$. Of course, the latter state has a higher degree of regularity than the former. However, as shown by Burgess & Scott (Reference Burgess and Scott2017, Reference Burgess and Scott2018), the former has also a quite high degree of regularity, which reflects itself in the distribution of vortices with respect to radius, which can be predicted using scaling arguments.
These observations indicate that the condensate is not necessarily reached through a local cascade as predicted by Kraichnan (Reference Kraichnan1967) but may also be reached through a non-local cascade with a flow field dominated by coherent vortices and a steep energy spectrum, which slowly builds up in the direction from small to large wavenumbers, in contrast to the steep spectrum in the simulation by Chertkov et al. (Reference Chertkov, Connaughton, Kolokolov and Lebedev2007), which was built up in the opposite way. The high degree of regularity and the steep energy spectrum of the final state suggest that it is dominated by highly non-local interactions, which is also in line with the fact that the number of local triads decreases with decreasing wavenumber. Consider the nonlinear transfer to the lowest-order interacting mode, 
$l= 2$, in spherical harmonic space, which is also equal to the negative energy flux through the mode 
$l=3$,
\begin{equation} T(2) ={-} \varPi_{E}(3) = \sum_{n=3}^{\infty} G(2,n,n+1) . \end{equation}
If we use the threshold 
$n/2 \le 5$ for locality, there are eight local triads, while the rest are non-local. It would be surprising if local triad interactions were dominant.
7. Applications
In Earth's atmosphere, there is a relatively localised kinetic energy source at 
$l \in [10, 20]$ (see Augier & Lindborg Reference Augier and Lindborg2013), where available potential energy (Lorenz Reference Lorenz1955) is transferred to kinetic energy in a process involving baroclinic instability (see Vallis Reference Vallis2006). The subsequent energy transfer in spherical harmonic space is complicated by the fact that the velocity field is not strictly rotational (or non-divergent). Although the rotational field is totally dominant in terms of magnitude, the divergent field is transferring a non-negligible part of the energy in the downscale direction. On the other hand, the rotational field is transferring energy only in the upscale direction, which has been shown in a number of studies (e.g. Boer & Sheperd Reference Boer and Sheperd1983; Augier & Lindborg Reference Augier and Lindborg2013; Brune & Becker Reference Brune and Becker2013; Burgess, Erler & Shepherd Reference Burgess, Erler and Shepherd2013). Energy is transferred upscale at a rate 
${\sim }1 \, {\textrm {W}}\,{\textrm {m}}^{-2}$, integrated over the height of the troposphere. This is the mechanism by which the general circulation of the atmosphere is maintained, as also predicted by Lorenz (Reference Lorenz1955), based on an analysis using a mean-field/eddy-field decomposition. There can be little doubt that this transfer is caused by first type interactions as predicted by Fjørtoft (Reference Fjørtoft1953). As shown by analyses of Boer & Sheperd (Reference Boer and Sheperd1983) and Burgess et al. (Reference Burgess, Erler and Shepherd2013), the rotational field is also transferring a substantial amount of enstrophy downscale, which is the other signature of first type interactions.
The fact that the source is localised at such low wavenumbers, 
$l \sim 10-20$, implies that there is no room for a broad upscale energy cascade range to develop. As for the enstrophy transfer, the approximate 
$l^{-3}$-spectrum which is observed in a quite narrow range in analyses of general circulation models, may be a sign of something like an enstrophy cascade range, although the strong non-locality of enstrophy transfer makes such an interpretation questionable. In any case, three-dimensional effects are becoming important already at 
$l \sim 50$, so for Earth's atmosphere the theory of two-dimensional turbulence cannot be used to make much more predictions than determining the energy and enstrophy transfer directions for the large scale rotational field. As for large scale vortical structures, such as cyclones, they do not resemble the very long lived coherent vortices of two-dimensional turbulence, but are highly three-dimensional structures carrying helicity (Lilly Reference Lilly1986).
As for the atmosphere of the giant gas planets, Saturn and Jupiter, there is growing evidence from observations as well as simulations that the theory of two-dimensional turbulence is highly relevant. Evidence of an upscale energy cascade range has been extracted from analyses of Jupiter's tropospheric winds using Cassini near-infrared imaging (Choi & Showman Reference Choi and Showman2011; Galperin et al. Reference Galperin, Young, Sukoriansky, Dikovskaya, Read, Lancester and Armstrong2014; Young & Read Reference Young and Read2017). Since the dynamics is highly anisotropic and also involves Rossby waves, the assumption of isotropy which is often used in studies of two-dimensional turbulence, is not applicable. An extended concept of ‘zonostrophic turbulence’ accounting for Rossby waves and the formation of jets has been developed by Sukoriansky, Galperin & Dikovskaya (Reference Sukoriansky, Galperin and Dikovskaya2002). According to this theory, the zonal mean spectrum, which is the spectrum associated with all spherical harmonic modes with zonal wavenumber zero, should scale as 
$l^{-5}$, while the spectrum made up by all other modes conforms to the predictions of the Kraichnan (Reference Kraichnan1967) theory. Recently, it has been shown (Cabanes, Spiga & Young Reference Cabanes, Spiga and Young2020) that the dynamics can be successfully simulated and the spectra as well as the spectral transfers show all signs of two-dimensional turbulence, including both an upscale energy cascade and a downscale enstrophy cascade. However, the prediction of a 
$l^{-5/3}$-spectrum in the upscale cascade range is only partially fulfilled. In the low wavenumber end of the spectrum, there is indeed a range with 
$l^{-5/3}$-scaling, but for somewhat larger wavenumbers there is a local maximum (see figure 6 in Cabanes et al. Reference Cabanes, Spiga and Young2020) of the same type as seen in simulations where a multitude of vortices are formed in the upscale cascade range. It remains to be investigated if the 
$l^{-5/3}$-spectrum can be reproduced in different set-ups and to what degree the observed upscale energy cascade is consistent with the theory of Kraichnan (Reference Kraichnan1967) of a local cascade.
8. Conclusion
As mentioned in the introduction there are several unresolved issues of the double cascade theory which are of general interest from a broad perspective of physics. In particular, this is true for all problems connected to the upscale energy cascade – locality of triad interactions in spectral space, formation of coherent vortices and formation of a condensate when the cascade reaches the scale of the confinement. The reason why these problems are so interesting from a fundamental point of view, is that the upscale cascade somehow seems to defy the trend from more to less order which is ubiquitous to systems with many degrees of freedom, as expressed in the second law of thermodynamics but also in the Kolmogorov theory of three-dimensional turbulence. Isotropisation is an essential element of such a trend. As energy is successively transferred from large to small scales in the three-dimensional energy cascade, information is lost and motions become more isotropic and less coherent with decreasing scale. In the upscale energy cascade there is an opposite trend. If a random force generates small scale incoherent motions with no preferred direction, the upscale cascade will transform these into large scale coherent and anisotropic motions that will finally condensate into a globally ordered state.
As we have shown, the spherical harmonic spectral energy equation has a similar form as the corresponding spectral energy equation in Fourier space, with the important difference that it can be derived without invoking any assumption of isotropy. It describes the nonlinear transfer between all scales of motion, represented by a discrete set of integer modes, 
$l = 1, 2, 3, \ldots$, with associated vector fields 
$\boldsymbol {u}_l$, which are coordinate invariant. Needless to say, no artificial boundary conditions have to be implemented on the sphere. These features make the spherical geometry the ideal test ground for exploration of two-dimensional turbulence by means of simulations. It is the hope of the authors that the present work will stimulate such research.
Acknowledgements
The authors would like to thank three anonymous reviewers for positive and constructive criticism.
Declaration of interest
The authors report no conflict of interest.
Appendix A. Operational rules
In this appendix we list the rules that permit us to operate with 
$\{\boldsymbol {\nabla }, \star, \boldsymbol {\cdot } \}$ on scalars and vectors on a sphere. We let 
$\psi$ be a scalar field and 
$\boldsymbol {A}$ and 
$\boldsymbol {B}$ two vector fields defined on a sphere with radius 
$a$.
(i) Gradient of a scalar, divergence of a vector and divergence of a rotated vector
(A 1)
$$\begin{gather} \boldsymbol{\nabla} \psi , \end{gather}$$(A 2)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{A} , \end{gather}$$(A 3)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{A} . \end{gather}$$(ii) Divergence of a rotated gradient
(A 4)\begin{equation} \boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{\nabla} \psi = 0 . \end{equation}(iii) Scalar and vector Laplace (see Appendix B)
(A 5)$$\begin{gather} \Delta \psi = \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\nabla} \psi , \end{gather}$$(A 6)$$\begin{gather}\Delta \boldsymbol{A} = \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{A}) -{\star} \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{A}) + \frac{\boldsymbol{A}}{a^2} . \end{gather}$$(iv) Commutation and distribution
(A 7)$$\begin{gather} \star{\star} \boldsymbol{A} ={-} \boldsymbol{A} , \end{gather}$$(A 8)$$\begin{gather}\star ( \psi \boldsymbol{A}) = \psi \star \boldsymbol{A} , \end{gather}$$(A 9)$$\begin{gather}\boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{B} = \boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{A} , \end{gather}$$(A 10)$$\begin{gather}\star \boldsymbol{A} \boldsymbol{\cdot} \star \boldsymbol{B} = \boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{B} , \end{gather}$$(A 11)$$\begin{gather}\star \boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{B} ={-} \boldsymbol{A} \boldsymbol{\cdot} \star \boldsymbol{B} , \end{gather}$$(A 12)$$\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} (\psi \boldsymbol{A}) = \boldsymbol{\nabla} \psi \boldsymbol{\cdot} \boldsymbol{A} + \psi \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{A} , \end{gather}$$(A 13)$$\begin{gather}\boldsymbol{\nabla} ( \boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{B}) = (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{B}) \star \boldsymbol{A} + (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{A}) \star \boldsymbol{B} + (\boldsymbol{A} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{B} + (\boldsymbol{B} \boldsymbol{\cdot} \boldsymbol{\nabla} ) \boldsymbol{A} . \end{gather}$$
Appendix B. The viscous term and the vector Laplacian
Here, we will establish the form of the Navier–Stokes equations (3.1) and the expression for the vector Laplacian (A6). We will do this using tensor notation, and translate to or from vector notation when needed. Thus we introduce any suitable system of two coordinates on the sphere, with metric tensor 
$g_{ab}$.
B.1. The viscous term
The equation of motion for an incompressible fluid of constant density 
$\rho$ is
\begin{equation} \rho\left(\frac{\partial u^a}{\partial t}+u^bu^a_{;b}\right)=\tau^{ab}_{;b}+F^a, \end{equation}
where 
$F$ is an external forcing, and for a linear viscous fluid the symmetric stress tensor is
\begin{equation} \tau^{ab}=2\mu e^{ab}-g^{ab}p, \end{equation}
where 
$\mu$ is the constant dynamic viscosity, 
$p$ is the pressure and the symmetric strain rate tensor is
\begin{equation} e^{ab}=\frac{1}{2}\left(g^{bc}u^a_{;c}+g^{ac}u^b_{;c}\right). \end{equation}In particular, the divergence of the viscous stress is
\begin{equation} \mu\left(g^{bc}u^a_{;cb}+g^{ac}u^b_{;cb}\right). \end{equation}We recognise the first term in the parenthesis as the vector Laplacian. The second term can be simplified by changing the order of the covariant derivatives, but since the sphere is a curved manifold, covariant derivatives do not commute, so the change of order involves the curvature tensor
\begin{equation} u^b_{;cb}=\left(u^b_{;b}\right)_{;c}-R^b_{\ dcb}u^d={-}R_{dc}u^d, \end{equation}
utilising that the fluid is incompressible so 
$u^b_{;b}=0$. For a two-dimensional manifold, the Ricci tensor can always be expressed in terms of the curvature scalar and the metric as
\begin{equation} R_{dc}=\frac{R}{2}g_{dc}, \end{equation}
and for a sphere of radius 
$a$, the curvature scalar is
\begin{equation} R={-}\frac{2}{a^2}, \end{equation}so the divergence of the viscous stress is
\begin{equation} \mu\left(g^{bc}u^a_{;cb}+g^{ac}\frac{1}{a^2}g_{dc}u^d\right)=\mu\left(g^{cb}u^a_{;cb}+\frac{u^a}{a^2}\right) . \end{equation}
Dividing by 
$\rho$, using 
$\mu /\rho = \nu$ and shifting to vector notation we have
\begin{equation} \nu\left(\Delta\boldsymbol{u}+\frac{\boldsymbol{u}}{a^2}\right), \end{equation}which is the viscous term found in (3.1).
B.2. The vector Laplacian
We will show that
\begin{equation} \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{A}) -{\star} \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{A})={-}\frac{\boldsymbol{A}}{a^2}+\Delta \boldsymbol{A} , \end{equation}
by translating the expression first into tensor notation and then back. This is straightforward, except maybe for the rotation operation on an oriented two-dimensional manifold. The rotation operation on contravariant vector components 
$A^a$ can be written in two equivalent ways using the Levi-Civita symbols
\begin{equation} \frac{1}{\sqrt{g}}\epsilon^{ba}g_{bc}A^c=\sqrt{g}g^{ab}\epsilon_{cb}A^c. \end{equation}Consider the negative of the second term of the left-hand side of (B10). In tensor notation it can be written as
$$\begin{gather} \frac{1}{\sqrt{g}}\epsilon^{ab}g_{bc}g^{cd}\left(\left(\sqrt{g}g^{ef}\epsilon_{fg}A^g\right)_{;e}\right)_{;d}=\epsilon^{ab}\epsilon_{fg}g^{ef}A^g_{;eb}\nonumber\\ =\left(\delta^a_f\delta^b_g-\delta^a_g\delta^b_f\right)g^{ef}A^g_{;eb}=g^{ea}A^b_{;eb}-g^{eb}A^a_{;eb}, \end{gather}$$where we recognise the second term as the negative vector Laplacian. The whole left-hand side of (B10) is then
$$\begin{gather} g^{ae}\left(A^b_{;b}\right)_{;e}-g^{ea}A^b_{;eb}+g^{eb}A^a_{;eb}=g^{ae}\left(A^b_{;be}-A^b_{;eb}\right)+g^{eb}A^a_{;eb}\nonumber\\ =g^{ae}R_{ec}A^c+g^{eb}A^a_{;eb}={-}\frac{A^a}{a^2}+g^{eb}A^a_{;eb}, \end{gather}$$again using the curvature of the sphere. Translating back to vector notation, we have the right-hand side of (B10).
Appendix C. Derivation of (3.1) from the three-dimensional Navier–Stokes equation
At the request of one of the reviewers we here derive the expression of the viscous term found in the two-dimensional incompressible Navier–Stokes equation (3.1) from the three-dimensional incompressible equation. To do this, we consider a thin layer of fluid within a spherical shell with boundaries at 
$r = a$ and 
$r = a(1+ \epsilon )$ and write the velocity field within the layer as 
${\boldsymbol v}(r,\theta, \phi )$. Boundary conditions of no radial velocity, 
$v_r=0$, and no shear stress, 
$\tau _{\theta r}=\tau _{\phi r}=0$, are applied on both boundaries. For a given tangential length scale the radial velocity must then tend to zero throughout the layer as 
$\epsilon \rightarrow 0$. Using 
$v_r = 0$, the radial shear stresses associated with 
$\boldsymbol {v}$ are (see e.g. Kundu, Cohen & Dowling Reference Kundu, Cohen and Dowling2016, Appendix B)
$$\begin{gather} \tau_{r \theta} = \tau_{\theta r} ={\mu} \left ( \frac{\partial v_{\theta}}{\partial r} - \frac{v_{\theta}}{r} \right )=\mu r\frac{\partial\left(v_\theta/r\right)}{\partial r} , \end{gather}$$
$$\begin{gather}\tau_{r \phi} = \tau_{\phi r} = {\mu} \left ( \frac{\partial v_{\phi}}{\partial r} - \frac{v_{\phi}}{r} \right )= \mu r\frac{\partial\left(v_\phi/r\right)}{\partial r} . \end{gather}$$Applying the stress free boundary conditions, and using the same reasoning as above, in the thin limit the condition
\begin{equation} \frac{\partial\left(v_\theta/r\right)}{\partial r}=\frac{\partial\left(v_\phi/r\right)}{\partial r}=0 \end{equation}must hold throughout the layer. This means we can write the velocity field within the layer as
\begin{equation} {\boldsymbol v}(r,\theta, \phi) = (r/a) \boldsymbol{u}(\theta, \phi) , \end{equation}
where 
$\boldsymbol {u}(\theta, \phi )$ is the purely tangential velocity at 
$r = a$. Using the textbook polar coordinate expression of the viscous term in the three-dimensional incompressible Navier–Stokes (see e.g. Kundu et al. Reference Kundu, Cohen and Dowling2016, Appendix B) we obtain
\begin{align} \nu \Delta^{(3)} {\boldsymbol v} |_{r=a} & = \nu \left (\nabla^2 u_\theta + \left ( 2 -\frac{1} {\sin^2\theta} \right ) \frac{u_\theta}{a^2} - \frac{2 \cos \theta}{a^2 \sin^2 \theta} \frac{\partial u_\phi}{\partial \phi} \right ) \boldsymbol{e}_\theta\nonumber\\ &\quad + \nu \left ( \nabla^2 u_\phi + \left ( 2 -\frac{1} {\sin^2\theta} \right ) \frac{u_\phi}{a^2} + \frac{2 \cos \theta}{a^2 \sin^2 \theta} \frac{\partial u_\theta}{\partial \phi} \right ) \boldsymbol{e}_\phi , \end{align}
where 
$\Delta ^{(3)}$ is the three-dimensional vector Laplace operator and
\begin{equation} \nabla^2 u_\theta = \frac{1}{a^2 \sin \theta} \frac{\partial }{\partial \theta} \left ( \sin \theta \frac{\partial u_{\theta}}{\partial \theta} \right ) + \frac{1}{a^2 \sin^2 \theta} \frac{\partial ^2 u_{\theta}}{\partial \phi^2} , \end{equation}
with a corresponding expression for 
$\nabla ^2 u_{\phi }$. The radial component of (C5) vanishes by 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {u} = 0$. The term
\begin{equation} 2 \nu \frac{\boldsymbol{u}} {a^2} = \frac{2 \nu}{a^2} (u_\theta \boldsymbol{e}_\theta + u_\phi \boldsymbol{e}_\phi) , \end{equation}
emanates from the radial derivative of 
$\boldsymbol v(r, \theta, \phi )$ in the three-dimensional Laplace operator. If this term is not taken into account (as in Fjørtoft Reference Fjørtoft1953; Tang & Orszag Reference Tang and Orszag1978; Boer Reference Boer1983) the viscous stress tensor will include a component acting as an external drag on the sphere, implying that angular momentum is not conserved and that a solid body rotation will decay.
A straightforward but somewhat lengthy calculation of the polar coordinate expression of the viscous term in (3.1) gives
\begin{align} \nu \left ( \Delta \boldsymbol{u} + \frac{\boldsymbol{u}}{a^2} \right )& = \nu \left ( \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}) -{\star} \boldsymbol{\nabla} (\boldsymbol{\nabla} \boldsymbol{\cdot} \star \boldsymbol{u}) + \frac{2\boldsymbol{u}}{a^2} \right ) \nonumber\\ & = \nu \left ( \boldsymbol{e}_\theta \frac{1}{a} \frac{\partial }{\partial \theta} + \boldsymbol{e}_{\phi} \frac{1}{a \sin \theta} \frac{\partial }{\partial \phi} \right) \left ( \frac{1}{a \sin \theta} \frac{\partial }{\partial \theta} ( \sin \theta u_\theta) + \frac{1}{a \sin \theta} \frac{\partial u_\phi}{\partial \phi} \right ) \nonumber\\ &\quad - \nu \left ( -\boldsymbol{e}_\theta \frac{1}{a\sin \theta} \frac{\partial }{\partial \phi} + \boldsymbol{e}_{\phi} \frac{1}{a } \frac{\partial }{\partial \theta} \right) \left ( -\frac{1}{a \sin \theta} \frac{\partial }{\partial \theta} ( \sin \theta u_\phi) + \frac{1}{a \sin \theta} \frac{\partial u_\theta}{\partial \phi} \right ) \nonumber\\ &\quad + \nu \frac{2} {a^2} ( u_{\theta} \boldsymbol{e}_\theta + u_\phi \boldsymbol{e}_\phi ) \nonumber\\ & = \nu \left ( \nabla^2 u_\theta + \left ( 2 -\frac{1} {\sin^2\theta} \right ) \frac{u_\theta}{a^2} - \frac{2 \cos \theta}{a^2 \sin^2 \theta} \frac{\partial u_\phi}{\partial \phi} \right ) \boldsymbol{e}_\theta \nonumber\\ &\quad + \nu \left ( \nabla^2 u_\phi + \left ( 2 -\frac{1} {\sin^2\theta} \right ) \frac{u_\phi}{a^2} + \frac{2 \cos \theta}{a^2 \sin^2 \theta} \frac{\partial u_\theta}{\partial \phi} \right ) \boldsymbol{e}_\phi , \end{align}
where we have skipped some obvious steps where we have taken derivatives of 
$\sin \theta$, 
$\cos \theta$ and 
$1/\sin \theta$ and noticed that terms including mixed derivatives cancel. Comparing (C5) and (C8) we can thus conclude
\begin{equation} \nu \Delta^{(3)} {\boldsymbol v} |_{r=a} = \nu \left ( \Delta \boldsymbol{u} + \frac{\boldsymbol{u}}{a^2} \right ) , \end{equation}showing that the viscous term of two-dimensional equation (3.1) can be derived from the viscous term of the three-dimensional equation.
Appendix D. Derivation of energy flux relation
In this appendix we derive the energy flux relation (5.9). Using the symmetry 
$G(l, n, s) = G(l, s, n)$ and (4.15) we can write the energy flux (5.7) as
$$\begin{gather} \varPi_{E}(k) = \frac{1}{2} \sum_{l=k}^{\infty} \sum_{n = 2}^{\infty} \sum_{s=2}^{\infty} G(l,n,s) = \frac{1}{2} \sum_{l=k}^{\infty} \left ( \sum_{n = 2}^{k-1} + \sum_{n=k}^{\infty} \right) \left ( \sum_{s=2}^{k-1} + \sum_{s=k}^{\infty} \right) G(l,n,s) \nonumber\\ = \frac{1}{2} \sum_{l=k}^{\infty} \sum_{n = 2}^{k-1}\sum_{s=2}^{k-1} G(l,n,s) + \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(l,n,s) . \end{gather}$$Now we have
$$\begin{gather} \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(l,n,s) ={-} \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(s,n,l) -\sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(n,s,l)\nonumber\\ ={-} \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(l,n,s) -\sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(n,s,l) , \end{gather}$$so that
\begin{equation} \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(l,n,s) ={-}\frac{1}{2} \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(n,s,l) . \end{equation}Substituting (D3) into (D1) and relabelling the names of the indices, we obtain
$$\begin{gather} \varPi_{E}(k) = \frac{1}{2} \sum_{l=k}^{\infty} \sum_{n = 2}^{k-1}\sum_{s=2}^{k-1} G(l,n,s) -\frac{1}{2} \sum_{l=k}^{\infty} \sum_{n= 2}^{k-1} \sum_{s=k}^{\infty} G(n,s,l) \nonumber\\ = \sum_{s=k}^{\infty} \sum_{n = 2}^{k-1}\sum_{l=2}^{n-1} G(s,n,l) - \sum_{l= 2}^{k-1}\sum_{s=k}^{\infty} \sum_{n=k}^{s-1} G(l,n,s) . \end{gather}$$
All terms in the two sums which do not satisfy 
$n < s < l + n$ and 
$l+n+s = {\textrm {odd}}$, are zero. We conclude that the energy flux can be written as in (5.9). The expression (5.10) for the enstrophy flux can be derived in the same way.
Open access
