Statistical Equilibrium Principles in 2D Fluid Flow: From Geophysical Fluids to the Solar Tachocline

An overview is presented of several diverse branches of work in the area of effectively 2D fluid equilibria which have in common that they are constrained by an infinite number of conservation laws. Broad concepts, and the enormous variety of physical phenomena that can be explored, are highlighted. These span, roughly in order of increasing complexity, Euler flow, nonlinear Rossby waves, 3D axisymmetric flow, shallow water dynamics, and 2D magnetohydrodynamics. The classical field theories describing these systems bear some resemblance to perhaps more familiar fluctuating membrane and continuous spin models, but the fluid physics drives these models into unconventional regimes exhibiting large scale jet and eddy structures. From a dynamical point of view these structures are the end result of various conserved variable forward and inverse cascades. The resulting balance between large scale structure and small scale fluctuations is controlled by the competition between energy and entropy in the system free energy, in turn highly tunable through setting the values of the conserved integrals. Although the statistical mechanical description of such systems is fully self-consistent, with remarkable mathematical structure and diversity of solutions, great care must be taken because the underlying assumptions, especially ergodicity, can be violated or at minimum lead to exceedingly long equilibration times. Generalization of the theory to include weak driving and dissipation (e.g., non-equilibrium statistical mechanics and associated linear response formalism) could provide additional insights, but has yet to be properly explored.


Introduction
Remarkable progress has been made over the past 30 years or so applying rigorous statistical equilibrium principles to classical fluid systems with increasing degrees of complexity [1]. The essential idea is that a freely decaying, strongly turbulent initial condition at late time is often observed to relax into a macroscopically smooth steady state (illustrated below in Figure 3). These ideas are especially interesting in two dimensions where inverse cascades can generate nontrivial macroscopic features, such as system-spanning eddies or jets, from purely small scale, but significantly nonlinear fluctuations. Moreover, additional strong constraints, that forbid 2D eddies from "turning over" and effectively self-canceling, lead to an infinite number of additional conserved integrals of the motion, known as Casimirs. Acting together, all of these lead to a similarly infinite number of possible late-time flow geometries. These are exemplified, e.g., by Jupiter's Great Red Spot, gas giant latitudinal band structure, polar vortices [2][3][4], and other planetary flows.
Some, but by no means all, of these near-steady state long-lived, structures might be considered as weakly driven, balanced by weak dissipation. It then becomes interesting to seek quantitative and qualitative insights using models in an idealized zero driving, zero dissipation limit. For realistic comparisons, these models may additionally require nontrivial multilayer vertical structure. Here we consider only the simplest models with dr dr ω(r) q(r) T G h (r, r ) ω(r ) q(r ) compression field q, surface height h, planar coordinate r = (x, y)] Continuous spin, long-range interacting Ising-type field ω, with spin + dr 1 2 gh(r) 2 − h(r)µ[ω(r)/h(r)] weighting µ(ω/h), tensor-coupled nonlinearly to Gaussian fields q, h 2D magnetohydrodynamics [Section 9; stream function ψ, K[A, ψ] = dr 1 2 |∇A(r)| 2 + 1 2 |∇ψ(r)| 2 magnetic vector potential A, planar section r = (x, y) orthogonal to electric current density J = −∇ 2 A alongẑ] −μ [A(r)]∇A(r) · ∇ψ(r) + µ[A(r)] Model is equivalent to that of a pair of gradient-coupled elastic membranes, external confining and coupling potentials µ(A),μ (A) Predictions for the late time equilibrium state, assuming that it is reached, are based only on certain macroscopic features of the initial condition, namely the values of the conserved integrals, including total energy, linear or angular momentum, and the Casimirs. Although insensitive to the details of the turbulent decay that gives rise to these states, such predictions, beyond their intrinsic interest, could provide useful consistency checks on results from late time direct numerical simulations. Conversely, lack of consistency, if indeed robustly borne out by the numerics, could point to existence interesting equilibration barriers and metastable behaviors. There is already significant evidence that such barriers are much more common in such highly constrained 2D flows than in, e.g., conventional particle systems, through a variety of mechanisms [19][20][21][22][23][24][25][26][27][28][29].
Finally, it is worth mentioning that characterization of stable steady state flows may be mathematically motivated by considerations other than convergence to equilibrium. Thus, for example, "maximally mixed" steady state flows that extremalize a single Casimir (at fixed energy) are investigated in Ref. [30], motivated by earlier ideas in Ref. [31]. The corresponding variational functionals are rather different from those emerging from statistical mechanics which simultaneously control the values of all Casimirs. It could well be that these maximally mixed flows are also equilbrium states (for some to-be-determined values of the thermodynamic parameters) but this possible connection has not yet been investigated.

Outline
The remainder of this paper is summarized as follows. We begin by presenting a fairly detailed derivation of the statistical equilibrium theory for the simplest possible model, the 2D Euler equation, which is fully described by the scalar vorticity. The ingredients of this theory follow a logical chain that is repeated, or extended as necessary, for the more complicated systems. In Section 2, the equations of motion are introduced and their reduction to the vorticity field dynamics demonstrated. The usual energy and momentum conservation laws are exhibited, followed by the Casimir constraints.
General equilibrium concepts are introduced in Section 3 in terms of invariant (steady state) measures over the phase space of all vorticity configurations. Identifying such measures relies on the Liouville theorem, which establishes a type of phase space incompressibility condition. Once proven, the allowed measures are constructed from the fluid conserved integrals themselves, and the exact choice corresponds to what is known as a statistical ensemble. The thermodynamic entropy, free energy, etc., follow from the logarithm of the global phase space integral (partition function) in the usual way. The grand canonical ensemble for the Euler equation is introduced as providing the most convenient mathematical framework.
The general statistical formalism is applied to the 2D Euler equation in Section 4. Perhaps surprisingly, given the infinite number of constraints, the system free energy may actually be derived exactly as an explicit variational equation-the long range Coulomb-like vortex interactions enable an exact mean field-type approximation [11][12][13]. The minima describe the various possible equilibrium states, whose large scale flow pattern varies with the specified conserved integral values. Critically, the Casimir constraints permit both positive and negative temperature equilibria, with the latter encouraging compact eddy structures reminiscent of Jupiter's Red Spot. There is again a very interesting competition between energy and entropy that controls the amplitude and size of such structures. Simple two-level system models are introduced that allow convenient exploration of these phenomena.
A brief discussion of some of the limitations of the statistical equilibrium hypothesis is presented in Section 5. Vortex mixing dynamics in 2D is clearly far more constrained than particle dynamics underlying conventional systems (though microscale viscosity, neglected here, in a sense bridges the two regimes). It should therefore not be too surprising that significant barriers to equilibration can occur [1]. Some of these barriers can actually be understood as local rather than global minima of the free energy functional. Examples include separated compact eddies that orbit each other, failing to merge (as would be entropically favored) above a critical separation [21]. Detailed numerical simulations show evidence for different levels of equilibration in different spatial regions, depending on the strength of local mixing dynamics [22]. Others are somewhat more mysterious: equilibration on the surface of a sphere (rather than in a flat bounded domain) is found to fail much more catastrophically, with a macroscopically fluctuating chaotic vorticity field surviving for all achievable computation times [26].
In Section 6, we discuss the most straightforward generalization of the Euler results to a more general class of single scalar field systems whose canonical structure automatically ensures an infinite set of Casimirs. Under reasonable conditions, the mean field approximation is again exact, and the free energy functional emerges from a Legendre transformation of the energy. An important example is the quasigeostrophic (QG) equation, a scalar field approximation to the shallow water Equations [32,33]. This system also has an additional approximate adiabatic invariant [19] that is completely separate from the standard conservation laws, and provides another possible equilibration barrier example.
In Section 7, we consider 3D axisymmetric flow in which azimuthal symmetry is imposed on flows confined to a cylinder (Taylor-Couette geometry). The equations of motion now reduce to a coupled pair of scalar equations describing coupled toroidal and poloidal flow, with only the former experiencing the Casimir constraints [34]. However it is the poloidal flow, within each range-height slice, that is most directly analogous to the Euler equation vorticity. The fact that it is now only indirectly influenced by the Casimirs drastically changes the character of the equilibrium state [35][36][37][38][39]. The poloidal vorticity exhibits no large scale structure, though the velocity field does maintain strong microscale fluctuations. The toroidal velocity field exhibits relatively simple radial band-like structure controlled by the Casimirs [39].
In Section 8, we consider the full shallow water equations, which may be reduced to three coupled scalar equations, with again only one of them, the potential vorticity equation, possessing Casimir constraints. The statistical fluctuations of both the compressional part of the velocity and the surface height remain very strong in equilibrium, and these drive similarly strong fluctuations in the vortex interactions [40], playing the role of an unbounded heat sink that precludes the existence of negative temperature eddy-like states [40,41]. This raises very interesting questions, which cannot be answered by an equilibrium theory alone, regarding the rate at which wave-eddy interactions dissipate such structures if they are created in the initial state, and how they might be maintained (as seen in planetary atmospheres and in experiments) outside of equilibrium. Most optimistically, there may be mechanisms by which additional weak dissipation processes, such as wave breaking, can act to differentially suppress the waves, maintaining the eddies as formally metastable near-equilibria. We exhibit a possible variational formalism, a fairly straightforward generalization of that describing Euler and QG equilibria, that might be used to approximately describe these [42,43]. This system also has a separate adiabatic invariant [20].
In Section 9, we consider magnetohydrodynamic flow of perfectly conducting fluids, which couple mass and electrical current flow through the Maxwell equations. This model has been used to model the solar tachocline [44,45] which marks the very thin 2D boundary between the rigidly rotating radiative interior and the differentially rotating exterior convective zone. The results here are significantly different than all previous examples because the Casimir constraints are tied to the magnetic vector potential instead of the vorticity [34,46,47]. The model that emerges maps onto a pair of interacting elastic membranes in an external confining potential controlled by the Casimirs [48]. The microscale fluctuations are purely Gaussian, and this allows a formally exact derivation of the free energy functional whose minima again determine the large scale structure of the magnetic and flow fields. In the solar context, the structure of these fields has implications for the transport of angular momentum between the two zones.
The paper is concluded in Section 10. It is remarkable how much physical structure the equilibrium theories contain, and how different this structure is for each of the examples treated. There are a number of other well known systems with Casimir constraints [34] that can still be explored. Near-equilibrium generalizations are also of great interest.

Two-Dimensional Euler Equation
It is useful to consider first the simplest system, the 2D Euler Equation [11][12][13] defined by the equation of motion (1) in some 2D domain D. The pressure p is determined by the incompressibility constraint (2)

Vorticity and Stream Function
The constraint (1) permits the stream function representation By taking the curl of both sides of (1), one obtains the vorticity equation which physically states that ω is freely advected by its own induced velocity field v, constructed below. From (3) follows the relation with formal solution ψ(r) = D dr G(r, r )ω(r ), (6) in which the Laplace Green function is the solution to − ∇ 2 G(r, r ) = δ(r − r ) together with the same boundary conditions, on both r and r , satisfied by ψ. Energy conservation requires free slip boundary conditions, equivalent to constant ψ (Dirichlet boundary conditions). If there are multiple boundaries Γ n , n = 1, 2, . . . , n ∂ , e.g., an annular (see Figure 1) or more general multi-holed domain, then ψ = ψ 0,n may be assigned separate values on each boundary and are also constants of the motion. The circulation about each boundary is also conserved. However, since the constants ψ n,0 uniquely define ψ, it follows that the γ n are not independently conserved, but are (linearly) related to the former. Equations (3) and (5) together uniquely determine v in terms of ω, so that the first line of (4) indeed represents a (scalar) closed evolution equation. . Strip and annular (or disc if R 2 = 0) geometries for which, respectively, a conserved linear momentum (10) or angular momentum (11) exists. The strip has periodic boundary conditions along x and Dirichlet boundary conditions on the lower and upper boundaries Γ 1,2 . The annulus has Dirichlet boundary conditions on both boundaries. The latter lead to two independent circulation integrals (8) for each domain (which are seen to actually have the same topology).

Conservation Laws
The conserved energy is just the kinetic energy dr ω(r)G(r, r )ω(r ) (9) in which the boundary conditions ensure absence of boundary terms in the integration by parts used to obtain the second line, and (6) has then been substituted to obtain the last line. If the domain is translation invariant along some directionl (infinite or periodic strip geometry, illustrated on the left in Figure 1) then the corresponding component of the linear is conserved. If the domain is rotation invariant (disc or annular geometry, illustrated on the right in Figure 1) then the vertical component of the angular momentum is conserved. Both of these can be written in the linear form with the choice α(r) =l × r or 1 2 r 2 , depending on the domain. Note that on a true spherical domain, the full vector angular momentum L is conserved and (11) is generalized appropriately.
The "self-advection" Equation (4) implies that any (1D) function of the vorticity is conserved. These may be conveniently summarized by conservation of the function for any value of σ, in terms of which These are very often exhibited in terms of the powers F(σ) = σ n , which are seen to generate the moments of g(σ).

Statistical Equilibrium Concepts
We now summarize the key statistical equilibrium concepts underlying the thermodynamic fluid treatment, especially the key role of microscale entropy. These concepts will serve to define the mathematical basis for computing thermodynamic functions and using them to characterize large scale steady state flows and other quantities of physical interest. A notional picture is illustrated in Figure 2, going back to the original ideas of Onsager [5]. Conventional positive temperature bound "molecule" (left) and unbound plasma states (middle) exhibit no large scale vorticity or flow structure. In this picture, the physically interesting fluid equilibria correspond to much higher energy flows (right) in which the charges are forced to segregate, effectively like attracting like. We will see that such states indeed emerge as negative temperature equilibria.  [5]. The low energy state (1) on the left corresponds to a molecular dipole state with strongly bound charges. The middle state (2) corresponds to a higher energy plasma-like state with unbounded charges but that continue to obey local charge neutrality. The state (3) on the right exhibits large scale structure obtained by increasing the energy even further, forcing the charges to segregate into separate non-neutral regions. This negative temperature state is accessible in fluid dynamics because the charges are not conventional momentum and kinetic energy carrying particles. In the vortex field description, charges carry only potential energy of interaction.
The standard underlying assumption, known as the ergodic hypothesis, is that very long time averages beginning from some given initial condition are equivalent to certain phase space averages over all field configurations consistent with the conservation laws. Figure 3 schematically illustrates this idea for Euler flow, in which the turbulent mixing process eventually produces a smooth looking steady state with the original discrete vorticity levels hidden at the finest scales.
This section will detail the phase space averaging process under the ergodic hypothesis. Ergodicity is almost never provable from first principles and it can indeed be violated even in conventional particle systems. As discussed in Section 5, and hinted at in other sections, violations are known to occur in fluid systems as well, through mechanisms that are understood to varying degrees [1,21,22]. This remains an open area of research. Highly schematic illustration of the turbulent mixing process that begins here with a well defined though irregular region of finite, fixed vorticity ω = q 0 , surrounded by a vorticity free (potential flow) region, ω = 0. Over time the vortex region stretches and folds to give rise as t → ∞ to a fully mixed smoothly varying macroscale steady state. However, the macro-view obscures the continuing microscale dynamics (illustrated in Figure 4) where restriction to values ω = 0, q is preserved, consistent with the Casimir constraints. Within each l-cell one may define the local vorticity distribution n 0 (r l , σ) which has a well defined continuum limit a, l → 0 but in such a way that l/a → ∞. Its first moment defines the equilibrium vorticity (58) and its area integral is constrained by the Casimir function (59). This illustrates the formal limiting process by which, e.g., a discrete set of (a-scale) vorticity levels controlled by the Casimirs produces a smooth (l-scale) average.

Phase Space Measure and the Liouville Theorem
At the purely mathematical level, the statistical equilibrium approach is based on characterizing invariant measures on the phase space Γ of all possible functions ω(r). Phase space integrals with respect to such a measure are therefore time independent and are used to construct physical equilibrium averages.
To be more specific, a probability density functional ρ[ω, t], which here assigns a positive real number to any given field realization ω = {ω(r)} r∈D , evolves according to the conservation law (16) in which ∇ ω is the (infinite dimensional) phase space gradient, and V[ω] is the phase space velocity whose vector components are defined by each point r ∈ D: derived from the equation of motion (4). The linear functional v[ω](r) ≡ v(r) is given by the curl of (6). The form (16) ensures conservation of probability for any phase space volume co-moving with the phase space flow. Now, an equilibrium probability density ρ eq [ω] allows one to define equilibrium averages in the form The functional integral is defined here by a limiting process in which r i , i = 1, 2, 3, . . . , N a = A D /a 2 , with A D the area of D, runs over a uniform grid (e.g., square lattice) with elements of area a 2 → 0. For all such averages to be timeindependent, ρ eq must be as well and hence obey On the other hand, the equation of motion for any functional I[ω, t], defined by takes the phase space advective form In particular, if I is a conserved integral then it must obey The key observation is that if the phase space flow obeys the "phase space incompressibility condition" then the equilibrium measure condition (19) reduces to Comparing (23), this corresponds to the requirement that ρ eq be a conserved integral. This is the content of the Liouville theorem.
Liouville Theorem for the Euler Equation The most transparent way to verify the Liouville theorem for the Euler equation, avoiding continuum functional derivatives, is to represent ω as a discrete orthogonal mode expansion on the finite domain D. We consider an expansion of the stream function in Laplacian eigenmodes: with positive eigenvalues λ l > 0 for a finite domain. The φ l obey the same (Dirichlet or periodic) boundary conditions that ψ does and may be taken to be real and orthonormal. It follows from (5) that and the equation of motion for ω l may be derived in the forṁ with coefficients W lmn = D drφ n (r)∇φ l (r) × ∇φ m (r).
These are totally antisymmetric with the third one obtained via integration by parts, and making use of the free slip boundary condition to eliminate the boundary term. Using this representation, one obtains However, the coefficients W mnn = 0 all vanish by virtue of the antisymmetry result. Thus, ω n does not actually appear on the right hand side of (29), trivially verifying the Liouville condition.
For completeness, some hints to a real space derivation may be provided as follows. The direct functional derivative produces This result, in the limit r → r, is quite singular, due to both the delta function and the logarithmic singularity in G. However, one may make sense of it by recognizing that in free space a single point vortex remains stationary. Thus, the logarithmic singularity G F = −(2π) −1 ln |r − r | in G does not contribute to self-induced motion, only the boundary-induced image correction G(r, r ) = G F (r, r ) + Φ(r, r ). Removing the free space contributions, the integral (23) may then be reduced to the boundary integral of n · ∇ × Φ(r, r), which vanishes in the limit because the flow due to the single, opposite sign, image vortex infinitesimally on the other side of the boundary obeys the free slip condition at r.

Microcanonical Ensemble
The choice of equilibrium measure goes by the name of statistical ensemble. Perhaps the most transparent choice is the microcanonical ensemble, in which one constrains a particular value c γ to each conserved integral C γ [ω] and we use the shorthand c = {c γ }. For the Euler equation this clearly involves an infinite product, which will be further characterized below. The energy is separated out explicitly for convenience. Equilibrium averages (18) by construction limit the support of the phase space integral to vorticity fields constrained by the specified values ε, c.
The partition function serves to normalize ρ µ as a probability density, but also defines the entropy function through the Boltzmann relation in which the factor 1/N a yields a finite, well defined result in the continuum limit (20), here seen to play the role of the thermodynamic (infinite volume) limit in conventional systems. Explicit examples will be given below. All thermodynamic quantities follow from the entropy function in the usual way. Most critically the inverse temperature is obtained from the energy derivative, and more generally the derivative defines the thermodynamic field µ γ conjugate to c γ .

Grand Canonical Ensemble
It is generally extremely difficult to compute delta function constrained integrals such as (35). Instead one seeks to make use of the thermodynamic analogue of Lagrange multipliers by switching to a smoother probability distribution. Thus, the grand canonical ensemble is defined by the Laplace transform with partition function and statistical functional now including fields µ = {µ γ }. These now replace the conserved integrals c as the fundamental thermodynamic variables. The subscript a on β allows for the fact that the inverse temperature β a = 1/T a might need be scaled nontrivially in order to obtain a consistent thermodynamic description in the continuum limit. It will in fact be shown below that the scaling is required, with finite values of β = 1/T smoothly controlling the equilibrium state. This scaling is essentially required to control a nontrivial balance between energy and entropy (fluctuation) effects. Roughly speaking, equilibrium flows have lower temperature, T a = Ta 2 → 0, than that of any conventional thermodynamic system! The physical meaning of this will be discussed below.

Thermodynamic Free Energy
The partition function is now related to the thermodynamic free energy by and, with these scalings, is also finite and well defined in the continuum limit. Note that β a = (β/A D )N a so that this actually involves the same a-scaling as the entropy (36). From the definition (39) the derivatives produce the thermodynamic averages of the conserved integrals, defined here bȳ The standard equivalence of ensembles in the thermodynamic limit, which requires showing that the averagesc γ [β, µ] are in fact infinitely sharply peaked about a single unique value of C γ [ω], can be shown to follow here from the continuum limit N a → ∞ (with peak width scaling as 1/ √ N a ). Other well known ensembles correspond to partial Laplace transforms over a subset of the conserved integrals. The canonical ensemble corresponds to transforming only the energy, resulting in statistical weight e −β a E[ω] multiplying the remaining delta functions. In a number of conventional systems the energy is actually the only conserved integral. There may be cases of "ensemble inequivalence" where dealing with the delta functions provides a more physically consistent approach [1]. However, even in such cases it is generally much simpler to apply the grand canonical approach and then use physical arguments to adapt it after the fact to more broadly enforce equivalence.

Grand Canonical Formulation of the Euler Equation
For the Euler equation the index γ includes the continuous index σ appearing in (14), and one obtains the more explicit form in which the conserved momentum (12), when it exists, enters with conjugate field µ P . The 1D field function µ(σ) is conjugate to the conserved function g(σ) ≡ g[ω; σ], promoted here to a functional of the vorticity. Inserting this form into (41) produces the field theory displayed in the first row of Table 1.

Thermodynamics of the Euler Equation: Exact Solution
We will now show, quite remarkably, that the Euler equation free energy (43) may be computed exactly. More specifically, the evaluation of the full phase space integral (40) may be reduced to a variational equation for the free energy from which the equilibrium vorticity function is obtained as a solution to a (highly nonlinear) PDE generated by the corresponding Euler-Lagrange equation. The derivation here will be physically motivated rather than rigorous-full details may be found in Ref. [13]. Such variational approaches often emerge as approximate "mean field" descriptions of conventional thermodynamic systems. Here the mean field form is in fact exact due to the long range (Coulomb-like) interactions (9) between vortices.

Mean Field Approach
The key property of the energy function (9) is that it is dominated by the long range nature of G(r, r ): in the macroscopic coherent flow regime of interest here the stream function (and therefore the advection velocity field) is dominated by the global integral (6) over the entire domain. In contrast to systems with local interactions, the contribution from a small area l 2 about r here scales as l 2 ln(l) → 0. It follows that if one considers a fluctuation ω(r) − ω 0 (r) about the equilibrium field one may accurately replace in which is the equilibrium stream function. The inverse relationship then also follows. In conventional particle systems G is typically a short ranged microscale interaction, ψ 0 is therefore dominated by local fluctuations on the same scale as ω, and (48) is at best approximate. Here the Casimirs strongly bound the fluctuations of ω, the Green function effectively performs a self-averaging operation so that ψ(r) − ψ 0 (r) → 0 in the continuum limit with probability one, and (48) becomes exact. This statement can be made rigorous by letting both the intermediate scale b and the grid scale a → 0, but with b/a → ∞ [13]. In this way the local integral of ω − ω 0 against the smooth function ψ 0 in (48) scales as a/b → 0. Thus, the Casimir constraints ensure that although ω is discontinuous from grid point to grid point on the microscale a, its local fluctuations are bounded by the support of g(σ). It follows also that the velocity v(r) is continuous and the stream function ψ is continuously differentiable. The form of ψ 0 must now be determined self consistently by using (48) to compute the free energy. Substituting (48) and (46) into (41) one obtains with 2D function This form is now purely local in the fluctuating field ω(r). The temperature scaling (42) is now seen to be chosen to enable the replacement and the partition function then follows in the product form in which we define the 1D function essentially the Laplace transform of e βµ(σ) . Explicit forms for W obtained from simple model forms for µ(σ) will be discussed below.
Taking the logarithm of (54) and restoring continuum notation, the final free energy functional (43) takes the form with the dependence on the thermodynamic fields β, µ now highlighted explicitly. The scaling (42) is again confirmed to yield a well defined finite result. A self consistent equation for ψ 0 is obtained by generalizing the free energy calculation to compute equilibrium averages (45). The fundamental quantity needed is the vorticity distribution function (illustrated in Figure 4) This simple result follows from the cancellation of the integrals over all other ω i = ω(r) between the numerator and denominator of (45). This function quantifies the fluctuations of the vortex field in the microscopic neighborhood of any given point r (defined by the l-cells in Figure 4). In particular, the mean vorticity is derived in the form The right hand side is a local function of ψ 0 (r), so that we have produced a type of nonlinear Poisson equation for ψ 0 .
In addition, the Casimirs (14) are recovered from the area integral which allows one, in principle, to invert for µ(σ) for specified g(σ). The identical result may be shown to follow from the functional derivative This derivative is performed only with respect to the explicit µ dependence in (56), keeping ψ 0 fixed. This works because the self-consistency condition (58) is equivalent to the extremum condition which zeros out the δψ 0 /δµ(σ) contribution to (60). The equilibrium momentum may similarly be derived either from ω 0 or from the free energy derivative. The mean fluid kinetic energy follows as well either by substituting ω 0 (r) = ∇ × v 0 (r) into (9) or from the β derivative exhibited in (44).

Microscale Entropy
The distribution function (57) also allows one to introduce the important concept of the microscale fluid entropy. The equilibrium flow defined by ω 0 and ψ 0 is smooth, in general infinitely differentiable on any finite physical length scale. The equilibration process may be thought of as the completion of the inverse cascade of energy, which serves to create the inhomogeneous flow on the domain scale A D , and the forward cascade of enstrophy (and all other Casimirs) to infinitesimal scales that render the microscale fluctuations invisible. Of course, additional physical dissipation processes such as viscosity will eventually smooth out these microscales, but this not necessary to make sense of the idealized fluid equilibria considered here.
Using (57) the equilibrium entropy may be expressed in the classic information theoretic form This precisely captures the information lost in going from the exact microscale specification of the finely mixed vorticity field ( Figure 4) at any given instant of time to the timeindependent equilibrium average, in which only ω 0 is specified. For any given distribution n 0 , not necessarily equilibrium, one may derive (65) from the Boltzmann formula which may be compared to the microcanonical expression (36). The derivation proceeds via the previously described limiting process in which one counts the total number of ways N[n 0 ] to distribute the (l/a) 2 vorticity levels contained in the intermediate scale area l 2 , with level populations constrained by n 0 (essentially an a-cell permutation count repeated over all l-cells). In fact, an alternative rigorous microcanonical approach to deriving the free energy functional (56) is to maximize S[n 0 ] subject to the all of the conserved integral constraints [13]. The maximal solution for n 0 is recovered precisely in the form (57).

Rotating Fluids and Generalization to the Beta Plane
Before turning to explicit examples and further generalization of the theory, it is worth treating the simplest extension to rotating fluids. The beta plane approximation incorporates planetary rotation through the generalization in which f (r) = 2ẑ · Ω = 2Ω sin(θ L ) is the Coriolis function derived from the local vertical projection of the angular rotation vector Ω corresponding to latitude θ L (r). The curl of this equation leads to self-advection of the potential vorticity exhibiting the sum of local and frame of reference rotation rates. The kinetic energy (9) remains unchanged, but is now expressed in terms of ω P by substituting ω = ω P − f . Similarly, for the momenta, which are now conserved only if f (r) possesses the required invariance-constant latitude (east-west) periodic strip, or disc or annulus surrounding the pole. The equilibrium free energy follows in a form identical to (56), but with in which F is the solution to Poisson equation For linear f = βy on a strip, or f = βr on a disc or annulus (beta plane linear approximation), one obtains the cubic form F = − 1 6 βy 3 or F = − 1 9 βr 3 . The result is the combination acting as an "external potential" ψ 0 − Ψ P inside the W function in (56). Since the two functional forms are different [linear or quadratic-see (12)

More General Curvilinear Domains
More generally the Euler equation on a 2D curved (in particular spherical) surface, with and without rotation, may developed as well [26]. The vorticity and stream function may be defined by adopting appropriate curvilinear coordinates, and the generalization of the self advection dynamics (4) for the vorticity then follows. The conserved Casimir area integrals then follow immediately as well, as does the statistical theory leading to a free energy functional in a form very similar to (56).

Simplified Model Examples
The equilibrium Equation (58) looks quite complicated, but some very interesting, physically meaningful results may be derived by specializing to few-parameter models. We will focus on the two-level system in which the vorticity field is constrained to take values 0 or q only (illustrated in Figure 3). Since for given domain area A D there is only a single degree of freedom, one may normalize in which the single conjugate field µ q is used to adjust the relative areas of the vortex "charges." Substituting into (58) one obtains the equilibrium equation with a Fermi-like distribution function on the right hand side, and in which for simplicity we set the momentum to zero (if it exists) by taking µ P = 0. For large β → ∞ (T → 0 + ) the solution is ω 0 = 0 on the region where ψ 0 < µ q /q, and ω 0 = q on the compliment, so that the equilibrium solution is also two-level. This solution corresponds to the lowest possible energy state, and by Gauss's law spreads the vorticity out as much as possible (equal-signed charges repel), distributing it up against the boundary of D. On the other hand, for large β → −∞ (T → 0 − ), which is perfectly allowed in this system, the two regions switch roles, with ω 0 = q on the region where ψ 0 > µ q /q and ω 0 = 0 on the compliment. The solution corresponds to the highest possible energy (equal-signed charges now effectively attract), and the result is a single compact vortex somewhere in the interior of D. Varying µ q varies the position of the vortex boundary, hence size of the vortex. As one varies −∞ < β < ∞ the vortex edge will be smeared out on the scale |T| = 1/|β| and the solution will continuously interpolate between these two extremes. Figure 5 illustrates these results for a unit disc domain. The solutions for this simple case are azimuthally symmetric, functions of the radius r alone. This behavior of the solution as one varies −∞ < β = ∂S/∂E < ∞ in accompanied by a very interesting picture of the energy dependence of the entropy S(E), illustrated in Figure 6. The entropy vanishes for |β| → ∞ (T → 0 ± ) corresponding to the minimum energy E min (vorticity compacted against the boundary) and maximum energy E max (vorticity compacted at the center). The maximum entropy occurs for β = 0 (maximally disordered uniform vorticity state at T → ±∞) but at some intermediate value of the energy. For conventional particle systems, the particle momenta are permitted to grow without bound and S(E) diverges with E → ∞-the curve never turns over and negative temperatures are forbidden. Figure 5. Example equilibrium vorticity profiles ω 0 (r) for the two level system on the unit disk for a sequence of inverse temperatures −∞ ≤ β ≤ ∞, obtained by numerically solving the nonlinear Laplace Equation (74). Vorticity level q = 1 occupies fractional area α = 0.2, hence total vorticity Ω 0 = πα. For each temperature, the Lagrange multiplier µ q (β) must be determined iteratively to satisfy this constraint. As seen, the β = −∞ (T = 0 − ) maximum energy solution gathers all vorticity near the disc center, while the β = +∞ (T = 0 + ) minimum energy solution compacts all vorticity against the disc boundary. The β = 0 (T = ±∞) maximum entropy solution distributes the vorticity uniformly (center panel).
More interesting behaviors may observed in annular domains [13,15,16,21] where the azimuthal symmetry may be broken (a form of second order phase transition). Dynamically, a zonal jet (symmetric vortex ring in this case) becomes unstable and at late time forms a simply connected Red Spot-like vortex. Within the equilibrium theory, the energy advantage of a more compact shape leads to spontaneous azimuthal symmetry breaking for decreasing negative temperature.
Other interesting behaviors may be explored using the three level system with A 0 + 2A q = A D , and equilibrium equation The model here is simplified by enforcing symmetry between charges ±q. High energy, negative temperature equilibria, for example, with two separated, opposite-signed vortex blobs may be constructed. On the other hand, low energy states correspond to fine-scale intermixing of the two charges, generating a (conventional) featureless, neutral system with no macroscale flow structure.
Breaking the symmetry between the charges, µ q = µ −q , allows one to separately control the relative size of these blobs, and eliminate full cancellation at positive temperatures. Figure 6. Schematic illustration of the entropy function S(E) associated with the two level system (72), and also of the point vortex system pictured in Figure 2. As described in the text, the Casimir constraints on the vorticity allow for both positive and negative temperatures, and corresponding entropy limited to a finite energy interval, vanishing with infinite slope at both ends. This general picture will hold for any g(σ) with bounded support. The dashed line corresponds to conventional particle systems in which the momentum degree of freedom can absorb unbounded energy.

Metastable Steady States
Ergocity is a statement that long time averages are equivalent to an equilibrium average, namely that any initial condition ω(r, 0) will explore essentially all of the phase space permitted by the basic conservation laws. There are a number of cases where this assumption can fail, and the highly constrained nature of 2D flows exacerbates this (see Section 3.1.2 of Ref. [1] for some discussion on this point, as well as the more recent review [24]). This is in contrast, for example, to conventional gases with their randomly moving atoms and molecules and interpenetrating motion trajectories. Thus, equilibration of 2D flows is perhaps more closely analogous to dense, glassy system dynamics with strong barriers to individual particle motions.
Some of these barriers can actually be understood as local rather than global minima of the free energy functional. Examples include separated compact eddies that orbit each other, failing to merge (as would be entropically favored) above a critical separation [21]. There are numerous examples of more conventional systems that show analogous behaviors, including decay of superflow (which requires, e.g., nucleation of an eventually systemspanning vortex ring), and metastability of certain crystal structures such as that of carbon's diamond. Detailed 2D Euler equation numerical simulations also show evidence for different levels of equilibration in different spatial regions, depending on the strength of local mixing dynamics [22].

Viscosity Effects
The varying roles of microscale viscosity should also be pointed out. Of course, viscous effects standing in for thermal exchange between the micro-and macro-levels (as well as bottom and other forms of friction with the world outside the idealized 2D domain), will eventually lead to strong violations of the vast majority of conservation laws, producing decay to an essentially trivial flow. The statistical mechanics approach can at best be valid on intermediate time scales where such effects can be ignored.
At a more subtle level, the l-cell picture in Figure 4 will be the first to be violated. The individual a-cells will diffusively mix to form a uniform, average vorticity essentially coinciding with the local average ω 0 (r) defined by the second line of (58). It is important to emphasize that this by itself does not violate the statistical mechanics predictions: the macro-scale flow is insensitive to such microscale averaging. Although the Casimir function g(σ) is lost in this process, replaced by the "diffusion-mixed" form it can be shown that ω 0 (r) can be consistently derived from g d (σ) as the extremum energy solution of the equilibrium Equation (58) (|T| → 0 or |β| → ∞, depending the sign of β prior to the action of the viscosity). This is formally established as the "dressed-vorticity" corollary in Section V-D of Ref. [13].
A more nuanced definition of the intermediate time scale is therefore that it be small enough that large scale flows are not significantly affected by viscosity, but not so small that it unnecessarily forbids the occurrence of simultaneous fine scale diffusive mixing during the course of the late time turbulent cascade.

Strongly Fluctuating Long-Lived States
Finally we note that equilibration dynamics can be strongly affected by the topology of the 2D domain. Thus, equilibration on the surface of a sphere [49] (rather than in a flat bounded or doubly periodic plane) is found to fail much more catastrophically, with a macroscopically fluctuating chaotic vorticity field surviving for all achievable computation times [25][26][27][28]. Conservation of the full angular momentum vector in the spherical geometry (rather than just a vertical component) ensures that the vorticity cannot condense into a dipole pattern if the initial state has zero total angular momentum. Numerical experiments show that most (but not all) of the vorticity condenses instead into a quadrupole with two positive vortices and two negative vortices; small satellite vortices also persist [26]. The quadrupolar configuration oscillates, likely chaotically, at long times, consistent with the dynamics of the much simpler problem of four point vortices on the surface of sphere [29].

General Statistical Theory of Single-Field Systems
We next consider generalization of the statistical approach to other 2D systems characterized by an infinite number of conserved Casimir-type integrals constraining the dynamics of a single scalar field [33]. The former requires the existence of a self-advecting field q(r, t) with equation of motion generalizing (4) in the sense that the relation between q and the velocity field v may be more general. A convenient way to constrain this relationship, and simultaneously ensure existence of a consistent statistical theory, is to demand that (78) be derived from a Hamiltonian equation of motion
In particular, if one defines the stream function by then (79) takes the form which is exactly (78) with the usual stream function relation (3), and from which incompressibility of v also follows immediately. Conservation of the Casimirs (13)-(15), with q replacing ω, follows immediately as well.

Statistical Mechanics
The Liouville theorem dr δq(r) δq(r) = 0 (83) follows directly from the Hamiltonian structure. Specifically, using a mode representation (26) for q, one obtainsq with coefficients defined by (30), and it follows that in which the vanishing follows because W lmn is antisymmetric in l, m [see (31)] while the mixed partial is even. Note that there is no assumption here that E[q] is quadratic in q, though for many of the standard examples it is.
The key consequence is that the phase space measure is defined simply by replacing ω by q in (19). The statistical ensembles (34) and (39), and the form (46) continues to define the Lagrange multiplier function µ(σ). The momentum term (12) is also obtained by simply substituting q for ω. Formally this follows from the identity {q, L} = ∂ ξ q, where ξ is the symmetry coordinate, which shows that L is the generator of translations along ξ.
The variational result for the free energy proceeds by following the steps (48)-(56), but with the replacement in which L[ψ 0 ] (given by the domain integral of − 1 2 |∇ψ 0 | 2 for the Euler equation) is the Legendre transform of E[q 0 ], obtained by inverting the relation ψ 0 [q 0 ] [generalizing (6)] to obtain q 0 [ψ 0 ] [generalizing (50)] and substituting the result into the first line of (86). This inverse relationship is also encoded in L via the general Legendre transform relation The free energy now generalizes to in which the Lagrange multiplier-Laplace transform W(τ) continues to be defined by (52) and (55). The local q-distribution function continues to take the from (57), and the equilibrium Equation (58) generalizes to

Quasi-Geostrophic Flow and Nonlinear Rossby Waves
Quasi-geostrophic (QG) flow, including the Coriolis term f described in Section 4.3, is defined by in which R 0 = 1/k R = c/ f is the Rossby radius of deformation, the length scale beyond which Coriolis effects begin to dominate gravitational/hydrostatic effects on the fluid dynamics, with c the speed of internal gravity waves. On Earth one may estimate this speed as follows: write c = g eff H eff where H eff is the effective fluid layer depth and g eff the effective acceleration due to gravity. For single layer shallow water theory, g eff = g and H eff = H are the "bare" physical values. For internal waves in a density stratified medium g eff /g ∝ δρ/ρ is reduced by the density constrast between layers, and H eff is the scale height, namely the effective height of the water column actually taking part in the motion (e.g., thermocline depth). In the Earth's oceans, δρ/ρ ∼ 10 −2 and the Kelvin wave speed c ∼ 2 m/s is therefore O(100) times smaller than the "bare" shallow water wave speed. The QG model, whose large scale wave excitations are known as Rossby waves, emerges from the shallow water equations, discussed in Section 8, in the limit where the surface height adiabatically follows the eddy motion via quasi-hydrostatic balance. Higher frequency traveling surface wave excitations are neglected. The energy function is in which the Green function now obeys the Poisson equation [compare (7)] In free space one obtains the modified Bessel function form which maintains the Euler equation logarithmic singularity near the origin, but decays exponentially ∼ e −|r−r |/R 0 on the scale of the Rossby radius (which depends strongly on latitude, but is on the order of 50 km at mid-latitudes on Earth). The free surface motions therefore act to screen the vortex charge at larger distances. With this adjustment of G, the statistical functional (41) continues to take the general form of the field theory displayed in the first row of Table 1 (though f has been dropped there for simplicity). The Legendre transform operation yields the form Again, the momentum functionals are identical to those of the Euler equation, with the same function α(r) as appearing in (10)- (12). The equilibria of this system have been explored by a number of authors [1,32,33]. An interesting aspect of the vortex screening, and resulting finite range interactions, is that eddies with size much larger than R 0 have identical physics as fluid droplets with finite surface tension. Thus, the transition between interior and exterior of the eddy occurs over length scale R 0 , and a surface energy per unit length Σ 0 (R 0 , β, µ) (see, e.g., Figure 1 in Ref. [33]) may be assigned to this interface. The shape of the eddy is obtained by minimizing the total surface energy subject to the effective external forces provided by the Coriolis and angular momentum effects. In particular, the Coriolis term is analogous to an external gravitational field and the vorticity is analogous to a mass density. It follows that the equilibrium state will tend to organize with "lighter" regions of lower vorticity floating on (northwards of) "heavier" regions of higher vorticity. This provides a partial explanation for the ubiquity of "zonal jet" structures, with compact eddies requiring a less commonly occurring balance of forces.

Adiabatic Conservation Laws and Slow Equilibration
The QG equation has an added complication that, in addition to conservation laws treated so far, it also has an approximate adiabatic invariant B[q] [19,20]. Like the energy, B is quadratic in q, is insensitive to the microscale fluctuations, hence dominated by the large scale flow. Its conservation improves as the flow becomes more weakly nonlinear, hence very often as the turbulent state relaxes and the flow equilibrates. It is likely, therefore, that this invariant acts as (yet another) barrier to full equilibration-its approximately conserved value, computed from the initial state, will typically be different from that computed from the equilibrium state based on energy and Casimir conservation alone.
It may be argued, for example, that preserving B constrains the inverse cascade to focus energy on wavevectors close to the y-axis (axis of rotation), hence (through the curl relation) enhancing the formation of zonal flows (organized normal to that axis). The equilibrium solutions (ignoring B) also often yield zonal flows, but the north-south geometry will in general be different. The full consequences of this competition deserve to be more fully explored.

Generalized Surface Quasigeostrophic Equations
The generalized quasigeostrophic (GSQG) equation replaces the vorticity-stream function relation with in which (−∆) α is multiplication by |k| 2α in Fourier space and α = 1 corresponds to the conventional Euler equation. The equation of motion for ω retains the self-advection form (4), but now with velocity The Casimirs (13) and (4) remain unchanged, but the energy function is now where G α = (−∆) −α is the Green function of the (−∆) α operator on D. At small separations it has the power law form ∼ |r − r | 2(α−1) , replacing the logarithm for α = 1. The generalized equation is useful since it allows one to study the mathematical properties of the flow as the |r − r | → 0 singularity is varied [50]. The derivation of the mean field Equation (58) follows exactly as before, except that −∇ 2 ψ 0 is replaced by (−∆) α ψ 0 in the first line. Care should be taken here because the mean field approximation (48) is valid only if the interaction is sufficiently long-ranged. This will fail for sufficiently small (perhaps negative) α, and exploring this would be an interesting question for future investigation. For example, the question of solution regularity remains open for α = 1 2 , and it would be interesting to see if this is reflected in the equilibrium flow properties.

3D Axisymmetric Flow
Over the next few sections we will briefly review applications of equilibrium ideas to yet more complicated fluid systems. More details may be found in the referenced literature. The steps outlined in the previous sections-Liouville theorem and statistical measure, choice of equilibrium ensemble, entropy and free energy functions-remain highly relevant, but the exact variational solution derived for the Euler equation is in general no longer available. Rather, it becomes an approximate tool, along lines similar to the use of mean field theories in conventional systems. Specifically, the equilibrium states, though still constrained by an infinite number of conserved integrals, now contain further degrees of freedom (such as a free surface height or other additional coupled field) that escape the constraints, and continue to exhibit fluctuations on finite scales.

Axisymmetric Equation of Motion
The case of 3D axisymmetric flow, illustrated in Figure 7, will be our first example of the impact of an additional degree of freedom, not constrained by Casimirs [35][36][37][38][39]. Under the constraint of azimuthal symmetry, and specializing to cylindrical coordinates, the full 3D Euler equation velocity field may be written in the form is the vertical component of the angular momentum density and characterizes "toroidal" flow around about the axis, while the 3D incompressibility condition allows one to express the "poloidal" flow components in terms of a stream function ψ(r, z) via v r = −(∂ z ψ)/r, v z = (∂ r ψ)/r. The latter is related to the poloidal vorticity ω θ =θ · ∇ × v via which serves to define a modified radial coordinate y = r 2 /2 and modified 2D Laplacian ∆ * . Defining the 2D coordinate ρ = (y, z), the formal inverse of the latter is obtained from the (Dirichlet) Green function relation generalizing (6) and (7). Defining the modified 2D gradient ∇ ρ = (∂ y , ∂ z ) and velocity w = ∇ ρ × ψ = (rv r , v z ), one obtains the incompressibility condition ∇ ρ · w = 0, and the Euler equation may be reduced to the coupled pair of scalar equations The first states that the toroidal velocity field is in essence a passive scalar that is freely advected by the poloidal velocity field generated by q and obtained from the curl of (101). The second states that the self-advection of the poloidal vorticity field is additionally forced by s, a type of Coriolis effect. In the absence of such forcing the q equation would be formally identical to the Euler Equation (4). The effects of this forcing play a critical role in the statistical equilibria, which therefore differ strongly from those of the Euler equation. Figure 7. Axisymmetric flow geometry confined to a cylinder of height H and inner and outer radii R 1 < R 2 . The pattern of flows is taken to be invariant under rotation about the cylinder axis, and is therefore specified by a toroidal flow field s about the axis, and a poloidal vorticity field q within any 2D radial planar section D.

Conservation Laws
In addition to conservation of total (kinetic) energy one obtains two classes of Casimir-type constraints. The first of (102) leads directly to conservation of the domain integral of any function F(s), which may be characterized by conservation of the function [compare (14)] For a strictly finite cylinder of height H, with Dirichlet boundary conditions on all surfaces, Equation (104) comprises all of the Casimirs-there is no constraint on q at all. However, in the case of periodic boundary conditions in z (termed a Taylor-Couette type geometry) with specified period H, it follows from the second of (102) that the domain integral of any combination of the form qF(s) is conserved as well, characterized by conservation of the functiong In essence, Dirichlet boundary conditions impose additional forces on the vertical motion of the fluid that destroy this constraint. In most of what follows we will assume periodic boundary conditions since it leads to more interesting results.
The additional constraint (105) implies that the mean value of q over each level set s(r) = σ is conserved, but may otherwise fluctuate arbitrarily. In particular, there is no control over the range of values that q may take, permitting unbounded fluctuations about this mean.

Axisymmetric Equilibria
Along similar lines to that derived in Section 3.1, the Liouville theorem leads to equilibrium measures that must take the form of a conserved integral, with phase space integral defined by the continuum limit of free integration over the s and q fields: Details may be found in App. A of Ref. [39]. The grand canonical ensemble is defined by with inverse temperature β a = β/a 2 again scaling with a [see (42)]. Lagrange multiplier functions µ(σ) andμ(σ) enforcing conservation of g(σ) andg(σ), respectively, are introduced through the grand canonical statistical functional This model, reproduced in the second row of Table 1, takes the form of a purely local field s, with no self interactions and µ(σ) playing the role of a local potential energy, linearly coupled to an unconstrained Gaussian field q. As such, its thermodynamic behavior bears little resemblance to that of the Euler Equation (41) with (9) and (46). We summarize here its basic properties-full details may again be found in Ref. [39]. The first observation is that the magnitude of q is controlled only by the positive definite quadratic form E G [q]. The linear term qμ(s) serves (by completing the square) only to shift the mean. This type of shift is exactly the degree of freedom required to enforce the second set of Casimirs (105). Being quadratic, the resulting Gaussian statistical averages over q are finite and well defined only for positive temperatures, β > 0. However, being Gaussian, arbitrarily high energy flows may be created at positive temperature, so all of the conservation laws continue to be satisfied. In contrast, for the 2D Euler equation negative temperatures may be required because the Casimir constraints also bound, through ω, the maximum energy of positive temperature states.
The end result is that because fluctuations about the local mean in q are uncontrolled, one obtains identically vanishing mean stream function ψ 0 (ρ) and poloidal vorticity q 0 (ρ) = −∆ * ψ 0 (ρ). Hidden from these are the finite averages of higher order quantities, such as the mean square velocity |∇ ρ ψ(ρ)| 2 ∝ T > 0 (an equipartition result). In this sense the equilibria are similar to those of conventional particle systems.
The second observation is that if the initial flow is such that |s| is bounded, then e βµ(s) will be as well. Thus, statistical averages over the field s are well defined irrespective of the value of β (either positive and negative, though as seen q requires β > 0). However the s(ρ) 2 energy contribution from s is purely local, and the long range Coulomb interaction effects seen in the 2D Euler case are absent here. Given the absence of any finite scale structure in q, the qμ(s) term may be shown to play a negligible role in the statistics of s, and one obtains the exact q-independent result for the local distribution of σ. In particular the local mean is derived in the form

Pure Poloidal Flows
If one imposes zero toroidal velocity, s = 0, the right hand side of the equation of motion (102) for q vanishes. Thus, q = ω/r becomes purely self-advecting and its Casimirs are now conserved. The equilibrium theory for q therefore becomes 2D Euler-like (rather than Gaussian), though now governed by the modified Laplacian (100) and Green function (101), which still retains its logarithmic character away from the symmetry axis.
This special case is of some interest because there have been recent results demonstrating finite time singularities for axial flows in the cylinder (R 1 = 0) obtained from a singular initial condition with q ∼ 1/r 1−α for r → 0 and sufficiently small α [51]. This initial flow singularity is also built into the equilibrium theory since it leads to the power law form g(σ) ∼ σ −1−2/(1−α) , σ → ∞, for the Casimir function (14) (with q here replacing ω). One expects then a corresponding power law form for the Lagrange multiplier function µ(σ). The effects of such singularities would certainly be interesting to explore, in particular whether the same values of α governing the finite time singularities play a special role. Of course, there are many initial conditions leading to the same conserved integral values, and a true finite time divergence may only occur for special choices. Thus, α may play a different role in the existence and structure of the equilibrium states.

Axisymmetric Flow Equilibration Issues
Just as for the Euler equation, there are significant questions regarding the convergence to equilibrium for Taylor-Couette flows of this type. In particular, experiments do appear to show very long lived negative temperature-type states, with q displaying large scale coherent structure (see Refs. [35,37,38] and references therein). Reasonable comparisons with experiments were obtained in Ref. [38] by artificially bounding |q| < M, with M remaining finite in the continuum limit a → 0, and applying Euler equation mean field ideas to obtain negative temperature states for the resulting altered model.
Elucidating the barriers to equilibration, limiting or greatly slowing the growth of |q| predicted by the model (108) as the forward cascade proceeds, remains an interesting open question. We will encounter very similar issues below in relation to the surface height field for the shallow water equations.

Shallow Water Dynamics and Wave-EDDY Interactions
The shallow water system, illustrated in Figure 8, is defined by the equations of motion in which v is the horizontal velocity, h is the (fluctuating) fluid free surface height over a flat bottom, z = 0, and for convenience we include the Coriolis parameter f from the outset. Comparing to the Euler form (67), the pressure gradient is provided by changes in surface height, and the second equation expresses incompressibility of the full 3D velocity by relating surface height change to the divergence of the mass current j = hv. These equations are derived from the 3D Euler equation in the formal asymptotic limit in which the length scale of horizontal variability (including the horizontal extent of the 2D domain D) is much larger than h, and v is approximated as independent of z. The vertical velocity is then v z = −z∇ · v, hence v z (r, h) = −h∇ · v, and it follows that the h equation may be equivalently written in the intuitive form Dh/Dt = v z (h).

Conservation Laws
It is straightforward to check that the potential vorticity is advectively conserved, DΩ/Dt = 0, which clearly reduces to (68) for fixed surface height. The corresponding conserved Casimir area integrals are Integrating this over σ, one obtains in particular conservation of the mean height Since the fluid is compressible, the additional compression field is required to fully reconstruct the velocity in the form v = ∇ × ψ − ∇φ (116) in which the stream function ψ and potential function φ are obtained by solving Free slip boundary conditions on v require as before Dirichlet boundary conditions on the stream function ψ, but Neumann boundary conditionsn · ∇φ| ∂D = 0, on the compression potential. Thus, one obtains with subscripts labeling the Green function boundary conditions. Both are long-ranged, with logarithmic singularities at the origin.
The conserved energy is a sum of kinetic and potential terms. By substituting (116) and (118), the kinetic term may be organized in the form in which the components of the 2 × 2 tensor Green function G h (r, r ) is an integral-product of h with gradients of G D and G N . The exact from is not needed in what follows since expressions in terms of ψ will reemerge as central in the statistical analysis. In the presence of translation or rotation symmetry, momentum conservation analogous to (10) or (11) also occurs, but will be neglected here for simplicity.

Liouville Theorem and Statistical Measures
Proving a Liouville theorem for this system is much more involved, and we only quote the result here-a full derivation may be found in App. A of Ref. [40]. The simplest approach, conceptually, is to treat the height h and the two components of the mass current j = hv as fundamental canonical variables. In terms of these it can be shown that the correct phase space integration measure, accompanying the conserved equilibrium density ρ eq [h, j], continues to be defined by approximating the domain D by a uniform mesh, with lattice parameter a, and freely integrating over the discretized fields: From this representation, using finite difference approximations to the gradients, one may change variables from v to (Ω, Q) to obtain The only slightly usual feature is the height measure h 4 dh coming from the various changes of variable.

Shallow Water Equilibria
The grand canonical form of the equilibrium measure ρ eq = Z −1 e −β a K is again obtained by introducing the Lagrange multiplier function µ(σ) to control the Casimir constraints. The statistical functional is also displayed, with a more compact notation, in the third row of Table 1. Like the Euler equation, and unlike for the axisymmetric flow model, the vortex degree of freedom Ω is directly constrained by the Casimirs. However, the additional (height and compression) degrees of freedom enter in a complicated way that makes this model very difficult to analyze. Height fluctuations, controlled only locally by the 1 2 gh 2 term, are strongly coupled to Ω, and forbid any simple reduction to a mean field type description.
In order to gain some intuition and make closest possible contact with the Euler equation, one may integrate out the Gaussian field Q to obtain the reduced functional in which the scalar Green function satisfies with Dirichlet boundary conditions. The resemblance to (41), together with (9) and (46), is clear. However, the presence of the rapidly varying, not necessarily low amplitude, height field, without any intrinsic correlations that might perhaps smooth it out, drastically effects G h . In particular, it is not smooth and hence strongly violates the conditions under which the mean field approximation described in Section 4.1 is derived. One may think of G h as generating a Coulomb-type interaction between vortices that retains a strong equilibrium fluctuation on finite length scales. Moreover, the 1 2 gh 2 term makes sense only at positive temperatures. Similar to the q field in the 3D axisymmetric model, height fluctuations absorb unbounded energy for increasing T, and are hence in principle capable of dissipating negative temperature-like vortex states and converting large scale vortex motion into height fluctuations.

Quasi-Hydrostatic Shallow Water Equilibria
There are, however, physical motivations, completely outside of equilibrium considerations, for seeking equilibria with smooth height fields. Thus, a forward-type cascade of high amplitude, small scale height fluctuations will eventually violate the long wavelength assumption entering the derivation of the shallow water equations. When these assumptions are violated the full 3D Euler equations will display shock wave formation, wave breaking, and other 3D motions that will serve to effectively dissipate strong wave motions without significantly impacting large scale eddy motions.
Interesting work for the future would be more careful investigations of the validity of such alternative routes to equilibrium. For now let us briefly explore the consequences. If the height is smooth on the scale of variation of Ω(r), then G h is smooth and, following steps analogous to the functional Taylor expansion (48), one obtains the mean field approximation in which the shallow water stream function Ψ associated with the mass current hv (which is indeed incompressible in equilibrium, and differs from the velocity stream function ψ introduced earlier), is defined by leading to Using (125) and (128) one may expresŝ The fully fluctuating field Ω now appears only in the final local term in (126), and one may now integrate it out to obtain the shallow water Free energy functional generalizing (56): in which, generalizing (54) and (55), we define where the h 2 prefactor comes from the phase space measure [the original h 4 i in (122) is reduced to h 2 i after performing the Q integral]. We again observe the required scaling β a = β/a 2 to obtain a finite result in the continuum limit.
The self-consistent equation for Ψ 0 (r) is obtained from the extremum condition δF /δΨ 0 (r) = 0, which yields Similar to Euler result (58), the self-consistency condition equates the mean vorticity derived from the equilibrium stream function [first line of (132)] with that computed from the local distribution function [second line of (132)], here emerging as a certain function W of Ψ 0 controlled by the Lagrange multipliers β, µ(σ).
The equation for h is obtained by applying the extremum condition δF /δh(r) = 0, |∇Ψ 0 (r)| 2 2h(r) 2 + gh(r) ≡ 1 2 |v 0 (r)| 2 + gh(r) which corresponds to the reasonable assumption that the dissipation process self consistently acts to minimize the free energy. This is formally correct for large β where height fluctuations are indeed small. Thus, more formally, the self-consistency requirement is that the dissipation process produces a new effectively low temperature system. The Lagrange multipliers µ(σ) will change as well so as to enforce approximately the same g(σ)-to the extent the large scale eddy degrees of freedom are unaffected by the high frequency wave suppression.
In the large β limit one can show that ∂ h W(h, Ψ 0 ) W 1 (Ψ 0 ) is independent of h. It follows then that (133) expresses the Bernoulli condition, namely that the sum of local kinetic energy and pressure is constant along stream lines (level curves of Ψ 0 ). This is indeed a rigorous requirement for steady flows. More generally, one may continue to apply (133) for moderate values of β as an approximate model in which some fluctuations in h are kept (and the Bernoulli condition is weakly violated).
Another interesting consequence is that, accepting (130) as an approximate free energy, negative temperatures are no longer precluded. Thus, W(h, τ) is perfectly well defined for β < 0 and solutions to (133) may be sought for both positive and negative β. As previously stated, negative temperature equilibria are formally unstable to leakage of energy into (positive-temperature) wave motions, but the physical coupling of large-scale flows to small-scale wave generation is extremely weak and it makes sense to develop a theory along these lines that neglects such effects. The key observation here is that compact eddy structures, such as Jupiter's Great Red Spot, having vorticity maxima confined away from the system boundaries, can only be interpreted as negative-temperature states. Such structures therefore lie outside the strict shallow water theory presented here and nonequilibrium dissipation arguments must therefore be invoked in order to make contact with the effective equilibrium descriptions ubiquitous in the literature [1].
The result (130) reduces to the Euler equation result (56) if one constrains h(r) = H 0 . As one relaxes this constraint the vorticity pattern will evolve somewhat to accommodate the sloping surface in response to quasi-hydrostatic force balance, as observed in [42,43]. However, one does not expect major changes from the 2D Euler result unless one drives the system to extremely high vorticity gradients, which is typically not of geophysical relevance.

2D Magnetohydrodynamics and the Solar Tachocline
Our final example is that of the ideal, perfectly electrically conducting fluid (relevant to the energy conserving limit) interacting with an external magnetic field. The effective 2D theory of interest here emerges as follows. We begin with the 3D magnetohydrodynamic (MHD) equations where Ω is the rotation vector. The right hand side of the first equation now includes, in addition to the pressure term, the Lorentz force acting on a parcel of fluid. The second equation is Faraday's law with electric field determined by the constraint E + v × B = 0, which zeroes out the net force on the charge imposed by the perfectly conducting limit. The equations are closed using Ampere's law J = ∇ × B, and the generalized pressure P (which includes also contributions from centrifugal force, gravity, etc.) continues to enforce the incompressibility condition (2). The constraint ∇ · B = 0 is automatically enforced by the second equation. The solar tachocline (an example simulation result for which is shown in Figure 9) is the observed sharp radial boundary between the solid body rotating radiative interior and differentially rotating outer convective zone. Here the current J = Jẑ passes normally through the surface, while v and B are in-plane. The incompressibility conditions allow one to define the stream function (3) together with the (z component of the) magnetic vector potential The vector Equation (134) now reduce to the pair of scalar equations with f = 2Ω sin(ϕ) defined by the solar latitude ϕ. The kinetic plus electromagnetic energy is conserved if both ψ and A obey Dirichlet (free slip) boundary conditions. For annular or periodic strip geometries the angular momentum (12) is conserved, and can alternatively be written in the form P = − D dr∇α · ∇ψ.
(138) Figure 9. Example numerically generated long-time (near-equilibrium) behavior of freely decaying 2D magnetohydrodynamics on the sphere. The zonal velocity field (above) and zonal magnetic field (below) undergo coupled dynamics according to (134), reducing to (135) and (136) in the 2D solar tachocline model [34,44,45]. In the (ψ, A) representation (141) of the statistical functional, where v is defined by the level curves of ψ and B is defined by the level curves of A, the model is that of two gradient-coupled membranes in an external potential which, among other things, tends to preferentially align the two vector fields.
It is immediately evident that the Casimirs are completely different here since the potential vorticity ω + f is no longer advectively conserved. Instead it is the potential A that is conserved, which has the drastic effect of imposing no direct control on the second derivative J. In fact, similar to the axisymmetric case (105), there are two sets of Casimirs with the second following from the fact that JB is orthogonal to ∇A. Dynamically, if B · ∇J happens to be small, one may expect to observe gradual evolution from Euler-type large scale eddy states to the quite different equilibria based on the vector potential Casimirs. The latter in particular permit diverging small scale vorticity fluctuations, as exhibited below. The presence of even a weak magnetic field in 2D MHD simulations has indeed been found to destroy the conventional inverse cascade, breaking up large scale eddy flows [44]. However, as described below, the new set of Casimirs (139) are also capable of generating large scale flows, but based on significantly different initial states with imposed structure on A rather than on ω.
A Liouville theorem may straightforwardly be proven for the pair ω, A so that the equilibrium phase space measure is in which J a is the Jacobian associated with the change of variable ω → ψ. This simply adds a constant to the free energy and drops out of any statistical average. Defining corresponding Lagrange multipliers µ(σ) andμ(σ), we consider then the grand canonical statistical functional in which integration by parts has been used to express everything in terms of at most first order gradients of the fields. This form is also displayed in the fourth row of Table 1 (with f again dropped for simplicity).
The physical model associated with ρ eq = Z −1 GC e −β a K is that of two membranes with "heights" A(r), ψ(r) and unit surface tension (the coefficient of the gradient-squared terms), and additionally coupled through their gradients. The term µ(A) + f (r)μ(A) is a smoothly position-dependent external potential, confining A near its minimum. The ψ membrane is not directly confined, but the gradient coupling favors B parallel toμ (A)v + λ∇ × α.
Using the scaling β a = β/a 2 one sees that the membrane experiences local Brownianlike fluctuations, with neighboring height differences scaling as a/ β. It follows that A and ψ are continuous, but have randomly fluctuating gradient, so that v and B fluctuate from site to site with scale 1/ β. The membranes are therefore globally smooth but microscopically rough. In fact one may make use of this separation of scales to write in which A 0 , ψ 0 are the equilibrium averages, to be determined self-consistently below, and δA, δψ = O(a/ β) are fluctuation corrections. Substituting these into (141) one obtains in which δφ ± = (δA ± δψ)/ √ 2 are independent Gaussian fields. For smooth A 0 , ψ 0 , all other terms, including those linear in δA, δψ, vanish with a → 0. The major complication here is that the coefficient ν [A 0 (r)] is not only position dependent, but yet to be determined.
The free energy functional follows from (142) in the form in which the Gaussian correction is defined by and has a well defined continuum limit. The equilibrium equations, obtained from δK 0 /δψ 0 (r) = 0, δF /δA 0 (r) = 0 yield, respectively is the magnetic-velocity Gaussian correlation function. The first equation provides a direct relation between the equilibrium velocity and magnetic field, being collinear up to a momentum conservation-induced mean flow subtraction-this is the effect of the gradient coupling term in (141). By substituting the curl of this relation into the second equation, it is straightforward to derive a closed equation for A 0 alone. These equations look quite complicated, but have a straightforward physical interpretation. The functional K[A 0 , ψ 0 ] reflects a classical surface tension minimization problem in the presence of the external potentials µ,μ. The second derivative term ω 0 + ν (A 0 )J 0 in the second equation represents a surface tension restoring force in response to the forcing terms on the right. The γ term represents the fluctuation corrections to the surface tension due to the membrane roughening effect. The self-consistent dependence on A 0 , ψ 0 arises from such effects as regions of strongly stretched membrane having reduced amplitude fluctuations. Example solutions of these equations, displaying similar large scale vortex flow patterns as the Euler equations, are shown in Ref. [48].
All equilibrium conserved integrals are derived as usual through differentiation with respect to the Lagrange multipliers: in which ε(r) is another microscale Gaussian fluctuation correction that may also be written terms of pair correlation functions [48]. Note that due to continuity of A, g(σ) is a large scale quantity, i.e., its own equilibrium average. Hence the level sets of a given initial condition A(r, t = 0) are exactly preserved (though perhaps significantly contorted) in the equilibrium function A 0 (r). On the other hand, due to strong (unbounded) fluctuations of ω, a microscale correction tog(σ) is evident. The physically observable fields are the membrane gradients B and v. Depending on the initial condition, their fluctuations, though bounded from point to point, could still be large compared to their mean values. This is physically quite different from the Euler equation where the second derivative has bounded fluctuations and the gradients are smooth. This has implications for the effects of dissipation which could be much stronger in this system, quelling micro-fluctuations and perhaps more rapidly bleeding energy out the large scale flow. The appearance (or not) of macroscale magnetic structure in the solar tachocline has significant implications for angular momentum transport between the two zones that it separates [44].
The example simulation result shown in Figure 9 is not intended as an equilibrium theory comparison-this will require future more careful study. However, it does verify that large scale magnetic field structures can survive for a long time even as the vorticity structure becomes much more diffuse. For this particular case the magnetic field magnitude is only weakly changed from its initial condition (not shown) while lack of vorticity conservation allows the zonal velocity magnitude to drop by nearly an order of magnitude.

Conclusions
In this article we have discussed the application of statistical mechanics to the characterization of certain classes of large scale 2D steady state flows, following, for example, the free decay of an initial turbulent state (Figure 3), highlighting the role of the competition between flow energy and microscale entropy production. The thermodynamic formalism makes sense only for systems whose dynamics is governed by a conserved Hamiltonian. When applied to fluid dynamics this limits consideration to idealized flows in which all dissipative terms are dropped. This, at minimum, limits the applicability to high Reynolds number flows with a large separation of scales between outer scale inertial, energy conserving dynamics, and small scaling mixing that eventually encounters viscous dissipation. With this separation, one may propose that the idealized models may provide reasonable predictions over an intermediate range of time scales that include a sufficient degree of intermediate scale equilibration. This is especially interesting in two dimensional models, where one encounters an infinite number of conserved integrals of the motion (Casimirs) beyond the standard total energy and momentum. These strongly constrain the flow and in cases of interest lead to the phenomenon of an inverse cascade of energy to large scales, balanced by an "enstrophy cascade" to smaller scales, namely a fine-scale mixing of low energy eddies (Figure 4). In a finite domain, the inverse cascade "condenses" into a system scale steady state structure. The goal of the thermodynamic treatment is to predict such structures based only on the values of the conserved integrals imposed by the initial flow-the only quantities "visible" to the statistical formalism. Given the very large number of such integrals, there are potentially many different large scale flow patterns that might be accessed, exemplified by long lived eddies such as Jupiter's Great Red Spot, zonal jet features, etc.
Following the classic construction of the statistical formalism (dynamics in phase space, ergodic hypothesis, Liouville theorem, invariant measures, choice of ensemble), the problem may be reduced to the analysis of a classical field theory (Table 1), with analogies to continuous spin Ising models (perhaps interacting with additional Gaussian degrees of freedom), and interacting elastic membrane models, depending on the exact problem and the fluid degree of freedom to which the Casimirs are applied. The fluid physics, however, drives these models into unusual regimes, e.g., of very high energy (negative temperatures) that are not normally encountered in more conventional versions of these models (Figures 2 and 6). In these regimes we have seen that the statistical approach, in the form of a thermodynamic free energy variational principle, is indeed able to produce the desired macroscopic flows. The formalism additionally lends insight into the role of the various conserved integrals in controlling the geometry of these flows. Simple examples for the 2D Euler equation are shown in Figure 5.
Despite the mathematical elegance of the theory and its predictions, there remain numerous questions regarding the validity of the underlying assumptions, especially the ergodic hypothesis and the convergence to a true equilibrium state [1]. In comparison to conventional particle systems, there are many possible barriers to equilibration, including extra adiabatic invariants (Section 6.3), metastable equilibria [21], and very long-lived chaotic states (Section 7.5). Some of these are well understood, others deserve more careful study.
There are also systems for which the equilibrium theory apparently works too well! Thus, the inclusion of additional physical degrees of freedom intended to make the model more physically realistic, such as surface motions in the shallow water Equations (Section 8), in principle destabilizes negative temperatures states, leading to an ultraviolet catastrophe of surface waves despite the Casimir constraints. In fact, long-lived planetary eddies are much more in line with predictions of the much simpler Euler or quasigeostrophic Equation [1]. Similar issues are seen in axisymmetric flows (Section 7) where an ultraviolet catastrophe of poloidal vorticity predicts only rather trivial large scale toroidal flows. In both cases the catastrophic coupling of the new small scale fluctuations to existing large scale structures is likely very weak, and high frequency wave or poloidal vorticity generation might better be thought of as an additional weak dissipation mechanism that can also be ignored over time scales of interest. The resulting quasi-hydrostatic limit of the shallow water equations provides one possible route to formally maintaining negative temperature states while still treating the surface height in a consistent manner (Section 8.4).
The previous discussion motivates a number of future investigations into a more careful treatment of additional dynamical time-and length-scale separations that could either hinder or aid statistical equilibrium approximations, and how to properly define the effective conserved integrals entering a new idealized flow model, e.g., through an appropriate spatial filter.
In addition, it is clear that very long-lived eddies, such as Jupiter's Red Spot, require some sort of driving force to survive. The weak driving-weak dissipation limit could perhaps be formulated through convergence to a near-equilibrium state in which the conserved integrals come into detailed balance, e.g., through some kind of Onsager nonequilibrium linear response theory applied to the fluid Hamiltonian. On the other hand, it is known that weak stochastic forcing can occasionally lead to rare, sudden, catastrophic changes to the state [52,53] so some care must be taken in finding the correct regime in which to formulate the problem.