Hamiltonian formulation and analysis of a collisionless fluid reconnection model

# Hamiltonian formulation and analysis of a collisionless fluid reconnection model

## Abstract

The Hamiltonian formulation of a plasma four-field fluid model that describes collisionless reconnection is presented. The formulation is noncanonical with a corresponding Lie-Poisson bracket. The bracket is used to obtain new independent families of invariants, so-called Casimir invariants, three of which are directly related to Lagrangian invariants of the system. The Casimirs are used to obtain a variational principle for equilibrium equations that generalize the Grad-Shafranov equation to include flow. Dipole and homogeneous equilibria are constructed. The linear dynamics of the latter is treated in detail in a Hamiltonian context: canonically conjugate variables are obtained; the dispersion relation is analyzed and exact thresholds for spectral stability are obtained; the canonical transformation to normal form is described; an unambiguous definition of negative energy modes is given; and thresholds sufficient for energy-Casimir stability are obtained. The Hamiltonian formulation also is used to obtain an expression for the collisionless conductivity and it is further used to describe the linear growth and nonlinear saturation of the collisionless tearing mode.

## 1 Introduction

Due to the lower dimensionality of configuration space as compared to phase space, fluid models of the plasma have an intrinsic computational advantage over kinetic models. For this reason it is of great interest to develop and apply such models even in cases where the collisionality is too small to provide a firm justification for their use. In particular, fluid models have made important contributions to the understanding of magnetic reconnection, a phenomenon that plays a key role in events such as solar flares, magnetospheric substorms, and sawtooth oscillations in tokamaks. [1, 2, 3, 4] They have also made key contributions to the understanding of plasma turbulence in the core [5, 6, 7, 8] and edge[9, 10, 11, 12] of magnetic confinement experiments in collisionless as well as collisional regimes. More generally, they offer the promise of being able to perform simulations of multiscale phenomena that are beyond the reach of kinetic models even after accounting for foreseeable advances in computation speed.[13]

An important consideration when constructing new plasma fluid models is the existence of a Hamiltonian structure (see [14, 15, 16] for reviews). Because the fundamental laws governing charged particle dynamics are Hamiltonian, dissipative terms, which ultimately arise from simplifications, must be accompanied by phenomenological constants such as resistivity and viscosity. When such phenomenological quantities are neglected, it is desirable that the resulting model be Hamiltonian, as is the case for the most important kinetic and fluid models of plasma physics. The preservation of the Hamiltonian structure provides some confidence that the truncations that are used to derive the fluid model have not introduced unphysical sources of dissipation. The presence of the Hamiltonian structure has the additional benefit of providing important tools for calculations. For example, the magnetohydrodynamic (MHD) energy principle is a consequence of the Hamiltonian nature of MHD.

Since the discovery of the noncanonical Hamiltonian structure of MHD [17], many plasma fluid models have been shown to possess a Hamiltonian description in terms of noncanonical Poisson brackets, e.g. [4, 6, 7, 8, 18, 19, 20, 21, 22, 23]. In some cases, the requirement that the dynamics has Hamiltonian form has been used to guide the construction[6, 20] and has led to the identification of new and physically important terms.[20] In another case, the absence of a Hamiltonian structure for a given model was shown to lead to the violation of the solubility conditions for the equilibrium equations.[24] The Hamiltonian structure of fluid models has also been shown to be important for the consistent calculation of zonal flow dynamics.[6, 7, 8]

Several fluid models have been proposed to study electromagnetic plasma dynamics (see [12, 25, 26] for reviews). Some of these models have been instrumental in advancing our understanding of magnetic reconnection.[20, 21, 25, 27, 4] In particular, the model of Ref.[20] led to the discovery of fast (compared to Sweet-Parker) magnetic reconnection by Aydemir.[28, 29] This model included the effects of finite ion temperature, but neglected electron inertia. An alternative model that neglected ion temperature but did include electron inertia as well as curvature effects was proposed around the same time by Hazeltine and Meiss (HM). The HM model was originally used to provide a unified description of the formation of current channels in semi-collisional and collisionless regimes. Its Hamiltonian nature, however, has not been investigated until now.

Interest in the effects of electron inertia was recently revived by a controversy over its influence on the rate of collisionless magnetic reconnection (CMR). In CMR, the “frozen-in” condition of MHD is broken by the inclusion of electron inertia instead of resistivity.[30, 31] This led Schep and collaborators to study two models that may be viewed as limiting cases of the HM model (aside from the fact that they avoid the Boussinesq approximation, retaining instead the density in the form ).[4, 21, 27] These authors constructed a Hamiltonian formulation for these models, which were subsequently used to demonstrate the role of phase mixing of the Lagrangian invariants during fast reconnection[22] and to investigate the role of instabilities of the nonlinearly developed current sheet.[32] More recently, Fitzpatrick and Porcelli (FP) have considered another limiting form of the HM model that, compared to the model derived by Schep et al., is valid for a wider range of values of , the latter indicating the ratio between the plasma pressure and the magnetic pressure based on the toroidal guide field. The FP model also extends the models of Schep et al. by including the effects of parallel ion compressibility. It has subsequently been used to study two-fluid effects on the Taylor problem[33] and on the linear growth of tearing modes.[34] Rogers et al. have recently shown that the predictions of the FP model for the linear growth rate of the tearing mode are in good agreement with those obtained with the gyrokinetic code GS2 when the ion temperature is not too large.[35]

In the present paper, we investigate the Hamiltonian structure of the FP version of the model of HM.[25, 36] The paper is organized as follows. In Sec. 3 we present the noncanonical Poisson bracket and show that this bracket produces the equations of motion with the appropriate Hamiltonian. In Sec. 4, we use the noncanonical Poisson bracket to obtain four infinite families of new Casimir invariants, three of which suggest that specific combinations of the field variables are Lagrangian invariants. In terms of these variables, the equations of motion and the Hamiltonian structure achieve a much simplified form. A preliminary version of these results was announced in [37].

The remaining sections describe a variety of applications that rely on the Hamiltonian structure. In Sec. 5 we describe a variational principle for equilibria of the system and show that they are governed by a generalized Grad-Shafranov system of a pair of coupled elliptic equations. We treat two examples: dipole equilibria with Bessel function solutions and homogeneous equilibria that support wave motion. Sec. 6 explores the latter example in further detail by showing how to construct conventional canonical variables for the linear dynamics. We show that the system possesses Alfvén-like and drift-shear modes, obtain exact stability thresholds, and give a definition of negative energy modes. We note that the Hamiltonian form is indispensable for an unambiguous definition of negative energy modes. We also obtain energy-stability conditions, sufficient conditions for stability akin to the criterion of MHD. Section 7 contains a derivation of the collisionless conductivity that relies on the Jacobi identity of the Hamiltonian formulation, which we then use to obtain the tearing-layer parameter and the growth rate for the collisionless tearing mode. In Sec. 8 we use the conservation of a Casimir invariant to obtain the nonlinearly saturated current profile and compare it to that obtained by Rutherford.[38] In Sec. 9 we summarize and conclude.

## 2 Model equations

The model of [36] is given by the following equations:

 ∂(ψ−d2e∇2ψ)∂t+[φ,ψ−d2e∇2ψ]−dβ[ψ,Z]=0, (1) ∂Z∂t+[φ,Z]−cβ[v,ψ]−dβ[∇2ψ,ψ]=0, (2) ∂∇2φ∂t+[φ,∇2φ]+[∇2ψ,ψ]=0, (3) ∂v∂t+[φ,v]−cβ[Z,ψ]=0. (4)

Equation (1) is a reduced Ohm’s law where the presence of finite electron inertia, which makes it possible for MR to take place, is indicated by the terms proportional to the electron skin depth . Equations (2), (3) and (4) are obtained from the electron vorticity equation, the vorticity equation, and the parallel momentum equation, respectively.

Considering a Cartesian coordinate system and taking as an ignorable coordinate, the fields , , and are related to the magnetic field and to the velocity field by the relations and , respectively. Here is a constant guide field, whereas and with indicating the ion skin depth. For small , , the sonic Larmor radius. The ions are assumed to be cold, but electron pressure perturbations are taken into account and are given by , with a constant background pressure and coupled to the magnetic field via the relation . Notice, here the parameter is defined as , and above all the quantities are expressed in a dimensionless form according to the following normalization: , , , where is a typical scale length of the problem, is a reference value for the poloidal magnetic field, and is the Alfvén speed based on and on the constant density. Finally, , for generic fields and .

## 3 Hamiltonian formulation

A desirable property for fluid models of the plasma is that the non-dissipative part of their equation of motion should admit a noncanonical Hamiltonian formulation [14, 15, 16]. In short this means that it is possible to reformulate the ideal part of an -field model as

 ∂ξi∂t={ξi,H},i=1,⋯,n, (5)

where are suitable field variables, is the Hamiltonian functional and is the Poisson bracket consisting of an antisymmetric bilinear form satisfying the Jacobi identity.

The first task in the derivation of a noncanonical Hamiltonian formulation is to identify a conserved functional, usually the energy, that can serve as the Hamiltonian of the model. If one considers, for instance, a square domain in the plane with doubly periodic boundary conditions, the four-field model (1)–(4) admits the following constant of the motion:

 H=12∫Dd2x(d2eJ2+|∇ψ|2+|∇φ|2+v2+Z2) (6)

with indicating the parallel current density. The quantity represents the total energy of the system. The first term refers to the kinetic energy due to the relative motion of the electrons with respect to ions along the direction. The third and fourth terms account for the kinetic energy, whereas the second and last terms account for the magnetic energy.

Adopting , , , and as field variables, i.e. , and (6) as Hamiltonian, it is possible to show that the model can indeed be cast in a noncanonical Hamiltonian form with the following Lie-Poisson bracket:

 {F,G}=∫d2x(U[Fξ,Gξ]U+ψe[Fξ,Gξ]ψe+Z[Fξ,Gξ]Z+v[Fξ,Gξ]v), (7)

where

 [Fξ,Gξ]U = [FU,GU] [Fξ,Gξ]Z = [FZ,GU]+[FU,GZ]−dβde2[Fψe,Gψe] +cβde2([Fv,Gψe]+[Fψe,Gv])−α[FZ,GZ]−cβγ[Fv,Gv] [Fξ,Gξ]ψe = [Fψe,GU]+[FU,Gψe]−dβ([FZ,Gψe]+[Fψe,GZ])+cβ([Fv,GZ]+[FZ,Gv]) [Fξ,Gξ]v = [Fv,GU]+[FU,Gv]+cβde2([FZ,Gψe]+[Fψe,GZ]−cβγ([Fv,GZ]+[FZ,Gv],

with , , and subscripts indicate functional differentiation. It is straightforward to show that Eq. (5) with the Hamiltonian given in Eq. (6) and the above bracket reproduces equations (1)-(4) of [36]. The Jacobi identity for the above Poisson bracket is readily established by applying the method described in Refs. [14, 39].

Because the number of parameters is escalating, we record their definitions here for later referral:

 dβ=cβdi,d=√d2e+d2i,γ=de2/di,α=dβ+cβde2/di=cβd2/di. (9)

The basic parameters of the model are , and , while , and are useful shorthands.

## 4 Casimir invariants and bracket normal form

Lie-Poisson brackets for noncanonical Hamiltonian systems are accompanied by the presence of Casimir invariants. A Casimir invariant is a functional that annihilates the Lie-Poisson bracket when paired with any other functional. That is, a Casimir must satisfy

 {F,C}=0, (10)

for every functional . Thus, Casimir invariants constrain the nonlinear dynamics generated by the Poisson bracket for any choice of Hamiltonian.

In order to identify the Casimirs of the four-field model we proceed in the following way. First, multiplying Eq. (4) times and adding it to Eq. (1), we find

 ∂D∂t+[φ,D]=0, (11)

where is the ion canonical momentum. Equation (11) indicates that the field is a Lagrangian invariant that is advected by the flow generated by the stream-function . The presence of this Lagrangian invariant also suggests that using as one of the variables will simplify the Lie-Poisson bracket. Indeed, upon replacing with as field variable, Eq. (10) for the four-field model becomes

 {F,C}=∫d2x(FU[CU,U]+FD[CU,D]+FU[CD,D] (12) +cβFv[CZ,D]+cβFZ[Cv,D]+FZ[CU,Z]+FU[CZ,Z] −αFZ[CZ,Z]−cβγFv[Cv,Z]+Fv[CU,v]+FU[Cv,v] −αFv[CZ,v]−αFZ[Cv,v])=0,

a simpler bracket.

For Eq. (10) to be satisfied for any , it is necessary for the coefficients of each of the functional derivatives of in (12) to vanish separately. This leads to the following system of equations for :

 [CU,D] = 0, (13) [CU,U]+[CD,D]+[CZ,Z]+[Cv,v] = 0, (14) −cβ[Cv,D]−[CU,Z]+α([CZ,Z]+[Cv,v]) = 0, (15) cβ[CZ,D]−cβγ[Cv,Z]+[CU,v]−α[CZ,v] = 0. (16)

The problem of finding the complete set of Casimir invariants is thus equivalent to solving the set of four coupled equations (13)–(16).

Beginning with (13), a functional integration shows that is of the form

 C(U,D,Z,v)=∫d2x(UF(D)+g(D,Z,v)), (17)

where and represent arbitrary functions of their arguments. The problem is now reduced to finding the functions and . Considering next Eq. (14), we see that this equation is automatically satisfied for any choice of with an integrand that depends only upon the field variables and not their spatial derivatives, and therefore imposes no constraints. Substitution of (17) into (15) yields

 (18)

where indicates derivative with respect to argument and subscripts on denote partial derivatives. In (18) the coefficients multiplying the brackets ’’ must vanish independently. These two relations lead to

 cβgv+αgD=ZF′(D)+K(D), (19)

where an arbitrary function of , and integration of (19) by the method of characteristics gives

 g(Z,v,D)=ZαF(D)+K(D)+G(Z,v−cβD/α), (20)

where . Similarly, insertion of (17) into (14) yields

 (cβgZv−F′(D)+αgZD)[D,v]+cβ(gZZ+γgvD)[D,Z]+(cβγgvv−αgZZ)[v,Z]=0, (21)

which gives three equations, but only provides new information. Defining , we see that satisfies the wave equation

 GZZ−cβγαGVV=0, (22)

and thus has the general solution

 G=∑±G±(±12α√cβγα(D−αcβv)−Z2α). (23)

Therefore, in light of (17), (20), and (23) we have the following four independent families of Casimir invariants:

 C1 = ∫d2x(U+Zα)F(D), (24) C2 = ∫d2xK(D), (25) C± = ∫d2xG±(±d2i2cβd3deD∓di2cβddev−di2cβd2Z), (26)

where we have scaled the arguments of for a reason that will become apparent soon.

Knowledge of the functional dependence of the Casimirs suggests a simplification of the Lie-Poisson bracket will occur if it is written in terms of the new variables

 D = D, (27) ζ = U+Zα, (28) T± = ±di2cβd3de(diD−d2v∓ddeZ), (29)

which possess the inverse relations

 D = D, (30) U = ζ+T++T−, (31) Z = −α(T++T−), (32) v = did2D−cβdeddi(T+−T−). (33)

Indeed, in the new variables, the Lie-Poisson bracket reads

 {F,G}=∫d2x(ζ[Fζ,Gζ]+D([FD,Gζ]+[Fζ,GD]) (34) +T−[FT−,GT−]+T+[FT+,GT+]),

which is a bracket normal form that relies on the scaling used above. The bracket in terms of the new variables reveals its algebraic structure: it is identified as a sum of direct product and semi-direct product parts [19, 15, 39] and, consequently, the Jacobi identity follows from general theory. Making use of the variables suggested by the form of the Casimirs, the model equations can be rewritten in the compact form

 ∂D∂t = −[φ,D], (35) ∂ζ∂t = −[φ,ζ]+d−2[D,ψ], (36) ∂T±∂t = −[φ±,T±], (37)

where for convenience we have defined

 φ±:=φ±cβddeψ. (38)

Note, the variable plays the role of a “generalized” vorticity, and our development reveals the existence of the three Lagrangian invariants , and associated with the families of Casimirs , and , respectively. The existence of such invariants implies that the values of , and are constant on the contour lines of , , and , respectively. By choosing “top-hat” functions for the free functions in the Casimirs , and , it follows that the area enclosed by the contour lines of the Lagrangian invariants remains constant. Notice also that and in the limit and become proportional to the Lagrangian invariants of the two-field model derived in [4]. The family is of a different nature and one of the constraints imposed by it is that the total value of within an area enclosed by a contour line of remains constant.

## 5 Equilibria

Equations governing the equilibrium of the FP model are most easily obtained by setting the time derivatives equal to zero in Eqs. (35)-(37) and solving for the the fields , , and . Alternatively, it is possible to derive equilibrium equations from a variational principle, the existence of which is assured by the Hamiltonian nature of the equations.[15] The construction of the variational principle is more laborious than the direct approach but the extra work is richly rewarded by the well-known benefits of variational principles. In particular, the variational principle provides a basis for studying the stability as well as the equilibria of the system.

For a system with a collection of Casimirs, extrema of the free energy functional are equilibria of the system. Such extrema can be derived by setting the first variation equal to zero. If , denotes the field variables of the system, this amounts to solving the equations , where the subscript indicates functional derivative with respect to . One advantage of this variational approach is that for equilibria obtained as extrema of , the second variation of provides criteria sufficient for stability.

### 5.1 General equilibria

For the model here, the free energy functional reads

 F[ξ] = ∫Dd2x⎡⎣d2e(∇2Lψe)22+|∇Lψe|22+|∇∇−2U|22+v22+Z22+K(ψe+div) (39) +

where the linear operator is defined by . Notice that the arguments of the functions present in the Casimirs are of course much less compact when written in terms of the variables . On the other hand, in terms of the variables the free energy functional reads

 F[Υ] = ∫Dd2x[c2βd4d2i(T2++T2−)+D22d2−12(ζ+T++T−)∇−2(ζ+T++T−) (40) − 12(ded2D+cβd(T+−T−))L(ded2D+cβd(T+−T−)) + K(D)+ζF(D)+G+(T+)+G−(T−)].

Equation (40) shows that in terms of the variables that are “natural” for the Casimirs, the expression for the Hamiltonian becomes complicated. Unfortunately, there exists no preferred set of variables in terms of which both the Hamiltonian and the Casimirs take a simple form. In order to obtain equilibrium solutions by means of the variational principle, it is convenient to choose the variables that are natural for the Casimirs, calculate the required functional derivatives of the Casimirs with respect to these variables, and then use the functional chain rule to obtain the functional derivatives of in terms of the variables . More specifically, by setting , the resulting equilibrium equations are given by

 Fζ = Hζ+∑jCjζ=Hζ+F(D)=0, (41) FD = HD+∑jCjD=HD+ζF′(D)+K′(D)=0, (42) FT± = HT±+∑jCjT±=HT±+G′±(T±)=0. (43)

where the index ranges over the set . The functional chain rule then can be used to evaluate the functional derivatives of the Hamiltonian,

 HD = d2ed2Hψe+did2Hv, Hζ = HU, HT± = ±cβdedHψe+HU−αHZ∓αdedHv.

The functional derivatives , , and can themselves easily be obtained from the Hamiltonian written in the form (6). The equilibrium equations (41)–(43) are then given by

 −φ+F(D) = 0, (44) −d2ed2∇2ψ+divd2+ζF′(D)+K′(D) = 0, (45) ∓cβded∇2ψ−φ−αZ∓αdevd+G′±(T±) = 0. (46)

These equations are expressed in a mixture of the and variables. In order to simplify them we eliminate from the last two equations by using

 d2e∇2ψ=ψ−D+div.

Using this in Eq. (45)-(46), we find

 −d−2ψ+ζF′(D)+^K′(D) = 0, (47) −φ±±cβddedi(diD−d2v∓ddeZ)+G′±(T±) = 0. (48)

where . Equation (48) may be simplified further by expressing and in terms of the variables using (32)-(33). This leads to the following complete system of equilibrium equations:

 −φ+F(D) = 0, (49) −φ±+^G′±(T±) = 0, (50) −d−2ψ+ζF′(D)+^K′(D) = 0, (51)

where . One easily verifies that Eqs. (49)-(51) describe equilibrium states by substituting them into Eqs. (35)-(37).

Continuing our present approach of eliminating the variables would now lead us to express the , , and in terms of integral operators acting on the variables. Clearly this is undesirable. Instead, we note that the above four equations express a dependency between the six quantities , , , , and . It is thus possible in principle to use these equations to express four these quantities in terms of the remaining two. If we choose and as the independent fields we will obtain a closed system of equilibrium equations of the form

 ∇2ψ = S(ψ,φ), (52) ∇2φ = P(ψ,φ), (53)

Equation (52) is a generalized version of the Grad-Shafranov equation, whereas (53) is an analogous equation that determines the equilibrium polarization.

In order to calculate the form of the functions and we invert Eq. (49)-(50):

 D = a(φ), (54) T± = t±(φ±), (55) ζ = d−2ψa′(φ)+b(φ), (56)

where and are the inverses of and respectively and . Solving these equations for and then yields (52)-(53) with

 S(ψ,φ) := ψd2e−a(φ)d2−cβdde[t+(φ+)−t−(φ−)], (57) P(ψ,φ) := b(φ)+t+(φ+)+t−(φ−)+a′(φ)ψd2. (58)

We complete the system by expressing the two remaining unknowns and in terms of and . From Eqs. (32) and (33), using (54), (55) there follows

 v(ψ,φ) = dia(φ)/d2−cβded(t+(φ+)−t−(φ−))/di; (59) Z(ψ,φ) = −α(t+(φ+)+t−(φ−)). (60)

We have thus shown that solving the equilibrium equations amounts to solving the coupled system of (52) and (53) for the unknowns and . This requires making choices for the free functions , , and . If one is only interested in solving the equilibrium problem, one may determine these free functions directly from physical considerations and the relationship between these functions and those appearing in the variational principle may be ignored. These relationships become important, however, if one wishes to use the variational principle either to solve the equilibrium or the stability problem (the variational principle has well-known advantages both for numerical and analytic applications). In this case the functions , and , appearing in the variational form may be determined in terms of , , and as described above. Variational treatments of two-fluid equilibria have been given by [40, 41, 42], and applications of variational principles to stability are discussed by [40, 41, 43, 44].

We conclude this section by noting two difficulties with the equilibrium equations (52)-(53). The first is that these equations may become hyperbolic in the presence of strong flows. Recent analyses of this problem can be found in [42, 45]. The second difficulty is that the right-hand sides are singular in the limits and , where is a macroscopic scale length.[40, 41, 46] That is, for macroscopic equilibria the derivatives in Eqs. (52)-(53) are multiplied by a small parameter, so that these equations form a stiff system. This has led to considerable grief, in particular for Field-Reversed Configuration (FRC) devices where a similar set of equations is encountered. Steinhauer has proposed a method for dealing with this problem that he named the “nearby fluid” approximation.[46, 47] In the following section we outline a similar approach to solving Eqs. (52)-(53).

### 5.2 Perturbative solution for macroscopic equilibria

In order to deal with the singular nature of the equilibrium Eqs. (52)-(53), we expand the fields in powers of the small parameters and and solve term by term. A byproduct of this procedure is the clarification of the physical meaning of the profile functions.

We begin by considering the limit . For convenience we define

 ^t±(±ψ+deφcβd):=decβdt±(φ±)

and expand this and other profile functions in powers of according to

 ^t±(ϕ)=^t(0)±(ϕ)+de^t(1)±(ϕ)+d2e^t(2)±(ϕ)+….

We then consider the equilibrium equations order by order in . From (53) and (58) we obtain at lowest order ; using this result in (52) and (57) we obtain to lowest order . To next order, (52) and (57) yield ; substituting this result in (53) and (58), we find the following equation for the vorticity

 ∇2φ=^b(φ)+a′(φ)ψ/di2−h(ψ)/dβ,

where is a profile function for the vorticity and .

From the terms of order unity in (52) and (57) we next find

 ∇2ψ=−a(φ)/d2i+h′(ψ)φ/dβ+I(ψ), (61)

where . Equation (61) is the Grad-Shafranov equation, where the term proportional to is the polarization current and describes the inductive current.

The parallel velocity may be obtained from (54): to order , , yields

 div+ψ=a(φ). (62)

Expanding (60) to order , gives

 Z+φ/dβ=h(ψ). (63)

The sum represents the electron stream-function. The fact that the electron stream-function is constant on surfaces of constant flux, as expressed by Eq. (63), is a statement of the frozen-in property.

We may carry out the limit in a similar way. From the ion momentum conservation Eq. (62), we obtain showing that the electrostatic potential must be a flux function to lowest order. It is convenient to introduce to denote the inverse of . We also define the Alfvénic Mach number . Note that Eq. (63) specifies that to lowest order, . In terms of these quantities, the vorticity equation shows that to lowest order,

 ^b(φ)=a(φ)a′(φ)/d2i−φ/c2β.

In order to eliminate from the Grad-Shafranov equation it is necessary to calculate the correction to the electrostatic potential. This is given by the vorticity equation,

 M′(∇ψ)2+M∇2ψ=(M−2−c−2β)φ(2)−h(1)(ψ).

Note that the potential exhibits a singularity for corresponding to the sound-wave resonance. Lastly, after eliminating the electrostatic potential from Eq. (61) we recover the MHD version of the Grad Shafranov equation,

 (1−M2)∇2ψ−MM′(∇ψ)2=^I(ψ).

In the following sections we present some explicit solutions of the equilibrium equations for simple profile functions.

The case of quadratic Casimir invariants is easily tractable. Choosing

and following the steps of Sec. 5.1 leads to

Upon inserting (65) into (52) and (53), we obtain

 ∇2ψ=S1ψ+S2φand∇2φ=P1ψ+P2φ, (66)

where and are constants that depend on , and the parameters of the system. Note, and are arbitrary except that . Consequently, Eqs. (66) have a variety of solutions that are closely related to the double-Beltrami flows investigated by Yoshida et al.[40, 41] Specifically, Refs. [40, 41] neglect electron inertial effects that are kept here, but they allow for more general geometry. In general, (66) can be diagonalized resulting in two decoupled equations of the form

 ∇2χi=−λiχi,i=1,2, (67)

where and the ’s are linear combinations of and . If a solution of this system is found, then one obtains and as particular linear combinations of and as described in Sec. 5.1.

Rather than describe the general solution, we give an example representative of the kinds of solutions that are possible. We assume a circular domain of unit radius, adopt polar coordinates , and adjusting the parameters , and so that are zeros (possibly distinct) of the first order Bessel function, i.e. . With these assumptions, we obtain the solution

 χi(r,θ)=AiJ1(√λir)cosθ, (68)

where the ’s are constants and each of the ’s has a dipolar structure like that depicted in Fig. 1.

### 5.4 Homogeneous equilibria

The quadratic Casimirs of (64) also yield homogeneous equilibria, i.e. equilibria for which the linear dynamics has constant coefficients. For this choice, the free energy functional of (40) can be written as follows:

 F=12∫Dd2x(ξT^Hξ+ΥT^AΥ), (69)

where

 ^H=⎛⎜ ⎜ ⎜⎝−L∇20000−∇−20000100001⎞⎟ ⎟ ⎟⎠and^A=⎛⎜ ⎜ ⎜ ⎜⎝ADAζ00Aζ00000A+0000A−⎞⎟ ⎟ ⎟ ⎟⎠. (70)

Recall and . Equations (27)-(29) amount to , where the matrix is given by

 T=⎛⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜⎝100di01dicβd20d2i2cβded30−di2cβd2−dedi2cβd3−d2i2cβded30−di2cβd2dedi2cβd3⎞⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟⎠. (71)

The free energy functional can then be written as a quadratic form:

 F=12∫Dd2x[ξT(^H+TT^AT)ξ], (72)

whence the equilibrium equations are obtained upon setting the functional derivatives to zero,

 (^H+TT^AT)ξ=0. (73)

Here we have assumed the formal self-adjointness of the operators and . Equation (73) is a linear homogeneous system of four equations with the four unknowns , , , and . The equilibria treated in Sec. 5.3 are solutions of this system, but of interest here are the homogeneous equilibria

 ψ0