This paper studies the nonlinear evolution of magnetic field turbulence in proximity of steady ideal Magnetohydrodynamics (MHD) configurations characterized by a small electric current, a small plasma flow, and approximate flux surfaces, a physical setting that is relevant for plasma confinement in stellarators. The aim is to gather insight on magnetic field dynamics, to elucidate accessibility and stability of three-dimensional MHD equilibria, as well as to formulate practical methods to compute them. Starting from the ideal MHD equations, a reduced dynamical system of two coupled nonlinear partial differential equations for the flux function and the angle variable associated with the Clebsch representation of the magnetic field is obtained. It is shown that under suitable boundary and gauge conditions such reduced system preserves magnetic energy, magnetic helicity, and total magnetic flux. The noncanonical Hamiltonian structure of the reduced system is identified, and used to show the nonlinear stability of steady solutions against perturbations involving only one Clebsch potential. The Hamiltonian structure is also applied to construct a dissipative dynamical system through the method of double brackets. This dissipative system enables the computation of MHD equilibria by minimizing energy until a critical point of the Hamiltonian is reached. Finally, an iterative scheme based on the alternate solution of the two steady equations in the reduced system is proposed as a further method to compute MHD equilibria. A theorem is proven which states that the iterative scheme converges to a nontrivial MHD equilbrium as long as solutions exist at each step of the iteration.

This study is concerned with the dynamics of the magnetic field around a magnetohydrodynamics (MHD) equilibrium1 
(1)
In this equation Bx denotes the equilibrium magnetic field in a smooth bounded domain ΩR3, x=x,y,z Cartesian coordinates, μ0 the vacuum permeability, and Px the equilibrium pressure field. For the purpose of the present paper, the dynamics around (1) is governed by the ideal MHD equations in Ω,
(2a)
(2b)
(2c)
(2d)
Here, ux,t, Bx,t, ρx,t, Px,t, and Ex,t denote the time-dependent velocity field, magnetic field, plasma mass density, pressure field, and electric field respectively. The closure of system (2) can be obtained by choosing an equation of state relating P to the other fields, and by determining E from the electron fluid momentum equation. These aspects will be discussed in detail later. The MHD equilibrium (1) can be obtained from system (2) by setting time derivatives to zero and taking a plasma flow u = 0. If the electric field is determined through the ideal Ohm’s law E + u × B = 0, the corresponding equilibrium electric field is simply E = 0.

Despite its relevance for the confinement of magnetized plasmas and the development of nuclear fusion reactors known as stellarators, a general theory concerning the existence of regular solutions B,P of the MHD equilibrium equations (1) is not available at present.2 This is because the characteristic surfaces associated with the first-order system of partial differential equations (PDEs) (1) depend on the unknown B.3,4 The existence of regular solutions is known for the special cases in which the pressure is constant or the fields B,P are invariant under some combination of continuous Euclidean isometries of R3 (rotations and translations). In the first case the magnetic field is a Beltrami field.5 In the second case, Eq. (1) reduces to the Grad–Shafranov equation, a nonlinear second order elliptic PDE for the flux function.6–8 Both cases are however not relevant for stellarators,9 which consist of toroidal vessels without trivial symmetries surrounded by coils with complex shapes whose purpose is to generate the field line twist required to minimize particle losses caused by cross-field drifts. In principle, stellarators achieve steady plasma confinement mostly through an externally produced vacuum magnetic field, and are therefore more suitable for continued operation compared to an axially symmetric tokamak where the field line twist is obtained by driving an electric current within the plasma. However, the lack of axial symmetry results in the breaking of conservation of vertical angular momentum, a fact that deteriorates the quality of plasma confinement. For this reason, in addition to (1) the equilibrium magnetic field within a stellarator must satisfy additional conditions, such as quasisymmetry, a property that ensures particle confinement by constraining particle orbits close to a given flux surface.10–12 

In this context, the aim of the present paper is (i) to obtain a closed set of reduced equations preserving the Hamiltonian structure13,14 of ideal MHD and describing the nonlinear evolution of the magnetic field in proximity of MHD equilibria (1) and in a physical regime relevant for stellarator plasmas, (ii) to use the derived equations to elucidate the stability properties of such equilibria, and (iii) to apply the derived equations to formulate dissipative and iterative schemes to construct nontrivial MHD equilibria (1) with nested flux surfaces and a non-vanishing pressure gradient in toroidal domains of arbitrary shape. In addition to providing a toy model of magnetic field turbulence in a physically relevant setting, we conjecture that the derived results may serve as a starting point for a mathematical proof of existence of nontrivial MHD equilibria (1) in toroidal domains without Euclidean symmetries.

The present paper is organized as follows. In Sec. II the ideal MHD equations (2) are reduced according to an ordering in which the plasma flow u, the electric current ∇ × B, and time derivatives /∂t are small, and the pressure field P is related to the mass density ρ and the velocity field u by a generalized Bernoulli principle. This ordering also implies the existence of approximate flux surfaces Ψ for the magnetic field. This fact is used to enforce at leading order a Clebsch representation15,16 of the magnetic field B = ∇Ψ × ∇Θ with Ψ the flux function and Θ a multi-valued (angle) variable. In the reduced system, the dynamics of the magnetic field is thus described by a pair of coupled equations for the two Clebsch potentials Ψ and Θ.

In Sec. III it is shown that under suitable boundary conditions and gauge conditions for the magnetic vector potential the reduced dynamics preserves magnetic energy, magnetic helicity, and total magnetic flux. These conservation laws are then applied to describe steady states in terms of critical points of a functional of Ψ and Θ given by a linear combination of magnetic energy and total magnetic flux. This variational formulation can be physically interpreted in analogy with Taylor relaxation18,19 in which magnetic energy is minimized under the constraint of magnetic helicity.17 Here, functionals involving higher order derivatives of the dynamical variables are dissipated at a faster rate by non-ideal (dissipative) mechanisms.20 

In Sec. IV the noncanonical Hamiltonian structure of the reduced equations is identified in terms of a Poisson bracket and a Hamiltonian functional.21,22 In particular, it is shown that the Poisson bracket satisfies all the Poisson bracket axioms, including the Jacobi identity.23 

In Sec. V the Hamiltonian structure obtained in Sec. IV is used to prove that steady solutions of the reduced dynamics are nonlinearly stable24–27 against turbulent fluctuations involving only one of the two Clebsch potentials. Here, we recall that a positive second variation of the Hamiltonian is not sufficient to guarantee nonlinear stability. In general, a norm on the space of solutions must be found such that the deviation of the perturbed solution at a given time is bound by the discrepancy of initial conditions in the prescribed topology. The type of nonlinear stability shown here effectively constrains the deviation of the perturbed Clebsch potential from initial conditions on the level sets of the other unperturbed Clebsch potential.

In Sec. VI the method of double brackets28,29 is used to formulate a dissipative dynamical system for the Clebsch potentials Ψ and Θ with the property that, instead of being constant, the Hamiltonian is progressively dissipated. This pair of equations has the structure of coupled diffusion equations. Double bracket dynamics is obtained by applying twice the Poisson bracket and represents an effective tool to compute energy minima while preserving the Casimir invariants that span the kernel of the Poisson bracket.30,31 The derived dissipative dynamical system may therefore be applied to compute MHD equilibria (1) corresponding to critical points of the Hamiltonian in a dynamical fashion.

In Sec. VII a second iterative scheme based on the alternate solution of the two steady equations of the reduced dynamical system obtained in Sec. II is discussed. Here, a theorem is proven which states that the iterative scheme converges to a nontrivial MHD equilbrium (with a non-vanishing pressure gradient) as long as solutions exist at each step of the iteration. The dissipative and iterative schemes obtained in Secs. VI and VII may be regarded as alternative methods to the constrained minimization of the plasma energy ΩB22μ0+Pγ1dV, where γ ≥ 0 is the adiabatic index, often used to numerically compute MHD equilibria.32,33

Concluding remarks are given in Sec. VIII.

The aim of this section is to derive a reduction of the ideal MHD equations that determines the nonlinear evolution of magnetic field turbulence within a plasma characterized by approximate flux surfaces. This system should be appropriate to describe the confinement regime of a tokamak or a stellarator provided that flux surfaces exist to some degree throughout time evolution. This statement will be made quantitatively precise later.

We start by considering the ideal MHD equations (2) in a smooth bounded domain ΩR3. The relationship between the pressure field P and the other field variables will be given later in the form of a generalized Bernoulli principle. The relation between the electric field E and the other fields is given by a generalized Ohm’s law following from the electron fluid momentum equation
(3)
Here, me, −e, nex,t, uex,t, and Pex,t denote the electron mass, charge, number density, fluid velocity, and pressure field, and we defined d/dt = /∂t + ue · ∇. Now recall that
(4)
where mi, δ = me/mi, ni, uix,t are the ion mass, the electron to ion mass ratio, the ion number density, and the ion fluid velocity, and we have assumed quasineutrality ni = ne = n. The electron momentum equation (3) therefore gives the generalized Ohm’s law
(5)
The closure of the MHD system (2) can thus be obtained by neglecting the last term involving the electron inertia, and by a proper ansatz on the electron pressure Pe. Indeed, assuming a barotropic equation of state for the electron pressure, Pe=Peρ, we obtain
(6)
with κ = mi/0 a physical constant associated with the Hall effect. We emphasize that the Hall effect, which scales with the ion Larmor radius, is absent in ideal MHD and represents a higher-order correction in Hall MHD.34 Therefore, it can be neglected, resulting in the ideal Ohm’s law ∇ × E = −∇ × (u × B). Remarkably, retaining the Hall term does not affect the analysis of the reduced equations derived below, allowing for a generalization within the Hall MHD framework. In the following, we shall also demand that the vector fields u, B, and ∇ × B, and the pressure P satisfy the boundary conditions
(7)
where Ω denotes the boundary of Ω and n the unit outward normal to Ω. Notice in particular that the boundary condition for ∇ × B implies that there is no net current flow across the bounding surface Ω, and it is expected to hold true as long as both ue and ui are tangent to Ω [recall Eq. (4)]. We also emphasize that the boundary conditions (7), which represent the physical constraints expected at the interface with a perfect conductor, are not applied to the ideal MHD equations (2), but to the reduced equations derived below. More details on how boundary conditions are applied to each set of reduced equations will be given later.
Let ϵ > 0 denote a small ordering parameter, L ∼ Ω1/3 the characteristic size of the system (e.g. the linear size of a stellarator), and T a reference time scale (for example, a small fraction of the confinement time scale). Assuming ρ > 0, we start by considering the following ordering conditions
(8a)
(8b)
(8c)
Physically, these requirements describe a magnetized plasma regime with small flow and small electric current and where the velocity field and the plasma density evolve slowly in response to the turbulent evolution of the magnetic field, which is driven by the non-vanishing of ∇ × E in Eq. (2b). This regime is relevant for example within stellarators, which are designed to minimize internal flows and currents. A few words are in order regarding the consistency of the ordering (8) with the plasma regimes commonly encountered in ideal MHD (see, e.g., Refs. 35 and 36). Although the pressure gradient ∇P is a second-order term in (8), this does not imply that the plasma beta (the ratio between mechanical pressure and magnetic pressure) is small. Indeed,
(9)
as long as P=P0+px,t, where the spatial constant P0 satisfies T2P0/ρL21 while the perturbation p obeys T2p/ρL2ϵ2.
Alternative orderings leading to the same governing equations will be described at the end of this section. Later we will also see that the smallness of ∇ × E implies that the magnetic field is endowed with approximate flux surfaces. Note that B/μ0ρ is the Alfvén speed, and that the ordering (8) does not involve conditions on the size of the typical particle number ρL3/mi. We now make the following ordering prescription,
(10)
When the density ρ=ρcR is constant this relationship can be satisfied through the Bernoulli principle P+ρcu2/2+h=0, where h is any function such that u · ∇h = 0. Hence, Eq. (10), which plays a role analogous to an equation of state, can be interpreted as a generalized Bernoulli principle. Notice also that Eq. (10) arises naturally from the plasma momentum equation (2a) when the system is steady and the vorticity ∇ × u and the electric current ∇ × B are small (a regime that is particularly relevant for stellarators). Equations (6) and (10) and the ordering conditions (8) can then be used to reduce system (2) to the leading order system
(11a)
(11b)
(11c)
(11d)
(11e)
Observe that in the limit u0 steady solutions of system (11) are described by the MHD equilibrium equations (1) and a barotropic equation of state P=Pρ. We also emphasize that, in addition to the induction Eq. (11c) and the divergence-free condition (11e) for the magnetic field, system (11) comprises five other equations [momentum equation (11a), equation of state (11b), and continuity equation (11d)] involving five additional variables u,ρ,P. In particular, given the magnetic field Bx,t at some instant t, the five variables u,ρ,P are determined by the five Eqs. (11a), (11b), and d(11d), which represent a first order nonlinear system of PDEs. Finally, we note that the Hall term, proportional to κ, in the induction equation (11c), is not a leading-order term according to the ordering specified in Eq. (8), and can therefore be omitted. However, as previously mentioned, this term would emerge in reductions starting from Hall MHD under an alternative ordering, such as the one given in Eq. (34) below. Importantly, including this effect does not impact the analysis of the reduced equations.
In the following, we are concerned with non-vacuum configurations such that the electric current and the magnetic field are not collinear, i.e., ×B×B0. By dotting Eq. (11a) with u and using Eq. (11b) we thus see that
(12)
where αx,t and βx,t are functions determined by system (11), provided that it admits solutions. Equation (12) implies that the dominant contributions to the flow velocity u are either in the direction of the electric current, or along the magnetic field itself.
Next, consider the ordering condition involving B/∂t = −∇ × E in (8) when the electric field is given by the generalized Ohm’s law (6) and the velocity field has expression (12). We know that uϵ, so that Eϵ. On the other hand, (8) prescribes that Bt×Eϵ2. Now recall that E = (κ/ρα)(∇ × B) × B at leading order. We thus see that, for ∇ × Eϵ2 to hold, the electric field must be of the type E = (κ/ρα)(∇ × B) × B = −ϵ∇Ψ + ϵ2a for some function Ψx,t and vector field ax,t. Note also that Bt=ϵ2×a, where the vector field a carries the rotational component of E. More precisely, for all t ≥ 0 there exists a function Ψx,t such that
(13)
where the factor ϵ/T in front of Ψ emphasizes the fact that Ψ scales as BL−2Ψ and that it has dimensions of magnetic flux, and ρ0R, ρ0 > 0, denotes a reference mass density. The function Ψ should be single-valued, because it corresponds to the electrostatic potential. Equation (13) implies that the magnetic field possesses approximate flux surfaces defined by level sets of Ψ. Recalling the boundary conditions (7), from (13) we must also have
(14)
where in this notation the sign ± depends on the orientation of ×B×B on Ω. Taking the square of both sides of this equation shows that Oϵ=0 on Ω. Hence,
(15)
This result also implies that Ψ is constant on Ω.
Since the magnetic field is solenoidal, Eq. (13) is sufficient to infer that there exists a single valued function Θ such that
(16)
in any sufficiently small neighborhood U ⊆ Ω provided that B is sufficiently regular (Lie-Darboux theorem37). Note that the normalization factor Lμ0ρ0/T must be kept because in general both B and ρ0 can be large [only the size of the ratio TB/Lμ0ρ is controlled by the ordering (8)]. The term proportional to Oϵ allows for small disturbances around the base field ∇Ψ × ∇Θ. The local nature of the Clebsch representation (16) can be overcome by promoting Θ to a multi-valued potential (angle-type variable). In the following, we thus assume a global representation for the magnetic field given by the leading order parametrization
(17)
It turns out38,39 that if a non-vanishing magnetic field B in a toroidal volume Ω is endowed with toroidal nested flux surfaces ζ, it can be written as B = ∇ζ × ∇Θζ where the multivalued potential Θζ is given by
(18)
Here, μ and ν denote toroidal and poloidal angle variables, ιζ the so-called rotational transform (the number of poloidal transits per toroidal transit of magnetic fields lines), and χ a single valued function. Note that however we kept a distinction between Ψ and ζ. Indeed, in our analysis sometimes it will be useful to introduce time-independent coordinates ζx, μx, and νx, and regard the potentials Ψ and Θ in Eq. (17) as functions of these coordinates and time t, with the time dependence of
(19)
carried by the rotational transform ιΨ and the single valued part χx,t.
Substituting the Clebsch representation (17) into the third equation in system (11) for magnetic field induction, defining
(20)
and using standard vector identities leads to
(21)
This equation can be rearranged as
(22)
If B0, the vector fields ∇Ψ and ∇Θ are linearly independent. Therefore, by dotting this equation by ∇Ψ and ∇Θ one finds
(23a)
(23b)
where fΨ,Θ and gΨ,Θ are functions of Ψ and Θ such that ∂f/Ψ = −∂g/Θ. Next, observe that the first (momentum) equation in the MHD system (11) implies that equilibria without flow satisfy
(24)
with P0=P0Ψ,Θ the pressure at the instant t = t0 in which the system is at equilibrium. Similarly, when /∂t = 0 the induction equation (11c) implies that
(25)
for some function Φ and where Q0 is the value of the function
(26)
at the instant t = t0. Equation (25) implies that Q0 is a function of P0. Comparing these results with steady states of (23), we find that at equilibrium
(27)
In order to fulfill the constraint ∂f0/Ψ = −∂g0/Θ, we demand that P0=P0Ψ so that f0 = 0 and g0=μ0Q0P0dP0/dΨ. The induction equation for the magnetic field can thus be written as
(28a)
(28b)
Note that solutions of system (28) produce exact time-dependent solutions of system (11) such that steady states without flow have pressure P0Ψ and Q0=Q0P0. Using Eq. (12), the vector field ξ can be written as
(29)
Recalling the Clebsch representation (17), system (28) can be equivalently expressed as
(30a)
(30b)
The two equations appearing in (30) can be regarded as a dynamical system which determines the nonlinear evolution of the magnetic field. Here, the function Q is evaluated through α, ρ, and P, which are determined from the solution of the Eqs. (11a), (11b), and d(11d) for the variables u,ρ,P.
A simple closure of system (30) can be obtained through the following reasoning. Suppose that we are interested in knowing the evolution of the magnetic field around some equilibrum (1) that we have obtained at some instant t = t0, for example, within a stellarator. At t = t0 the electric field is irrotational since 0 = tB = −∇ × E. Recalling (6) we have (25) so that Q is a function of P. Furthermore, we may identify the flux function with the equilibrium pressure field, P0 = λΨ with λR a constant bearing units of pressure over magnetic flux, without loss of generality. When the system is perturbed at some t = t1, we may consider a regime in which the fields u, ρ, and P react passively to changes in the magnetic field so that the functional form of α, ρ, and P is preserved for tt1. This amounts to assuming relations of the type
(31)
Then, system (30) reduces to an independent nonlinear system of two coupled PDEs for the variables Ψ and Θ,
(32a)
(32b)
A possible choice of boundary conditions for this closed dynamical system is
(33)
Notice that these boundary conditions are compatible with B · n = 0 and P = constant on Ω. Furthermore, since Ψ/∂t = −Q0∇ × B · ∇Ψ and Ψ is constant on Ω, they also imply ∇ × B · n = 0 on Ω. We also observe that since in general Θ is a multi-valued (angle) variable, for computational purposes it may be convenient to perform a change of variables. For example, if Ω is a toroidal volume spanned by coordinates ζ,μ,ν where level sets of ζ define nested toroidal surfaces within Ω, and μ, ν are toroidal and poloidal angles, one may set Θ = ιμ + ν + χ, and consider system (32) in terms of Ψ and χ, with ι′ = /dΨ = 0 on Ω and where χx,t is a single-valued function satisfying χn=ιμ+νn on Ω.
In the following, we shall focus our attention on the derived dynamical systems (11) and (30) and, in particular, on the model system (32). Finally, we observe that while the ordering (8) was considered for its relevance in stellarator applications, the same governing equations (11), (30), and (32) can be obtained under more general orderings. For example,
(34a)
(34b)

In this section we first discuss the invariants of systems (11), (30), and (32). Then, these invariants are used to construct a variational principle describing steady configurations of system (32). These steady states correspond to MHD equilibria (1) and can be understood as the result of a constrained relaxation process in which the weakest invariant is dissipated while the others are kept constant.

System (11) is endowed with invariants. To see this, first observe that the induction equation therein can be written as
(35)
where we used (12) and (26). Let us now make some considerations on the appropriate boundary conditions for the magnetic field B in Eq. (35). Remember that, in general, Q may depend on B according to (11). Hence, Eq. (35) involves at least second-order spatial derivatives of the magnetic field. By analogy with the magnetic diffusion equation tB = ΔB with ∇ · B = 0, the boundary condition B · n = 0 alone is insufficient to uniquely determine the solution. Therefore, at least one additional boundary condition, such as the requirement (∇ × B) · n = 0 in Eq. (7), is necessary to ensure a well-posed problem.
The magnetic energy of the system is given by
(36)
From (35), it follows that
(37)
where we have used the boundary conditions (7). This shows that the magnetic energy MΩ is an invariant of system (11). Let us verify that MΩ is also an invariant of the reduced systems (30) and (32) under the first boundary condition in (33), i.e. Ψ = constant on Ω. Noting that Ψ/∂t = −Q∇ × B · ∇Ψ = 0 as well as n × ∇Ψ = 0 on Ω, and considering the case in which Θ/∂t is single-valued, we have
(38)
Observe that for Θ = ιμ + ν + χ with ιR, the quantity Θ/∂t = ∂χ/∂t remains single-valued. However, if ι=ιΨ, Θ/∂t is not single valued. In this case conservation of magnetic energy still holds,
(39)
Here, we used the same boundary conditions Ψ/∂t = 0 and n × ∇Ψ = 0 on Ω.
Next, consider the magnetic helicity
(40)
where Ax,t is a single-valued vector potential such that B = ∇ × A. We have
(41)
On the other hand, the induction equation (35) implies that
(42)
where qx,tkercurl, i.e. ∇ × q = 0. Due to gauge freedom, the vector field q can be absorbed in the definition of A without loss of generality so that at the boundary we have
(43)
where we used the boundary conditions (7). Hence, the magnetic helicity KΩ is an invariant of system (11).
The magnetic helicity KΩ is also an invariant of (30) and (32) under the boundary condition Ψ = constant on Ω. However, it degenerates to a trivial invariant KΩ = 0 when Ω is connected.40 To see this, first assume that Θ/∂t and ∇Θ are single-valued, and consider the single-valued vector potential
(44)
with q0xkercurl a time-independent gauge vector field such that ∇ × q0 = 0. We have
(45)
as well as
(46)
where in the last passage we used the fact that Ψ/∂t = −Q∇ × B · ∇Ψ = 0 on Ω due to the boundary conditions (33), the fact that n × ∇Ψ = 0 on Ω, and the fact that Θ/∂t is single-valued. If ∇Θ and Θ/∂t are not single-valued, with Θ given by (19), conservation of magnetic helicity still holds. Indeed, setting q0=μΨdιdΨdΨ (note that now q0 is multi-valued), and considering the single-valued vector potential
(47)
one has
(48)
where we defined ΞΨ=ξdΨ with ξΨ=ΨdιdΨdΨ. Then, the rate of change in helicity is
(49)
This integral can be written as a surface integral, which vanishes due to the boundary conditions Ψ/∂t = 0 and ∇Ψ × n = 0 on Ω.
These results show that the magnetic helicity KΩ defined in Eqs. (45) and (48) is an invariant of systems (30) and (32). However, if the boundary Ω is a connected surface (e.g. the boundary of a solid toroidal volume), magnetic helicity can be written as
(50)
where Ψ|Ω=ΨeR and Ξ|Ω=ΞeR are the boundary values of Ψ and Ξ. This shows that the magnetic helicity KΩ becomes a trivial invariant whenever Ω defines a connected surface. Nevertheless, KΩ is nontrivial when Ω is not a connected surface (such as when Ω is the boundary of a hollow toroidal volume).
Now consider the functional
(51)
where f is any function of Ψ. When f = Ψ, FΩ is the total magnetic flux
(52)
On the other hand, recall that system (30) corresponds to the induction equation of system (11) under the assumption (17) regarding the existence of approximate flux surfaces. Therefore, the rate of change of FΩ following from system (30) is given by
(53)
When Q=Q0Ψ this integral identically vanishes. This shows that the functional FΩ is an invariant of (32). It should be emphasized that the conservation of FΩ is, in general, a property that is favorable to plasma confinement. Indeed, if we regard the density ρρΨ as a function of the magnetic flux, conservation of total particle number ΩρΨdV implies, for example, that large migrations of particles from regions of large Ψ to regions of low Ψ cannot occur without breaking the constancy of the density weighted magnetic flux ΩρΨΨdV.

The invariants of systems (11) and (32) are summarized in Tables I and II respectively.

TABLE I.

Invariants of system (11). Field conditions for the conservation of magnetic helicity KΩ specify the gauge q/∂t of the vector potential A.

InvariantExpressionField conditionsBoundary conditions
Magnetic energy MΩ 12μ0ΩB2dV  B · n = 0, ∇ × B · n = 0 
Magnetic helicity KΩ 12ΩABdV qt=AtQ×B×B=0 B · n = 0, ∇ × B · n = 0 
InvariantExpressionField conditionsBoundary conditions
Magnetic energy MΩ 12μ0ΩB2dV  B · n = 0, ∇ × B · n = 0 
Magnetic helicity KΩ 12ΩABdV qt=AtQ×B×B=0 B · n = 0, ∇ × B · n = 0 
TABLE II.

Invariants of system (32). The magnetic helicity KΩ degenerates to a trivial invariant KΩ = 0 when Ω is a connected surface.

InvariantExpressionField conditionsBoundary conditions
Magnetic energy MΩ 12μ0ΩB2dV B = ∇Ψ × ∇Θ Ψ = constant 
Magnetic helicity KΩ 12ΩABdV B = ∇Ψ × ∇Θ, A = q0 + Ψ∇Θ Ψ = constant, 
Ω not connected 
Magnetic flux FΩ ΩfΨdV B = ∇Ψ × ∇Θ Ψ = constant 
InvariantExpressionField conditionsBoundary conditions
Magnetic energy MΩ 12μ0ΩB2dV B = ∇Ψ × ∇Θ Ψ = constant 
Magnetic helicity KΩ 12ΩABdV B = ∇Ψ × ∇Θ, A = q0 + Ψ∇Θ Ψ = constant, 
Ω not connected 
Magnetic flux FΩ ΩfΨdV B = ∇Ψ × ∇Θ Ψ = constant 
In the remaining part of this section we restrict our attention to the variational formulation of steady states associated with the model system (32). From a physical standpoint one expects equilibrium states to correspond to critical points of an energy functional. In the present setting, the energy involved is that associated with the magnetic field B = ∇Ψ × ∇Θ. Here, we consider the target functional
(54)
Note that WΩ comprises the magnetic energy MΩ and the functional of the magnetic flux FΩ, which are invariants of system (32). In practice, non-ideal processes involving dissipation result in faster violation of ideal invariants that include higher order derivatives of the dynamical variables.20 Hence, we expect equilibrium states to be the result of the minimization by dissipation of MΩ under the constraint of preserved magnetic flux FΩ. The constant λ thus plays the role of a Langrange multiplier. The relaxation scenario described above is analogous to so-called Taylor relaxation in which a Beltrami state is produced as a result of a dissipation process in which magnetic energy is minimized under the constraint of magnetic helicity.18 Notice also that Taylor relaxation can be expected in the context of system (11), whose invariants include the magnetic energy MΩ and the magnetic helicity KΩ.
To simplify the analysis, in the reminder of this study we shall restrict ourselves to the case in which both Θ/∂t and ∇Θ are single-valued [if Θ is given by (19), this means ιR]. We remark that however the case in which /dΨ ≠ 0 can be handled in a similar fashion by working with Ψ and the single-valued component χ. Then, one can verify41,42 that setting to zero the first variation of the functional WΩ with respect to Ψ and Θ under the assumption that δΨ and δΘ identically vanish on the boundary, or δΨ = 0 and Ψ = constant on the boundary, gives the system of equations
(55a)
(55b)
which are equivalent to the MHD equilibrium equations (1) with pressure P = λf upon substitution of B = ∇Ψ × ∇Θ. Notice also that solutions of (55) give steady solutions of (32) when the choice f = Ψ is made.

The aim of this section is to show that the model system (32) is endowed with a Hamiltonian structure. This property will later be used to discuss the nonlinear stability of steady solutions.

Let δF/δΨ denote the functional derivative of the functional F:XR on the state space X with respect to ΨX. We have

Proposition 1.
System (32) is a Hamiltonian system with Poisson bracket
(56)
and Hamiltonian
(57)

To verify Proposition 1, first recall that in a Hamiltonian system the evolution of a physical observable FX* satisfies the equation of motion
(58)
where X is a vector space, X* the set of differentiable functionals F:XR, and HΩX* the Hamiltonian, and the Poisson bracket ,:X*×X*X* satisfies the axioms of bilinearity, alternativity, antisymmetry, Leibniz rule, and Jacobi identity,
(59a)
(59b)
(59c)
(59d)
(59e)
a,bR and F,G,HX*. For completeness, we remark that antisymmetry (59c) follows from bilinearity (59a) and alternativity (59b) by evaluation of F+G,F+G. The Leibniz rule (59d) ensures that the Poisson bracket acts as a differential operator, while the Jacobi identity (59e) assigns the Lie-algebra structure.
In the present setting, the state space X denotes the function space to which the dynamical variables Ψ and Θ of system (32) belong, while X* represents the vector space of differentiable functionals of Ψ and Θ. Below, we show that system (32) can be cast in the noncanonical Hamiltonian form
(60a)
(60b)
To see this, let us first verify that the bracket (56) correctly generates system (32). We have
(61)
where we have assumed that variations δΘ vanish on the boundary Ω when evaluating δHΩ/δΘ. The same result can be obtained by setting Ψ = constant on Ω instead of δΘ = 0 on Ω. Similarly,
(62)
where we have assumed that variations δΨ vanish on the boundary Ω when evaluating δHΩ/δΨ. This shows that the bracket (56) generates system (32) according to (60).
We are now left with the task of verifying that the bracket (56) satisfies the Poisson bracket axioms (59). The verification of bilinearity (59a), alternativity (59b), antisymmetry (59c), and Leibniz rule (59d) is immediate. Denoting with ↻ summation of even permutations, and introducing the simplified notation FΨ = δF/δΨ for functional derivatives, the Jacobi identity (59e) can be evaluated as
(63)
where we used the fact that terms involving second order functional derivatives21–23 of F, G, and H and terms containing dQ0/dΨ cancel upon summation of even permutations. Equation (63) shows that the bracket (56) also satisfies the Jacobi identity, and thus it is a Poisson bracket.
Let χ=Ψ,ΘX denote a point in the space of solutions X of system (32). The Hamiltonian nature of system (32) implies that critical points of the Hamiltonian HΩ, i.e. points χ0=Ψ0,Θ0X such that the first variation of HΩ vanishes,
(64)
correspond to steady states of system (32). Indeed, at χ0 both δHΩ/δΨ and δHΩ/δΘ vanish, implying that Ψ/∂t = Θ/∂t = 0 [recall Eq. (60)]. Information on the stability of a critical point χ0 can be gained with the aid of conservation of energy if norms 1:XR and 2:XR and positive real constants C,C can be found such that
(65)
Here, the simplified notation χt=χx,t has been used. Note that if Eq. (65) holds, then a solution χ of system (32) that initially differs from the critical point χ0 by the amount δχ0=χ0χ0 remains close to χ0 at all later times t ≥ 0 as prescribed by the norms 1 and 2. Usually, 1 and 2 satisfy the norm axioms, and assign a topology to the state space X. However, this requirement can be relaxed by demanding 1 and 2 to be positive real valued functions that provide some measure of the distance between two points in some subspace of the state space X. Furthermore, it is common to have a single distance function =1=2. Using the expression (57) for the Hamiltonian HΩ and the boundary conditions δΨ = δΘ = 0 on Ω or δΨ = 0 and Ψ = constant on Ω, one finds that
(66)
Due to the presence of second order and third order terms in the variations δΨ and δΘ without a definite sign, the difference (66) cannot be used to define the functions 1 and 2 in general. However, the situation is different if we consider perturbations that involve only one of the Clebsch potentials Ψ and Θ, i.e., if either δΨ or δΘ vanishes. In particular, we have the following

Proposition 2.
Critical points χ0=Ψ0,Θ0 of system (32) are nonlinearly stable against perturbations of Ψ0 in the distance ΨΘ2=12ΩΨ×Θ2dV. In particular, for all t ≥ 0
(67)

The Proof of Proposition 2 follows by evaluating the difference (66) for δΘ = 0. Similarly,

Proposition 3.
Critical points χ0=Ψ0,Θ0 of system (32) are nonlinearly stable against perturbations of Θ0 in the distance ΘΨ2=12ΩΨ×Θ2dV. In particular, for all t ≥ 0
(68)

Again, the Proof of Proposition 3 follows by evaluation of the difference (66) for δΨ = 0.

We conclude this section by observing that the distance functions ΨΘ2 and ΘΨ2 behave as seminorms due to the degeneracy brought by the cross product. For example, the distance function ΨΘ2 is degenerate since ΨΘ2=Ψ+gΘΘ2 for any function gΘ. Such degeneracy can be removed by restricting the state space X to include only those functions Ψ such that ΣΘΨdS=0 where ΣΘ=xΩ:Θx=cR is a level set of Θ with surface element dS. Then, ΨΘ2=0 if and only if Ψ=gΘ. On the other hand, ΣΘΨdS=gΘΣΘ=0 if and only if g = 0. Alternatively, one may simply interpret Propositions 2 and 3 as constraints on the size of ΨtΨ0 and ΘtΘ0 across ∇Θ0 and ∇Ψ0 respectively (the gradients of the variations are constrained on the submanifolds defined by level sets of Θ0 and Ψ0).

In this section we are concerned with the following system of two nonlinear PDEs in Ω,
(69a)
(69b)
subject to Q0Ψ>0 and the boundary conditions (33) and where the positive constants γ, σ bear physical units. System (69) has the following properties:

Proposition 4.
Steady states of system (6) correspond to MHD equilibria
(70)
with B = ∇Ψ × ∇Θ. Furthermore, the energy HΩ is progressively dissipated
(71)

Hence, if solutions Ψ,Θ exist in the limit t → +, they are nontrivial critical points of HΩ. We therefore suggest that the dissipative system (69) can be applied to numerically compute MHD equilibria with nested flux surfaces.

To prove Proposition 4, let us first explain how system (69) is derived. Consider an n-dimensional Hamiltonian system
(72)
where z=z1,,zn are the phase space coordinates, H=Hz the Hamiltonian function (energy) of the system, and Jij=Jji, i, j = 1, …, n, the components of the Poisson tensor J. The antisymmetry of J ensures conservation of energy according to Ḣ=Hiżi=JijHiHj=0. Here, the notation Hi = ∂H/∂zi was used. Applying the Poisson tensor twice to the Hamiltonian H has the opposite effect: the dynamical system
(73)
where gjk = gkj are the components of a symmetric positive semi-definite covariant 2-tensor, dissipates the energy H according to
(74)
where the inequality follows from the positive semi-definiteness of the tensor g. The type of relaxation described by system (73) is called “double bracket dissipation,”28,29 because it is obtained by repeated application of the Poisson tensor J defining the Poisson bracket f,g=fiJijgj. Double bracket dissipation can be used as an efficient method to compute nontrivial steady states of ideal dynamical systems such as the Euler equations because the resulting equations preserve the Casimir invariants spanning the kernel of the Poisson tensor.30,31
Now recall that system (32) has the Hamiltonian form (60), which can be equivalently written as
(75)
where
(76)
is the contravariant symplectic Poisson 2-tensor. Double bracket dissipation for system (75) can therefore be obtained according to
(77)
where the constant diagonal covariant 2-tensor
(78)
which plays the role of g in Eq. (73), serves the purpose of keeping the consistency of physical units. One can verify that evaluation of (77) gives the anticipated system (69).
Let us now examine the main properties of system (69). Setting Ψ/∂t = 0 and Θ/∂t = 0 it is immediately clear that the resulting solution Ψ,Θ is such that the magnetic field B = ∇Ψ × ∇Θ solves Eq. (70). Now consider the rate of energy change,
(79)
The boundary condition Ψ = constant on Ω implies Ψ/∂t = 0 there. Furthermore, recall that the unit outward normal to Ω satisfies n × ∇Ψ = 0. Using (69) we thus find that
(80)
In the last passage we used the hypothesis Q0Ψ,γ,σ>0. This inequality shows that if the limit
(81)
exists, the corresponding solutions Ψ,Θ=limt+Ψ,Θ of the dissipative system (69) are nontrivial solutions of (70). System (69) therefore provides a dynamical method to search for MHD equilibria that correspond to critical points δHΩΨ,Θ=0 of HΩ. Note that the key difference between the Hamiltonian system (32) and the dissipative system (69) is that they allow the search for critical points in the state space X along different directions. More precisely, in the former case if a steady solution is found it corresponds to a critical point having the same value of HΩ as that associated with initial conditions due to conservation of energy, while in the latter case steady solutions will exhibit a lower value of HΩ as a consequence of (80).
In this section we propose an iterative scheme to construct nontrivial MHD equilibria of the type (70) in Ω. The iterative scheme is based on the observation that nontrivial steady solutions of system (32) are given by the following system of two coupled PDEs for the Clebsch potentials Ψ,Θ in Ω:
(82a)
(82b)
In general, taking two arbitrary functions Ψ0,Θ0X they will not solve system (82). The idea is to introduce new Clebsch potentials Θ0, Ψ1, Θ1, Ψ2, Θ2, …, Θi−1, Ψi, … by iteratively solving Eqs. (82a) and (82b) in the volume Ω starting from some function Ψ0x, that is
(83a)
(83b)
(83c)
(83d)
(83e)
(83f)
(83g)
(83h)
(83i)
and so on, so that limiΘi1,Ψi=Θ,Ψ hopefully converges to a regular solution Θ,Ψ of system (82). Here, we observe that the first equation in the iteration (83a) serves the purpose of determining Θ0x from the “initial condition” Ψ0. In this context, a possible choice of boundary conditions is (33). The following theorem states that, if at each step of the iteration solutions of a prescribed regularity exist, then the iteration converges to a solution of (82) [and (70)] with the same regularity. More precisely we have

Theorem 1.
Assume μ0λ ≠ 0 and consider an iterative scheme in which Eqs. (82a) and (82b) are solved alternately in Ω
(84a)
(84b)
starting from a given pair Θ0x,Ψ0xX such that
(85)
with ∇Ψ0 × ∇Θ00. Suppose that during the iteration solutions exist and are nontrivial, i.e., ∇Ψi × ∇Θi0 for i ≥ 1. Further assume that the limit
(86)
exists. Then, the pair Θ,Ψ solves Eq. (82). Furthermore, the vector field B = ∇Ψ × ∇Θ defines a nontrivial MHD equilibrium solving Eq. (70).

Proof.
We must show that the iteration procedure described above converges to the desired solution Θ,Ψ of (82). We begin by noting that when a solution Ψi of Eq. (84a) is found, the quantity
(87)
is also a solution of (84a). Therefore, we can restrict the space of solutions X to those pairs Θi1,Ψi such that
(88)
Next, we define
(89)
where Ψi,ΘjX, i, j = 0, 1, 2, 3, …, are determined in Ω iteratively according to Eq. (84), and the last inequality follows from the property (88). Next, observe that setting Ψ = Ψ1, Θ = Θ0, δΨ = Ψ0 − Ψ1, and δΘ = 0 from Eq. (66) one has
(90)
Suppose that H00H10 = 0. Then, Ψ1=Ψ0+fΘ0 for some function f of Θ0. This implies that we have found a nontrivial solution of system (82) given by Θ = Θ0 and Ψ = Ψ0. This solution also defines an MHD equilibrium (70) with B = ∇Ψ0 × ∇Θ0. Indeed, Θ0 is, by construction, a solution of (85), while
(91)
We may therefore restrict our attention to the case H00 > H10. In a similar manner, one finds
(92)
with δΘ = Θ0 − Θ1. Again, the case H10 = H11 gives a nontrivial solution B = ∇Ψ1 × ∇Θ0 of system (82) since Θ1=Θ0+gΨ1 for some smooth function gΨ1. Indeed,
(93a)
(93b)
Hence, either one finds a solution, or H00 > H10 > H11. Repeating this procedure one may therefore construct a decreasing sequence
(94)
where the last inequality follows from the fact at each step of the iteration solutions are nontrivial (in particular ∇Ψi × ∇Θi0 by hypothesis while ∇Ψi × ∇Θi−10 since μ0λ ≠ 0 for all i ≥ 1). As explained above, the decreasing sequence may be interrupted if two contiguous steps possess the same energy, implying that a solution has been found after a finite number of iterations.
Next, consider the pair Θi1,Ψi. We may quantify the degree at which Θi1,Ψi fails to be a solution of system (82) through the functional
(95)
where δΘi = Θi−1 − Θi. Indeed, when ΔHi = 0, one has δΘi=gΨi for some smooth function gΨi. Hence,
(96a)
(96b)
Now suppose that many iterations (84) are carried out, and that at each step i the functional ΔHi is evaluated. We claim that
(97)
To see this, first note that ΔHi < H00 < + by construction. Next, let ΔHn ∈ (0, H00) denote the value of ΔHi at the nth iteration. Evidently, the number of times ΔHi can be equal to or exceed ΔHn is limited to be at most the natural number N such that N < H00HnN + 1. Hence, after a sufficiently large number m of iterations we must have ΔHi < ΔHn for all im > n. Now suppose that although the sequence ΔHi is decreasing, its limit inferior does not reach zero, i.e.,
(98)
for some δR. However, this is a contradiction, since one can always find a natural number M such that > H00. We must therefore conclude that
(99)
In particular, this implies that there exists some sufficiently large number of iterations s such that 0 ≤ ΔHsϵ for any arbitrarily small ϵ > 0. By the same argument as above, if ΔHs > 0 the sequence ΔHi can equal or exceed ΔHs a finite number of times after which
(100)
The limit of the sequence ΔHi therefore belongs to [0, ϵ] for any arbitrarily small ϵ > 0. Equation (97) thus follows. This result implies
(101)
Recalling (86), it follows that limiδΘi=gΨ. We therefore find
(102a)
(102b)
Hence, the pair
(103)
defines a solution of system (82). Furthermore, μ0λ ≠ 0 in (82) implies that the vector field B = ∇Ψ × ∇Θ is non-vanishing and that it is a solution of the MHD equilibrium equations (70).

Finally, we observe that if a solution Θn1,Ψn of system (82) is found after a finite number of iterations n, successive iterations will simply return the same solution, e.g., Θn,Ψn+1=Θn1,Ψn. Hence, in this case Θ,Ψ=Θn1,Ψn.□

Remark 1.

If Ω is a hollow toroidal volume with boundary Ω corresponding to two distinct level sets of a smooth function Ψ0CΩ, with ∇Ψ00 in Ω, and if level sets of Ψ0 foliate Ω with nested toroidal surfaces, Theorem 1 of Ref. 42 ensures that Eq. (85) always has a nontrivial solution Θ0 such that ∇Ψ0 × ∇Θ00. Furthermore, the angle variable Θ0 is not unique, but solutions exist in the form Θ0 = + + χ0, where μ, ν are toroidal and poloidal angle variables, the integers M, N determine the rotational transform of the vector field ∇Ψ0 × ∇Θ0, and the function χ0x is single-valued. The same result applies when solving for Θi at any step (84b) of the iteration provided that Ψi satisfies the same properties listed above for Ψ0 in Ω.

Remark 2.

An argument analogous to that used in the Proof of Theorem 1 in Ref. 42 shows that for a given angle variable Θi−1 in (84a), a solution Ψi can be obtained by reducing Eq. (84a) to a two-dimensional elliptic equation on each level set of Θi−1 and by joining solutions corresponding to adjacent level sets.

Remark 3.

In light of Remarks 1 and 2 above, if one could show that at each step of the iteration the solutions Θi−1 and Ψi, i ≥ 1 preserve their properties (in particular, Θi remains an angle variable and Ψi foliates the domain with nested toroidal surfaces) then, combining this result with Theorem 1 proved in this section, one would have obtained a proof of the existence of MHD equilibria (70) in hollow toroidal volumes of arbitrary shape. In such construction, although no control is available on the form of the flux surfaces Ψ within Ω, one can conjecture that, if they exist, solutions B = ∇Ψ × ∇Θ with different rotational transforms can be obtained by appropriate choice of the integers M, N mentioned in Remark 1.

In this study, starting from the ideal MHD equations (2) and with the aid of Clebsch potentials, we derived a reduced set of Eqs. (11) and (30), as well as (32) describing the nonlinear evolution of magnetic field turbulence in proximity of MHD equilibria (1). The ordering (8) used to arrive at these equations is appropriate for a plasma with small flow, small electric current, approximate flux surfaces, and slow time variation. This setting is expected to be relevant for stellarator plasmas. The same governing equations can be obtained under the more general ordering (34) in which both the plasma flow and the electric current are not small. We showed that the reduced equations possess invariants. In particular, the closed system (32) preserves magnetic energy, magnetic helicity, and total magnetic flux under suitable boundary and gauge conditions. Furthermore, it exhibits a noncanonical Hamiltonian structure with Poisson bracket (56) and Hamiltonian (57) (Proposition 1 of Sec. IV). Such Hamiltonian structure can be used to examine the stability properties of steady solutions: we found that MHD equilibria (70) are nonlinearly stable against perturbations involving a single Clebsch potential in the sense of Propositions 2 and 3 of Sec. V. The Hamiltonian structure can also be applied to obtain a dissipative dynamical system (69) with the property that the Hamiltonian of the system is progressively dissipated as described by Proposition 4 of Sec. VI. System (69), which comprises two coupled diffusion equations for the Clebsch potentials, thus provides a dynamical method to compute nontrivial MHD equilibria (70) by minimizing the Hamiltonian (57). We further proposed a second scheme to compute MHD equilibria (70) based on the iterative solution of the two coupled Eq. (82). Here, Theorem 1 shows that, if solutions exists at each step of the iteration, the process must converge toward a solution of system (82) and thus to a nontrivial MHD equilibrium of the type (70).

The reduced equations derived in the present paper can be regarded as a toy model of turbulence that can be useful to assess dynamical accessibility and stability of MHD equilibria in physically relevant regimes. Furthermore, they provide two practical approaches (a dissipative one and an iterative one) to numerically compute MHD equilibria. Finally, as outlined in the remarks at the end of Sec. VII, we conjecture that the iterative scheme of Sec. VII may represent the basis for a mathematical proof of the existence of MHD equilibria with a non-vanishing pressure gradient in hollow tori of arbitrary shape (that is, configurations in which the boundary Ω is not invariant under some combination of Euclidean isometries).

N.S. would like to thank Z. Yoshida for useful discussion.

The research of N.S. was partially supported by JSPS KAKENHI Grant Nos. 21K13851. and 22H04936.

The authors have no conflicts to disclose.

Naoki Sato: Conceptualization (lead); Formal analysis (lead); Funding acquisition (lead); Methodology (lead); Supervision (lead); Writing – original draft (lead); Writing – review & editing (lead). Michio Yamada: Conceptualization (supporting); Formal analysis (supporting); Methodology (supporting); Supervision (supporting); Writing – original draft (supporting); Writing – review & editing (supporting).

Data sharing is not applicable to this article as no new data were created or analyzed in this study.

1.
M. D.
Kruskal
and
R. M.
Kulsrud
, “
Equilibrium of a magnetically confined plasma in a toroid
,”
Phys. Fluids
1
,
265
(
1958
).
2.
C. L.
Surdo
, “
Global magnetofluidostatic fields (an unsolved PDE problem)
,”
Int. J. Math. Math. Sci.
9
,
123
130
(
1986
).
3.
Z.
Yoshida
and
H.
Yamada
, “
Structurally-unstable electrostatic potentials in plasmas
,”
Prog. Theor. Phys.
84
,
203
(
1990
).
4.
H.
Grad
, “
Reducible problems in magneto-fluid dynamic steady flows
,”
Rev. Mod. Phys.
32
,
830
(
1960
).
5.
Z.
Yoshida
and
Y.
Giga
, “
Remarks on spectra of operator rot
,”
Math. Z.
204
,
235
245
(
1990
).
6.
H.
Grad
, “
Toroidal containment of a plasma
,”
Phys. Fluids
10
,
137
(
1967
).
7.
J. W.
Edenstrasser
, “
Unified treatment of symmetric MHD equilibria
,”
J. Plasma Phys.
24
,
299
313
(
1980
).
8.
J. W.
Edenstrasser
, “
The only three classes of symmetric MHD equilibria
,”
J. Plasma Phys.
24
,
515
518
(
1980
).
9.
P.
Helander
, “
Theory of plasma confinement in non-axisymmetric magnetic fields
,”
Rep. Prog. Phys.
77
,
087001
(
2014
).
10.
E.
Rodriguez
,
P.
Helander
, and
A.
Bhattacharjee
, “
Necessary and sufficient conditions for quasisymmetry
,”
Phys. Plasmas
27
,
062501
(
2020
).
11.
M.
Landreman
and
E.
Paul
, “
Magnetic fields with precise quasisymmetry for plasma confinement
,”
Phys. Rev. Lett.
128
,
035001
(
2022
).
12.
N.
Sato
, “
Existence of weakly quasisymmetric magnetic fields without rotational transform in asymmetric toroidal domains
,”
Sci. Rep.
12
,
11322
(
2022
).
13.
P. J.
Morrison
, “
Hamiltonian description of the ideal fluid
,”
Rev. Mod. Phys.
70
,
467
521
(
1998
).
14.
H. M.
Abdelhamid
,
Y.
Kawazura
, and
Z.
Yoshida
, “
Hamiltonian formalism of extended magnetohydrodynamics
,”
J. Phys. A: Math. Theor.
48
(
23
),
235502
(
2015
).
15.
Z.
Yoshida
, “
Clebsch parameterization: Basic properties and remarks on its applications
,”
J. Math. Phys.
50
,
113101
(
2009
).
16.
Z.
Yoshida
and
P. J.
Morrison
, “
Epi-two-dimensional fluid flow: A new topological paradigm for dimensionality
,”
Phys. Rev. Lett.
119
,
244501
(
2017
).
17.
L.
Woltjer
, “
A theorem on force-free magnetic fields
,”
Proc. Natl. Acad. Sci. U. S. A.
44
,
489
(
1958
).
18.
J. B.
Taylor
, “
Relaxation of toroidal plasma and generation of reverse magnetic fields
,”
Phys. Rev. Lett.
33
,
1139
(
1974
).
19.
J. B.
Taylor
, “
Relaxation and magnetic reconnection in plasmas
,”
Rev. Mod. Phys.
58
,
741
(
1986
).
20.
Z.
Yoshida
and
S. M.
Mahajan
, “
Variational principles and self-organization in two-fluid plasmas
,”
Phys. Rev. Lett.
88
,
095001
(
2002
).
21.
P. J.
Morrison
, “
Poisson brackets for fluids and plasmas
,”
AIP Conf. Proc.
88
,
13
46
(
1982
).
22.
R.
Littlejohn
, “
Singular Poisson tensors
,”
AIP Conf. Proc.
88
,
47
66
(
1982
).
23.
P. J.
Olver
, “
The Jacobi identity
,” in
Applications of Lie Groups to Differential Equations
, 2nd ed. (
Springer-Verlag
,
New York
,
1993
), pp.
436
445
.
24.
D. D.
Holm
,
J. E.
Marsden
,
T.
Ratiu
, and
A.
Weinstein
, “
Nonlinear stability of fluid and plasma equilibria
,”
Phys. Rep.
123
,
1
116
(
1985
).
25.
C.
Tronci
,
E.
Tassi
, and
P. J.
Morrison
, “
Energy-Casimir stability of hybrid Vlasov-MHD models
,”
J. Phys. A: Math. Theor.
48
,
185501
(
2015
).
26.
G.
Rein
, “
Non-linear stability for the Vlasov–Poisson system—the energy-Casimir method
,”
Math. Methods Appl. Sci.
17
,
1129
1140
(
1994
).
27.
V. I.
Arnold
and
B. A.
Khesin
, “
Stability criteria for steady flows
,” in
Topological Methods in Hydrodynamics
(
Springer
,
1998
), pp.
89
96
.
28.
P. J.
Morrison
, “
A paradigm for joined Hamiltonian and dissipative systems
,”
Physica D
18
,
410
419
(
1986
).
29.
P. J.
Morrison
, “
Thoughts on brackets and dissipation: Old and new
,”
J. Phys.: Conf. Ser.
169
,
012006
(
2009
).
30.
M.
Furukawa
,
T.
Watanabe
,
P. J.
Morrison
, and
K.
Ichiguchi
, “
Calculation of large-aspect-ratio tokamak and toroidally-averaged stellarator equilibria of high-beta reduced magnetohydrodynamics via simulated annealing
,”
Phys. Plasmas
25
,
082506
(
2018
).
31.
G. K.
Vallis
,
G. F.
Carnevale
, and
W. R.
Young
, “
Extremal energy properties and construction of stable solutions of the Euler equations
,”
J. Fluid Mech.
207
,
133
152
(
1989
).
32.
R.
Chodura
and
A.
Schlüter
, “
A 3D code for MHD equilibrium and stability
,”
J. Comput. Phys.
41
,
68
88
(
1981
).
33.
S. P.
Hirshman
and
J. C.
Whitson
, “
Steepest-descent moment method for three-dimensional magnetohydrodynamic equilibria
,”
Phys. Fluids
26
,
3553
(
1983
).
34.
M.
Acheritogaray
,
P.
Degond
,
A.
Frouvelle
, and
J. G.
Liu
, “
Kinetic formulation and global existence for the Hall-magneto-hydrodynamics system
,”
Kinet. Relat. Models
4
(
4
),
901
918
(
2011
).
35.
R.
Fitzpatrick
, “
MHD equations
,” in
Plasma Physics: An Introduction
(
CRC Press
,
2015
), pp.
107
108
.
36.
J. P.
Freidberg
, “
The ideal MHD model
,” in
Ideal MHD
(
Cambridge University Press
,
Cambridge
,
2014
), pp.
7
35
.
37.
M.
de León
,
Methods of Differential Geometry in Analytical Mechanics
(
Elsevier
,
New York
,
1989
), pp.
250
253
.
38.
A. H.
Boozer
, “
Guiding center drift equations
,”
Phys. Fluids
23
,
904
908
(
1980
).
39.
A. H.
Boozer
, “
Establishment of magnetic coordinates for a given magnetic field
,”
Phys. Fluids
25
,
520
521
(
1982
).
40.
D.
Pfefferlé
,
L.
Noakes
, and
D.
Perrella
, “
Gauge freedom in magnetostatics and the effect on helicity in toroidal volumes
,”
J. Math. Phys.
62
,
093505
(
2021
).
41.
H.
Grad
and
H.
Rubin
, “
Hydromagnetic equilibria and force free fields
,” in
Proceedings of the 2nd United Nations International Conference on the Peaceful Uses of Atomic Energy
(
United Nations Publication
,
1958
), Vol.
31
, pp.
190
197
.
42.
N.
Sato
and
M.
Yamada
, “
Nested invariant tori foliating a vector field and its curl: Toward MHD equilibria and steady Euler flows in toroidal domains without continuous Euclidean isometries
,”
J. Math. Phys.
64
,
081505
(
2023
).