2.1. Helicity conservation and Taylor relaxation

The conservation of magnetic fluxes threading cuts through ${\it\Omega}$ (or ${\it\Omega}_{i}$ in MRxMHD) that leave it topologically connected follows simply from $\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{B}=0$ and the boundary condition (1.1) everywhere on $\partial {\it\Omega}$ (because of the no-gaps assumption mentioned in the Introduction). Thus, while these fluxes are IMHD invariants, they must also be invariants in any physical model; flux conservation is said to be completely robust.

As it is conserved even under reconnection (in the limit as resistivity approaches zero, see e.g. equation (20) ff. of Jensen & Chu (Reference Jensen and Chu1984)) and has a robust topological interpretation (see e.g. Berger (Reference Berger1999) for a heuristic review or Arnold & Khesin (Reference Arnold and Khesin1998) for a more mathematical treatment) the most robust of the remaining IMHD invariants is widely accepted to be the magnetic helicity $2{\rm\mu}_{0}K_{{\it\Omega}}$ , where, Bhattacharjee & Dewar (Reference Bhattacharjee and Dewar1982), we define the invariant $K_{{\it\Omega}}$ as

with $\boldsymbol{A}$ a vector potential giving $\boldsymbol{B}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{A}$ and $\text{d}V=\text{d}^{3}x$ the volume element. Because $K$ has one less gradient of $\boldsymbol{A}$ than the magnetic energy,

it can also be argued, see e.g. § I.C of Taylor (Reference Taylor1986), that, in a weakly resistive plasma with small-scale turbulent fluctuations, $K$ decays slower with time than $W^{B}$ .

Other general ideal invariants can be related by Noether’s theorems to some symmetries in appropriate parameterizations of field variables in the IMHD action (for example, the cross helicity $\int _{{\it\Omega}}\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}V$ pertains to a relabelling symmetry in the Lagrangian representation of the fields), Salmon (Reference Salmon1988), Padhye & Morrison (Reference Padhye and Morrison1996a ,Reference Padhye and Morrison b ), Webb & Zank (Reference Webb and Zank2007), Webb et al. (Reference Webb, Dasgupta, McKenzie, Hu and Zank2014a ,Reference Webb, Dasgupta, McKenzie, Hu and Zank b ), Araki (Reference Araki2015). However, in the spirit of Taylor relaxation, we choose the minimal set required to obtain a non-trivial solution and thus keep only $K_{{\it\Omega}}$ as the only non-holonomic constraint involving $\boldsymbol{v}$ or $\boldsymbol{B}$ .

It is readily shown, using the tangential- $\boldsymbol{B}$ condition (1.1), that $K_{{\it\Omega}}$ is invariant under gauge transformations $\boldsymbol{A}\mapsto \boldsymbol{A}+\boldsymbol{{\rm\nabla}}{\it\chi}$ as long as ${\it\chi}$ is single valued (implying conservation of line integrals $\oint _{\partial {\it\Omega}}\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ around loops on the boundary, which, by Stokes’ theorem, is equivalent to the above-mentioned conservation of magnetic fluxes). As tangential $\boldsymbol{B}$ is to be a holonomic constraint rather than a natural boundary condition, we do not treat $\boldsymbol{A}$ on the boundary as freely variable and can constrain ${\it\chi}$ . Thus we do not need to use either the Bevir–Gray (subtraction of products of toroidal and poloidal loop integrals, Bevir & Gray (Reference Bevir, Gray and Lewis1982)) or relative helicity (subtraction of vacuum-field helicity, Jensen & Chu (Reference Jensen and Chu1984)) modifications of the helicity, the latter fact also implying there is no physical necessity to decompose the magnetic field into a vacuum (harmonic) and a plasma-current-generated component (though it may still be useful conceptually and mathematically, Yoshida & Giga (Reference Yoshida and Giga1990), Yoshida & Dewar (Reference Yoshida and Dewar2012)).

The Woltjer–Taylor variational principle (originally proposed, though with less physical motivation, by Woltjer (Reference Woltjer1958)) states that the final relaxed state is that which minimizes the magnetic energy, (2.2) (the negative of which occurs in (1.2)), under the magnetic helicity constraint, implemented by minimizing $W_{{\it\Omega}}^{B}-{\it\mu}K_{{\it\Omega}}$ under variations of the magnetic vector potential $\boldsymbol{A}$ , ${\it\mu}$ being a Lagrange multiplierFootnote ² . The resulting Euler–Lagrange equation is the Beltrami equation,

This describes a force-free field, i.e. one with the current $\boldsymbol{j}=\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}/{\rm\mu}_{0}$ parallel to $\boldsymbol{B}$ , implying $\boldsymbol{{\rm\nabla}}p=0$ globally in an equilibrium plasma. Thus this single-region relaxation principle describes only plasmas with no thermal confinement.

2.2. Generalization of Taylor relaxation

As the Woltjer–Taylor variational principle is not explicitly based on knowledge of the detailed sub-macroscopic physics leading to relaxation and self-organization, its applicability to modelling a given system can only be justified empirically.Footnote ³ We have already remarked in the Introduction that in fact Taylor relaxation theory in its original form is too simple to apply to modern fusion devices, but that it can be extended naturally by supplementing magnetic helicity with further global ideal invariants, thus preserving much of its simplicity but increasing its flexibility in applications to modelling fusion plasmas.

The fundamental basis of our generalized MHD relaxation principles is the requirement that the states they describe be a subset of the states allowed within ideal MHD. This is ensured, Bhattacharjee & Dewar (Reference Bhattacharjee and Dewar1982), Dewar et al. (Reference Dewar, Hole, McGann, Mills and Hudson2008), by using only constraints from a subset of those implied by IMHD. We take this formal criterion as the paramount principle for constructing consistent modifications of IMHD, regarding considerations of possible subscale physics that might lead to breaking of some ideal invariants and not others only as an heuristic guide in choosing an appropriate subset of IMHD constraints. For instance, in choosing relaxation subregions in which to apply the Woltjer–Taylor variational principle, we do not necessarily assume the magnetic field is wholly or partially chaotic, though Beltrami solutions can accommodate such cases, Dombre et al. (Reference Dombre, Frisch, Greene, Hénon, Mehr and Soward1986). Justification for the choice of constraints must ultimately be empirical, by comparison either with experiment or ab initio simulations.

A generalization of the Taylor relaxation idea by increasing the number of constraints was proposed by Bhattacharjee & Dewar (Reference Bhattacharjee and Dewar1982), but the smooth IMHD invariants chosen then are not well defined in a non-integrable magnetic field with islands and chaotic regions. More recently, Hudson, Hole & Dewar (Reference Hudson, Hole and Dewar2007), Dewar et al. (Reference Dewar, Hole, McGann, Mills and Hudson2008), generalizations of Taylor relaxation theory were proposed based on the assumption that Taylor-relaxed plasma can coexist with current sheets that act as transport barriers partitioning ${\it\Omega}$ into multiple regions ${\it\Omega}_{i}$ , invariant under field-line flow, Hudson et al. (Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012). To describe this approach we have introduced the terminologies RxMHD when ${\it\Omega}$ is not partitioned, and MRxMHD when it is, the D (for ‘dynamics’) being justified below. The ${\it\delta}$ -function currents in MRxMHD are compatible with IMHD so they may be regarded as singular alternatives to the smooth IMHD constraints of Bhattacharjee & Dewar (Reference Bhattacharjee and Dewar1982). As in Bhattacharjee & Dewar (Reference Bhattacharjee and Dewar1982), we also introduce entropy constraints to allow a non-trivial pressure profile and retain the non-singular magnetic helicity invariant(s), (2.1), but separately conserved in each MRxMHD subregion.

Our development of MRxMHD is implicitly based on the multiple time scale scenario sketched below (a conceptual framework motivating the formal development; whether there are situations where it approximates physical reality remains to be explored in further work):

1. (1) A fast relaxation time scale during which all but a finite number of IMHD constraints are broken through thermal diffusion and micro-reconnection events (associated with unspecified mechanisms like micro-tearing turbulence, high-order resonant structuresFootnote ⁴ and field-line chaos). On this time scale, the system self-organizes into multiple Taylor states in disjoint subregions ${\it\Omega}_{i}$ with non-disjoint boundaries $\partial {\it\Omega}_{i}$ (geometrically fixed on this time scale due to plasma inertia) supporting current sheets on their common interfaces ${\it\Gamma}_{i,j}$ (recent simulations by Smiet et al. (Reference Smiet, Candelaresi, Thompson, Swearngin, Dalhuizen and Bouwmeester2015) give some support for this scenario). The tangential- $\boldsymbol{B}$ boundary condition, (1.1) is satisfied on both sides of these interfaces, but in general $\boldsymbol{B}$ suffers a tangential discontinuity across them. (To represent discontinuities across an interface ${\it\Gamma}_{i,j}$ , we distinguish its inward and outward faces ${\it\Gamma}_{i,j}^{\mp }$ by the superscripts $-$ and $+$ , respectively, following McGann et al. (Reference McGann, Hudson, Dewar and von Nessi2010).)

2. (2) An intermediate dynamical time scale (the time scale treated in this paper) during which the plasma, including a number of embedded current sheets separating subregions within which magnetic helicity, magnetic fluxes, mass and entropy are conserved, evolves adiabatically with respect to the relaxation time scale as the geometric shapes of the boundaries $\partial {\it\Omega}_{i}$ evolve dynamically from their initial conditions, and possibly in response to external forcing from the ‘switching on’ of boundary ripple, Hahm & Kulsrud (Reference Hahm and Kulsrud1985), Dewar et al. (Reference Dewar, Bhattacharjee, Kulsrud and Wright2013), Comisso et al. (Reference Comisso, Grasso and Waelbroeck2015a ,Reference Comisso, Grasso and Waelbroeck b ). Low-order resonant structures within the plasma that are excited, Boozer & Pomphrey (Reference Boozer and Pomphrey2010), White (Reference White2013) by geometric change resist the formation of magnetic islands by developing shielding current sheets.

3. (3) A long diffusion (possibly anomalous) time scale on which plasma and magnetic flux leak and mix between subregions through weak spots in the current sheets ${\it\Gamma}_{i,j}$ , violating the mass and flux isolation of the subregions assumed in MRxMHD and also violating entropy conservation. New subregions may form, changing the topological structure of the system. Phenomena on this time scale are not treated in this paper.

Figure 2.

Four sections of the regions shown in figure 1 at toroidal angles ${\it\zeta}=0$ , ${\rm\pi}/2$ , ${\rm\pi}$ and $3{\rm\pi}/2$ , showing helical rotation of the islands ${\it\Omega}_{1a}$ and ${\it\Omega}_{1b}$ about the central region ${\it\Omega}_{0}$ . The outer region, ${\it\Omega}_{2}$ , is an annular genus-2 toroid bounded by ${\it\Gamma}_{1a,2}^{+}$ , ${\it\Gamma}_{1b,2}^{+}$ and ${\it\Gamma}_{2,v}^{-}$ , the plasma-vacuum interface (green torus).

To illustrate these concepts, a multi-region case of interest is shown schematically in figures 1 and 2: an island chain $\{{\it\Omega}_{1a},{\it\Omega}_{1b}\}$ twisting, with linking number 2, Berger (Reference Berger1999), around a simple toroid ${\it\Omega}_{0}$ , both contained within an annular toroid ${\it\Omega}_{2}$ . This illustrates the flexibility of the MRxMHD partition; it is not limited to simply nested tori as it can include island or plasmoid chains. (This requires a more complicated labelling system for the relaxation regions and current sheets, but if we limit ourselves to primary island chains, we can still use a sequential numbering system and indicate the component islands and separatrices of the chain using the lettering scheme illustrated here.)

In single-region RxMHD, the topology of ${\it\Omega}$ determines the number, ${\it\nu}$ Footnote ⁵ , of independent fluxes ${\it\Phi}^{l}$ , $1\leqslant l\leqslant {\it\nu}$ , that must be specified for uniqueness of a harmonic (vacuum) field $\boldsymbol{B}_{H}$ solution such that $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}_{H}=0$ in ${\it\Omega}$ , $\boldsymbol{n}\boldsymbol{\cdot }\boldsymbol{B}_{H}=0$ on $\partial {\it\Omega}$ . For ${\it\mu}$ not an eigenvalueFootnote ⁶ of the Beltrami equation, (2.3), with homogenous boundary conditions (see e.g. § IV of Taylor (Reference Taylor1986)), specifying the ${\it\nu}$ fluxes ${\it\Phi}^{l}$ also specifies the Beltrami field uniquely, Yoshida & Giga (Reference Yoshida and Giga1990). Similarly, in MRxMHD we need to determine the genus ${\it\nu}_{i}$ of each relaxation subregion and specify its fluxes ${\it\Phi}_{i}^{l}$ , which are invariant under relaxation.

In the case of an annular toroid ${\it\Omega}_{i}$ (e.g. ${\it\Omega}_{2}$ in figures 1 and 2), $\partial {\it\Omega}_{i}$ consists of two disjoint tori (e.g. ${\it\Gamma}_{1a}^{+}\cup {\it\Gamma}_{1b}^{+}$ and ${\it\Gamma}_{2}^{-}$ ). A standard single-handled torus (i.e. with one hole) can be covered by a single coordinate chart, typically using a poloidal angle ${\it\theta}$ and a toroidal angle ${\it\zeta}$ . Assuming its boundaries to be two such standard tori, the genus of an annular toroid is ${\it\nu}=2$ , as both toroidal and poloidal cuts are required to make it simply connected.

Higher-genus cases might arise in toroidal confinement when treating doublet/multi-pinch, Taylor (Reference Taylor1986), or bundle divertor, Stott, Wilson & Gibson (Reference Stott, Wilson and Gibson1977), configurations. In the bundle divertor case, the plasma-vacuum boundary $\partial {\it\Omega}$ is a two handled torus, which cannot be described by a single toroidal–poloidal coordinate system but must instead be partitioned into two separate patches (coordinate charts). We shall not consider such exotic cases further in this paper.

The conservation of the ${\it\Phi}_{i}^{l}$ implies boundary constraints on the vector potential $\boldsymbol{A}$ : From Stokes’ theorem

where $\boldsymbol{n}$ is the unit normal at a point on the $l$ th topologically distinct surface of section ${\it\sigma}_{i}^{l}$ cutting ${\it\Omega}_{i}$ and ${\it\gamma}_{i}^{l}=\partial {\it\sigma}_{i}^{l}\in \partial {\it\Omega}_{i}$ is a loop around the boundary of ${\it\sigma}_{i}^{l}$ , with direction with respect to that of $\boldsymbol{n}$ given by the right-hand rule.

As the loops ${\it\gamma}_{i}^{l}$ lie on the boundary $\partial {\it\Omega}_{i}$ , which is composed of the current-sheet interfaces ${\it\Gamma}_{i,j}$ where there are ${\it\delta}$ -function currents causing tangential discontinuities in $\boldsymbol{B}$ , one might think the values of the loop integrals would depend on whether the loops traverse the inner or outer faces of the ${\it\Gamma}_{i,j}$ . However, it does not matter which faces are used as $\boldsymbol{B}$ remains finite within the current sheet; hence, being of infinitesimal width, a current sheet contains only infinitesimal flux. However, this continuity of its loop integrals does not necessarily mean $\boldsymbol{A}$ itself is continuous across current sheets, as a discontinuous gauge term $\boldsymbol{{\rm\nabla}}{\it\chi}$ does not affect the loop integrals (provided ${\it\chi}$ is single valued). This freedom allows a coordinate-dependent gauge to be used in each relaxation region, as in the SPEC code, Hudson et al. (Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012).

It is important to recognize that the fluxes ${\it\Phi}_{i}^{l}$ depend only on magnetic fields within ${\it\Omega}_{i}$ , so their individual conservation constrains line integrals of $\boldsymbol{A}$ only around loops that enclose the plasma within ${\it\Omega}_{i}$ . However, taking into account the conservation of all the ${\it\Phi}_{i}^{l}$ leaves only the toroidal line integral $\oint _{\partial {\it\Omega}}^{tor}\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ on the plasma–vacuum interface unconstrained. This represents the external poloidal flux threading the hole in the torus $\partial {\it\Omega}$ . While this flux is arbitrary as far as the physics of the plasma within ${\it\Omega}$ is concerned, it is still conserved because we are assuming the wall acts as a superconducting shell which traps the external poloidal flux threading it.

3.1. Lagrangian and Eulerian fluid kinematics

Central to the Lagrangian approach to fluid mechanics is the concept of fluid elements, whose motions with respect to time $t$ through a three-dimensional Cartesian frame (points in which we designate by the vector $\boldsymbol{x}\equiv x\boldsymbol{e}_{x}+y\boldsymbol{e}_{y}+z\boldsymbol{e}_{z}$ ), are described by a family of feasible trajectories (pathlines respecting the constraints) $\boldsymbol{x}^{t}=\boldsymbol{r}^{t}(\boldsymbol{x}_{0})$ , labelled by $\boldsymbol{x}_{0}$ , the initial positions of fluid elements at an arbitrary time $t=t_{0}$ . The fluid elements are advected by the Eulerian velocity vector field $\boldsymbol{v}(\boldsymbol{x},t)$ to their positions at arbitrary time $t$ through the time evolution function $\boldsymbol{r}^{t}(\boldsymbol{x})$ defined as the integral of the following equation and initial condition

for all $\boldsymbol{x}$ in the domain of interest. (Note that this makes $\boldsymbol{r}^{t}$ implicitly a function of $t_{0}$ , which we can make explicit when needed using the notation $\boldsymbol{r}^{t}(\boldsymbol{x}\mid t_{0})$ , constant parameters such as $t_{0}$ being listed after the vertical bar $|$ .)

Where an argument (other than $t$ ) is not specified, we treat $\boldsymbol{r}^{t}$ as a map $\mathbb{R}^{3}\rightarrow \mathbb{R}^{3}$ (i.e. a 3-vector function of 3-vectors), but often it is necessary to recognize that it is also a functional of $\boldsymbol{v}(\boldsymbol{x},t)$ , which will be indicated explicitly when required using the notation $\boldsymbol{r}^{t}[\boldsymbol{v}](\boldsymbol{x})$ . This provides a very flexible notation that can be adapted for generating other maps.

Suppose $\boldsymbol{v}$ is also a function of some time-independent parameter (say $s$ , denoted by $\boldsymbol{v}(\boldsymbol{x},t\mid s)$ ) then $\boldsymbol{r}^{t}$ will also be a function of $s$ , denoted $\boldsymbol{r}^{t}(\boldsymbol{x}\mid s,t_{0})$ . Suppose further that $\boldsymbol{x}$ lies on an arbitrary curve $\boldsymbol{x}=\boldsymbol{f}(s)$ at time $t_{0}$ and denote the resulting family of trajectories by $\boldsymbol{R}^{t}(s)\equiv \boldsymbol{r}^{t}(\boldsymbol{f}(s)\mid s,t_{0})$ . Differentiating both sides of the equation of motion in (3.1) with respect to $s$ , we have

where $\partial _{s}$ means the partial derivative with respect to $s$ and $\boldsymbol{{\rm\nabla}}\boldsymbol{v}$ denotes $\boldsymbol{{\rm\nabla}}_{\boldsymbol{x}}\boldsymbol{v}(\boldsymbol{x},t)$ .

We now use the special case of the above result where $\boldsymbol{v}$ is independent of $s$ to build up some useful differential-geometric evolution results. First, denoting an infinitesimal line element advected by the fluid by $\text{d}\boldsymbol{l}^{t}\equiv \text{d}\boldsymbol{R}^{t}$ , (3.2) immediately gives

When acting on all points in a region ${\it\Omega}$ , the evolution function $\boldsymbol{r}^{t}$ defines the Lagrangian map, mapping an initial region ${\it\Omega}_{0}$ onto its image ${\it\Omega}^{t}$ . From (3.3) applied to the sides of an infinitesimal rhomboid within ${\it\Omega}^{t}$ , the advection equation for infinitesimal volumes $\text{d}V^{t}$ is found to be

which is equivalent to the evolution equation, $(\text{d}/\text{d}t)J(t)=(\boldsymbol{{\rm\nabla}}\boldsymbol{\cdot }\boldsymbol{v})J(t)$ , for the Jacobian, $J(t)=\text{d}V^{t}/\text{d}V_{0}\equiv \partial (x^{t},y^{t},z^{t})/\partial (x_{0},y_{0},z_{0})$ , of the transformation from initial coordinates of fluid elements to the corresponding coordinates at time $t$ .

Likewise, the boundary $\partial {\it\Omega}_{0}$ maps onto $\partial {\it\Omega}^{t}$ , within which area elements $\text{d}\boldsymbol{S}^{t}\equiv \boldsymbol{n}^{t}\,\text{d}S^{t}$ advect according to

Representing a Lagrangian map as a dynamical flow induced by an Eulerian velocity field allows connection to be established with the modern Lie algebra approach to fluid dynamics, Arnold & Khesin (Reference Arnold and Khesin1998). However, although Lie operator methods are useful in Hamiltonian perturbation theory, Dewar (Reference Dewar1976), the simpler Lagrangian-based variational approach used in this paper avoids the need for most of this abstract machinery. (Likewise for abstract differential geometry.) Nevertheless, we shall find the extension of our Eulerian–Lagrangian mapping notation to include flows other than time evolution makes for a compact notation, and, being semi-Eulerian, leads to a more familiar form for the perturbation expansion of the Lagrangian than the strictly Lagrangian approach, Dewar (Reference Dewar1970).

Application of Hamilton’s Principle requires us to vary trial fluid element pathlines to find the Euler–Lagrange equations that determine which such pathlines are actually physical. Thus we introduce a new flow that maps the position vectors of fluid elements from their unvaried positions at each time $t$ to their varied positions at the same time by using a flow analogous to the Lagrangian map defined in (3.1), but with $t$ replaced by a dimensionless variation parameter ${\it\epsilon}$ (typically small), and with the variational map generated by a variational velocity ${\bf\nu}$ . Thus, in (3.1)–(3.5) replace $t$ with ${\it\epsilon}$ , the initial time $t_{0}$ with 0, the time evolution flow $\boldsymbol{r}^{t}[\boldsymbol{v}]$ with the ${\it\epsilon}$ -flow $\boldsymbol{r}^{{\it\epsilon}}[{\bf\nu}]$ and the velocity $\boldsymbol{v}(\boldsymbol{x},t)$ with ${\bf\nu}(\boldsymbol{x},{\it\epsilon}\mid t)$ , to give

with the parametric time-dependence $(\ldots \mid t)$ now left implicit.

The varied Lagrangian map $\boldsymbol{r}^{t}[\boldsymbol{v}_{{\sim}}]$ is now found by composing the unvaried Lagrangian map with the variational map,

where $\,\circ \,$ denotes composition of functions: $f\,\circ \,g(x)\equiv f(g(x))$ and subscript ${\sim}$ denotes a varied quantity. This implicitly defines the varied velocity field $\boldsymbol{v}_{{\sim}}$ , which is to be found in terms of the unvaried position of a representative fluid element $\boldsymbol{x}^{t}\equiv \boldsymbol{r}^{t}(\boldsymbol{x}_{0})$ as the total time derivative of the varied position, $\boldsymbol{x}_{{\sim}}^{t}\equiv \boldsymbol{r}^{t}[\boldsymbol{v}_{{\sim}}](\boldsymbol{x}_{0})\equiv \boldsymbol{r}^{{\it\epsilon}}(\boldsymbol{x}^{t}\mid t)$ , as

where $\boldsymbol{x}_{{\sim}}$ here denotes the Eulerian representation of the varied position, $\boldsymbol{x}_{{\sim}}(\boldsymbol{x},t)\equiv \boldsymbol{r}^{{\it\epsilon}}(\boldsymbol{x}\mid t)$ , and $\text{D}_{t}\equiv \partial _{t}+\boldsymbol{v}(\boldsymbol{x},t)\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ is the advective derivative.

In MRxMHD, we also need to consider the case of fluid elements on the common interfaces (current sheets) ${\it\Gamma}_{i,j}=\partial {\it\Omega}_{i}\cap \partial {\it\Omega}_{j}$ separating subregions ${\it\Omega}_{i}$ and ${\it\Omega}_{j}$ . The shape of the interface, which we represent as the level surface $f_{(i,j)}=0$ of an appropriate smooth function $f(\boldsymbol{x})$ changing monotonically across the surface, is not known a priori so must be subject to variation in applying Hamilton’s Principle. Thus we must introduce the variation parameter ${\it\epsilon}$ in representing ${\it\Gamma}_{i,j}$ geometrically: $f_{(i,j)}(\boldsymbol{x}\mid t,{\it\epsilon})=\pm 0$ . We use the notation $\boldsymbol{x}_{{\sim}\pm }^{t}$ to distinguish which side of the interface a varied fluid element is on: $f_{(i,j)}(\boldsymbol{x}_{{\sim}\pm }^{t}\mid t,{\it\epsilon})=\pm 0$ . Taking the total derivative of left and right-hand sides of this expression with respect to $t$ , on both sides of the interface, we find $\boldsymbol{n}\boldsymbol{\cdot }\unicode[STIX]{x27E6}\boldsymbol{v}_{{\sim}}\unicode[STIX]{x27E7}=0$ , $\boldsymbol{n}\equiv \boldsymbol{{\rm\nabla}}f_{(i,j)}/|\boldsymbol{{\rm\nabla}}f_{(i,j)}|$ being the unit normal and $\unicode[STIX]{x27E6}\cdot \unicode[STIX]{x27E7}$ denoting the jump in a quantity as the evaluation point crosses the interface (so $\unicode[STIX]{x27E6}\partial _{t}f_{(i,j)}(\boldsymbol{x}\mid t,{\it\epsilon})\unicode[STIX]{x27E7}=0$ , as $f_{(i,j)}$ is assumed smooth). Similarly, total differentiation with respect to ${\it\epsilon}$ gives

This states that the normal component of ${\bf\nu}$ is constrained to be continuous across the interface, but otherwise it is unconstrained. Rather, ${\it\Gamma}_{i,j}$ is advected with ${\bf\nu}$ during variations at constant $t$ , just as it is under time evolution.

An exception is the case where the plasma is confined by a prescribed, though possibly time dependent, boundary ‘wall’ (w). This is the special case that $f_{(i,j=w)}$ is not a function of ${\it\epsilon}$ , so one obtains the constraint $\boldsymbol{n}\boldsymbol{\cdot }{\bf\nu}=0$ at constant $t$ on ${\it\Gamma}_{(i,w)}^{-}$ , with the tangential components unconstrained.

So far we have treated ${\it\epsilon}$ as a finite parameter, on a par with $t$ . However, in this paper we use the variational transformation only for calculating the first variation of the action in Hamilton’s Principle, so we need the variational map only to linear order: $\boldsymbol{r}^{{\it\epsilon}}[{\bf\nu}](\boldsymbol{x})=\boldsymbol{x}+{\it\epsilon}{\bf\nu}(\boldsymbol{x},0\mid t)+O({\it\epsilon}^{2})$ . Thus, defining the Lagrangian variation in position ${\rm\Delta}\boldsymbol{x}$ through $\boldsymbol{x}_{{\sim}}=\boldsymbol{x}+{\it\epsilon}{\rm\Delta}\boldsymbol{x}(\boldsymbol{x},t)+O({\it\epsilon}^{2})$ , we have ${\rm\Delta}\boldsymbol{x}(\boldsymbol{x},t)={\bf\nu}(\boldsymbol{x},0\mid t)$ . Similarly, defining the Lagrangian variation in velocity, ${\rm\Delta}\boldsymbol{v}$ , to be such that $\boldsymbol{v}_{{\sim}}(\boldsymbol{x}_{{\sim}},t)=\boldsymbol{v}(\boldsymbol{x},t)+{\it\epsilon}{\rm\Delta}\boldsymbol{v}(\boldsymbol{x},t)+O({\it\epsilon}^{2})$ , we have, from (3.8),

In the above we have followed Newcomb (Reference Newcomb1962) in using ${\rm\Delta}f\equiv \lim _{{\it\epsilon}\rightarrow 0}[f_{{\sim}}(\boldsymbol{x}_{{\sim}},{\it\epsilon}\mid t)-f(\boldsymbol{x},{\it\epsilon}\mid t)]/{\it\epsilon}$ to denote the Lagrangian variation in an arbitrary field $f$ , while the corresponding Eulerian variation ${\it\delta}f$ is defined by ${\it\delta}f\equiv \lim _{{\it\epsilon}\rightarrow 0}[f_{{\sim}}(\boldsymbol{x},{\it\epsilon}\mid t)-f(\boldsymbol{x},{\it\epsilon}\mid t)]/{\it\epsilon}$ , i.e. with both the varied and unvaried field evaluated at the same, unvaried, position $\boldsymbol{x}$ . Thus ${\it\Delta}$ may be regarded as the operator $\lim _{{\it\epsilon}\rightarrow 0}\text{d}/\text{d}{\it\epsilon}$ while ${\it\delta}$ is the operator $\lim _{{\it\epsilon}\rightarrow 0}\partial /\partial {\it\epsilon}$ .

By definition, the two operators are related by ${\it\Delta}={\it\delta}+{\rm\Delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ . Applying both sides to $\boldsymbol{x}$ , it is easily verified as a consistency check that ${\it\delta}\boldsymbol{x}=0$ . Also, both ${\it\delta}$ and ${\it\Delta}$ being differential operators, the product rule, e.g. ${\rm\Delta}(fg)=({\rm\Delta}f)g+g{\rm\Delta}f$ , and commutation relation

apply, and correspondingly for ${\it\delta}$ , where $f(\boldsymbol{x},{\it\epsilon})$ and $g(\boldsymbol{x},{\it\epsilon})$ are arbitrary.

3.2. Holonomically constrained and free variations

In this subsection we extend the general formalism developed above for fluid kinematics to treat variation and perturbation of fields, specifically mass density ${\it\rho}(\boldsymbol{x},t)$ , pressure $p(\boldsymbol{x},t)$ and magnetic vector potential $\boldsymbol{A}(\boldsymbol{x},t)$ (and hence magnetic field $\boldsymbol{B}\equiv \boldsymbol{{\rm\nabla}}\times \boldsymbol{A}$ ).

The density (mass conservation) equation is the lowest and most robust one in the hierarchy of moment equations used in deriving fluid models from kinetic theory. In fact, freezing mass into fluid elements seems fundamental to any fluid theory, so we build the holonomic mass conservation constraints ${\it\rho}_{{\sim}}(\boldsymbol{x}_{{\sim}}^{t},t\mid {\it\epsilon})\,\text{d}V_{{\sim}}^{t}={\it\rho}(\boldsymbol{x}^{t},t)\,\text{d}V^{t}={\it\rho}(\boldsymbol{x}_{0},t_{0})\,\text{d}V_{0}$ into both IMHD (as in Newcomb (Reference Newcomb1962), Dewar (Reference Dewar1970)) and into our new MRxMHD formulation. From (3.4) and its ${\it\epsilon}$ -flow analogue, these imply

or, decomposing $\text{d}/\text{d}t$ as $\partial _{t}+\boldsymbol{v}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ , $\text{d}/\text{d}{\it\epsilon}$ as $\partial _{{\it\epsilon}}+{\bf\nu}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ ,

As explained above, the Lagrangian variation in density is ${\rm\Delta}{\it\rho}=\text{d}{\it\rho}_{{\sim}}/\text{d}{\it\epsilon}\mid _{{\it\epsilon}=0}$ , and the Eulerian variation is ${\it\delta}{\it\rho}=\partial {\it\rho}_{{\sim}}/\partial {\it\epsilon}\mid _{{\it\epsilon}=0}$ . That is, from (3.12) and (3.13),

In MRxMHD (unlike IMHD, Newcomb (Reference Newcomb1962), Dewar (Reference Dewar1970)) pressure and magnetic field $\boldsymbol{{\rm\nabla}}\times \boldsymbol{A}$ are not holonomically constrained by the freezing of entropy and flux microscopically into each fluid element but are freely variable, expressed by writing pressure and vector potential as $p(\boldsymbol{x},t|{\it\epsilon})$ and $\boldsymbol{A}(\boldsymbol{x},t|{\it\epsilon})$ , with their ${\it\epsilon}$ -derivatives ${\it\delta}p$ and ${\it\delta}\boldsymbol{A}$ being arbitrary variations in Hamilton’s Principle (as is ${\rm\Delta}\boldsymbol{x}$ ). (This is not of course to say Lagrangian variations of free, i.e. not holonomically constrained, fields do not exist, but rather that their Eulerian variations ${\it\delta}$ are primary, with their Lagrangian variations being defined by the rule ${\it\Delta}={\it\delta}+{\rm\Delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}$ above.) The relation of interface kinematics to the holonomic tangential- $\boldsymbol{B}$ constraint, (1.1), is developed in appendix B.

Consider an integral of the form $L=\int _{{\it\Omega}}{\mathcal{L}}\text{d}V$ , where ${\mathcal{L}}$ is a function of various fields to be varied in Hamilton’s Principle (in our case ${\it\rho}$ , $\boldsymbol{v}$ , $p$ and $\boldsymbol{A}$ ), which are thus functions of ${\it\epsilon}$ . By differentiating $L$ with respect to ${\it\epsilon}$ , using the ${\it\epsilon}$ -flow analogue of (3.4) (after changing variables to $x^{0},y^{0},z^{0}$ , where, as in Frieman & Rotenberg (Reference Frieman and Rotenberg1960), superscript $0$ refers to ${\it\epsilon}=0$ ) and integration by parts (Gauss’ theorem) we find the convenient identity

which allows us to use the Eulerian variations ${\it\delta}=\partial _{{\it\epsilon}}$ defined above to restrict the occurrence of ${\rm\Delta}\boldsymbol{x}$ to those fields ( ${\it\rho}$ and $\boldsymbol{v}$ ) that are constrained to vary with fluid motions. (An alternative, but equivalent, approach is to transform to $x_{0},y_{0},z_{0}$ at the outset, using Lagrangian variations and the ${\it\epsilon}$ -flow analogue of (3.5) to obtain surface terms.)

In applying (3.15) to compute the action variation ${\it\delta}{\mathcal{S}}=\int \text{d}t\int _{{\it\Omega}^{t}}\text{d}V{\mathcal{L}}$ , further surface terms will arise after straightforward integration by parts to remove spatial derivatives of ${\rm\Delta}\boldsymbol{x}$ . However, the integration by parts arising from terms containing time derivatives is more subtle because it involves a change of variables from $x,y,z$ to $x_{0},y_{0},z_{0}$ to enable integration by parts, and then a change back to $x,y,z$ . In our case, the only such term arises from the kinetic energy term, a function of $\boldsymbol{v}$ , the time derivative arising from the Eulerian variation ${\it\delta}\boldsymbol{v}={\rm\Delta}\boldsymbol{v}-{\rm\Delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v}=\text{D}_{t}{\rm\Delta}\boldsymbol{x}-{\rm\Delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{{\rm\nabla}}\boldsymbol{v}$ by (3.10).

Consider a Lagrangian density of the form ${\mathcal{L}}(\boldsymbol{v})$ . Then, noting that ${\rm\Delta}\boldsymbol{x}$ is always taken to vanish at the endpoints of the time integration in Hamilton’s Principle, and using (3.4) we find the contribution to the action from ${\it\delta}\boldsymbol{v}$ ,

which has no surface term.

3.3. Lagrangian formulation of MRxMHD

To derive the dynamics of multi-region relaxed plasmas we use the same Lagrangian density ${\mathcal{L}}^{\mathit{MHD}}$ as in (1.2), integrating over each subvolume ${\it\Omega}_{i}$ and augmenting with appropriate Lagrange multiplier terms to form the Lagrangian in ${\it\Omega}_{i}$ , $L_{i}$ . We then sum to form the total effective Lagrangian

where ${\mathcal{R}}$ denotes the set of plasma relaxation regions, with their Lagrangians $L_{i}$ being given by

where $L_{i}^{\mathit{MHD}}$ is given by the integral of ${\mathcal{L}}^{\mathit{MHD}}$ , defined in (1.2), over ${\it\Omega}_{i}$ , the $S_{i}$ are the entropy invariants ( $mS_{en}$ in the notation of appendix A) given by integrating the right-hand side of (A 4) over ${\it\Omega}_{i}$ ,

and the magnetic helicity invariants $K_{i}$ are as in (2.1), evaluated over ${\it\Omega}_{i}$ . (The constant ${\it\kappa}$ , required to make the argument of $\ln$ dimensionless, is arbitrary for our purposes but is identified physically in appendix A.)

The Lagrange multipliers ${\it\tau}_{i}$ and ${\it\mu}_{i}$ are constant during variation in Hamilton’s Principle (i.e. are independent of ${\it\epsilon}$ ) but may change with $t$ , as they must be chosen to make $S_{i}=S_{i0}$ and $K_{i}=K_{i0}$ during evolution under the Euler–Lagrange equations. The constant reference values $S_{i0}$ and $K_{i0}$ are the respective initial values at $t=t_{0}$ evaluated over ${\it\Omega}_{i0}$ , making $L_{i}=L_{i}^{\mathit{MHD}}$ when the ${\it\tau}_{i}$ and ${\it\mu}_{i}$ are adjusted to satisfy the conservation conditions above. Subtracting off the constant $K_{i0}$ also makes $L_{i}$ independent of the constant ${\it\kappa}$ (because of the holonomic conservation of $\int {\it\rho}\,\text{d}V$ ), so the curious fact that the physical value of ${\it\kappa}$ , (A 5), involves Planck’s constant is irrelevant.

If there is a vacuum region, ${\it\Omega}_{v}$ , between a plasma–vacuum interface and the wall, then this may be treated similarly, but with ${\it\rho}$ set to zero and with the entropy and helicity constraints deleted, leaving the Lagrangian density

After applying (3.15)–(3.16) and appropriate integration by parts, the variation of the action must be of the general form

where ${\mathcal{R}}^{+}$ denotes the set of plasma relaxation regions plus the vacuum region (though note that only the variational derivative ${\it\delta}{\mathcal{S}}_{v}/{\it\delta}\boldsymbol{A}$ is non-zero in ${\it\Omega}_{v}$ ). The interval over which the time integral is taken does not need to be specified as variations at the endpoints are taken to vanish in Hamilton’s Principle.

Hamilton’s Variational Principle is the statement that the Euler–Lagrange equations following from requiring ${\it\delta}{\mathcal{S}}=0$ for all variations of the independent fields determine these fields physically. In the above variational derivative notation, these are the Euler–Lagrange equations ${\it\delta}{\mathcal{S}}_{i}/{\it\delta}\boldsymbol{A}=0$ , ${\it\delta}{\mathcal{S}}_{i}/{\it\delta}p=0$ , ${\it\delta}{\mathcal{S}}_{i}/{\it\delta}\boldsymbol{x}=0$ within the volumes ${\it\Omega}_{i}$ . To find the surface Euler–Lagrange equations the constraint (3.9) needs to be considered, the consequences of which will be discussed in § 3.5.

3.4. Volume variations

Inserting (3.18) in (3.17), varying ${\mathcal{S}}=\int L\,\text{d}t$ , integrating by parts and comparing with (3.21), we identify the variational derivative with respect to the vector potential as

giving as the corresponding Euler–Lagrange equation, ${\it\delta}{\mathcal{S}}_{i}/{\it\delta}\boldsymbol{A}=0$ , the Beltrami equation (2.3). In the case of the vacuum Lagrangian $L_{v}$ , ${\it\mu}_{i}$ is set to zero, giving the statement that the vacuum field is harmonic, $\boldsymbol{{\rm\nabla}}\times \boldsymbol{B}=0$ . (While the Euler–Lagrange equation in the vacuum region is the same as the Beltrami equation with ${\it\mu}$ set to zero, note that we have actually deleted the constraint of constant magnetic helicity in the vacuum region, so the vacuum is not completely equivalent to a currentless plasma.)

The variational derivative with respect to the pressure is

the corresponding Euler–Lagrange equation being

with ${\it\tau}_{i}$ identified (see (A 2)) as the specific temperature $T_{i}/m$ in ${\it\Omega}_{i}$ , where $m$ is the effective ion mass $m_{i}/Z_{eff}$ ( ${\it\tau}_{i}$ is also $C_{i}^{2}$ , where $C_{i}$ is the ion sound speed).

Note that, despite using internal energy and entropy densities consistent with a microscopically isentropic equation of state, we have only enforced macroscopic entropy conservation over a whole subregion, leading to the microscopically isothermal equation of state (3.24) ( ${\it\tau}_{i}$ being spatially constant). Rapid equilibration of temperature within ${\it\Omega}_{i}$ is compatible with magnetic surfaces being destroyed within a relaxation region, hence poor local thermal confinement, but the bounding interfaces are assumed to be thermal transport barriers so the temperature can jump across each interface. Thus, while temperature profiles are restricted to being piece-wise constant, they are otherwise arbitrary, enabling the use of MRxMHD to model hot, magnetically confined plasmas (see e.g. Hudson et al. Reference Hudson, Dewar, Dennis, Hole, McGann, von Nessi and Lazerson2012).

As already remarked, the Lagrange multipliers ${\it\mu}_{i}$ and ${\it\tau}_{i}$ may change with time to maintain the constancy of their respective constraints. Now that the ${\it\tau}_{i}$ have been identified as temperatures, a simple thought experiment makes it physically clear that this must be so: suppose the ${\it\rho}_{i}$ , and thus the $p_{i}$ , are spatially constant within each ${\it\Omega}_{i}$ , then, just as for adiabatically deformed bags of ideal gas, ${\it\tau}_{i}V_{i}^{{\it\gamma}-1}=\text{const.}$ and any change in volume $V_{i}$ leads to a change in temperature ${\it\tau}_{i}$ .

Using (3.14) and (3.16) we then find the variational derivative with respect to fluid element positions

where the last line follows from the isothermal equation of state (3.24). The corresponding Euler–Lagrange equation ${\it\delta}{\mathcal{S}}_{i}/{\it\delta}\boldsymbol{x}=0$ is the equation of motion for a compressible Euler fluid in momentum conservation form. Note the extraordinary simplicity of this result, with $\boldsymbol{v}$ being decoupled from $\boldsymbol{B}$ as the Beltrami equation (2.3) implies the Lorentz force $\boldsymbol{j}\times \boldsymbol{B}$ is zero, showing that MRxMHD supports only steady flows and sound waves, with phase velocity $(T/m)^{1/2}$ , within the relaxation regions. (However surface waves on the interfaces can involve perturbations of $\boldsymbol{B}$ .)

3.5. Surface variations

It is shown in appendix B that tangential $\boldsymbol{B}$ at the interfaces ${\it\Gamma}_{i,j}$ implies a holonomic constraint on variations of the tangential component, $\boldsymbol{A}_{tgt}\equiv (\unicode[STIX]{x1D644}-\boldsymbol{n}\boldsymbol{n})\boldsymbol{\cdot }\boldsymbol{A}$ , of the vector potential,

for all $\boldsymbol{x}$ on the interface or boundary, where ${\it\delta}{\it\chi}$ is an arbitrary single-valued gauge potential. The normal component $\boldsymbol{n}\boldsymbol{\cdot }{\it\delta}\boldsymbol{A}$ is unconstrained. It is also shown in appendix B that this constraint implies invariance of line integrals of $\boldsymbol{A}$ around loops on these surfaces, and thus conservation of flux does not need to be imposed as an extra constraint.

Taking into account § 3.5, the surface variational derivative in (3.21) is found to be

which is the velocity-independent part of the stress tensor, Dewar (Reference Dewar1970) dotted with  $\boldsymbol{n}_{i}$ .

However, before we can apply this result, we need to take into account the fact that a boundary $\partial {\it\Omega}_{i}$ is made up of interfaces between ${\it\Omega}_{i}$ and neighbouring regions ${\it\Omega}_{j}$ (say) across which, by (3.9), the normal components of ${\rm\Delta}\boldsymbol{x}$ are continuous. Taking into account the outward normals of contiguous regions being oppositely directed, this implies the constraint ${\rm\Delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{n}_{i}=-{\rm\Delta}\boldsymbol{x}\boldsymbol{\cdot }\boldsymbol{n}_{j}$ on the common interfaces ${\it\Gamma}_{i,j}$ . This coupling gives the surface Euler–Lagrange equation, the natural boundary condition between interfaces, as the continuity condition

which is the same as the jump condition for advected discontinuities in ideal MHD (ee e.g. § 5.12 of Hosking & Dewar (Reference Hosking and Dewar2015)). (The tangential components ${\rm\Delta}\boldsymbol{x}_{tgt}$ are separately variable, but give no natural boundary conditions because they do not appear in ${\it\delta}{\mathcal{S}}$ .)