Three-dimensional magnetohydrodynamic equilibria with continuous magnetic fields

Abstract

A brief critique is presented of some different classes of magnetohydrodynamic equilibrium solutions based on their continuity properties and whether the magnetic field is integrable or not. A generalized energy functional is introduced that is comprised of alternating ideal regions, with nested flux surfaces with an irrational rotational transform, and Taylor-relaxed regions, possibly with magnetic islands and chaos. The equilibrium states have globally continuous magnetic fields, and may be constructed for arbitrary three-dimensional plasma boundaries and appropriately prescribed pressure and rotational-transform profiles.

1 Introduction

A fundamental requirement for magnetically confining plasmas for fusion research is to construct configurations for which the macroscopic forces acting on the plasma are balanced. The simplest, non-trivial equilibrium model considers only the pressure gradient and Lorentz forces, and force balance is described by

(1.1) $$\begin{eqnarray}\unicode[STIX]{x1D735}p=\boldsymbol{j}\times \boldsymbol{B},\end{eqnarray}$$

where $\unicode[STIX]{x1D735}p$ is the pressure gradient, $\boldsymbol{j}=\unicode[STIX]{x1D735}\times \boldsymbol{B}$ is the current density and $\boldsymbol{B}$ is the magnetic field. This equation is sometimes referred to as the ideal force balance equation, and it can be derived as the Euler–Lagrange equation for states that minimize the plasma energy functional under ideal variations (Bernstein et al. 1958; Kruskal & Kulsrud 1958; Hirshman, Sanchez & Cook 2011). The energy functional and its variations will be described below.

Despite the dramatic oversimplification of plasma dynamics, this equation is widely used to define the equilibrium. Indeed, it is used because of the simplicity: accurate numerical evaluations for simple plasma models are, understandably, faster than those of more complicated models. It therefore becomes practical to compute the hundreds of thousands of equilibria required for experimental design optimization, equilibrium reconstruction and so on, in strongly shaped, three-dimensional (3-D) geometries. Furthermore, if the macroscopic forces acting on the plasma are not at least approximately balanced, then there is little point in considering the microscopic forces.

Preferably, exact solutions should be elucidated that can be approximated with standard numerical discretizations consistent with the mathematical structure of the solutions, so the numerical error will reliably and predictably decrease with increasing numerical resolution. As with all differential equations, boundary conditions must be supplied to obtain a unique solution (for the sake of simplicity, this paper will ignore the possibility of bifurcations, for which two distinct solutions may be found for the same boundary conditions). In fact, the correct choice of boundary conditions is crucially important in guaranteeing the existence of well-defined solutions.

There are fundamental mathematical problems with (1.1) that are associated with its elliptic and hyperbolic characteristics (Grad 1960; Shiraishi, Ohsaki & Yoshida 2005), which this paper will not address. The mixed ideal–relaxed equilibrium model introduced below will, in the ‘ideal regions’, avoid the difficulties associated with the real characteristics by following Betancourt & Garabedian (1985) in assuming the existence of nested toroidal flux surfaces, which allows the equation $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p=0$ to be immediately solved by $p=p(\unicode[STIX]{x1D713})$ , where $\unicode[STIX]{x1D713}$ labels the enclosed toroidal flux. In the ‘relaxed’ regions, attention will be restricted to a subset of solutions of (1.1), namely linear force-free fields that satisfy $\unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}\boldsymbol{B}$ for constant $\unicode[STIX]{x1D707}$ , and the assumption of nested surfaces is not required.

This paper shall restrict attention to the so-called fixed-boundary case, for which the plasma boundary is prescribed, herein assumed to be smooth, and for which $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ , where $\boldsymbol{n}$ is normal. It is, however, simple to generalize the following to the free-boundary case, for which a supporting ‘vacuum’ field generated by currents external to the plasma must be provided.

This paper shall also adopt what may be called the ‘equilibrium’ approach: the pressure, $p(\unicode[STIX]{x1D713})$ , is to be provided and is required to not change during the calculation. Depending on the particular class of equilibrium to be constructed, at least one other profile function must usually be provided, such as the parallel current density, $\unicode[STIX]{x1D707}(\unicode[STIX]{x1D713})$ , or the rotational transform, $\text{-}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})$ . The equilibrium calculation is then to determine the magnetic field that satisfies force balance and is consistent with the given plasma boundary and the given profiles. Note that, typically, if the parallel current density is specified a priori, then the rotational transform is only known a posteriori, and vice versa.

The equilibrium approach is in contrast to, for example, what may be called the ‘transport’ approach, whereby an initial pressure and magnetic field both evolve dynamically in time (or iteratively) according to, for example, the resistive, extended magnetohydrodynamic (MHD) equations (Sovinec et al. 2003; Jardin, Breslau & Ferraro 2007) towards what might be called a resistive, or ‘Ohmic’, steady state (Park et al. 1986; Suzuki et al. 2006; Schlutt et al. 2012, 2013). For example, the pressure might be allowed to evolve according to an anisotropic diffusion law, which is effectively a transport equation. The transport approach certainly has merit and can include additional, non-ideal physics; however, it does not easily lend itself to constructing an equilibrium state with a given pressure. (See also the simulated annealing method advanced by Furukawa & Morrison (2016), which advances an initial state according to a modified set of equations derived from reduced MHD with constrained Casimirs.)

Equation (1.1) implies $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p=0$ , so the pressure is constant along each magnetic field line. This constraint has important consequences: the pressure, which is an ‘input’, is intimately related to the magnetic field, which is an ‘output’ of the numerical calculation. A necessary feature of equilibrium codes is to appropriately constrain the magnetic field to ensure that intact magnetic flux surfaces coincide with the prescribed pressure gradients: an equilibrium code that solves ideal force balance must constrain the topology of the field to be consistent with the given pressure.

1.1 Different classes of solution

By restricting attention to axisymmetric configurations with a rotational symmetry, $\unicode[STIX]{x1D735}p=\boldsymbol{j}\times \boldsymbol{B}$ reduces to the Grad–Shafranov equation (Grad & Rubin 1958; Shafranov 1966; Boozer 2005). The ignorable coordinate guarantees the existence of solutions with integrable magnetic fields. Here, the word ‘integrable’ is used in the dynamical systems context (Goldstein 1980) to refer to magnetic fields with a continuously nested family of ‘flux’ surfaces that remain invariant under the magnetic field-line flow. Arbitrary smooth functions for the pressure and current-density profiles, for example, may be admitted.

Hereafter, this paper will consider the ‘3-D’ case, for which the plasma boundary does not have a continuous symmetry or an ignorable coordinate, and for which the magnetic field may or may not be integrable, depending on whether $\unicode[STIX]{x1D6FF}$ -function current-densities (i.e. sheet-currents) are admitted or not. Identifying computationally tractable, physically acceptable solutions is much more complicated than in the two-dimensional case. Since the early days of research into magnetically confined plasma it was recognized that 3-D MHD equilibrium states may be ‘pathological’ (Grad 1967).

There are several problems that must be addressed, depending on the class of solution that one seeks.

Solutions can be categorized as having either continuous or discontinuous pressure and magnetic field, $p$ and $\boldsymbol{B}$ , with either continuous or discontinuous derivatives, $\unicode[STIX]{x1D735}p$ and $\boldsymbol{j}$ , and with either integrable or non-integrable magnetic fields. Hereafter, we will use the word ‘smooth’ to describe a function that is continuous and has continuous first derivatives.

Identification of the continuity properties of the solution is crucial as this determines which numerical discretizations may be employed. The continuity properties of the solution to a differential equation are partly determined by the continuity properties of the supplied boundary conditions. To obtain continuous solutions, the pressure and rotational transform must also be continuous, but this is not sufficient: it is also required to ensure that any singularities that may be present in the differential equation are avoided.

1.1.1 Smooth pressure, smooth non-integrable field

It seems reasonable to seek 3-D solutions with a smooth pressure and a smooth magnetic field. Being analogous to $1(1/2)$ dimensional Hamiltonian systems (Cary & Littlejohn 1983), continuous smooth 3-D magnetic field-line flows with shear are typically non-integrable (Lichtenberg & Lieberman 1992; Meiss 1992), possessing a fractal mix of (i) invariant surfaces known as Kolmogorov, Arnold Moser (KAM) surfaces (Moser 1973; Arnold 1978), which have ‘sufficiently irrational’ rotational transform, (ii) magnetic islands, which appear where the rotational transform is rational and (iii) chaotic ‘irregular’ field lines, which are associated with the unstable manifolds of the periodic field lines and ergodically fill a highly non-trivial volume. (Note that a magnetic vector field may be a smooth function of position, $\boldsymbol{B}(\boldsymbol{x}+\unicode[STIX]{x1D6FF}\boldsymbol{x})\approx \boldsymbol{B}(\boldsymbol{x})+\unicode[STIX]{x1D735}\boldsymbol{B}(\boldsymbol{x})\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{x}$ , but the magnetic field lines may be chaotic/irregular.) From $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p=0$ , it follows that any non-trivial, continuous pressure consistent with such a field must also be fractal, with $\unicode[STIX]{x1D735}p=0$ across the chaotic volumes and with non-zero, finite pressure gradients at a non-zero measure of KAM surfaces. The KAM surfaces nowhere densely fill a finite volume, and thus an uncountable infinity of discontinuities in the pressure gradient must arise (Kraus & Hudson 2017). Solutions with an infinity of discontinuities are intractable from a numerical perspective. Discontinuities in the pressure gradient drive discontinuities in the perpendicular current density, $\boldsymbol{j}_{\bot }=\boldsymbol{B}\times \unicode[STIX]{x1D735}p/B^{2}$ , and the magnetic field is not smooth.

Given an arbitrary, non-integrable magnetic field, it is a highly non-trivial problem to determine the fractal topological structure of the magnetic field lines. Which irrational surfaces survive 3-D perturbations depends in part on how ‘irrational’ the rotational transform is and how the system is perturbed from integrability. Individual KAM surfaces can be identified (with significant computational cost) using Greene’s residue criterion (Greene 1979); however, no one has yet, to the authors’ knowledge, described how to determine the measure of phase-space that is occupied with KAM surfaces for a given, non-integrable field.

It is the inverse of this task that is required for the equilibrium approach: one must first provide a continuous pressure profile with a fractally discontinuous gradient, and then appropriately constrain the representation of the non-integrable magnetic field to be topologically consistent with this given profile, i.e. to ensure that the flux surfaces coincide with the pressure gradients.

It is quite difficult to work with explicitly fractal functions. For example, consider the pressure-gradient profile defined by the Diophantine condition, which plays a prominent role in KAM theory and thus also in determining the structure of non-integrable magnetic fields,

(1.2) $$\begin{eqnarray}p^{\prime }(x)=\left\{\begin{array}{@{}ll@{}}-1,\quad & \text{if }|x-n/m|>d/m^{k},\forall n,m,\\ 0,\quad & \text{otherwise},\end{array}\right.\end{eqnarray}$$

where $d>0$ and $k\geqslant 2$ . The pressure gradient is zero in a non-zero neighbourhood of all rationals, $x=n/m$ . This function is not Riemannian-integrable. A standard discretization to compute the pressure on axis, with $p(1)=0$ , given by $p(0)=-\sum _{i=1}^{N}p^{\prime }(x_{i})\unicode[STIX]{x0394}x$ , where $x_{i}=i/N$ and $\unicode[STIX]{x0394}x=1/N$ , fails spectacularly, as do higher-order quadratures that are based on regular grids (Kraus & Hudson 2017).

To approximate such ‘fractal’ equilibria with non-integrable magnetic fields, a more reliable approach is to first provide well-defined, non-fractal pressure and rotational-transform profiles, which in turn provide a well-defined, non-fractal equilibrium that can be approximated with standard numerical discretizations to arbitrary accuracy, and then to consider the limiting properties of a sequence of such equilibria as the pressure and rotational-transform profiles approach fractals. We shall return to this idea later.

1.1.2 Smooth pressure, smooth integrable field

Instead of admitting equilibria with non-integrable fields, an alternative is to seek solutions with a smooth pressure and smooth, integrable magnetic field (Bauer, Betancourt & Garabedian 1982; Hirshman & Whitson 1983; Hirshman, van Rij & Merkel 1986; Taylor 1994). Such fields, having continuously nested flux surfaces, presumably are consistent with smooth pressure and transform profiles; however, unphysical currents arise near the rational rotational-transform surfaces.

The perpendicular current density consistent with (1.1) is $\boldsymbol{j}_{\bot }=\boldsymbol{B}\times \unicode[STIX]{x1D735}p/B^{2}$ . By enforcing $\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}=0$ , with $\boldsymbol{j}=j_{\Vert }\boldsymbol{B}+\boldsymbol{j}_{\bot }$ , a magnetic differential equation then determines the parallel current, $\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}j_{\Vert }=-\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}_{\bot }$ . Magnetic differential equations are densely singular, and thus are intractable numerically. For integrable fields, straight field-line coordinates, $\boldsymbol{x}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ , can be constructed and the magnetic field can be written $\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\text{-}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ . The Fourier harmonics of $j_{\Vert }$ must satisfy (Bhattacharjee et al. 1995)

(1.3) $$\begin{eqnarray}j_{\Vert ,m,n}=\frac{i(\sqrt{g}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}_{\bot })_{m,n}}{x}+\unicode[STIX]{x1D6E5}_{m,n}\unicode[STIX]{x1D6FF}(x),\end{eqnarray}$$

where $\unicode[STIX]{x1D6E5}_{m,n}$ is an as-yet undetermined constant and $x(\unicode[STIX]{x1D713})\equiv m\text{-}\unicode[STIX]{x1D704}(\unicode[STIX]{x1D713})-n$ . The Jacobian satisfies $1/\sqrt{g}=\boldsymbol{B}\boldsymbol{\cdot }\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ .

The $\unicode[STIX]{x1D6FF}$ -function current density is just a mathematical approximation of localized currents, and is acceptable in a macroscopic, perfectly conducting ideal-MHD model. (For example, the current-density associated with a finite current passing along a very thin strand of super-conducting wire is extremely well approximated by a $\unicode[STIX]{x1D6FF}$ -function.) Including $\unicode[STIX]{x1D6FF}$ -functions in the current density will result in a non-smooth magnetic field.

The $1/x$ singularity is far more problematic. For a special choice of straight field-line angles, namely Boozer coordinates (Boozer 1982; D’haeseleer et al. 1991), the magnetic field may be written $\boldsymbol{B}=\unicode[STIX]{x1D6FD}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}+I(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+G(\unicode[STIX]{x1D713})\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}$ , so that $1/B^{2}=\sqrt{g}/(G+\text{-}\unicode[STIX]{x1D704}I)$ , and

(1.4) $$\begin{eqnarray}(\sqrt{g}\unicode[STIX]{x1D735}\boldsymbol{\cdot }\boldsymbol{j}_{\bot })_{m,n}=\frac{p^{\prime }\sqrt{g}_{m,n}(nI-mG)}{G+\text{-}\unicode[STIX]{x1D704}I}.\end{eqnarray}$$

The magnitude of $\sqrt{g}_{m,n}$ may be considered to be an ‘output’ quantity: it is determined by the geometry of, and the tangential magnetic field on, the rational surfaces, both of which are determined by the magnetic field. For an arbitrary boundary, there is no apparent a priori control over the geometry of the internal flux surfaces.

Assuming the pressure satisfies $p(x)\approx p+p^{\prime }x+p^{\prime \prime }x^{2}/2+\ldots$ , the current through a cross-sectional surface bounded by $x=\unicode[STIX]{x1D716}$ and $x=\unicode[STIX]{x1D6FF}$ , and $\unicode[STIX]{x1D703}=0$ and $\unicode[STIX]{x1D703}=\unicode[STIX]{x03C0}/m$ , associated with the resonant harmonic of the parallel current density described by (1.3) is

(1.5) $$\begin{eqnarray}-\!\frac{2}{m}\frac{i(nI-mG)}{(G+\text{-}\unicode[STIX]{x1D704}I)}\frac{p^{\prime }\sqrt{g}_{m,n}}{\text{-}\unicode[STIX]{x1D704}^{\prime }}(\ln \unicode[STIX]{x1D6FF}-\ln \unicode[STIX]{x1D716}),\end{eqnarray}$$

where all terms are evaluated at the rational surface. This approaches infinity as $\unicode[STIX]{x1D716}$ approaches zero.

This shows that there are cross-sectional surfaces close to every rational surface through which the total current is infinite, and this is unphysical. To guarantee such problems are avoided, and assuming that there are no restrictions on $\sqrt{g}_{m,n}$ , the pressure gradient must be zero on each rational surface. The next order term for the current through the cross-sectional surface is proportional to $p^{\prime \prime }(\unicode[STIX]{x1D6FF}-\unicode[STIX]{x1D716})$ , and so we must require that $p^{\prime \prime }<\infty$ . For any system with shear, the rational surfaces densely fill space, and so either the pressure profile is trivial, with $p^{\prime }=0$ everywhere, or the pressure gradient must be discontinuous.

There is another possibility: rather than flattening the pressure to avoid the logarithmic infinities in the parallel current, one may restrict attention to so-called ‘healed’ configurations, for which the resonant harmonic of the Jacobian, $\sqrt{g}_{m,n}$ , vanishes at each resonant surface (Weitzner 2014; Zakharov 2015; Weitzner 2016). Such a condition could only be satisfied for a restricted class of 3-D plasma boundaries.

There is another problem with ideal-MHD equilibria with integrable magnetic fields and rational surfaces, which is frequently overlooked: the solutions are not analytic functions of the boundary. The equation describing the first-order plasma displacement, under the constraints of ideal-MHD, induced by a small deformation to the boundary is ${\mathcal{L}}_{0}[\unicode[STIX]{x1D743}]\equiv \unicode[STIX]{x1D6FF}\boldsymbol{j}[\unicode[STIX]{x1D743}]\times \boldsymbol{B}+\boldsymbol{j}\times \unicode[STIX]{x1D6FF}\boldsymbol{B}[\unicode[STIX]{x1D743}]-\unicode[STIX]{x1D735}\unicode[STIX]{x1D6FF}p[\unicode[STIX]{x1D743}]=0$ . (Expressions relating the perturbed field, $\unicode[STIX]{x1D6FF}\boldsymbol{B}$ , and pressure, $\unicode[STIX]{x1D6FF}p$ , to ideal plasma displacements, $\unicode[STIX]{x1D743}$ , are given below.) As discussed by Rosenbluth, Dagazian & Rutherford (1973), this is a singular equation, and the perturbed surfaces overlap and perturbation theory breaks down. The problem of non-analyticity led Rosenbluth et al. (1973) to consider a nonlinear treatment of 3-D ‘kink’ states, and this analysis has recently been revisited in the context of understanding the effect of resonant magnetic perturbations (RMPs) in tokamak plasmas (Loizu & Helander 2017).

1.1.3 Discontinuous pressure, discontinuous non-integrable magnetic field

Discontinuous and non-smooth solutions to differential equations are not a problem per se. Well-defined equilibrium solutions with a finite number of discontinuities have been introduced. In 1996, stepped-pressure equilibrium states were introduced by Bruno & Laurence (1996), and theorems were provided that guarantee the existence of such equilibria, provided the 3D deviation from axisymmetry was sufficiently small. These configurations were recognized as extrema of the multi-region, relaxed MHD (MRxMHD) energy functional that was later introduced by Dewar and co-workers (Hole, Hudson & Dewar 2006; Hudson, Hole & Dewar 2007; Dewar et al. 2008; Hudson et al. 2012).

Stepped-pressure equilibria can be thought of as being comprised of a finite number of nested Taylor states (Taylor 1974, 1986), in each of which the pressure is flat and the field satisfies a Beltrami equation, $\unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}\boldsymbol{B}$ with constant $\unicode[STIX]{x1D707}$ . The constraints of ideal-MHD are not continuously enforced, and this eliminates the problem of non-analyticity at the rational surfaces. The magnetic field may reconnect, i.e. the topology is not constrained, and magnetic islands will generally open at resonances. Where islands overlap, field-line chaos can emerge. For such ‘irregular’ field lines, the rotational transform is not well defined.

The discontinuities in the pressure coincide with a finite set of ‘ideal interfaces’, ${\mathcal{I}}_{i}$ , with strongly irrational rotational transform, that separate adjacent Taylor states. Strongly irrational numbers may, for example (Hudson et al. 2012), be simply expressed as $\text{-}\unicode[STIX]{x1D704}=(p_{1}+\unicode[STIX]{x1D6FE}\,p_{2})/(q_{1}+\unicode[STIX]{x1D6FE}\,q_{2})$ , where $\unicode[STIX]{x1D6FE}=(1+\sqrt{5})/2$ is the golden mean and $p_{1}/q_{1}$ and $p_{2}/q_{2}$ are neighbouring rationals (Meiss 1992). On these interfaces, the magnetic field is constrained to remain tangential, and the discontinuities in the pressure are balanced by discontinuities in the field strength, so that the ‘total pressure,’ $P\equiv p+B^{2}/2$ , is continuous across the ${\mathcal{I}}_{i}$ . The existence of tangential discontinuities in $\boldsymbol{B}$ implies the existence of sheet-currents. Stepped-pressure states, or MRxMHD states as they are also called, are almost everywhere relaxed but include a discrete set of (zero-volume) ideal interfaces. Example pressure and rotational-transform profiles are shown in figure 1.

Figure 1. Piecewise-constant, discontinuous pressure profile (left), and the discretely defined, strongly irrational rotational-transform profile, with some low-order island chains for illustration.

1.1.4 Continuous pressure, discontinuous integrable magnetic field

Another class of solutions with discontinuous magnetic fields, which are globally ideal, was introduced recently by Loizu et al. (2015b ), namely stepped-transform equilibria: equilibria with continuously nested flux surfaces with a discontinuous rotational transform. These were introduced after investigations (Loizu et al. 2015a ) into the $1/x$ and $\unicode[STIX]{x1D6FF}$ -function current-densities in ideal-MHD equilibria with integrable fields revealed the necessity to enforce infinite shear, $\text{-}\unicode[STIX]{x1D704}^{\prime }=\infty$ , at the rational surfaces in order to obtain consistent solutions. Effectively, the rational surface is removed from the equilibrium, and the non-integrable current-densities are avoided. Stepped-transform states can self-consistently support globally smooth, arbitrary pressure profiles. Removing the rational surfaces also removes the problem of non-analyticity, provided the discontinuities in the rotational transform across the rationals exceeds a minimum value – the sine qua non condition (Loizu et al. 2015b ) – for which analytic estimates were provided. The discontinuities in the rotational transform imply discontinuities in the tangential magnetic field, and so sheet-currents must also exist in these solutions.

The original investigation (Loizu et al. 2015b ) of these stepped-transform states was restricted to cylindrical geometry, with only one resonant deformation, and so only one rational surface was of concern, and so only one discontinuity in the rotational-transform profile was required to eliminate the pathologies. In the general case with an arbitrary 3-D boundary, every rational surface would generally result in unphysical currents. It is easy to generalize the concept to define equilibria with piecewise-constant rotational transform, for which the rotational transform is everywhere strongly irrational, and for which there is a finite collection of discontinuities/sheet-currents. Example profiles are shown in figure 2.

Both the stepped-pressure and the stepped-transform classes of equilibria possess sheet-currents and discontinuous magnetic fields. This is acceptable within a macroscopic, ideal-MHD context, as well as from a mathematical perspective, as a finite set of discontinuities is easy to accommodate numerically. The discontinuities in the magnetic field may create difficulties for subsequent calculations, for example, gyrokinetic calculations of transport.

In this paper, a new class of well-defined, numerically tractable, non-fractal equilibria is introduced that allows for non-integrable magnetic fields that are continuous, i.e. for which there are no sheet-currents. These states are a combination of the piecewise-constant rotational-transform equilibria with nested flux surfaces and smooth pressure profiles, and the piecewise-constant pressure equilibria with, generally, magnetic islands and chaotic field lines.

Figure 2. Arbitrary, smooth pressure profile (left), and the piecewise-constant, strongly irrational, discontinuous rotational-transform profile (below). No island chains are admitted.

2 Combined ideal–relaxed energy functional

The new equilibrium states are comprised of alternating ideal and relaxed regions and are extrema of the mixed ideal–relaxed energy functional, which will now be described. Restricting attention to toroidal configurations, the plasma volume is partitioned into $N$ subregions, ${\mathcal{R}}_{i}$ , $i=1,\ldots ,N$ , and we denote the toroidal boundaries separating the subregions by ${\mathcal{I}}_{i}$ . The magnetic axis (or axes) lies in ${\mathcal{R}}_{1}$ , which is a toroid and is bounded by ${\mathcal{I}}_{1}$ . For $i=2,\ldots ,N$ the ${\mathcal{R}}_{i}$ are annular, and $\unicode[STIX]{x2202}{\mathcal{R}}_{i}={\mathcal{I}}_{i-1}\cup {\mathcal{I}}_{i}$ . The outermost boundary, ${\mathcal{I}}_{N}$ , is coincident with the plasma boundary. On each of the ${\mathcal{I}}_{i}$ the magnetic field is constrained to be tangential, $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ . In each ${\mathcal{R}}_{i}$ , the plasma energy (Kruskal & Kulsrud 1958) is

(2.1) $$\begin{eqnarray}W_{i}\equiv \int _{{\mathcal{R}}_{i}}\left(\frac{p}{\unicode[STIX]{x1D6FE}-1}+\frac{B^{2}}{2}\right)\,\text{d}v.\end{eqnarray}$$

The equilibrium states minimize $W_{i}$ in each volume with respect to variations in the pressure and the magnetic field, but with suitable constraints imposed so as to avoid trivial solutions, and with respect to deformations in the internal boundaries, i.e. the ${\mathcal{I}}_{i}$ for $i=1,N-1$ .

In the ideal regions we restrict attention to integrable magnetic fields, with nested flux surfaces, which may be labelled by the enclosed toroidal flux. The equation of state, $d_{t}(p/\unicode[STIX]{x1D70C}^{\unicode[STIX]{x1D6FE}})=0$ , where $d_{t}\equiv \unicode[STIX]{x2202}_{t}+\boldsymbol{v}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$ and $\boldsymbol{v}$ is the ‘velocity’ of an assumed plasma displacement, $\boldsymbol{v}=\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D743}$ , may be combined with mass conservation, $\unicode[STIX]{x2202}_{t}\unicode[STIX]{x1D70C}+\unicode[STIX]{x1D735}\boldsymbol{\cdot }(\unicode[STIX]{x1D70C}\boldsymbol{v})=0$ , to obtain an equation that relates the ideal variation in the pressure to the plasma displacement, $\unicode[STIX]{x1D6FF}p=(\unicode[STIX]{x1D6FE}-1)\unicode[STIX]{x1D743}\boldsymbol{\cdot }\unicode[STIX]{x1D735}p-\unicode[STIX]{x1D6FE}\unicode[STIX]{x1D735}\boldsymbol{\cdot }(p\,\unicode[STIX]{x1D743})$ . Variations in the magnetic field are related to $\unicode[STIX]{x1D743}$ by Faraday’s law, $\unicode[STIX]{x2202}_{t}\boldsymbol{B}=\unicode[STIX]{x1D735}\times \boldsymbol{E}$ , and the ideal Ohm’s law, $\boldsymbol{E}+\boldsymbol{v}\times \boldsymbol{B}=0$ , where $\boldsymbol{E}$ is the electric field, and we write $\unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D735}\times (\unicode[STIX]{x1D743}\times \boldsymbol{B})$ . Note that this last constraint does not allow the topology of the field to change. The first variation of $W_{i}$ is

(2.2) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}W_{i}=\int _{{\mathcal{R}}_{i}}(\unicode[STIX]{x1D735}p-\boldsymbol{j}\times \boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D743}\,\text{d}v-\int _{\unicode[STIX]{x2202}{\mathcal{R}}_{i}}(p+B^{2}/2)\,\unicode[STIX]{x1D743}\boldsymbol{\cdot }\,\text{d}\boldsymbol{s}.\end{eqnarray}$$

In the Taylor-relaxed regions, the variations in the field and pressure are not related to (internal) plasma displacements. The mass and entropy constraints do not apply to individual fluid elements but instead to the entire volume, and the constraint on the pressure is $p_{i}V_{i}^{\unicode[STIX]{x1D6FE}}=a_{i}$ , where $V_{i}$ is the volume of ${\mathcal{R}}_{i}$ and $a_{i}$ is a constant. The internal energy in ${\mathcal{R}}_{i}$ is $\int _{{\mathcal{R}}_{i}}p_{i}/(\unicode[STIX]{x1D6FE}-1)\,\text{d}v$ $=$ $a_{i}V_{i}^{(1-\unicode[STIX]{x1D6FE})}/(\unicode[STIX]{x1D6FE}-1)$ , and the first variation of this due to a deformation, $\unicode[STIX]{x1D743}$ , of the boundary is $-p\int _{\unicode[STIX]{x2202}{\mathcal{R}}_{i}}\unicode[STIX]{x1D743}\boldsymbol{\cdot }\text{d}\boldsymbol{s}$ . The variation of the magnetic field is arbitrary, $\unicode[STIX]{x1D6FF}\boldsymbol{B}=\unicode[STIX]{x1D735}\times \unicode[STIX]{x1D6FF}\boldsymbol{A}$ , except for (i) constraints on the enclosed toroidal and poloidal fluxes, $\unicode[STIX]{x1D6F9}_{t,i}\equiv \int _{{\mathcal{P}}}\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ and $\unicode[STIX]{x1D6F9}_{p,i}\equiv \int _{{\mathcal{T}}}\boldsymbol{A}\boldsymbol{\cdot }\text{d}\boldsymbol{l}$ , where ${\mathcal{P}}$ and ${\mathcal{T}}$ are suitable poloidal and toroidal loops, (ii) conservation of the global helicity in each relaxed region,

(2.3) $$\begin{eqnarray}H_{i}\equiv \int _{{\mathcal{R}}_{i}}\boldsymbol{A}\boldsymbol{\cdot }\boldsymbol{B}\,\text{d}v,\end{eqnarray}$$

and (iii) the constraint that $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ on $\unicode[STIX]{x2202}{\mathcal{R}}_{i}$ . Much can be said about the helicity constraint (Woltjer 1958; Taylor 1974, 1986; Berger 1999), and we refer the interested reader to the recent paper by Moffat (2015).

The flux constraints can be enforced by constraining the representation for the vector potential, and the helicity constraint can be enforced by introducing a Lagrange multiplier, $\unicode[STIX]{x1D707}$ . The constrained energy functional in the relaxed regions is

(2.4) $$\begin{eqnarray}{\mathcal{F}}_{i}\equiv W_{i}-\frac{\unicode[STIX]{x1D707}}{2}(H_{i}-H_{i,0}).\end{eqnarray}$$

Note that if ${\mathcal{R}}_{i}$ is the innermost, toroidal region, the poloidal flux is not defined and only the constraints on the helicity and toroidal flux are required.

The first variation is

(2.5) $$\begin{eqnarray}\unicode[STIX]{x1D6FF}F_{i}=\int _{{\mathcal{R}}_{i}}(\unicode[STIX]{x1D735}\times \boldsymbol{B}-\unicode[STIX]{x1D707}\boldsymbol{B})\boldsymbol{\cdot }\unicode[STIX]{x1D6FF}\boldsymbol{A}\,\text{d}v-\int _{\unicode[STIX]{x2202}{\mathcal{R}}_{i}}(p+B^{2}/2)\,\unicode[STIX]{x1D743}\boldsymbol{\cdot }\,\text{d}\boldsymbol{s},\end{eqnarray}$$

where $\boldsymbol{A}=\unicode[STIX]{x1D743}\times \boldsymbol{B}$ has been used on the ${\mathcal{I}}_{i}$ .

The total constrained energy functional for the ideal–relaxed plasma is

(2.6) $$\begin{eqnarray}{\mathcal{F}}\equiv \mathop{\sum }_{i\in I}W_{i}+\mathop{\sum }_{j\in J}F_{j},\end{eqnarray}$$

where, for example, $I\equiv \{1,3,5,\ldots \}$ and $J\equiv \{2,4,6,\ldots \}$ , which makes the innermost volume an ideal region. Alternatively, a relaxed region may be assumed for the innermost volume, in which case $I\equiv \{2,4,6,\ldots \}$ and $J\equiv \{1,3,5,\ldots \}$ .

The Euler–Lagrange equations for extremizing states are as follows: in the ideal regions we have $\unicode[STIX]{x1D735}p=\boldsymbol{j}\times \boldsymbol{B}$ , in the relaxed regions we have $p=\text{const.}$ and $\unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}\boldsymbol{B}$ , and across the ${\mathcal{I}}_{i}$ we have $[[p+B^{2}/2]]=0$ . Note that fields that satisfy $\unicode[STIX]{x1D735}\times \boldsymbol{B}=\unicode[STIX]{x1D707}\boldsymbol{B}$ also satisfy $\unicode[STIX]{x1D735}p=\boldsymbol{j}\times \boldsymbol{B}$ , somewhat trivially, with $\unicode[STIX]{x1D735}p=0$ , so these mixed ideal–relaxed states globally satisfy $\unicode[STIX]{x1D735}p=\boldsymbol{j}\times \boldsymbol{B}$ .

Having presented a combined ideal–relaxed energy functional and derived the Euler–Lagrange equations governing extremal states, there are some subtleties concerning the prescribed pressure and rotational transform that must be addressed to eliminate the formation of sheet-currents. We seek solutions that have smooth pressure and continuous magnetic fields, so the pressure and the pressure gradient in each ideal region at each ${\mathcal{I}}_{i}$ must match that in the adjacent relaxed regions, where the pressure gradient is zero. To avoid the non-integrable current-densities described above, rational surfaces must be avoided in the ideal regions, so in the ideal regions we restrict attention to magnetic fields of the form $\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\text{-}\unicode[STIX]{x1D704}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ , where $\text{-}\unicode[STIX]{x1D704}_{i}$ is a strongly irrational constant.

Because of the possibility of reconnection and the formation of islands and irregular field lines, in the relaxed regions the rotational transform may not be globally defined. It is well defined on the ${\mathcal{I}}_{i}$ , which, because of the constraint $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ , remain as intact flux surfaces. However, if the Beltrami field is to be defined by prescribing the enclosed toroidal and poloidal fluxes and the helicity, the rotational transform on the ${\mathcal{I}}_{i}$ is a priori unknown, and must be computed a posteriori. We cannot a priori guarantee that an initial selection for $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}_{t,i}$ , $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}_{p,i}$ and $H_{i}$ is consistent with the existence of continuous rotational transform across the ${\mathcal{I}}_{i}$ . It will generally be required to iterate on the parallel current density – more formally, to iterate on $\unicode[STIX]{x0394}\unicode[STIX]{x1D713}_{p,i}$ and $H_{i,0}$ – in the relaxed regions to obtain the desired (single-valued) rotational-transform profile on the adjacent ${\mathcal{I}}_{i}$ .

We thus have described an equilibrium with a globally smooth pressure profile with ‘flattening’ across the rational surfaces, and with a piecewise-flat, piecewise-a priori-unknown rotational-transform profile. Smooth pressure gradients are supported in the ideal regions, which are filled with flux surfaces with a constant, strongly irrational rotational transform. Magnetic islands and chaotic field lines are allowed in the relaxed regions, in which the pressure gradient is zero, the rotational transform may or may not be defined, and $\boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{B}/B^{2}=\unicode[STIX]{x1D707}_{i}$ is a constant. Example pressure and transform profiles are shown in figures 3(a) and 3(b).

Figure 3. An example, in cylindrical geometry, of a mixed ideal–relaxed equilibrium with four regions. (a) The pressure profile. (b) The rotational-transform profile. (c) The ‘toroidal’ and poloidal components of the magnetic field, $B^{z}$ and $B^{\unicode[STIX]{x1D703}}$ . (d) The pressure gradient. (e) The poloidal component of the perpendicular current density, $j_{\bot }^{\unicode[STIX]{x1D703}}$ . (f) The parallel current density, $j_{\Vert }\equiv \boldsymbol{j}\boldsymbol{\cdot }\boldsymbol{B}/B^{2}$ .

Figure 3 also shows the magnetic field and current density for an example mixed ideal–relaxed equilibrium in a cylinder with ‘major radius’ $R=1$ and minor radius $a=0.1$ . The magnetic field is written $\boldsymbol{B}=rB^{\unicode[STIX]{x1D703}}(r)\hat{\unicode[STIX]{x1D703}}+B^{z}(r)\hat{z}$ . The current density is written $\boldsymbol{j}=\unicode[STIX]{x1D735}\times \boldsymbol{B}$ , which is manifestly divergence-free. In the ideal regions, the equilibrium equation reduces to

(2.7) $$\begin{eqnarray}\frac{\text{d}B^{z}}{\text{d}r}=\frac{-(p^{\prime }+2r\text{-}\unicode[STIX]{x1D704}^{2}B^{z2}+r^{2}\text{-}\unicode[STIX]{x1D704}^{\prime }\text{-}\unicode[STIX]{x1D704}B^{z2})}{B^{z}(1+r^{2}\text{-}\unicode[STIX]{x1D704}^{2})},\end{eqnarray}$$

where $p(r)$ and $\text{-}\unicode[STIX]{x1D704}(r)=\text{-}\unicode[STIX]{x1D704}_{i}=$ strongly irrational are assumed given, and assuming $B^{z}(0)=1$ . In the relaxed regions, the magnetic field is determined by continuing the integration according to

(2.8a,b ) $$\begin{eqnarray}\frac{\text{d}B^{\unicode[STIX]{x1D703}}}{\text{d}r}=\frac{\unicode[STIX]{x1D707}B^{z}-2B^{\unicode[STIX]{x1D703}}}{r},\quad \frac{\text{d}B^{z}}{\text{d}r}=-\unicode[STIX]{x1D707}rB^{\unicode[STIX]{x1D703}},\end{eqnarray}$$

where $\unicode[STIX]{x1D707}$ must be provided. The rotational transform will vary across the relaxed region in a fashion that is initially unknown. Thus, an iteration over $\unicode[STIX]{x1D707}$ is required to ensure that the rotational transform is equal to a prescribed, strongly irrational value, $\text{-}\unicode[STIX]{x1D704}_{i+2}$ , at the boundary of the next ideal region.

It is interesting to note that the perpendicular current density, $\boldsymbol{j}_{\bot }$ , is continuous. In the ideal regions, $\boldsymbol{j}_{\bot }=\boldsymbol{B}\times \unicode[STIX]{x1D735}p/B^{2}$ , and so if $\boldsymbol{B}$ and $\unicode[STIX]{x1D735}p$ are continuous, as they are, then $\boldsymbol{j}_{\bot }$ is continuous. The pressure gradient has been carefully chosen to ensure that $\unicode[STIX]{x1D735}p=0$ at the boundary separating the ideal and relaxed regions, and this ensures that $\boldsymbol{j}_{\bot }=0$ at these boundaries. In the relaxed regions, $\boldsymbol{j}_{\bot }$ is zero by construction.

The parallel current density is not continuous. In the ideal regions, the parallel current density depends on the given pressure and rotational-transform profiles, and is determined a posteriori as part of the equilibrium solution. In the relaxed regions, the parallel current density must be determined iteratively to ensure that the rotational-transform profile is continuous at the boundaries between the ideal and relaxed regions. There is insufficient freedom to ensure that the rotational transform is also smooth, and at the boundaries the parallel current density is generally discontinuous.

For the cylindrical example shown in figure 3, the cylindrical symmetry guarantees the existence of nested flux surfaces in both the ideal and the relaxed regions. In general geometry, this will not be the case. In the ideal regions, with strongly irrational rotational transform, the existence of nested flux surfaces is enforced by construction. The rational rotational-transform surfaces must lie in the relaxed regions, where relaxation and island formation is allowed. For the chaotic (irregular) field lines that may emerge, the rotational transform is not well defined.

We make some brief comments regarding a possible numerical construction for a general, 3-D equilibrium that is a combination of the algorithms already implemented in the VMEC (Hirshman & Whitson 1983; Hirshman et al. 1986) and SPEC (Hudson et al. 2012; Loizu, Hudson & Nührenberg 2016) codes. In the ideal regions, given the representation $\boldsymbol{B}=\unicode[STIX]{x1D735}\unicode[STIX]{x1D713}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+\text{-}\unicode[STIX]{x1D704}_{i}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\times \unicode[STIX]{x1D735}\unicode[STIX]{x1D713}$ , the numerical task amounts to finding the coordinate interpolation, $\boldsymbol{x}(\unicode[STIX]{x1D713},\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ , between the ${\mathcal{I}}_{i}$ that minimizes $W_{i}$ . This, essentially, is the approach adopted in VMEC. In the relaxed regions, by using a suitable gauge for the magnetic vector potential the magnetic can be represented as $\boldsymbol{B}=\unicode[STIX]{x1D735}\times \left(A_{\unicode[STIX]{x1D703}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D703}+A_{\unicode[STIX]{x1D701}}\unicode[STIX]{x1D735}\unicode[STIX]{x1D701}\right)$ , and the numerical task amounts to finding the functions $A_{\unicode[STIX]{x1D703}}(s,\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ and $A_{\unicode[STIX]{x1D701}}(s,\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ that extremize ${\mathcal{F}}_{i}$ , with suitable constraints imposed to enforce the boundary conditions that $\boldsymbol{B}\boldsymbol{\cdot }\boldsymbol{n}=0$ on the ${\mathcal{I}}_{i}$ and the flux constraints, and where $\boldsymbol{x}(s,\unicode[STIX]{x1D703},\unicode[STIX]{x1D701})$ is an arbitrary coordinate interpolation between the ${\mathcal{I}}_{i}$ . This is the approach adopted in SPEC. After computing the magnetic fields in each ${\mathcal{R}}_{i}$ , the geometry of the ${\mathcal{I}}_{i}$ must be adjusted (and the fields in each region recomputed) to satisfy continuity of the total pressure, $P\equiv p+B^{2}/2$ , across the ${\mathcal{I}}_{i}$ .

The new class of continuous solutions introduced herein can recover both classes of discontinuous solutions, namely those with discontinuous pressure and those with a discontinuous rotational transform, simply by letting the volume of the ideal or relaxed regions reduce to zero as desired. This can be enforced by constraining the toroidal flux in the appropriate regions.

Also, the number of volumes can become arbitrarily large. In practice, any acceptable pressure and transform profiles can be well approximated. Examples of what appear to be ‘fractal’ profiles are shown in figure 4.

Figure 4. Similar to figure 3, but with 64 regions.

This paper does not consider whether solutions with smooth pressure and continuous magnetic fields are preferable to the discontinuous solutions with sheet-currents. Ultimately, the question of which class of equilibria best models the observations may only be answered by validation. Towards this goal, it is certainly interesting to note that the pressure profile shown in figure 3(a) has a similar form to pressure profiles constructed by Ichiguchi et al. (2001) and Ichiguchi & Carreras (2011), who demonstrated that equilibria with flattened pressure across the rational surfaces seem to account for some experimental observations in the Large Helical Device (LHD) experiment.

On a similar note, we remark that the approach adopted in this paper cannot be trivially adapted to compute equilibria with a prescribed current-density profile. By taking the defining profiles to be the pressure and the current density, the rotational-transform profile is only known a posteriori. If, by chance, rational rotational-transform surfaces coincide with pressure gradients, then the pathologies associated with infinite currents and non-analyticity discussed in § 1.1.2 are not avoided. This suggests, for a given 3-D boundary, that not all pressure and current-density profiles are consistent with well-defined equilibrium solutions.

This point illustrates an implicit motivation for the approach adopted in this paper. The pathologies with 3-D equilibria are associated with rational rotational-transform surfaces. By specifying the rotational-transform profile a priori, we can ensure that these ‘surfaces’ are located in the relaxed regions, where the pressure gradient is zero and reconnection and island formation is allowed.

Even though a prescribed current density cannot be directly enforced, a given current-density profile can presumably be achieved, at least approximately, by iterating on the rotational transform.

We may expect that there will be a minimum allowed value for the jumps in the rotational transform across the relaxed volumes that are similar to the sine qua non condition described by Loizu et al. (2015b ). This condition is required to ensure that linear perturbation theory does not result in overlapping geometry, i.e. that the solutions are analytic functions of the 3-D boundary. We intend to explore this in future work.