Computation of multi-region, relaxed magnetohydrodynamic equilibria with prescribed toroidal current profile

The stepped-pressure equilibrium code (SPEC) (Hudson et al., Phys. Plasmas, vol. 19, issue 11, 2012, 112502) is extended to allow the computation of multi-region, relaxed magnetohydrodynamics (MRxMHD) equilibria at prescribed toroidal current profile. Toroidal currents are expressed in the framework of the MRxMHD theory, exhibiting spatial separation between pressure driven and externally driven currents. Additionally, analytical force balance derivatives at constant toroidal current are deployed in order to maintain SPEC's advantageous speed. The newly implemented capability is verified in screw pinch and classical stellarator geometries, and is applied to obtain the equilibrium $\beta$-limit of a classical stellarator without net toroidal currents. This new capability opens the possibility to study the effect of toroidal current on three-dimensional equilibria with the SPEC code.


Introduction
The computation of three-dimensional (3-D) magnetohydrodynamic (MHD) equilibria is of major importance for magnetic confinement devices such as tokamaks, reversed-field pinches or stellarators. In particular, the computation of equilibria at fixed toroidal current profile is crucial for basic physics studies (Loizu et al. 2017;Suzuki, Watanabe & Sakakibara 2020), equilibrium reconstruction (Lao et al. 1985;Hanson et al. 2009) and stellarator optimization (Geiger et al. 2010(Geiger et al. , 2015. Three-dimensional MHD equilibria consist of a complex mixture of nested flux surfaces, magnetic islands and chaotic field lines (Helander 2014;Hudson & Kraus 2017), hence their computation represents a great challenge. In fact, there is still no consensus in the community on how to best approach the computation of 3-D MHD equilibria (Hudson & Nakajima 2010). The multi-region, relaxed magnetohydrodynamic (MRxMHD) theory (Dewar et al. 2015;Hole, Hudson & Dewar 2006) has been developed to address this question. The MRxMHD minimizes the MHD energy functional (Kruskal & Kulsrud 2 A. Baillod, J. Loizu, Z.S. Qu, A. Kumar and J.P. Graves  1958) while keeping the magnetic helicity and the magnetic fluxes constant (Dewar et al. 2015) in a finite set of N vol nested volumes V l at constant pressure (see figure 1) but otherwise allowing arbitrary, non-ideal variations in the magnetic field. The interfaces I l separating the volumes are ideal flux surfaces which, therefore, cannot undergo magnetic reconnection, effectively constraining discretely the magnetic field topology. In the limit of a single volume, N vol = 1, Taylor's relaxation theory (Taylor 1974(Taylor , 1986) is recovered, while in the limit of an infinite number of volumes, N vol → ∞, it has been proved that MRxMHD recovers ideal MHD (Dennis et al. 2013a) in the limit of continuously nested flux surfaces, thereby bridging the gap between both theories.
The stepped-pressure equilibrium code (SPEC) (Hudson et al. 2012b,a) computes MRxMHD equilibria numerically given discrete input profiles. The SPEC code has been rigorously verified in stellarator geometry (Loizu, Hudson & Nührenberg 2016b), and has been successfully applied to study current sheets at rational surfaces (Loizu et al. 2015a,b), linear MHD instabilities (Kumar et al. 2021), tearing mode stability (Loizu & Hudson 2019) and nonlinear tearing saturation , equilibrium β-limits in a classical stellarator (Loizu et al. 2017), the penetration and amplification of resonant magnetic perturbations in the ideal limit (Loizu et al. 2016b) and relaxation phenomenon in reversed field pinches such as the formation of helical states (Dennis et al. 2013b) or the relaxation of flow during sawteeth (Dennis et al. 2014;Qu et al. 2020b). The MRxMHD has also been extended theoretically to study two-fluid effects (Lingam, Abdelhamid current on equilibrium beta limits, or the study of the sensitivity of a given equilibrium to toroidal current fluctuations. In this paper we describe the implementation of the new capability for SPEC, that allows MRxMHD equilibria to be calculated at prescribed toroidal current profiles. We first provide, in § 2, a brief reminder of the MRxMHD theory and a description of how plasma currents are evaluated in this framework. In § 3, the implementation of the new capability and its parallelization is described, and analytical derivatives that speed up the SPEC calculations are derived. In § 4, the new capability is verified against analytical results in cylindrical geometry and numerical results obtained with an older, verified (Loizu et al. 2016b) version of SPEC at fixed rotational transform in toroidal geometry. In § 5, we apply the new numerical tool to study the effect of pressure on the magnetic topology of a classical stellarator without net toroidal current, and compare the obtained ideal equilibrium β-limit to the value predicted by the high beta stellarator (HBS) theory (Wakatani 1998;Freidberg 2014). We conclude with a discussion in § 6.

Currents in MRxMHD
The plasma is divided into N vol nested volumes V l , l ∈ {1, . . . , N vol }, so that the MHD energy W l (Kruskal & Kulsrud 1958) local to each volume can be written as where p l is the pressure, B = |B| is the magnetic field strength, μ 0 is the vacuum permeability, γ is the adiabatic constant and dv is an infinitesimal volume element. The MRxMHD energy functional is (Hudson et al. 2012a) where μ l is a Lagrange multiplier, K l is the magnetic helicity in volume l and K l,0 the magnetic helicity constraint. The magnetic helicity is defined as where A l is the vector potential of the magnetic field B l , i.e. B l = ∇ × A l . In each volume V l , l ∈ {1, . . . , N vol }, the magnetic field B l is varied while keeping the toroidal magnetic flux, ψ t,l , and poloidal magnetic flux, ψ p,l , constant, until the MRxMHD energy (see (2.2)) is minimized. The corresponding Euler-Lagrange equations (Hudson et al. 2012a) describe a force-free magnetic field B l satisfying a Beltrami equation, In addition, the total pressure (plasma and magnetic pressure) is required to be continuous across each volume interface I l , In each volume V l , the solution to (2.4) is completely determined by three scalars (for example {μ l , ψ t,l , ψ p,l }), the geometry of the interfaces bounding the volume and a boundary condition for the magnetic field normal to the interfaces B l ·n k , k = {l − 1, l}, withn k a unit vector perpendicular to the interface k. In MRxMHD, the interfaces are imposed to be flux surfaces, with B l ·n k = 0. In the innermost volume, which is topologically different from the others, only two scalars are required in addition to the geometrical degrees of freedom and the condition B 1 ·n 1 = 0.
In MRxMHD, two spatially distinct net toroidal current profiles coexist, namely currents flowing in the volumes, {I v l,φ } l={1,...,N vol } , and surface currents flowing at the volumes' interfaces, {I s l,φ } l={1,...,N vol −1} (current sheets), where the subscript φ refers to the toroidal angle. The volume current I v l,φ in volume V l is easily evaluated using (2.4) and Ampere's law, where S l,φ is a constant-φ surface in volume V l and dS l,φ is the differential surface element normal to S l,φ . Volume currents include externally driven currents such as electron cyclotron current drive (ECCD), neutral beam current drive (known as NBCD) or Ohmic current. Equation (2.6) might be surprising since toroidal currents are usually expressed in terms of functions of the poloidal fluxes and not the toroidal fluxes. In essence, the poloidal flux dependence is contained in μ l , which is related to the parallel current density, as μ l = μ 0 j l · B l /B 2 l , with j l the current density in volume V l . The surface current I s l,φ at interface I l can be evaluated using Ampere's law, where Γ l is a closed curve following the interface I l poloidally andB θ is the m = n = 0 Fourier mode of the covariant component of the poloidal magnetic field. In (2.7), the poloidal and toroidal angles, θ and φ, are as-of-yet arbitrary. However, the surface currents I s l,φ are, as expected, independent of these angles choice, since the surface currents only depend on the m = n = 0 mode of the field. Surface currents represent all equilibrium pressure-driven currents, such as diamagnetic, Pfirsch-Schlüter and bootstrap currents, as well as shielding currents arising when an ideal interface is positioned on a resonance (Loizu et al. 2015a). As a side note, we remark that while ideal MHD equilibria are defined by two free functions (e.g. the pressure and the rotational transform profiles, p(ψ t ) and ι -(ψ t ), or the pressure and the current profiles, p(ψ t ) and I φ (ψ t )), MRxMHD requires two scalars to determine the solution in a volume V l , in addition to the pressure and toroidal flux. This can be considered as three independent discrete profiles that are required to determine an equilibrium. Examples are {p l , μ l , ψ p,l } l=1,...,N vol , {p l , μ l , K l } l=1,...,N vol or {p l , ι -− l , ι -+ l } l=1,...,N vol , with ι -± l the rotational transform on the inner and outer side of the interface I l , or {p l , I v l,φ , I s l,φ } l=1,...,N vol , as functions of {ψ t,l } l=1,...,N vol .
2.1. Current discretization Typically, continuous current profiles are provided by analytical models or after equilibrium reconstruction using experimental data. We now discuss how these profiles can be represented in the framework of MRxMHD. Consider an externally driven current profile, e.g. ECCD, provided as the enclosed toroidal current as a function of the toroidal magnetic flux, i.e. I φ,ECCD (ψ t ), and a pressure-driven current profile, e.g. the bootstrap FIGURE 2. Sketch of a pressure profile as a function of the toroidal flux. Blue, continuous pressure profile obtained via experiment or analytical model; red, SPEC discretized pressure profile; black dashed lines, volume interfaces. current, provided similarly as the enclosed toroidal current as a function of the toroidal flux, I φ,BS (ψ t ). We also assume that the pressure profile, p(ψ t ), the number of volumes, N vol , and their enclosed toroidal fluxes, {ψ t,l } l=1,...,N vol , are given (see figure 2). The question of how many volumes and where to position their interfaces to best represent a given pressure profile is not addressed in this paper.
A proposed representation of these current density profiles in MRxMHD is achieved as follows. The ECCD current is an externally driven, parallel current and is thus represented as a volume current since it flows parallel to the field lines; on the other hand, the bootstrap current is a pressure-driven, self-generated current and is represented as a surface current, since it is localized at the pressure gradients. The volume currents are obtained by integrating the externally driven current density in each volume (figure 3), which is simply given by the difference and the surface currents are obtained by integrating the pressure driven current density around each interface (figure 4), which is expressed as 2 , otherwise, (2.11)  with ψ a the total toroidal magnetic flux enclosed by the plasma. In (2.10)-(2.11), care has been taken for the first and last interfaces, where the surface of integration has been extended to include the current density from the magnetic axis and up to the plasma boundary. Note that this difference in the definition of the first and last surface currents vanishes as the number of volumes N vol is increased. Equations (2.8)-(2.11) are only one possible discretization of the continuous current profiles, proposed by the authors for illustration. Advantages of this particular representation are: (i) that the total toroidal current is always preserved; (ii) that the currents are approximately localized at the same location in the discretized than in the continuous case. The following sections do not depend on the particular choice of discretization of the currents.

Implementation of the current constraint in SPEC
3.1. The SPEC code The SPEC code can run in three different geometries, namely in slab (Loizu & Hudson 2019), cylindrical (Loizu et al. 2016a) and toroidal geometry (Loizu et al. 2016b). The coordinates used to describe position are in slab geometry, R cos θî + R sin θĵ + φk, in cylindrical geometry, R cos φî + R sin φĵ + Zk, in toroidal geometry, where {î,ĵ,k} is the unitary Cartesian basis, θ is a poloidal angle and φ is the usual toroidal angle. At each interface I l , the geometry of the surface is described by decomposing the coordinates R and Z on a Fourier basis, where N p is the number of field periods, M pol and N tor are the poloidal and toroidal mode numbers above which Fourier series are truncated, i.e. m = {0, . . . , M pol }, n = {−N tor , . . . , N tor }, and stellarator symmetry has been assumed for simplicity. Between interfaces l and l + 1, coordinates are constructed by linear interpolation of (R, Z) l and (R, Z) l+1 using a radial-like coordinate s. Geometrical degrees of freedom are {R l,m,n , Z l,m,n } ≡ {x i } i=1,...,N , with, e.g. in the case of a fixed-boundary, stellarator symmetric equilibrium in toroidal geometry, The magnetic field is represented by the magnetic vector potential, B l = ∇ × A l , which is described by using a Chebyshev-Fourier basis, where L rad is the radial resolution, T k is the Chebyshev polynomial of order k and i = {θ, φ}. Note that gauge freedom is used to set A l,s = 0 ∀l. The coefficients A l,i,k,m,n are the vector potential degrees of freedom, packed in a single array a l per volume. Fixed-boundary SPEC can solve the MRxMHD system defined by (2.4) and (2.5) for a given plasma boundary geometry, pressure profile {p l } l=1,...,N vol and profiles {ψ t,l , ψ p,l , μ l } l=1,...,N vol . In this case, the Beltrami equation in volume l, (2.4), can be cast into a linear system G l · a l = C l , (3.6) where G l and C l are matrices that depend only on the geometrical degrees of freedom ..,N , and linearly on ψ t,l , ψ p,l and μ l .
The constrained quantities {μ l , ψ t,l , ψ p,l } l=1,...,N vol can be replaced by other independent triplets such as {ψ t,l , ι -+ l , ι -− l } l=1,...,N vol , with ι -± l the rotational transform evaluated on the inner and outer side of interface I l , respectively. In this case, SPEC iterates on {ψ p,l , μ l } l=1,...,N vol until the solution satisfies the input rotational transform profile. In what follows, this method will be referred to as the 'rotational transform constraint'. Finally, SPEC uses a Newton method to iterate on {x i } i=1,...,N until the force balance, (2.5), is satisfied. This is achieved by expanding the force imbalance on a Fourier basis, and packing the Fourier components F l,m,n in a single array F of size N. Force balance is satisfied once all the Fourier modes F l,m,n , herein {F j } j=1,...,N , are sufficiently small.
It is worth noting that even though the magnetic helicity is not conserved when SPEC iterates on {x i } i=1,...,N to find an equilibrium that matches a given input ..,N vol , the final equilibrium satisfies the MRxMHD equilibrium equations ((2.4)-(2.5)). There is a magnetic helicity profile {K l } l=1,...,N vol corresponding to this equilibrium which is unknown a priori. Thus, the same equilibrium could have been found by minimizing the MRxMHD energy functional while keeping the magnetic helicity profile constant if the initial state had the same magnetic helicity profile (bifurcations are not considered in this paper). This capability is also available in SPEC, and details can be found in the literature (Hudson et al. 2012a;Dennis et al. 2013b).
Lastly, SPEC has been recently extended to allow free-boundary calculations , where the plasma is surrounded by a vacuum region bounded by a computational boundary. The solution in the vacuum region, which is a special instance of a Beltrami field, with μ = 0, is determined by a different set of scalars, namely the total poloidal current passing through the torus hole, I coil , and the total toroidal plasma current, I pl . Using Ampere's law, these currents can be written as withB + V,φ the m = n = 0 Fourier mode of the covariant toroidal magnetic field on the computational boundary andB − V,θ the m = n = 0 Fourier mode of the covariant poloidal magnetic field on the outer side of the plasma boundary. Additionally, the computational boundary is not necessarily a flux surface, thus B V ·n c = 0 with B V the magnetic field in the vacuum region andn c a unit vector normal to the computational boundary. The solution satisfying the constraints I coil and I pl is obtained by iteration over the poloidal and toroidal fluxes {ψ p,V , ψ t,V } in the vacuum volume.

Constraining the toroidal current profiles in SPEC
The SPEC code has been extended to allow the triplet {ψ t,l , I v l,φ , I s l,φ } l=1,...,N vol as a constraint, herein 'current constraint'. In the case of the rotational transform constraint, SPEC finds the solution to the linear system (2.4) volume by volume and iterates on {μ l , ψ p,l } until the field has the desired rotational transform at the volume interfaces. In the case of the current constraint, the constants {μ l } l=1,...,N vol are determined using (2.6), without the need for iterations, and this directly constrains the value of volume currents {I v l,φ } l=1,...,N vol . Regarding the poloidal fluxes, it can be shown (Appendix A) that the surface currents depend linearly on the poloidal fluxes, i.e.
( 3.10) where I s and ψ p are arrays containing all {I s l } l=1,...,N vol −1 and all {ψ p,l } l=2,...,N vol , respectively, and the matrix M and the array Q depend only on the geometry of the interfaces {x i } i=1,...,N . In this section, we consider the geometry, toroidal fluxes and the constants {μ l } l=1,...,N vol to be fixed and seek how the poloidal flux profile ψ p has to be constrained in order to obtain a surface current profile matching the input profile I s .
The unknown Q is eliminated by subtracting (3.10) evaluated at two different values of ψ p , i.e. evaluated once at ψ p and once at ψ p , where ψ p is an arbitrary choice of poloidal fluxes, and I s is the surface current profile calculated from the Beltrami fields {a l } l=1,...,N vol obtained when the poloidal fluxes are constrained to the values ψ p . The matrix M is evaluated by taking the derivatives of (3.10) with respect to the poloidal fluxes, i.e.
which leads to the following bidiagonal matrix: (3.13) Derivatives ofB ± θ,l with respect to the poloidal flux can be easily obtained by applying matrix perturbation theory on the linear system (3.6) (Hudson et al. 2012a), (3.14) Due to the linear nature of (3.10), the coefficients of the matrix M, i.e. (3.12), are independent of ψ p and thus can be evaluated once at any arbitrary value of ψ p . We thus conveniently evaluate them at ψ p . Equation (3.11) is then solved to obtain the poloidal flux profile ψ p . Finally, instead of solving a second time the Beltrami equation, (2.4), at ψ p , we take advantage of the linear dependence of a on ψ p , (i.e. 3.6), and solve where A l,i is one element of a l . The algorithm flow is summarized in figure 5.
In the case of a free-boundary computation, the toroidal flux in the vacuum region is varied to satisfy the poloidal linking current I coil . This slightly modifies the linear system (3.11). Details can be found in Appendix B.
Constraining the toroidal current profile takes away the control of other profiles such as the rotational transform or the magnetic helicity. However, as for the case of the rotational transform constraint, the equilibrium can be accessed by a relaxation process at constant magnetic helicity if the final magnetic helicity is known a priori. The MRxMHD equations are thus still satisfied by an equilibrium obtained by constraining the toroidal current profiles.
The new current constraint has been parallelized with MPI in a similar fashion to the other constraints. Each volume is associated with one central processing unit (CPU); since the solution to the Beltrami equation (2.4) in a volume is independent from other volumes, each CPU can solve the linear system (3.6) in parallel. Finally, the master CPU gathers all required derivatives to construct the matrix M and solves the linear system (3.11), before broadcasting the values of {ψ p,l } l=2,...,N vol and {a l } l=1,...,N vol to all CPUs.

Force gradient
The Newton algorithm used in SPEC to iterate on the interfaces' geometry uses analytic derivatives, which is faster than finite differentiation. To keep good performance while using the current constraint, derivatives of the force Fourier coefficients, {F j } j=1,...,N , with respect to the geometrical degrees of freedom, {x i } i={1,...,N} , at constant {ψ t,l , I v l,φ , I s l,φ } l=1,...,N vol are provided. Derivatives are first evaluated in real space and then Fourier transformed. Using the chain rule, where B − l , B + l are the magnetic field strength on the inner and outer side of volume l, respectively, and the pressure, p l , is considered constant in each volume with respect to variations in the geometry, μ l and ψ p,l . Note that all derivatives are taken at constant toroidal flux, volume current and surface current. Enforcing dψ t,l /dx i = 0 and dI v l,φ /dx i = 0 leads to dμ l /dx i = 0 using (2.6). The surface current constraint, dI s l,φ /dx i = 0, leads to a system of coupled equations using (2.7), which can be written as a linear system using the matrix M defined in (3.13), Derivatives ofB θ,l with respect to ψ p,l and x i can be obtained by applying matrix perturbation theory to the Beltrami system (3.6). The solution of (3.19), together with (3.17) provides the required derivatives of the force with respect to the geometry. The derivatives of the Fourier components of the force, dF j /dx i , are obtained by taking the Fourier transform of (3.17) and are packed in a matrix ∇F of size N 2 , henceforth named 'force gradient'. Appendix B provides details on the free-boundary case.

Verification of the current constraint
In this section we present a rigorous verification of the new capability of SPEC against analytical solutions in a screw pinch geometry and against a reference SPEC solution obtained with the rotational transform constraint in a classical stellarator geometry. All results presented in this paper were obtained with SPEC version 2.10.

Verification in cylindrical geometry
We consider a fixed-boundary screw pinch MRxMHD equilibrium that only depends on the radius R and whose solutions can be written analytically (Appendix C). We choose a set of somewhat arbitrary input parameters, i.e. a cylinder of minor radius a = 1 and length L = 2π, N vol = 3, p l = 0 ∀l ∈ {1, 2, 3}, ψ t = {1/9, 4/9, 1} Tm 2 , μ 0 I v = {0.2, 0.2, 0.4} Tm and μ 0 I s = {−0.4, 0.5} Tm, which uniquely define the analytical solution. The SPEC code is then run with the same input parameters and the solutions are compared (figure 6). Very good agreement between the analytical solution and the SPEC solution is obtained. Note that by constraining the toroidal current profiles, we lose control on the rotational transform profile, since only two profiles can be constrained in addition to the pressure profile. Hence discontinuities in the magnetic field components arise at the volume interfaces, even when there is no pressure, unless the input parameters are carefully selected.
The force gradient can also be expressed in terms of Bessel function integrals (see Appendix C). Figure 7 shows the normalized maximum absolute error between the force gradient obtained with SPEC and that obtained analytically as a function of the radial resolution L rad . As L rad is increased, exponential convergence is observed up until 10 −13 , where the error in the evaluation of the Bessel integrals starts to dominate.

Verification in toroidal geometry
A verification is proposed here in the more complex case of a free-boundary, rotating ellipse (also called classical stellarator) equilibrium with five field periods (N p = 5), multiple poloidal and toroidal modes (M pol = 4 and N tor = 2) and seven plasma volumes (N vol = 7). The pressure is set to zero and the computational boundary is defined by with R 00 = 10 m, Z 00 = 0 m, R 10 = −Z 10 = 1 m, R 11 = Z 11 = 0.25 m and N p = 5 (see figure 8). We suppose here that some hypothetical coils with a total current of μ 0 I coil = 42.87 Tm are able to generate a vacuum field without normal component to the computational boundary. The total toroidal magnetic flux in the plasma is set to ψ a = 0.61 Tm 2 and the toroidal magnetic flux in the vacuum region is set to ψ t,V = 1.39 Tm 2 , adding up to a total toroidal magnetic flux enclosed by the computational boundary of ψ t,V + ψ a = 2 Tm 2 .
To the authors' knowledge, no analytical solution to (2.4) and (2.5) exists in this geometry. The verification is thus carried out as follows. First, a rotational transform constraint case is run with an input ι --profile that is chosen to be 10 % larger than the vacuum rotational transform ι vac , i.e. ι -= 1.10 × ι vac , so that there is a non-zero contribution from the current to the rotational transform. The volume and surface currents are evaluated from the obtained equilibrium and used to run a current constraint calculation to obtain a second equilibrium. The same initial guess for the geometry and the interfaces is used for both calculations. The rotational transform profileῑ -is then extracted from the second equilibrium and compared with the reference ι --profile. The vacuum rotational transform profile, as well as the profiles ι -andῑ -are shown in figure 9(a). The toroidal current enclosed by the plasma is mostly contained in the volumes and adds up to a total of ∼2.7 kA, see figure 9(b). As expected the surfaces currents I s φ,l remain small (< 10 −2 kA), since there are no pressure gradients to drive them. The constraint on the rotational transform ι -is enforced on each side of the volumes' interfaces, indicated by grey dashed lines on figure 9(a). The value of ι -at the computational boundary is not constrained. Agreement between ι -andῑ -up to a relative error of max(|ι -−ῑ -|/|ι -|) ∼ 10 −5 is observed, showing that the same equilibrium can be obtained using either constraint. The maximum error between both profiles decreases as the numerical resolution is increased (data not shown). To verify the force gradient components, we use a fourth-order centred finite difference formula (Fornberg 1988), to obtain ∇F FD , i.e. a finite-difference estimate of the force gradient, and compare it with ∇F, i.e. the force gradient calculated in SPEC by using analytical derivatives. The finite difference estimate is evaluated by perturbing each geometrical degree of freedom {x i } l=1,...,N by a constant value R. Convergence as R → 0 is shown in figure 10. A convergence of order O( R 4 ) is observed down to ∼ 10 −11 for R ∼ 10 −4 . For lower values of R, the finite difference approximation error is dominated by round-off error. This shows that the analytical derivatives (the force gradient) are correctly implemented in SPEC.

Ideal equilibrium β-limit in a classical stellarator
In this section, it is shown that the SPEC current constraint can be used to recover the HBS theory prediction of the classical stellarator ideal equilibrium β-limit (Freidberg 2014;Wakatani 1998), when zero net toroidal current is considered.
A previous study (Loizu et al. 2017) showed remarkable agreement between SPEC calculations and the HBS theory in a simplified case, with the pressure profile approximated by a single step. An additional limitation, enforced by the SPEC version at that time, was that the vacuum region had to be approximated by a large plasma volume where the pressure and currents were set to zero, and a fixed-boundary would be applied. This approach is equivalent to assuming a free-boundary calculation with an infinitely conducting wall at the computational boundary. At large β, however, a strong Shafranov shift is present, and this approximation could have an impact on the result.
To understand the effect of these assumptions we consider here both fixed-and free-boundary calculations with a stepped-pressure profile approximating a Solov'ev's profile p = p 0 (1 − ψ t /ψ a ) where ψ a is the total toroidal flux enclosed by the plasma. The In panel (c), the ideal equilibrium β-limit has been exceeded and a central island opened outside the plasma. A large value of β has been selected for illustration purposes. Here (a) ( β = 0 %, ι a ≈ 0.28); (b) ( β ≈ 0.31 %, ι a ≈ 0.13); (c) ( β ≈ 0.62 %, ι a = 0.0). computational boundary is defined by (4.1)-(4.2). Fixed-boundary calculations were made first, with seven plasma volumes surrounded by an eighth large, vacuum-like volume so that it is similar to a free-boundary calculation. The toroidal flux profile has been chosen such that ψ a = 7 l=1 ψ t,l = 0.25 Tm 2 and ψ t,8 = 0.75 Tm 2 , leading to a total toroidal flux enclosed by the plasma and the vacuum-like region of 8 l=1 ψ t,l = 1 Tm 2 . Free-boundary input files were generated to replicate the same equilibrium in vacuum, using μ 0 I coil = 21.43 Tm and ψ a = 0.25 Tm 2 , which correspond to an equilibrium with ψ t,V = 0.75 Tm 2 . In both fixed-and free-boundary calculations, both the surface currents and the volume currents were set to zero, i.e. I vol φ,l = I surf φ,l = 0 for all l. The only control parameter remaining is p 0 , which was increased in order to increase the plasma average β until a central m = 1, n = 0 island opens. The island emerges when the rotational transform at the plasma edge, ι a , i.e. the rotational transform on the outer side of the plasma-vacuum interface, reaches zero (see figure 11). This is defined as the ideal equilibrium β-limit.
The physical mechanisms leading to the emergence of a separatrix are complex. In brief, the combination of non-zero poloidal harmonics of the Pfirsch-Schlüter and diamagnetic currents at the volumes' interfaces and the effect of the Shafranov shift perturbs sufficiently the poloidal magnetic field so that the rotational transform at the plasma edge decreases, until eventually it reaches zero. These results will be explained in more detail in a future publication.
The HBS theory, building on the pioneering work of Greene and Johnson's stellarator expansion theory (Greene & Johnson 1961), predicts that with ι v the rotational transform at the plasma edge in vacuum and ι -I a contribution from the plasma toroidal current, with I φ the total toroidal plasma current, B 0 the magnetic field strength on axis and a the effective minor radius at the plasma edge, i.e. a = ψ a /(ψ a + ψ t,V )r eff , with r eff the effective radius of the ellipse, given by the square root of the product between the ellipse major radius r max = |R 10 + R 11 | and the ellipse minor radius r min = |Z 10 + Z 11 |, i.e. r eff = √ r max r min . The parameter ν is defined as with β the volume averaged plasma β and a = a/R 00 . Setting ι a = 0 in (5.1) and solving for β leads to a prediction for the ideal equilibrium β-limit, In the case of interest, R 00 = 10 m, R 10 = −Z 10 = 1 m, R 11 = Z 11 = 0.25 m, leading to r eff ≈ 0.97 m and a ≈ 0.48 m. Using I φ = 0A and ι v ≈ 0.28, one obtains β lim,HBS ≈ 0.38 %. Figure 12 shows the rotational transform profile at different values of β (figure 12a) and the values of ι a obtained with SPEC as a function of β and compares them with the analytical prediction given by (5.1) (figure 12b). Good agreement is observed between the fixed-boundary and free-boundary calculations, showing that the fixed-boundary assumption made by Loizu et al. (2017) has only a small effect on the ideal equilibrium β-limit prediction. In addition, the free-boundary calculation and the HBS theory agree well, and predict approximately the same ideal equilibrium β-limit. The small but finite difference between the results is most likely due to the HBS theory, which employs an expansion in aspect ratio and assumes that the plasma boundary is circular. In fact, the number of volumes used in SPEC to approximate the continuous pressure profile does not significantly influence the result of this study (data not shown).
As a final remark, note that in our calculations, the plasma boundary is topologically constrained to be an ideal surface and cannot undergo reconnection. Another important constraint is that the pressure profile is fixed and cannot evolve. Other approaches without these constraints could lead to different results and conclusions. One natural question is then to ask about the physical validity of constraining the topology of a set of flux surfaces in MRxMHD. First studies about the existence of MRxMHD interfaces have been carried out (McGann et al. 2010;Qu et al. 2021), pointing towards certain existence criteria. The present study, however, aims at verifying the implementation of the toroidal current constraint in SPEC by retrieving well-established mathematical results. Future work will focus on code validation by using experimental data.

Conclusion
Toroidal currents in MRxMHD have been derived and expressed in terms of simple quantities. Two profiles have been identified: the volume current profile, flowing through the volumes; and the surface current profile, flowing at each volume interface. A physical interpretation has been given for each of the currents. Both profiles have been implemented as new constraints in the SPEC code, which can now compute MRxMHD equilibria for a given toroidal current profile. Analytical derivatives of the force on each volume interface with respect to the interfaces' geometry at fixed toroidal current have been derived and implemented in SPEC. These derivatives speed up substantially the Newton iterations on the interface geometries.
Both the new constraint and the force gradient implementation have been verified in slab, cylindrical and toroidal geometries. We presented in this paper only the latter two. In cylindrical geometry, we considered an axisymmetric screw pinch, where the obtained equilibria and force gradient could be compared with analytical solutions. In toroidal geometry, a classical stellarator geometry has been considered. The equilibrium has been verified to match the equilibrium obtained by constraining the rotational transform profile in SPEC, and the force gradient has been compared with a finite difference estimate.
Finally, the calculation of the ideal equilibrium β-limit in a classical stellarator with zero net toroidal current has been presented as a first application of the new SPEC capabilities. The free-boundary calculation showed very good agreement both with a previous, simpler study (Loizu et al. 2017) and with the HBS theory (Wakatani 1998;Freidberg 2014). In the future, the effect of bootstrap current on the equilibrium β-limit in a classical stellarator will be investigated. Furthermore, this new tool opens the possibility of comparing SPEC with previous results obtained with other equilibrium codes such as PIES (Drevlak et al. 2005;Hirsch et al. 2008), and perform similar studies on present experiments such as W7-X. In particular, it is envisaged to study W7-X equilibrium β-limits and its robustness to fluctuations in the toroidal current profile.

Editor Per Helander thanks the referees for their advice in evaluating this article.
Funding This work has been carried out within the framework of the EUROfusion Consortium and has received funding from the Euratom research and training program 2014-2018 and 2019-2020 under grant agreement no. 633053. The views and opinions expressed herein do not necessarily reflect those of the European Commission (A.B., J.L., J.P.G.). This work is partly funded by Australian ARC projects DP170102606 (Z.Q., A.K.). This work was also supported by a grant from the Simons Foundation/SFARI (560651, AB) (A.B., J.L., Z.Q., A.K.).

Declaration of interests
The authors report no conflict of interest.

Appendix A
In this appendix we show the linear relation between the surface currents (2.7) and the poloidal fluxes. We rewrite here (2.7) for convenience, We use general coordinates notation, with u i ≡ {s, θ, φ}, {∇s, ∇θ, ∇φ} the contravariant basis and e i ≡ {e s , e θ , e φ } the covariant basis. This derivation is local to a volume and we drop the subscript l everywhere for simplicity.
We first show that the surface currents depend linearly on the vector potential degrees of freedom a. The contravariant components of the magnetic field are obtained from ∇ × A = B, where ijk is the Levi-Civita tensor, √ g is the Jacobian, and the Einstein summation convention has been used. The covariant component of the magnetic field can then easily be expressed by B θ = g kθ B k . The m = n = 0 Fourier mode of B θ is theñ with S the total area of the flux surface, which depends only on geometrical quantities.
where the prime denotes the derivative with respect to the main argument. Equations (A3) and (A4)-(A6) combined show the linear dependence ofB θ on a. Finally, the Beltrami equation (3.6) provides a linear relation between a and {ψ p , ψ t , μ}. All relations being linear, this shows that the surface currents depend linearly on the poloidal and toroidal magnetic fluxes.

Appendix B
In this appendix, the main differences between a fixed-and free-boundary calculation with the new developed current constraint are outlined. The linear system (3.11) has to be rewritten by extending the arrays ψ and I with two new pairs of scalars (ψ p,V , ψ t,V ) and (I s N vol , I coil ), namely ψ ≡ (ψ p,2 , . . . , ψ N vol , ψ p,V , ψ t,V ) t and I ≡ (I s 1 , . . . , I s N vol , I coil ) t . Then, with the matrix M Fr , withB − φ,V the m = n = 0 Fourier mode of the covariant toroidal magnetic field on the plasma boundary outer side. Regarding (3.15), no changes are needed in the plasma volumes. In the vacuum region, however, the toroidal flux is not fixed and an additional term is needed, where the subscript V denotes the vacuum region. Regarding the force gradient, the derivative of the toroidal flux with respect to the geometry is non-zero in the vacuum region. This means that An additional equation is required for dψ t,V /dx i , and is provided by leading to Appendix C The solution to the Beltrami equation (2.4) in the lth volume of an axisymmetric cylinder is B l = c l,1 rJ 1 (μ l r) + c l,2 rY 1 (μ l r) ∇θ + c l,1 J 0 (μ l r) + c l,2 Y 0 (μ l r) ∇φ, where the usual (r, θ, φ) cylindrical coordinate system has been used, J i and Y i are the Bessel functions of the ith order of the first and second kind, respectively, and c l,1 , c l,2 are integration constants. Here ∇θ and ∇φ are the contravariant basis vectors. In addition, since B θ must vanish at the origin, we have that c 1,2 = 0. Indeed, the asymptotic expansion of Y 1 (x) close to x = 0 gives (Abramowitz & Stegun 1972) lim r→0 c 1,2 rY 1 (μ 1 r) ∼ lim r→0 −c 1,2 r 2 πr = − 2c 1,2 π , which is only zero if c 1,2 = 0.
We consider now the case of a screw pinch with three inner volumes, N vol = 3. The assumed constrained profiles are the toroidal flux {ψ t,l } l=1,2,3 , the volume current {I v l,φ } l=1,2,3 and the surface current {I s l,φ } l=1,2 . The constraint on the toroidal flux is ψ t,l = S l,φ B · ∇φ √ g dr dθ (C3) = R l R l−1 dr 2π 0 dθ c l,1 J 0 (μ l r)r + c l,2 Y 0 (μ l r)r (C4) where √ g = r is the Jacobian and S l,φ is a constant-φ surface in volume l. The Bessel function integrals have been renamed asJ l andȲ l , and R l is the radius of the lth interface. The constraints on the currents lead to Solving for {c l,1 , c l,2 } is equivalent to solving the linear system Derivatives of the force F l = [(B l+1 (R l )) 2 − (B l (R l )) 2 ]/2 can also be expressed analytically, leading to with l, j ∈ {1, 2}. Consider, e.g. the derivative of B l (R k ), k = {l − 1, l}, [B l (R k )] 2 = c l,1 J 1 (μ l R k ) + c l,2 Y 1 (μ l R k ) 2 + c l,1 J 0 (μ l R k ) + c l,2 Y 0 (μ l R k ) 2 , (C10) B l ∂B l ∂R k = (c l,1 J 1 + c l,2 Y 1 )(c l,1 J 1 + c l,1 μ l J 1 + c l,2 Y 1 + c l,2 μ l Y 1 ) +(c l,1 J 0 + c l,2 Y 0 )(c l,0 J 0 + c l,2 μ l J 0 + c l,2 Y 0 + c l,2 μ l Y 0 ), where the denotes a derivative with respect to the function argument, and all Bessel functions are evaluated at μ l R k . Finally, all derivatives must be taken at constant ψ t,l , I v l,φ and I s l,φ . In particular, the coefficients dc l,i /dR k are obtained from derivatives of (C8) with respect to R k .