Aurore Naso

Simon Thalabard

Gilles Collette

Pierre-Henri Chavanis

Bérengère Dubrulle

Statistical mechanics of Beltrami flows in axisymmetric geometry: Equilibria and bifurcations

Keywords:

We characterize the thermodynamical equilibrium states of axisymmetric Euler-Beltrami flows. They have the form of coherent structures presenting one or several cells. We find the relevant control parameters and derive the corresponding equations of state. We prove the coexistence of several equilibrium states for a given value of the control parameter like in 2D turbulence [Chavanis & Sommeria, J. Fluid Mech. 314, 267 (1996)]. We explore the stability of these equilibrium states and show that all states are saddle points of entropy and can, in principle, be destabilized by a perturbation with a larger wavenumber, resulting in a structure at the smallest available scale. This mechanism is therefore reminiscent of the 3D Richardson energy cascade towards smaller and smaller scales. Therefore, our system is truly intermediate between 2D turbulence (coherent structures) and 3D turbulence (energy cascade). We further explore numerically the robustness of the equilibrium states with respect to random perturbations using a relaxation algorithm in both canonical and microcanonical ensembles. We show that saddle points of entropy can be very robust and therefore play a role in the dynamics. We evidence differences in the robustness of the solutions in the canonical and microcanonical ensembles. A scenario of bifurcation between two different equilibria (with one or two cells) is proposed and discussed in connection with a recent observation of a turbulent bifurcation in a von Kármán experiment [Ravelet et al., Phys. Rev. Lett. 93, 164501 (2004)].

I. INTRODUCTION

The statistical mechanics of systems with long-range interactions has recently attracted a lot of attention [START_REF] Dauxois | Dynamics and Thermodynamics of Systems with Long-Range Interactions[END_REF]. Typical systems with long-range interactions include self-gravitating systems [2,3], two-dimensional vortices [START_REF] Chavanis | Statistical mechanics of two-dimensional vortices and stellar systems[END_REF], non-neutral plasmas [START_REF] Dubin | [END_REF], free electrons lasers [6] and toy models such as the Hamiltonian Mean Field (HMF) model [7,8]. Unusual properties of systems with long-range interaction such as negative specific heats or ensembles inequivalence have been evidenced and linked with lack of additivity [9]. In addition, a striking property of these systems is the rapid formation of quasi stationary self-organized states (coherent structures) such as galaxies in the universe [START_REF] Binney | Galactic Dynamics[END_REF], large scale vortices in geophysical and astrophysical flows [11][12][13][14][15] or quasi-stationary states in the HMF model [16][17][18]. These QSSs can be explained in terms of statistical mechanics using the theory developed by Lynden-Bell [19] for the Vlasov equation or by Miller [20] and Robert & Sommeria [21] for the 2D Euler equation.

Two-dimensional vortices interact via a logarithmic potential. Interaction of vortices in 3D turbulence is weaker than in 2D turbulence, but still long-range. Due to dissipative anomaly and vortex stretching, statistical mechanics of 3D turbulence has so far eluded theories. Recent progress was recently made considering 3D inviscid axisymmetric flows [22,23] that are intermediate between 2D and 3D flows: they are subject to vortex stretching like in 3D turbulence, but locally conserve a scalar quantity in the ideal limit, like in 2D turbulence. It is therefore interesting to study whether these systems obey the peculiarities observed in other systems with long-range interactions such as violent relaxation, existence of long-lived quasi-stationary states, negative specific heats and ensembles inequivalence.

The general study of the stability of axisymmetric flows, and the possible occurrence of phase transitions, is difficult due to the presence of an infinite number of Casimir invariants linked with the axisymmetry of the flow. In a previous paper [23], hereafter Paper I, we have considered a simplified axisymmetric Euler system characterized by only three conserved quantities: the fine-grained energy E f.g. , the helicity H and the angular momentum I. We have developed the corresponding statistical mechanics and shown that equilibrium states of this system have the form of Beltrami mean flows on which are superimposed Gaussian fluctuations. We have shown that the maximization of entropy S at fixed helicity H, angular momentum I and microscopic energy E f.g. (microcanonical ensemble) is equivalent to the maximization of free energy J = S -βE f.g. at fixed helicity H and angular momentum I (canonical ensemble). These variational principles are also equivalent to the minimization of macroscopic energy E c.g. at fixed helicity H and angular momentum I. This provides a justification of the minimum energy principle (selective decay) from statistical mechanics. We have furthermore discussed the analogy with the simplified thermodynamical approach of 2D turbulence developed in [24] based on only three conserved quantities: the fine-grained enstrophy Γ f.g.

2 , the energy E and the circulation Γ. We have shown that equilibrium states of this system have the form of Beltrami mean flows (linear vorticity-stream function relationship) on which are superimposed Gaussian fluctuations. We have shown that the maximization of entropy S at fixed energy E, circulation Γ and microscopic enstrophy Γ f.g.

2

(microcanonical ensemble) is equivalent to the maximization of grand potential S = Sα 2 Γ f.g.

2

at fixed energy E and circulation Γ (grand microcanonical ensemble). These variational principles are also equivalent to the minimization of macroscopic enstrophy Γ c.g.

2

at fixed energy E and circulation Γ. This provides a justification of a minimum enstrophy principle (selective decay) from statistical mechanics. In the analogy between 2D turbulence and 3D axisymmetric turbulence, the energy plays the role of the enstrophy.

In the present paper, we study more closely the equilibrium states of axisymmetric flows and explore their stability. We show that all critical points of macroscopic energy at fixed helicity and angular momentum are saddle points, so that they are unstable in a strict sense. Indeed, there is no minimum (macroscopic) energy state at fixed helicity and angular momentum (either globally or locally) because we can always decrease the energy by considering a perturbation at smaller scales. This is reminiscent of the Richardson energy cascade in 3D turbulence. Inversely, in 2D turbulence, there exists minimum enstrophy states that develop at large scales (inverse cascade). Therefore, our system is intermediate between 2D and 3D turbulence: there exists equilibrium states in the form of coherent structures (that are solutions of a mean field differential equation) like in 2D turbulence, but they are saddle points of macroscopic energy and are expected to cascade towards smaller and smaller scales like in 3D turbulence. However, we give arguments showing that saddle points can be robust in practice and play a role in the dynamics. Indeed, they are unstable only for some particular (optimal) perturbations and can persist for a long time if the system does not spontaneously generate these perturbations. Therefore, these large-scale coherent structures can play a role in the dynamics and they have indeed been observed in experiments of von Kármán flows [25]. In order to make this idea more precise, we have explored their stability numerically using phenomenological relaxation equations derived in [23]. We have found some domains of robustness in the parameter space. In particular, the one cell structure is highly robust for large values of the angular momentum I > I c and becomes weakly robust for low values of the angular momentum. In that case, we expect a phase transition (bifurcation) from the one-cell structure to the two-cells structure. We have also found that the value of the critical angular momentum I c changes depending whether we use relaxation equations associated with a canonical (fixed temperature) or microcanonical (fixed microscopic energy) description. At low temperatures T , we have evidence a new kind of "ensembles inequivalence" characterizing the robustness of saddle points with respect to random perturbations.

The paper is organized as follows: In Sec. II, we set-up the various notations and hypotheses we are going to use. The computation and characterization of equilibrium states is done in Sec. III. The stability analysis of these equilibrium states is performed in Sec. III G where we show analytically that all states are unstable with respect to large wavenumber perturbations. We evidence a process of energy condensation at small scales that is reminiscent of the Richardson cascade. We explore numerically the robustness of the equilibria in both canonical and microcanonical ensembles in Sec. IV. Our numerical method is probabilistic and rather involved. A discussion of our results is done in Sec. V where a bifurcation scenario relevant to the turbulent experimental von Kármán flow is suggested.

II. THEORETICAL SET-UP

A. Hypotheses and Notations

We consider a system with a cylindrical geometry enclosed in the volume delimited above and below by surfaces z = 0 and z = 2h, and radially by 0 ≤ r ≤ R. Like in Paper I, we consider an axisymmetric Euler-Beltrami system characterized by a velocity field u, with axisymmetric time averaged u. We furthermore assume that the only relevant invariants of the axisymmetric Euler equations for our problem are the averaged energy E = 1 2 u 2 dr, the averaged helicity H = u • ω dr and the averaged angular momentum I = σ dr where σ = ru θ . We introduce the potential vorticity ξ = ω θ /r and the stream function ψ such that u r = -r -1 ∂ z ψ and u z = r -1 ∂ r ψ. They are related to each other by the generalized Laplacian operator

∆ * ψ ≡ 1 r 2 ∂ 2 ψ ∂z 2 + 1 r ∂ ∂r 1 r ∂ψ ∂r = -ξ. (1)
In actual turbulent von Kármán experiments, we have been able to observe that the largest part of the kinetic energy is contained in the toroidal motions. It is therefore natural, as a first elementary step, to consider a model in which only toroidal fluctuations are considered, and suppose that the fluctuations in the other (poloidal) directions are simply frozen. With such an assumption, poloidal vorticity fluctuations are allowed, but toroidal vorticity fluctuations are excluded. We therefore only include a fraction of the vorticity fluctuations, that presumably become predominant at small scale, due to the existence of vortex stretching. As shown below and in the next paper [START_REF] Dubrulle | Paper III: Comparison with experiments[END_REF], this simplification however still allows for vortex stretching and energy cascades towards smaller scales, and leads to predictions that are in good agreement with experiments. Moreover, our hypotheses lead to a model that is self-contained and analytically tractable. According to our hypotheses, neither ξ nor ψ fluctuates in time: ξ = ξ and ψ = ψ. In that case, the conserved quantities can be rewritten

E f.g. = 1 2 ξψ + 1 2 σ 2 r 2 , (2)
H = ξ σ , (3)
I = σ , (4)
where f denotes the spatial average [32]

f ≡ 1 hR 2 R 0 2h 0 rdrdz f (r, z). (5)
The helicity and the angular momentum are robust constraints because they can be expressed in terms of coarsegrained quantities ξ and σ. By contrast, the energy is a fragile constraint because it cannot be expressed in terms of coarse-grained quantities. Indeed, it involves the fluctuations of angular momentum σ 2 . To emphasize that point, we have introduced the notation E ≡ E f.g. to designate the fine-grained (microscopic) energy. Splitting σ into a mean part σ and a fluctuating part δσ, we define the coarse-grained (macroscopic) energy by

E c.g. = 1 2 ξψ + 1 2 σ 2 r 2 . (6)
Then, the energy contained in the fluctuations is simply

E f luct ≡ E f.g. -E c.g. = 1 2 σ 2 r 2 , (7)
where

σ 2 ≡ σ 2 -σ 2 , (8)
is the local centered variance of angular momentum. We stress that the microscopic energy E = E f.g. is conserved while the macroscopic energy E c.g. is not conserved and is likely to decrease (see below).

B. The two statistical ensembles and the selective decay principle

In Paper I, we have developed a simplified thermodynamic approach of axisymmetric flows under the abovementioned hypothesis. Let ρ(r, η) denote the PDF of σ and let us recall the expression of the entropy

S [ρ] = -ρ ln ρ drdη. (9)
We have proven the equivalence between the microcanonical ensemble max ρ,ξ

{S[ρ] | E f.g. , H, I, ρdη = 1}, (10)
and the canonical ensemble max ρ,ξ

{J[ρ] = S -βE f.g. | H, I, ρdη = 1}. (11)
In each ensemble, the critical points are determined by the first order condition δS -βδE f.g. -µδH -αδI = 0. The equilibrium distribution is Gaussian

ρ(r, η) = β 2πr 2 1/2 e -β 2r 2 (η-σ) 2 , ( 12
)
the mean flow is a Beltrami state

σ = - β µ ψ, (13)
ξ = - βσ µr 2 - α µ , (14)
and the centered variance of angular momentum is

σ 2 = r 2 β . ( 15
)
These equations determine critical points of the variational problems [START_REF] Binney | Galactic Dynamics[END_REF] and ( 11) that cancel the first order variations of the thermodynamical potential. Clearly, [START_REF] Binney | Galactic Dynamics[END_REF] and (11) have the same critical points. Furthermore, it is shown in Paper I that [START_REF] Binney | Galactic Dynamics[END_REF] and (11) are equivalent for the maximization problem linked with the sign of the second order variations of the thermodynamical potential: a critical point determined by Eqs. ( 12)-( 15) is a maximum of S at fixed microscopic energy, helicity and angular momentum iff it is a maximum of J at fixed helicity and angular momentum. This equivalence is not generic. We always have the implication (11) ⇒ (10) but the reciprocal may be wrong. Here, the microcanonical and canonical ensembles are equivalent due to the quadratic nature of the microscopic energy ∼ σ 2 . We note that, according to Eq. ( 15), β is positive. In the canonical ensemble, β is prescribed. In the microcanonical ensemble, β is a Lagrange multiplier that must be related to the energy E = E f.g. . According to Eqs. ( 7) and ( 15), we find that β > 0 is determined by the condition

E f.g. -E c.g. = E f luct = 1 2β . ( 16
)
This relation shows that T = 1/β plays the role of a temperature associated with the fluctuations of angular momentum [33]. Finally, we have proven in Paper I that the two variational problems [START_REF] Binney | Galactic Dynamics[END_REF] and (11) are equivalent to max

σ,ξ {J[σ, ξ] = -βE c.g. | H, I}, (17)
or equivalently min

σ,ξ {E c.g. [σ, ξ] | H, I}, (18)
in the sense that the solution of [START_REF] Binney | Galactic Dynamics[END_REF] or (11) is given by Eq. ( 12) where (σ, ξ) are the solutions of ( 17) or (18). This justifies a selective decay principle from statistical mechanics. Indeed, it is often argued that an axisymmetric turbulent flow should evolve so as to minimize energy at fixed helicity and angular momentum. In general, this phenomenological principle is motivated by viscosity or other dissipative processes. In our approach, it is justified by the maximum entropy principle [START_REF] Binney | Galactic Dynamics[END_REF] of statistical mechanics when a coarse-graining is introduced. In the sequel, we shall study the maximization problem (18) since it is simpler than [START_REF] Binney | Galactic Dynamics[END_REF] or (11), albeit equivalent. Remark: although the variational problems [START_REF] Binney | Galactic Dynamics[END_REF] and (11) determining equilibrium states are equivalent, this does not mean that the relaxation equations associated with these variational problems are equivalent. To take an analogy, the Boltzmann (microcanonical) and the Kramers (canonical) equations have the same equilibrium states -the Maxwell distribution-but a different dynamics. In the following, we will show that the equilibrium variational problems [START_REF] Binney | Galactic Dynamics[END_REF] and (11) have no solution. Indeed, there is no maximum of entropy at fixed E f.g. , H and I and no minimum of free energy at fixed H and I. All the critical points of [START_REF] Binney | Galactic Dynamics[END_REF] and (11) are saddle points of the thermodynamical potentials. Then, the idea is to consider the out-of-equilibrium problem, introduce relaxation equations and study the robustness of saddle points with respect to random perturbations. For what concerns the out-of-equilibrium problem, the microcanonical and canonical ensembles may be inequivalent. We will see that they are indeed inequivalent.

III. MINIMUM ENERGY STATES

A. Critical points

In this section, we shall study the minimization problem min

σ,ξ {E c.g. [σ, ξ] | H, I}. (19)
The critical points of macroscopic energy at fixed helicity and angular momentum are determined by the condition

δE c.g. + µδH + αδI = 0, (20)
where µ (helical potential) and α (chemical potential) are Lagrange multipliers. Introducing the notations B = -1/µ and C = -α/µ, the variations on δξ and δσ lead to

σ = Bψ, (21)
ξ = B σ r 2 + C, (22)
which are equivalent to Eqs. ( 13)-( 14) up to a change of notations. In the following, it will be convenient to work with the new field φ = ψ/r. It is easy to check that

∆ * ψ = 1 r ∆φ - φ r 2 , ( 23
)
where ∆ is the usual Laplacian. Therefore, Eq. ( 1) becomes

-∆φ + φ r 2 = rξ, (24)
and the previous equations can be rewritten

σ = Brφ, (25)
ξ = B 2 r φ + C, (26)
where φ is solution of

-∆φ + φ r 2 = B 2 φ + Cr, (27)
with φ = 0 on the boundary. This is the fundamental differential equation of the problem. Note that a particular solution of this differential equation is

φ part = - C B 2 r, (28)
but it does not satisfy the boundary conditions. Using Eqs. (3) and (4), the helicity and the angular momentum are given by

H -CI = B 3 φ 2 , (29)
I = B φr . (30)
These equations are relationships between (B, C) and (H, I). Remark: we have not taken into account the conservation of circulation Γ = ξ because this would lead to a term AB/r in the r.h.s. of Eq. ( 27) that diverges as r → 0.

B. The different solutions

To construct the different solutions of Eq. ( 27) and study their stability, we shall follow the general procedure developed by Chavanis & Sommeria [26] for the 2D Euler equation. We first introduce an eigenmode decomposition to compute all critical points of (19). Then, we investigate their stability by determining whether they are (local) minima of macroscopic energy or saddle points.

The eigenmodes

We first assume that

C = 0. ( 31
)
In that case, the differential equation ( 27) becomes

-∆φ + φ r 2 = B 2 φ, (32)
with φ = 0 on the domain boundary. We introduce the eigenfunctions φ mn of the operator L ≡ -∆ + 1 r 2 . They are defined by

Lφ mn ≡ -∆φ mn + φ mn r 2 = B 2 mn φ mn , (33)
with φ mn = 0 on the domain boundary. It is easy to show that the eigenvalues Λ mn of L are positive (hence the notation Λ mn = B 2 mn ). Indeed, we have φ mn Lφ mn = Λ mn φ 2 mn and φ mn Lφ mn = (∇φ mn

) 2 + φ 2 mn r 2
≥ 0, which proves the result. It is also easy to show that the eigenfunctions are orthogonal with respect to the scalar product

f g ≡ 1 hR 2 R 0 2h 0 rdrdz f (r, z)g(r, z). (34)
Finally, we normalize them so that φ mn φ m ′ n ′ = δ mm ′ δ nn ′ . The eigenvalues and eigenfunctions of the operator L can be determined analytically. The differential equation (33) can be rewritten

∂ 2 φ ∂r 2 + 1 r ∂φ ∂r + ∂ 2 φ ∂z 2 - φ r 2 = -B 2 φ. (35)
We look for solutions in the form φ(r, z) = f (r)g(z). This yields

f ′′ f + 1 r f ′ f - 1 r 2 + B 2 = - g ′′ g ≡ κ 2 , (36)
where the sign of the constant has been chosen in order to satisfy the boundary condition φ = 0 in z = 0 and z = 2h. The differential equation for g is readily solved and we obtain

g(z) = sin(κz), (37)
with

κ = n π 2h , ( 38
)
where n is a strictly positive integer. On the other hand, the differential equation for f is

r 2 f ′′ + rf ′ -r 2 (κ 2 -B 2 )f -f = 0. ( 39
)
If we define λ 2 = B 2κ 2 and x = λr, the foregoing equation can be rewritten

x 2 f ′′ + xf ′ + x 2 f -f = 0. ( 40
)
This is a Bessel equation whose solution is

f = J 1 (x). (41)
Now, the boundary condition φ(R) = 0 implies f (λR) = 0 so that

λR = j 1m , (42)
where j 1m is the m-th zero of Bessel function J 1 . In conclusion, the eigenvalues are

B mn 2 = j 1m R 2 + nπ 2h 2 , (43)
and the eigenfunctions are

φ mn = N mn J 1 j 1m r R sin nπz 2h , (44)
with the normalization constant

N mn = 2 J 2 2 (j 1m ) . ( 45
)
The mode (m, n) corresponds to m cells in the r-direction and n cells in the z-direction. We shall distinguish two kinds of modes, according to their properties regarding the symmetry R with respect to the plane z = h. The odd eigenmodes denoted φ ′ mn are such that Rφ ′ mn = -φ ′ mn and correspond to n even. They have zero mean value in the z direction ( 2h 0 φ ′ mn dz = 0). For example, the mode (1, 2) is a two-cells solution in the vertical direction. The even eigenmodes denoted φ ′′ mn are such that Rφ ′′ mn = φ ′′ mn and correspond to n odd. They have non zero mean value in the vertical direction ( 2h 0 φ ′′ mn dz = 0). In particular, the mode (1, 1) is a one-cell solution. Returning to Eq. ( 32), this differential equation has solutions only for quantized values of B = B mn (eigenvalues) and the corresponding solutions (eigenfunctions) are

φ = H B 3 mn 1/2 φ mn , (46)
where we have used the helicity constraint (29) to determine the normalization constant. Note that Eq. ( 29) implies that B mn and H have the same sign, so that the square root is always defined. Substituting this result in Eq. ( 30), and introducing the control parameter

Λ = I 2 H , (47)
we find that these solutions exist only for Λ = Λ mn with

Λ mn = 1 B mn φ mn r 2 . ( 48
)
For the odd eigenmodes φ ′ mn , we have Λ = 0 and for the even eigenmodes φ ′′ mn , we have Λ ′′ mn = 1

B ′′ mn φ ′′ mn r 2 .

The continuum

We now assume that C = 0 and define

φ B = φ C . (49)
In that case, the fundamental differential equation ( 27) becomes

-∆φ B + φ B r 2 = B 2 φ B + r, (50)
with φ B = 0 on the domain boundary. We also assume that B = B mn . In that case, Eq. ( 50) admits a unique solution that can be obtained by expanding φ B on the eigenmodes. Using the identity

f = mn f φ mn φ mn , (51)
we get

φ B = mn φ ′′ mn r B ′′ mn 2 -B 2 φ ′′ mn . (52)
Of course, φ B can also be obtained by solving the differential equation ( 27) numerically. Note that this solution is even since only the even modes are "excited". Substituting Eq. (49) in Eq. ( 30), we obtain

C = I B φ B r . ( 53
)
Then, substituting Eqs. ( 49) and (53) in Eq. ( 29), we get

Λ = B φ B r 2 φ B r + B 2 φ 2 B . ( 54
)
This equation gives a relationship between B and Λ. Then, C is determined by Eq. ( 53). These equations can therefore be viewed as the equations of state of the system. They determine the branch formed by the solutions of the continuum. Using

φ B r = mn φ ′′ mn r 2 B ′′ mn 2 -B 2 , (55)
and

φ 2 B = mn φ ′′ mn r 2 (B ′′ mn 2 -B 2 ) 2 , (56)
we obtain

φ B r + B 2 φ 2 B = mn (B ′′ mn ) 2 φ ′′ mn r 2 (B ′′ mn 2 -B 2 ) 2 > 0. ( 57
)
This implies that B is of the same sign as Λ, hence H. Furthermore, Λ is an odd function of B. In the sequel, we shall consider only cases with H ≥ 0, i.e. B ≥ 0 and Λ ≥ 0 for illustration and figures. Note that Eq. (54) involves the important function

F (B) ≡ φ B r = mn φ ′′ mn r 2 B ′′ mn 2 -B 2 . ( 58
)
For Λ = 0, the inverse helical potential is

B = 0 or B = B (n) * where B (n) * is any zero of F , i.e. F (B (n) * ) = φ B (n) * r = 0. ( 59
)
For simplicity, we shall call B * = B

(1) * the first zero of F . This first zero is always between the first and the second even eigenmodes (see Appendix A). Its location with respect to the first odd eigenmode B ′ 12 depends on the aspect ratio of the cylinder: for h/R > 0.53, we have B * > B ′ 12 (case L-for Large aspect ratio) while for h/R < 0.53, B * < B ′ 12 (case S-for Small aspect ratio).

The mixed solutions

We now consider the case where C = 0 and B = B mn . For B = B ′′ mn , we recover the eigenfunction φ ′′ mn as a limit case. Therefore, the even eigenmodes are limit points of the main branch. On the other hand, for B = B ′ mn , the solution of Eq. ( 50) is not unique. Indeed, we can always add to the solution (52) an eigenmode χ mn φ ′ mn . This leads to the mixed solution

φ M = m ′ n ′ φ ′′ m ′ n ′ r (B ′′ m ′ n ′ ) 2 -(B ′ mn ) 2 φ ′′ m ′ n ′ + χ mn φ ′ mn . (60)
The "proportion" χ mn of the eigenmode present in the mixed solution is determined by the control parameter Λ.

Taking the norm of φ M and its scalar product with r, we get

φ 2 M = χ 2 mn + m ′ n ′ φ ′′ m ′ n ′ r 2 ((B ′′ m ′ n ′ ) 2 -(B ′ mn ) 2 ) 2 , = χ 2 mn + φ 2 B , φ M r = m ′ n ′ φ ′′ m ′ n ′ r 2 (B ′′ m ′ n ′ ) 2 -(B ′ mn ) 2 = φ B r . (61)
Substituting these results in Eqs. ( 29) and ( 30), we find that χ mn is determined by Λ according to

Λ = B φ B r 2 φ B r + B 2 (χ 2 mn + φ 2 B ) , (62)
with B = B ′ mn . These mixed solutions exist in the range 0 ≤ Λ ≤ Λ(B ′ mn ) and they form a plateau at constant B = B ′ mn . For χ mn → +∞, we recover the odd eigenmode φ ′ mn at Λ = 0 and for χ mn = 0, the plateau connects the branch of continuum solutions. The mixed solutions are therefore symmetry breaking solutions. They can be seen as a mixture of a continuum solution and an eigenmode solution, like in situations with different phase coexistence.

C. The helical potential curve

In this section, we plot B as a function of Λ. For given I = 0, this curve determines the inverse helical potential 1/µ as a function of the inverse helicity 1/H (conjugate variables). It is represented in Figs. 1 and4 for the cases L and S respectively. One sees that, for a given value of the control parameter Λ, there exists multiple solutions with different values of B. We will see in Sec. III D. that, for a given value of Λ, the macroscopic energy E c.g. decreases as B increases. Therefore, low values of B correspond to high energies states and high values of B correspond to low energies states.

Case L: In this case B * > B ′ 12 and the curve B(Λ) looks typically like in Fig. 1. For a given value of Λ, we have different solutions as represented in Fig. 2. The highest energy solution is a one cell solution (continuum branch), that we choose to call "direct monopole". The second one is a two vertical cells solution (mixed branch). The cells are symmetric for Λ = 0 but one of the two cells grows for increasing Λ. The third highest energy solution is another one-cell solution (continuum branch) rotating in a direction opposite to that of the highest energy solution. We therefore call it a "reversed monopole". We call these three respective branches of solutions "branch 1" and "branch 2" for the continuum solutions, and "mixed branch" for the mixed solutions. The branches 1 and 2 connect each other at Λ(B ′′ 11 ), the location of the first even eigenmode. A typical sequence of variation of the stream function with increasing Λ on these three branches is given in Fig. 3. One sees that, as we increase Λ on the mixed branch, the two cells solution, with a mixing layer at z = h continuously transforms itself into a one cell solution, via a continuous shift of the mixing layer towards the vertical boundary.

Case S: In this case B * < B ′ 12 and the curve B(Λ) looks typically like in Fig. 4. The highest energy solution is a one cell solution (continuum branch), (direct monopole). The second solution is another one-cell solution (continuum branch ) rotating in the opposite direction (reversed monopole). The third solution is a two horizontal cells solutions (continuum branch). Some stream functions are represented in Fig. 5.

Remark: there is a maximum value of Λ max above which there is no critical point of energy at fixed helicity and angular momentum. In that case, the system is expected to cascade towards smaller and smaller scales since there is no possibility to be blocked in a "saddle point". This is a bit similar to the Antonov instability in stellar dynamics due to the absence of critical point of entropy at fixed mass and energy below a critical value of energy [2,3,27,28]. In that case, the system is expected to collapse (gravothermal catastrophe). It is not yet clear whether a similar process can be achieved in experiments of turbulent axisymmetric flows. For Λ > Λ max , the system could become non-axisymmetric ruling out the theoretical analysis.

D. The coarse-grained energy

In the previous section, we have found several solutions with different values of B for each value of the control parameter Λ < Λ max . According to the variational principle (19), we should select the solution with the minimum macroscopic energy. Combining Eqs. (3), ( 4), ( 6), ( 21) and (22), we obtain the relation For the eigenmodes (C = 0), we find that

E c.g = 1 B H - 1 2 CI . ( 63
E c.g H = 1 B mn . ( 64
)
Let us consider the odd eigenmodes φ ′ mn that exist for Λ = 0 only. They are in competition with each other. We see that there is no minimum energy state since the energy decreases when (m, n) increase, i.e. when the eigenmodes develop smaller and smaller scales. Therefore, the minimum energy state corresponds to the structure concentrated at the smallest accessible scale. For the solutions of the continuum, using Eq. (53), the macroscopic energy is

E c.g H = 1 B 1 - Λ 2B φ B r . ( 65
)
We can easily plot it as a function of B (see Figs. 6 and7). Combining Figs. 1, 4, 6 and 7, we see that, for a given value of Λ, the solution with the smallest macroscopic energy corresponds to the highest B, i.e. to small-scale structures.

This is in complete opposition to what happens in 2D turbulence. In that case, the counterpart of the macroscopic energy E c.g. is the macroscopic enstrophy Γ c.g.

2

and the minimum enstrophy state corresponds to structures spreading at the largest scale. Strikingly, the bifurcation diagram in 2D turbulence [26] is reversed with respect to the present one. In conclusion, there is no global minimum of macroscopic energy at fixed helicity and angular momentum. We can always decrease the macroscopic energy by considering structures at smaller and smaller scales. Since [START_REF] Binney | Galactic Dynamics[END_REF], (11) and ( 18) are equivalent, we also conclude that there is no global maximum of entropy at fixed microscopic energy, helicity and angular momentum. We may note a similar fact in astrophysics. It is well-known that a stellar system has no global entropy maximum at fixed mass and energy [2,3,27,28]. This is associated to gravitational collapse (called the gravothermal catastrophe in the microcanonical ensemble) leading to the formation of binary stars. However, in the astrophysical problem, there exists local entropy maxima (metastable states) at fixed mass and energy if the energy is sufficiently high (above the Antonov energy). Similarly, we could investigate the existence of metastable states in the present problem. However, we will show in Sec. III G that there is no local minimum of macroscopic energy at fixed helicity and angular momentum. All the critical points ( 21

E. Chemical potential curve

In our system, the chemical potential is α = C/B. For given H, we have to plot α as a function of I (conjugate variables). The chemical potential is zero for the eigenmodes. Using the equation of state (53), we can express α for the continuum solutions as

α |H| = |Λ| B 2 φ B r , ( 66
)
where Λ is expressed as a function of B by Eq. (54). Therefore, Eq. (66) gives α/ |H| as a function of B. Eliminating B between Eqs. (66) and (54), we obtain α/ |H| as a function of Λ for the continuum. For the mixed solutions, we have

α |H| = |Λ| (B ′ mn ) 2 φ B ′ mn r , (67)

F. Caloric curve

If we come back to the initial variational problem [START_REF] Binney | Galactic Dynamics[END_REF], the caloric curve should give β as a function of the microscopic energy E = E f.g. (conjugate variables) for fixed values of H and I. Now, the temperature is determined by the expression

E f.g. = E c.g. + 1 2β . (68)
For given H and I, we can determine the discrete values of B (n) and the corresponding discrete values of E c.g.

(n) as explained previously. Then, for each discrete value, the temperature is related to the energy by Eq. (68). Therefore, the mean flow (Beltrami state) is fully determined by H and I and, for a given mean flow, the variance of the fluctuations (temperature) is determined by the energy E = E f.g. according to

1 2 T = E -E c.g. (n) . (69)
In conclusion, the caloric curve T (E), or more properly the series of equilibria, is formed by a a discrete number of straight lines with value at the origin -E c.g.

(n) and with constant specific heats C = dE/dT = 1/2. The specific heat is positive since the microcanonical and canonical ensembles are equivalent in our problem.

G. Stability analysis

In this section, we prove that the critical points of macroscopic energy at fixed helicity and angular momentum are all saddle points. A critical point of macroscopic energy at fixed helicity and angular momentum is a minimum (resp. maximum) iff the second order variations

δ 2 F ≡ δ 2 E c.g. - δ 2 H B = 1 2 rδξδφ + (δσ) 2 2r 2 - 1 B δξδσ rdrdz, (70)
are definite positive (resp. definite negative) for all perturbations that conserve helicity and angular momentum at first order, i.e. δH = ξδσ + σδξ = 0 and δI = δσ = 0. Adapting the procedure of Chavanis & Sommeria [26] to the present context, we shall determine sufficient conditions of instability.

(i) Let us prove that there is no local maximum of macroscopic energy at fixed angular momentum and helicity. Consider first the even solutions, including the continuum solutions and the even eigenmodes. We choose a perturbation such that δσ is odd and δξ = δφ = 0. For symmetry reason, this perturbation does not change I nor H at first order. On the other hand, for this perturbation δ 2 F = (δσ) 2 2r 2 rdrdz > 0. Consider now the odd eigenmodes. We choose a perturbation of the form δξ = δφ = 0 and δσ = B * rφ B * , where φ B * is the first continuum solution such that rφ B * = 0. For this perturbation, we have δσ = B * rφ B * = 0, σδξ = 0 and ξδσ = (B ′ mn ) 2 B * φ B * φ ′ mn = 0 since φ ′ mn is orthogonal to φ B * . Therefore, this perturbation does not change the helicity and the angular momentum at first order. On the other hand, for this perturbation δ 2 F = (δσ) 2 2r 2 rdrdz > 0. As a result, the critical points of macroscopic energy at fixed helicity and angular momentum cannot be energy maxima since we can always find particular perturbations that increase the energy while conserving the constraints.

(ii) Let us prove that there is no local minimum of macroscopic energy at fixed angular momentum and helicity. To that purpose, we consider perturbations of the form δσ = rφ ′ MN and δξ = B ′ MN φ ′ MN /r. The corresponding stream function is δφ = φ ′ MN /B ′ MN . Consider first the even solutions, including the continuum solutions and the even eigenmodes. In that case, we have δσ

= rφ ′ MN = 0, ξδσ = CB 2 φ B φ ′ MN + C rφ ′ MN = 0 and σδξ = CBB ′ MN φ B φ ′ MN = 0 since φ ′ MN is orthogonal to φ B .
The preceding relations remain valid for the odd eigenmodes (m, n) provided that (M, N ) = (m, n). Therefore, these perturbations do not change the helicity and the angular momentum at first order. On the other hand, for these perturbations, we have

δ 2 F = 1 - B ′ MN B . (71)
Thus, for given B and (M, N ) sufficiently large [START_REF]) = (m, n) if the critical point is an odd eigenmode[END_REF], i.e. B ′ MN > B, we have δ 2 F < 0. As a result, the critical points of macroscopic energy at fixed helicity and angular momentum cannot be energy minima since we can always find particular perturbations that decrease the energy while conserving the constraints.

In conclusion, the critical points of macroscopic energy at fixed helicity and angular momentum are saddle points since we can find perturbations making δ 2 F positive and perturbations making δ 2 F negative. This analysis shows that all Beltrami solutions are unstable. However, saddle points may be characterized by very long lifetimes as long as the system does not explore dangerous perturbations that destabilize them. This motivates the numerical stability analysis of Sec. IV.

Remark: Let us consider the odd eigenmode [START_REF] Dauxois | Dynamics and Thermodynamics of Systems with Long-Range Interactions[END_REF]2). We have seen that it can be destabilized by a perturbation φ ′ 14 or by a perturbation φ ′ MN at smaller scale. Let us now consider the effect of a perturbation of the form δξ = 1 + B 2 * φ B * /r and δσ = B * rφ B * , where φ B * is the first continuum mode such that rφ B * = 0. The corresponding stream function is δφ = φ B * . For this perturbation, we have δσ

= B * rφ B * = 0, ξδσ = (B ′ 12 ) 2 B * φ B * φ ′ 12 = 0 and σδξ = B ′ 12 rφ ′ 12 + B ′ 12 B 2 * φ ′ 12 φ B * = 0 since φ ′
12 is orthogonal to φ B * . Therefore, this perturbation does not change the helicity and the momentum at first order. For this perturbation, we have in addition

δ 2 F = B 2 * 1 - B * B ′ 12 φ 2 B * . ( 72
)
This quantity is negative when B ′ 12 < B * corresponding to case L. This implies that the eigenmode (1, 2) is also destabilized by the perturbation δφ = φ B * which is at larger scale than the perturbations φ ′ 14 .

IV. NUMERICAL STABILITY ANALYSIS: ROBUSTNESS OF SADDLE POINTS

The stability analysis performed in Sec. III G has shown that all the critical points of entropy at fixed microscopic energy, helicity and angular momentum are saddle points. We shall now investigate their robustness by using the relaxation equations derived in Paper I (for a review of relaxation equations in the context of 2D hydrodynamics, see [29]). These relaxation equations can serve as numerical algorithms to compute maximum entropy states or minimum energy states with relevant constraints. Their study is interesting in its own right since these equations constitute non trivial dynamical systems leading to rich bifurcations. Although these relaxation equations do not provide a parametrization of turbulence (we have no rigorous argument for that), they may however give an idea of the true dynamical evolution of the flow. In that respect, it would be interesting to compare these relaxation equations with Navier-Stokes simulations. This will, however, not be attempted in the present paper.

By construction, the relaxation equations monotonically increase entropy, or decrease energy, with relevant constraints. Different generic evolutions are possible: (i) they can relax towards a fully stable state (global maximum of entropy or global minimum of energy); (ii) they can relax towards a metastable state (local maximum of entropy or local minimum of energy); (iii) they do not relax towards a steady state and develop structures at smaller and smaller scales. In the present situation, we have seen that there are no stable and metastable states. Therefore, the stability analysis of Sec. III G predicts that the system should cascade towards smaller and smaller scales without limit (except the one fixed by the finite resolution of the simulations). This is a possible regime (see top of Fig. 10) but this is not what is generically observed in the experiments where long-lived structures at large scales are found (like at the bottom of Fig. 10). Here, we explore the possibility that these long-lived structures are saddle points of entropy or energy with relevant constraints. These saddle points are steady states of the relaxation equations. Although they are unstable (strictly speaking), we argue that these saddle points can be long-lived and relatively robust (this idea was previously developed for 2D flows in [24]). Indeed, they are unstable only for certain (dangerous) perturbations, but not for all perturbations. Therefore, they can be stable as long as the system does not explore dangerous perturbations that destabilize them. Of course, the rigorous characterization of this form of stability is extremely complex. In order to test this idea in a simple manner, we shall use the relaxation equations and study the robustness of the saddle points with respect to them.

A. Numerical method

Generalities

Our stability analysis is based on the numerical integration of the relaxation equations

∂ξ ∂t = -χ(βψ + µσ), (73)
∂σ ∂t = -D β σ r 2 + µξ + α , (74)
where D and χ are given functions of r and z, and β, µ and α evolve in time (see below) so as to guarantee the conservation of the invariants.

In the canonical ensemble, the temperature β is fixed and the conserved quantities are the helicity and the angular momentum. The Lagrange multipliers µ(t) and α(t) are computed at each time so as to guarantee the conservation of H and I. One may check that they are solutions of the system of algebraic equations (see Paper I)

Dξ α(t) + χσ 2 + Dξ 2 µ(t) = -β χψσ + Dξ σ r 2 (75) D α(t) + Dξ µ(t) = -β Dσ r 2 . (76)
These relaxation equations are associated with the maximization problem (11) provided that, at any given time, the distribution of angular momentum is given by Eq. ( 12) with constant β (see Paper I for details). By properly redefining the Lagrange multipliers, they are also associated with the minimization problem (18). In the microcanonical ensemble, the conserved quantities are E, H and I. In the sequel, it will be convenient to fix the time dependence of β by imposing β -1 (t) = 2 (E -E c.g (t)) at each time. Taking into account the two other invariants, one may check that β(t), α(t) and µ(t) are solution of the system of algebraic equations

Dξ α(t) + χσ 2 + Dξ 2 µ(t) + χψσ + D ξσ r 2 β(t) = 0, (77)
D α(t) + Dξ µ(t) + Dσ r 2 β(t) = 0, (78) 2
(E -E c.g ) β(t) = 1. (79)
These relaxation equations are associated with the maximization problem [START_REF] Binney | Galactic Dynamics[END_REF] provided that, at any given time, the distribution of angular momentum is given by Eq. ( 12) with β = β(t) (see Paper I for details).

In the sequel we focus on the special case D = D * r 2 and χ = χ * r -2 , where D * and χ * are constants, that allows a simple numerical treatement of the relaxation equations by projection along the Beltrami eigenmodes:

σ = Nm 1 rs n φ n , (80) ξ
= Nm 1 r -1 x n φ n , (81)
ψ = Nm 1 rp n φ n , (82)
where n = (m, n), 0 ≤ n ≤ N , 0 ≤ m ≤ M label the modes and N m = N × M is the number of modes. In that case, Eqs. ( 73) and (74) can be transformed into a set of N m ODEs:

ṡ n = -D * [βs n + µx n + α rφ n ] , (83) ẋn
= -χ * [βp n + µs n ] , (84)
p n = B -2 n x n , (85)
where B n is such that -∆φ n +φ n /r 2 = B 2 n φ n . Note that the constraints couple Eqs. ( 83)-(85) through the parameters β, α and µ. To investigate the robustness of a given stationary solution, we first perturb it with a suitable perturbation (see below), and then follow its dynamics thanks to the relaxation equations. Two typical time evolutions are provided in Fig 10 : if the solution is fragile with respect to the perturbation, it will cascade to another solution (usually the solution of smallest scale permitted by our resolution); if the solution is robust with respect to this perturbation, it will eventually return to its initial unperturbed state. To quantify the robustness of a given solution, we define a probabilistic stability criterion by computing the probability for the solution "to escape" from its basin of attraction. To that purpose, we select a threshold δ and compute at each time the probability of escape

p esc (t) = Prob r -2 (δσ(t)) 2 r -2 (δσ(0 + )) 2 > δ , (86)
using N p realizations with perturbations drawn at random at t = 0 from a suitable ensemble (see below). This allows us to define "statistically fragile" solutions as those for which p esc (t) → 1 when t → ∞, the others being referred to as "statistically robust". In practice, the limit t → ∞ is not accessible. We thus generalize this notion to a "finite time", by considering the asymptotic value of p esc reached at the largest time of the simulation, t max . In addition, the asymptotic value of p esc (t) provides a mean to quantify the degree of robustness of a solution. Examples are given in Fig. 11, for a fragile and for a robust solution. As can be seen, the fragile solution is fragile whatever the threshold δ. However, the degree of robustness of a solution depends on the threshold δ. Quite naturally, the larger the threshold, the more robust the solution. Remark: although the variational problems [START_REF] Binney | Galactic Dynamics[END_REF], (11) and ( 18) are equivalent, and all lead to the absence of stable equilibrium state, the corresponding relaxation equations described previously are different. Therefore, the robustness of the saddle points will be different in the canonical and microcanonical settings. This can be viewed as a form of "ensembles inequivalence" for an out-of-equilibrium situation.

Perturbations

The stability must be investigated using perturbations that rigorously conserve the integral constraints. This puts some conditions regarding the shape of the possible perturbations that we can use.

In the canonical ensemble, the integral constraints are H and I. Given an initial stationary solution (σ ⋆ , ξ ⋆ ), the perturbations (δσ, δξ) must obey ξ ⋆ δσ + δξσ ⋆ + δξδσ = 0, (87) δσ = 0.

(88)

One can check that this set of constraints is satisfied by any perturbation of the form

δσ = ǫr φ i1 - rφ i1 rφ i0 φ i0 (89) δξ = ǫr -1 s ⋆ i2 -1 rφ i1 rφ i0 x ⋆ i0 -x ⋆ i1 φ i2 (90) δψ = ǫrs ⋆ i2 -1 B -2 i2 rφ i1 rφ i0 x ⋆ i0 -x ⋆ i1 φ i2 ( 91
)
where ǫ is the amplitude of the perturbation, i 0 labels an even mode while i 1 and i 2 label two different modes different from i 0 such that σ ⋆ φ i2 = 0. Following Eqs. ( 80) and (81), we have set σ ⋆ = rs ⋆ i φ i and ξ ⋆ = r -1 x ⋆ i φ i . In the sequel, we fix the amplitude of the perturbation ǫ through the norm A 2 p = r -2 (δσ) 2 by imposing

ǫ = A p 1 + rφ i1 rφ i0 2 -1 2 .
The modes i 0 , i 1 and i 2 are chosen randomly according to the following procedure: i) we draw i 0 following a uniform law among the N × M even modes; ii) we draw i 2 following a uniform law among the N × M or N × M + 1 modes of the set of allowed σ ⋆ , excluding i 0 . This mode is therefore necessarily even for solution of continuum, and often even for mixed solutions; iii) we draw i 1 following a uniform law among the 2 × N × M even and odd modes, excluding i 0 and i 2 . This choice allows the generation of N p random perturbations with the same amplitude A p .

In the microcanonical ensemble, the relaxation equations conserve in addition the energy. To satisfy this additional constraint, we choose the perturbations according to the same procedure as in the canonical case, and then determine the initial value of the temperature β 0 = β(t = 0) in order to guarantee the conservation of the energy [35]. As explained previously, this amounts to taking β -1 0 = 2 (E -E c.g (σ ⋆ + δσ, ξ ⋆ + δξ)) . In the following, we shall group the perturbations into subclasses such that perturbations of the same class have the same temperature β 0 or, equivalently, the same macroscopic energy E c.g (σ ⋆ + δσ, ξ ⋆ + δξ). Note that the initial temperature of the perturbation differs from the temperature of the equilibrium state which is given by β -1 eq = 2 (E -E c.g (σ ⋆ , ξ ⋆ )) .

Parameters

In the sequel, we focus on the stability analysis in the case L, for the first three branches of solutions, relevant for comparison with experiments, see Paper III [START_REF] Dubrulle | Paper III: Comparison with experiments[END_REF]. Our parameters are as follows:

• The number of modes is N m = 2 N M with N = 10 (radial modes) and M = 12 (vertical modes) corresponding to 120 even modes and 120 odd modes. The radial and vertical lengths are R = 1.2 and h = 1.4.

• The amplitude of the perturbations is A p = 0.05. We consider N p = 200 realizations for each given stationary solution.

• The parameters D * and χ * are both taken equal to 1. The relaxation equations are integrated using an implicit Heun scheme. The time step is empirically chosen proportional to β(t = 0) -1 . For β(t = 0) = 1, the time step is 0.02. We have checked that this time step is small enough to guarantee the numerical conservation of I and H (canonical case) or E, I and H (microcanonical case).

Numerical protocol

For any value of Λ on a given branch of solutions, we proceed as follows:

Canonical ensemble i) We fix the value of the temperature β (it remains constant during the evolution). In the sequel, we focus on five arbitrary values, β = 1000, 100, 10, 1, 0.1, chosen so as to span a wide range.

ii) We compute the Beltrami solution (σ ⋆ , ξ ⋆ ) corresponding to a prescribed value of Λ on the given branch.

iii) We generate N p perturbed initial conditions leaving unchanged the helicity and the angular momentum of (σ ⋆ , ξ ⋆ ). iv) We evolve the perturbed initial conditions through Eqs. ( 83)-(85) and Eq. (75,76) for a certain amount of time t max .

Microcanonical ensemble

i) We fix the value of the energy E (it remains constant during the evolution). In the sequel, it is fixed after arbitrary choice of five values of the temperature, β 0 = 1000, 100, 10, 1, 0.1, chosen so as to span a wide range. Once β 0 has been fixed, the total energy is then fixed. It can vary from one realization to the other, but does not vary along the evolution.

ii) We compute the stationary Beltrami solution (σ ⋆ , ξ ⋆ ) corresponding to a prescribed value of Λ on the given branch.

iii) We generate N p perturbed initial conditions leaving unchanged the helicity and the angular momentum of (σ ⋆ , ξ ⋆ ). iv) We group together the perturbations that have the same initial temperature β -1 0 = 2(E -E c.g. ) measuring the initial energy of the fluctuations (equivalently, these perturbations have the same value of macroscopic energy E c.g. ). iv) We evolve the perturbed initial condition through Eqs. ( 83)-( 85) and (77-79) for a certain amount of time t max .

B. Numerical results

Robustness in the canonical ensemble

On the three branches, we computed the value p esc (∞) (computed at t max ) as a function of Λ for different temperatures T = β -1 . The results are displayed on Fig. 12. The mixed branch (vertical dipoles) and branch 2 (reversed monopoles) are found very robust for high threshold, and still retain a certain degree of robustness for a small threshold, with a 40 per cent probability of escape. There is no clear dependence on the temperature. This is natural, since temperature can be eliminated by a suitable rescaling of time (or of coefficients D and χ) and redefinition of Lagrange parameters. The behavior on branch 1 (direct monopoles) is more contrasted and provides a very clear transition around the critical value Λ c ≈ 0.1. For Λ < Λ c , the probability to escape is close to 1, meaning large fragility of the branch. For Λ > Λ c , the branch becomes much more robust, reaching a larger degree of robustness than the two other branches for high threshold, while reaching the same robustness for small threshold. The different behaviors are summarized on Fig. 13.

Note that the value Λ c = 0.1 is somewhat arbitrary. Indeed, increasing t max further, we observed the same qualitative scenario, with an increased value of Λ c . We also observed that over sufficiently long time, the branch 2 tends to become unstable, past a value of the order Λ = 0.25.

Robustness in the microcanonical ensemble

On the three branches, we computed the value p esc (∞) (computed at t max ) as a function of Λ for classes of perturbations with different initial temperature T 0 = β -1 0 . Note that the initial temperature fixes the amplitude of the velocity fluctuations. The results are displayed on Fig. 14. For the mixed branch (vertical dipoles) and branch 2 (reversed monopoles), the microcanonical results do not noticeably differ from the canonical results: the two branches are found very robust for large threshold, and still retain a certain degree of robustness for a small threshold, with a 40 per cent probability of escape. There is no clear dependence on the initial temperature. There is therefore no ensembles inequivalence for these two branches. This is not true anymore for branch 1 (direct monopoles). Indeed, one still observes a transition from robustness to fragility around a critical value Λ c but this quantity depends on the initial temperature T 0 : it takes a value Λ c ≈ 0.1 at large initial temperatures (large velocity fluctuations) and then decreases to 0 for small initial temperatures (small velocity fluctuations). The difference of robustness observed between the two ensembles may be seen as a kind of inequivalence of ensembles at small initial temperatures. The different behaviors are summarized on Fig. 15. Like in the canonical case, we checked that an increase of t max results in a larger fragility of the branch 1 and 2 towards small Λ, at a given temperature.

Remark: note that perturbations with small initial temperature have large macroscopic energies corresponding to perturbations at large scales. According to Sec. III G such perturbations are less destabilizing than perturbations at small scales (associated with small macroscopic energies hence large temperatures). This may explain the numerical results. We have studied the thermodynamics of axisymmetric Euler-Beltrami flows and proved the coexistence of several equilibrium states for the same values of the control parameters. All these states are saddle points of entropy but they can have very long lifetime as long as the system does not spontaneously develop dangerous perturbations. We have numerically explored the robustness of some of these states by using relaxation equations in the canonical and microcanonical ensembles. The dipoles (mixed branch) and the reversed monopoles (branch 2) were found to be rather robust in both ensembles. Furthermore, in the microcanonical ensemble there is no dependence on the initial temperature on these branches. By contrast, the direct monopoles (branch 1) display a sharp transition around a critical value Λ c . The value of Λ c increases with increasing integration time. In the microcanonical ensemble, this value also decreases with decreasing initial temperature, resulting in a difference of robustness in the canonical and microcanonical ensembles. This difference may be seen as a kind of "ensembles inequivalence". This is, however, a very unconventional terminology since it concerns here the robustness of saddle points with respect to random perturbations that keep the energy or the temperature fixed, over a finite amount of time.

B. Bifurcation scenario

The simulations have shown that the dipole (two-cells solution) is relatively robust for any value of the angular momentum. On the other hand, the direct monopole (one-cell solution) is very fragile at low angular momentum but becomes robust at high angular momentum. In that case, it is even more robust than the dipole. Therefore, increasing the total angular momentum of the flow, one expects to observe a transition from the two-cells solution (antisymmetric with respect to the middle plane) to the one-cell solution (symmetric with respect to the middle plane). This bifurcation scenario is sketched in Fig. 16. It is reminiscent of the turbulent transition reported in the von Kármán flow [START_REF] Ravelet | [END_REF] in which the initial two-cells flow observed at zero global rotation suddenly bifurcates when the rotation is large enough. Once the bifurcation has taken place, the level of fluctuation is experimentally observed to decrease strongly, resulting in a decrease of the statistical temperature. In our scenario, this means that the monopole branch is suddenly stabilized with respect to redecrease of the total angular momentum of the flow, resulting in a hysteresis that has also been observed experimentally. It would therefore be interesting to investigate more closely the relevance of our scenario to the experimental system. This is done in the next paper [START_REF] Dubrulle | Paper III: Comparison with experiments[END_REF].

C. Richardson energy cascade

An interesting outcome of our study lies in the fate of the solutions when they are destabilized by a dangerous perturbation: due to the energy minimization principle, the unstable solution tends to "cascade" towards a higher wavenumber solution in a way reminiscent to the Richardson energy cascade of 3D turbulence (see Fig. 10). The cascade stops when the largest available wavenumber is reached, since dangerous perturbations are necessarily at smaller scale than the achieved state. This form of energy condensation at the smallest scale may be seen as an interesting counterpart (in the opposite sense) of the large scale energy condensation observed in 2D turbulence via the inverse energy cascade process. This is the signature of the 2D and a half nature of our system, intermediate between 2D and 3D turbulence.

VI. CONCLUSION

We have characterized the thermodynamical equilibrium states of axisymmetric Euler-Beltrami flows and proved the coexistence of several equilibrium states for a given value of the control parameter like in 2D turbulence [26]. We further showed that all states are saddle points of entropy and can, in principle, be destabilized by a perturbation with a larger wavenumber, resulting in a structure at the smallest available scale. This mechanism is therefore reminiscent of the 3D Richardson energy cascade towards smaller and smaller scales. Therefore, our system is truly intermediate

FIG. 1 : 4 FIG. 2 :

 142 FIG. 1: Top: B as a function of Λ for case L (we have taken R = 1.4 and h = 1.2). For a given value of Λ (we have taken Λ = 0.05), the solutions of the continuum are denoted by red circles and the mixed solutions by green circles. The mixed solution branches are drawn using dotted lines. One observes multiplicity of solutions: at given Λ correspond several solutions with different B.

4 FIG. 3 :

 43 FIG. 3: Stream function φ along the three branches of solution at Λ = 0.004, 0.1 and 0.3 from left to right for each branch. Top= branch 1, direct monopole; Middle: mixed branch (vertical dipole); Bottom: branch 2,reversed monopole. For simplicity, we show at each point only one solution, corresponding to a given sign of I. The solution corresponding to opposite sign of I can be found by a change φ → -φ.

12 ×FIG. 4 : 1 FIG. 5 :

 12415 FIG. 4: Top: B as a function of Λ for case S (here R = 1.4 and h = 0.5). The solutions of the continuum are denoted by red circles and the mixed solutions by green circles. The mixed solution branches are drawn using dotted lines.

FIG. 6 :

 6 FIG.6: Macroscopic energy E c.g /H as a function of B for case L. The energy of the even eigenmodes are denoted by red circles and the energy of the odd eigenmodes by green circles.

FIG. 7 :

 7 FIG. 7:Macroscopic energy E c.g /H as a function of B for case S. The energy of the even eigenmodes are denoted by red circles and the energy of the odd eigenmodes by green circles.

11 FIG. 8 : 11 FIG. 9 :

 118119 FIG.8: Chemical potential versus angular momentum in case L. The chemical potential of the even eigenmodes are denoted by red circles and the chemical potential of the odd eigenmodes by green circles. The mixed solution branches are drawn using dotted lines.

FIG. 10 :

 10 FIG. 10: Time evolution of the stream function of a solution of branch 1 of the continuum. Top: in the fragile case; Bottom: in the robust case.

FIG. 11 :

 11 FIG. 11: Time evolution of pesc made with 200 perturbations around two Beltrami states. Top: fragile solution; Bottom: robust solution. Two different thresholds are used: δ = 1.1 (continuous line) δ = 20 (dotted line).

FIG. 12 :

 12 FIG. 12: Canonical ensemble. pesc(∞) as function of Λ on the three branches (left: branch 1, middle: branch 2, right: mixed branch) at different temperatures 1/β: 0.001(•), 0.01(•), 0.1(•), 1(•), and 10(•). Two different thresholds are used: δ = 1.1 (continuous line) = 20 (dotted line).

FIG. 13 :

 13 FIG.13: Robustness of the three branches in the canonical case computed at a given tmax. The lines are increasingly fat with increasing pesc(∞), i.e. robustness. Note that the value of Λc increases with increasing tmax.

FIG. 14 :

 14 FIG. 14: Microcanonical ensemble. pesc(∞) as function of Λ on the three branches (left: branch 1; middle: branch 2; right: mixed branch) at different initial temperatures 1/β0: 0.001(•), 0.01(•), 0.1(•), 1(•), and 10(•). Two different thresholds are used: δ = 1.1 (continuous line) δ = 20 (dotted line).

FIG. 15 :

 15 FIG.15: Robustness of the three branches in the microcanonical case computed at a given tmax. Left: for a low initial temperature; right: for a high initial temperature. The lines are increasingly fat with increasing pesc(∞), i.e. robustness. Note that the value of Λc increases with increasing tmax.

FIG. 16 :

 16 FIG.16: Bifurcation scenario. Left: the system starts at Λ = 0, on the mixed branch, in a two-cells topology (vertical dipole); Middle: increasing Λ up to Λc (red arrow), the system follows the mixed branch. One cell grows at the expense of the other, resulting in a shift of the mixing layer upwards; Right: At Λc, the system bifurcates towards branch 1 (more stable), resulting in a one-cell topology (direct monopole). In this new state, the fluctuations are much milder (empirical fact from experiments), thereby allowing a stabilization of the branch towards lower Λ. Therefore, the monopole subsists beyond this point, even after a redecreasing of Λ, depicted by the black arrow.

between 2D turbulence (coherent structures) and 3D turbulence (energy cascade). Through a numerical exploration of the robustness of the equilibrium states with respect to random perturbations using a relaxation algorithm in both canonical and microcanonical ensembles, we showed however that these saddle points of entropy can be very robust and therefore play a role in the dynamics. We evidenced differences in the robustness of the solutions in the canonical and microcanonical ensembles leading to a theoretical scenario of bifurcation between two different equilibria (with one or two cells) that resembles a recent observation of a turbulent bifurcation in a von Kármán experiment [31].

This work was supported by European Contract WALLTURB.

Appendix A

We show that the first zero of F (B), denoted B * , is always between the first B ′′ 1 and the second B ′′ 2 even eigenmode. To that purpose, we note that if B < B ′′ 1 then (B ′′ mn 2 -B 2 ) > 0 for any (m, n) so that F (B) > 0. There is no discontinuity of F in the interval [0, B ′′ 1 [ so that there is no zero in that interval. Consider now the interval ]B ′′ 1 , B ′′ 2 [. In that interval, F is also continuous and increasing since

2 ) -. Therefore, there exists a unique value of B * in the range ]B ′′ 1 , B ′′ 2 [, such that F (B * ) = 0. This shows that the first zero of F lies in between the first two even eigenmodes. This property remains true for the successive values of B (n) * and the successive even eigenmodes.