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.

## I. INTRODUCTION

^{1}

*μ*

_{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 Ω,

*P*to the other fields, and by determining

**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**

*E***=**

*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 structure^{13,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 ∇ ×

**, and time derivatives**

*B**∂*/

*∂t*are small, and the pressure field

*P*is related to the mass density

*ρ*and the velocity field

**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 representation**

*u*^{15,16}of the magnetic field

**= ∇Ψ × ∇Θ 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 Θ.**

*B*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 relaxation^{18,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 stable^{24–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 brackets^{28,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 $\u222b\Omega B22\mu 0+P\gamma \u22121dV$, where *γ* ≥ 0 is the adiabatic index, often used to numerically compute MHD equilibria.^{32,33}

Concluding remarks are given in Sec. VIII.

## II. A REDUCED IDEAL MHD SYSTEM FOR MAGNETIC FIELD TURBULENCE IN PLASMAS WITH SMALL CURRENTS AND APPROXIMATE FLUX SURFACES

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.

*P*and the other field variables will be given later in the form of a generalized Bernoulli principle. The relation between the electric field

**and the other fields is given by a generalized Ohm’s law following from the electron fluid momentum equation**

*E**m*

_{e}, −

*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*+

*u*_{e}· ∇. Now recall that

*m*

_{i},

*δ*=

*m*

_{e}/

*m*

_{i},

*n*

_{i}, $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

*n*

_{i}=

*n*

_{e}=

*n*. The electron momentum equation (3) therefore gives the generalized Ohm’s law

*P*

_{e}. Indeed, assuming a barotropic equation of state for the electron pressure, $Pe=Pe\rho $, we obtain

*κ*=

*m*

_{i}/

*eμ*

_{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***). 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**

*B***,**

*u***, and ∇ ×**

*B***, and the pressure**

*B**P*satisfy the boundary conditions

*∂*Ω denotes the boundary of Ω and

**the unit outward normal to**

*n**∂*Ω. Notice in particular that the boundary condition for ∇ ×

**implies that there is no net current flow across the bounding surface**

*B**∂*Ω, and it is expected to hold true as long as both

*u*_{e}and

*u*_{i}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.

*ϵ*> 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

**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 ∇**

*E**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,

*p*obeys

*T*

^{2}

*p*/

*ρL*

^{2}∼

*ϵ*

^{2}.

**implies that the magnetic field is endowed with approximate flux surfaces. Note that $B/\mu 0\rho $ is the Alfvén speed, and that the ordering (8) does not involve conditions on the size of the typical particle number**

*E**ρL*

^{3}/

*m*

_{i}. We now make the following ordering prescription,

*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 ∇ ×

**and the electric current ∇ ×**

*u***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**

*B***→**

*u***0**steady solutions of system (11) are described by the MHD equilibrium equations (1) and a barotropic equation of state $P=P\rho $. 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,\rho ,P$. In particular, given the magnetic field $Bx,t$ at some instant

*t*, the five variables $u,\rho ,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.

**and using Eq. (11b) we thus see that**

*u***are either in the direction of the electric current, or along the magnetic field itself.**

*u**∂*

**/**

*B**∂t*= −∇ ×

**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**

*E***∼**

*u**ϵ*, so that

**∼**

*E**ϵ*. On the other hand, (8) prescribes that $\u2202B\u2202t\u223c\u2207\xd7E\u223c\u03f52$. Now recall that

**= (**

*E**κ*/

*ρ*−

*α*)(∇ ×

**) ×**

*B***at leading order. We thus see that, for ∇ ×**

*B***∼**

*E**ϵ*

^{2}to hold, the electric field must be of the type

**= (**

*E**κ*/

*ρ*−

*α*)(∇ ×

**) ×**

*B***= −**

*B**ϵ*∇Ψ +

*ϵ*

^{2}

**for some function $\Psi x,t$ and vector field $ax,t$. Note also that $\u2202B\u2202t=\u2212\u03f52\u2207\xd7a$, where the vector field**

*a***carries the rotational component of**

*a***. More precisely, for all**

*E**t*≥ 0 there exists a function $\Psi x,t$ such that

*ϵ*/

*T*in front of Ψ emphasizes the fact that Ψ scales as

**∼**

*B**L*

^{−2}Ψ and that it has dimensions of magnetic flux, and $\rho 0\u2208R$,

*ρ*

_{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

*∂*Ω. Taking the square of both sides of this equation shows that $O\u03f5=0$ on

*∂*Ω. Hence,

*∂*Ω.

*U*⊆ Ω provided that

**is sufficiently regular (Lie-Darboux theorem**

*B*^{37}). Note that the normalization factor $L\mu 0\rho 0/T$ must be kept because in general both

**and**

*B**ρ*

_{0}can be large [only the size of the ratio $TB/L\mu 0\rho $ is controlled by the ordering (8)]. The term proportional to $O\u03f5$ 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

^{38,39}that if a non-vanishing magnetic field

**in a toroidal volume Ω is endowed with toroidal nested flux surfaces**

*B**ζ*, it can be written as

**= ∇**

*B**ζ*× ∇Θ

_{ζ}where the multivalued potential Θ

_{ζ}is given by

*μ*and

*ν*denote toroidal and poloidal angle variables, $\iota \zeta $ 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 $\zeta x$, $\mu x$, and $\nu x$, and regard the potentials Ψ and Θ in Eq. (17) as functions of these coordinates and time

*t*, with the time dependence of

**≠**

*B***0**, the vector fields ∇Ψ and ∇Θ are linearly independent. Therefore, by dotting this equation by ∇Ψ and ∇Θ one finds

*∂f*/

*∂*Ψ = −

*∂g*/

*∂*Θ. Next, observe that the first (momentum) equation in the MHD system (11) implies that equilibria without flow satisfy

*t*=

*t*

_{0}in which the system is at equilibrium. Similarly, when

*∂*/

*∂t*= 0 the induction equation (11c) implies that

*Q*

_{0}is the value of the function

*t*=

*t*

_{0}. Equation (25) implies that

*Q*

_{0}is a function of

*P*

_{0}. Comparing these results with steady states of (23), we find that at equilibrium

*∂f*

_{0}/

*∂*Ψ = −

*∂g*

_{0}/

*∂*Θ, we demand that $P0=P0\Psi $ so that

*f*

_{0}= 0 and $g0=\mu 0Q0P0dP0/d\Psi $. The induction equation for the magnetic field can thus be written as

**can be written as**

*ξ**Q*is evaluated through

*α*,

*ρ*, and

*P*, which are determined from the solution of the Eqs. (11a), (11b), and d(11d) for the variables $u,\rho ,P$.

*t*=

*t*

_{0}, for example, within a stellarator. At

*t*=

*t*

_{0}the electric field is irrotational since

**0**=

*∂*

_{t}

**= −∇ ×**

*B***. Recalling (6) we have (25) so that**

*E**Q*is a function of

*P*. Furthermore, we may identify the flux function with the equilibrium pressure field,

*P*

_{0}=

*λ*Ψ with $\lambda \u2208R$ a constant bearing units of pressure over magnetic flux, without loss of generality. When the system is perturbed at some

*t*=

*t*

_{1}, 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

*t*≥

*t*

_{1}. This amounts to assuming relations of the type

**·**

*B***= 0 and**

*n**P*= constant on

*∂*Ω. Furthermore, since

*∂*Ψ/

*∂t*= −

*Q*

_{0}∇ ×

**· ∇Ψ and Ψ is constant on**

*B**∂*Ω, they also imply ∇ ×

**·**

*B***= 0 on**

*n**∂*Ω. 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 $\zeta ,\mu ,\nu $ 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ι*/

*d*Ψ = 0 on

*∂*Ω and where $\chi x,t$ is a single-valued function satisfying $\u2207\chi \u22c5n=\u2212\iota \u2207\mu +\u2207\nu \u22c5n$ on

*∂*Ω.

## III. CONSERVATION LAWS AND RELAXED STATES

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.

### A. Conservation of magnetic energy, helicity, and flux

**in Eq. (35). Remember that, in general,**

*B**Q*may depend on

**according to (11). Hence, Eq. (35) involves at least second-order spatial derivatives of the magnetic field. By analogy with the magnetic diffusion equation**

*B**∂*

_{t}

**= Δ**

*B***with ∇ ·**

*B***= 0, the boundary condition**

*B***·**

*B***= 0 alone is insufficient to uniquely determine the solution. Therefore, at least one additional boundary condition, such as the requirement (∇ ×**

*n***) ·**

*B***= 0 in Eq. (7), is necessary to ensure a well-posed problem.**

*n**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*∇ ×

**· ∇Ψ = 0 as well as**

*B***× ∇Ψ =**

*n***0**on

*∂*Ω, and considering the case in which

*∂*Θ/

*∂t*is single-valued, we have

*ιμ*+

*ν*+

*χ*with $\iota \u2208R$, the quantity

*∂*Θ/

*∂t*=

*∂χ*/

*∂t*remains single-valued. However, if $\iota =\iota \Psi $,

*∂*Θ/

*∂t*is not single valued. In this case conservation of magnetic energy still holds,

*∂*Ψ/

*∂t*= 0 and

**× ∇Ψ =**

*n***0**on

*∂*Ω.

**= ∇ ×**

*B***. We have**

*A***=**

*q***0**. Due to gauge freedom, the vector field

**can be absorbed in the definition of**

*q***without loss of generality so that at the boundary we have**

*A**K*

_{Ω}is an invariant of system (11).

*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

*q*_{0}= 0. We have

*∂*Ψ/

*∂t*= −

*Q*∇ ×

**· ∇Ψ = 0 on**

*B**∂*Ω 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=\u2212\u2207\mu \u222b\Psi d\iota d\Psi d\Psi $ (note that now

*q*_{0}is multi-valued), and considering the single-valued vector potential

*∂*Ψ/

*∂t*= 0 and ∇Ψ ×

**=**

*n***0**on

*∂*Ω.

*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

*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).

*f*is any function of Ψ. When

*f*= Ψ,

*F*

_{Ω}is the total magnetic flux

*F*

_{Ω}following from system (30) is given by

*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 $\rho \u2248\rho \Psi $ as a function of the magnetic flux, conservation of total particle number $\u222b\Omega \rho \Psi 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 $\u222b\Omega \rho \Psi \Psi dV$.

Invariant . | Expression . | Field conditions . | Boundary conditions . |
---|---|---|---|

Magnetic energy M_{Ω} | $12\mu 0\u222b\Omega B2dV$ | · B = 0, ∇ × n · B = 0 n | |

Magnetic helicity K_{Ω} | $12\u222b\Omega A\u22c5BdV$ | $\u2202q\u2202t=\u2202A\u2202t\u2212Q\u2207\xd7B\xd7B=0$ | · B = 0, ∇ × n · B = 0 n |

Invariant . | Expression . | Field conditions . | Boundary conditions . |
---|---|---|---|

Magnetic energy M_{Ω} | $12\mu 0\u222b\Omega B2dV$ | · B = 0, ∇ × n · B = 0 n | |

Magnetic helicity K_{Ω} | $12\u222b\Omega A\u22c5BdV$ | $\u2202q\u2202t=\u2202A\u2202t\u2212Q\u2207\xd7B\xd7B=0$ | · B = 0, ∇ × n · B = 0 n |

Invariant . | Expression . | Field conditions . | Boundary conditions . |
---|---|---|---|

Magnetic energy M_{Ω} | $12\mu 0\u222b\Omega B2dV$ | = ∇Ψ × ∇Θ B | Ψ = constant |

Magnetic helicity K_{Ω} | $12\u222b\Omega A\u22c5BdV$ | = ∇Ψ × ∇Θ, B = Aq_{0} + Ψ∇Θ | Ψ = constant, |

∂Ω not connected | |||

Magnetic flux F_{Ω} | $\u222b\Omega f\Psi dV$ | = ∇Ψ × ∇Θ B | Ψ = constant |

Invariant . | Expression . | Field conditions . | Boundary conditions . |
---|---|---|---|

Magnetic energy M_{Ω} | $12\mu 0\u222b\Omega B2dV$ | = ∇Ψ × ∇Θ B | Ψ = constant |

Magnetic helicity K_{Ω} | $12\u222b\Omega A\u22c5BdV$ | = ∇Ψ × ∇Θ, B = Aq_{0} + Ψ∇Θ | Ψ = constant, |

∂Ω not connected | |||

Magnetic flux F_{Ω} | $\u222b\Omega f\Psi dV$ | = ∇Ψ × ∇Θ B | Ψ = constant |

### B. Steady states and relaxation

**= ∇Ψ × ∇Θ. Here, we consider the target functional**

*B**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*

_{Ω}.

*∂*Θ/

*∂t*and ∇Θ are single-valued [if Θ is given by (19), this means $\iota \u2208R$]. We remark that however the case in which

*dι*/

*d*Ψ ≠ 0 can be handled in a similar fashion by working with Ψ and the single-valued component

*χ*. Then, one can verify

^{41,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

*P*=

*λf*upon substitution of

**= ∇Ψ × ∇Θ. Notice also that solutions of (55) give steady solutions of (32) when the choice**

*B**f*= Ψ is made.

## IV. HAMILTONIAN STRUCTURE

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:X\u2192R$ on the state space $X$ with respect to $\Psi \u2208X$. We have

*System (32) is a Hamiltonian system with Poisson bracket*

*and Hamiltonian*

*δ*Θ vanish on the boundary

*∂*Ω when evaluating

*δH*

_{Ω}/

*δ*Θ. The same result can be obtained by setting Ψ = constant on

*∂*Ω instead of

*δ*Θ = 0 on

*∂*Ω. Similarly,

*δ*Ψ vanish on the boundary

*∂*Ω when evaluating

*δH*

_{Ω}/

*δ*Ψ. This shows that the bracket (56) generates system (32) according to (60).

*F*

_{Ψ}=

*δF*/

*δ*Ψ for functional derivatives, the Jacobi identity (59e) can be evaluated as

^{21–23}of

*F*,

*G*, and

*H*and terms containing

*dQ*

_{0}/

*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.

## V. REMARKS ON THE NONLINEAR STABILITY OF STEADY SOLUTIONS

*H*

_{Ω}, i.e. points $\chi 0=\Psi 0,\Theta 0\u2208X$ such that the first variation of

*H*

_{Ω}vanishes,

*χ*

_{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 $\u22c51:X\u2192R$ and $\u22c52:X\u2192R$ and positive real constants $C,C\u2032$ can be found such that

*χ*of system (32) that initially differs from the critical point

*χ*

_{0}by the amount $\delta \chi 0=\chi 0\u2212\chi 0$ remains close to

*χ*

_{0}at all later times

*t*≥ 0 as prescribed by the norms $\u22c51$ and $\u22c52$. Usually, $\u22c51$ and $\u22c52$ satisfy the norm axioms, and assign a topology to the state space $X$. However, this requirement can be relaxed by demanding $\u22c51$ and $\u22c52$ 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 $\u22c5=\u22c51=\u22c52$. Using the expression (57) for the Hamiltonian

*H*

_{Ω}and the boundary conditions

*δ*Ψ =

*δ*Θ = 0 on

*∂*Ω or

*δ*Ψ = 0 and Ψ = constant on

*∂*Ω, one finds that

*δ*Ψ and

*δ*Θ without a definite sign, the difference (66) cannot be used to define the functions $\u22c51$ and $\u22c52$ 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

*Critical points*$\chi 0=\Psi 0,\Theta 0$

*of system (32) are nonlinearly stable against perturbations of*Ψ

_{0}

*in the distance*$\Psi \Theta 2=12\u222b\Omega \u2207\Psi \xd7\u2207\Theta 2dV$

*. In particular, for all*

*t*≥ 0

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

*Critical points*$\chi 0=\Psi 0,\Theta 0$

*of system (32) are nonlinearly stable against perturbations of*Θ

_{0}

*in the distance*$\Theta \Psi 2=12\u222b\Omega \u2207\Psi \xd7\u2207\Theta 2dV$

*. In particular, for all*

*t*≥ 0

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 $\Psi \Theta 2$ and $\Theta \Psi 2$ behave as seminorms due to the degeneracy brought by the cross product. For example, the distance function $\Psi \Theta 2$ is degenerate since $\Psi \Theta 2=\Psi +g\Theta \Theta 2$ for any function $g\Theta $. Such degeneracy can be removed by restricting the state space $X$ to include only those functions Ψ such that $\u222b\Sigma \Theta \Psi dS=0$ where $\Sigma \Theta =x\u2208\Omega :\Theta x=c\u2208R$ is a level set of Θ with surface element *dS*. Then, $\Psi \Theta 2=0$ if and only if $\Psi =g\Theta $. On the other hand, $\u222b\Sigma \Theta \Psi dS=g\Theta \Sigma \Theta =0$ if and only if *g* = 0. Alternatively, one may simply interpret Propositions 2 and 3 as constraints on the size of $\u2207\Psi t\u2212\Psi 0$ and $\u2207\Theta t\u2212\Theta 0$ across ∇Θ_{0} and ∇Ψ_{0} respectively (the gradients of the variations are constrained on the submanifolds defined by level sets of Θ_{0} and Ψ_{0}).

## VI. CONSTRUCTION OF MHD EQUILIBRIA BY DOUBLE BRACKET DISSIPATION

*γ*,

*σ*bear physical units. System (69) has the following properties:

*Steady states of system (6) correspond to MHD equilibria*

*with*

**= ∇Ψ × ∇Θ**

*B**. Furthermore, the energy*

*H*

_{Ω}

*is progressively dissipated*

Hence, if solutions $\Psi ,\Theta $ 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.

*n*-dimensional Hamiltonian system

*i*,

*j*= 1, …,

*n*, the components of the Poisson tensor $J$. The antisymmetry of $J$ ensures conservation of energy according to $H\u0307=Hiz\u0307i=JijHiHj=0$. Here, the notation

*H*

_{i}=

*∂H*/

*∂z*

^{i}was used. Applying the Poisson tensor twice to the Hamiltonian

*H*has the opposite effect: the dynamical system

*g*

_{jk}=

*g*

_{kj}are the components of a symmetric positive semi-definite covariant 2-tensor, dissipates the energy

*H*according to

*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}

*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).

*∂*Ψ/

*∂t*= 0 and

*∂*Θ/

*∂t*= 0 it is immediately clear that the resulting solution $\Psi ,\Theta $ is such that the magnetic field

**= ∇Ψ × ∇Θ solves Eq. (70). Now consider the rate of energy change,**

*B**∂*Ω implies

*∂*Ψ/

*∂t*= 0 there. Furthermore, recall that the unit outward normal to

*∂*Ω satisfies

**× ∇Ψ =**

*n***0**. Using (69) we thus find that

*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).

## VII. CONSTRUCTION OF MHD EQUILIBRIA BY ITERATION

_{0}, Ψ

_{1}, Θ

_{1}, Ψ

_{2}, Θ

_{2}, …, Θ

_{i−1}, Ψ

_{i}, … by iteratively solving Eqs. (82a) and (82b) in the volume Ω starting from some function $\Psi 0x$, that is

_{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

*Assume*

*μ*

_{0}

*λ*≠ 0

*and consider an iterative scheme in which Eqs. (82a) and (82b) are solved alternately in*Ω

*starting from a given pair*$\Theta 0x,\Psi 0x\u2208X$

*such that*

*with*∇Ψ

_{0}× ∇Θ

_{0}≠

**0**

*. Suppose that during the iteration solutions exist and are nontrivial, i.e.,*∇Ψ

_{i}× ∇Θ

_{i}≠

**0**

*for*

*i*≥ 1

*. Further assume that the limit*

*exists. Then, the pair*$\Theta \u221e,\Psi \u221e$

*solves Eq. (82). Furthermore, the vector field*

**= ∇Ψ**

*B*_{∞}× ∇Θ

_{∞}

*defines a nontrivial MHD equilibrium solving Eq. (70).*

_{i}of Eq. (84a) is found, the quantity

*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

*H*

_{00}−

*H*

_{10}= 0. Then, $\Psi 1=\Psi 0+f\Theta 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

*H*

_{00}>

*H*

_{10}. In a similar manner, one finds

*δ*Θ = Θ

_{0}− Θ

_{1}. Again, the case

*H*

_{10}=

*H*

_{11}gives a nontrivial solution

**= ∇Ψ**

*B*_{1}× ∇Θ

_{0}of system (82) since $\Theta 1=\Theta 0+g\Psi 1$ for some smooth function $g\Psi 1$. Indeed,

*H*

_{00}>

*H*

_{10}>

*H*

_{11}. Repeating this procedure one may therefore construct a decreasing sequence

_{i}× ∇Θ

_{i}≠

**0**by hypothesis while ∇Ψ

_{i}× ∇Θ

_{i−1}≠

**0**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.

*δ*Θ

_{i}= Θ

_{i−1}− Θ

_{i}. Indeed, when Δ

*H*

_{i}= 0, one has $\delta \Theta i=g\Psi i$ for some smooth function $g\Psi i$. Hence,

*i*the functional Δ

*H*

_{i}is evaluated. We claim that

*H*

_{i}<

*H*

_{00}< +

*∞*by construction. Next, let Δ

*H*

_{n}∈ (0,

*H*

_{00}) denote the value of Δ

*H*

_{i}at the

*n*th iteration. Evidently, the number of times Δ

*H*

_{i}can be equal to or exceed Δ

*H*

_{n}is limited to be at most the natural number

*N*such that

*N*<

*H*

_{00}/Δ

*H*

_{n}≤

*N*+ 1. Hence, after a sufficiently large number

*m*of iterations we must have Δ

*H*

_{i}< Δ

*H*

_{n}for all

*i*≥

*m*>

*n*. Now suppose that although the sequence Δ

*H*

_{i}is decreasing, its limit inferior does not reach zero, i.e.,

*M*such that

*Mδ*>

*H*

_{00}. We must therefore conclude that

*s*such that 0 ≤ Δ

*H*

_{s}≤

*ϵ*for any arbitrarily small

*ϵ*> 0. By the same argument as above, if Δ

*H*

_{s}> 0 the sequence Δ

*H*

_{i}can equal or exceed Δ

*H*

_{s}a finite number of times after which

*H*

_{i}therefore belongs to [0,

*ϵ*] for any arbitrarily small

*ϵ*> 0. Equation (97) thus follows. This result implies

*μ*

_{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 $\Theta n\u22121,\Psi n$ of system (82) is found after a finite number of iterations *n*, successive iterations will simply return the same solution, e.g., $\Theta n,\Psi n+1=\Theta n\u22121,\Psi n$. Hence, in this case $\Theta \u221e,\Psi \u221e=\Theta n\u22121,\Psi n$.□

*If* Ω *is a hollow toroidal volume with boundary* *∂*Ω *corresponding to two distinct level sets of a smooth function* $\Psi 0\u2208C\u221e\Omega $*, with* ∇Ψ_{0} ≠ **0** *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} × ∇Θ_{0} ≠ **0***. Furthermore, the angle variable* Θ_{0} *is not unique, but solutions exist in the form* Θ_{0} = *Mμ* + *Nν* + *χ*_{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* $\chi 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* Ω.

*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.*

*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.*

## VIII. CONCLUDING REMARKS

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).

## ACKNOWLEDGMENTS

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.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

**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 AVAILABILITY

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

## REFERENCES

*Applications of Lie Groups to Differential Equations*

*Topological Methods in Hydrodynamics*

*Plasma Physics: An Introduction*

*Ideal MHD*

*Methods of Differential Geometry in Analytical Mechanics*