MHD equilibria with nonconstant pressure in nondegenerate toroidal domains

We prove the existence of piecewise smooth MHD equilibria in three-dimensional toroidal domains of $\mathbf{R}^3$ where the pressure is constant on the boundary but not in the interior. The pressure is piecewise constant and the plasma current exhibits an arbitrary number of current sheets. We also establish the existence of free boundary steady states surrounded by vacuum with an external surface current. The toroidal domains where these equilibria are shown to exist do not need to be small perturbations of an axisymmetric domain, and in fact they can have any knotted topology. The building blocks we use in our construction are analytic toroidal domains satisfying a certain nondegeneracy condition, which roughly states that there exists a force-free field that is ergodic on the surface of the domain. The proof involves three main ingredients: a gluing construction of piecewise smooth MHD equilibria, a Hamilton-Jacobi equation on the two-dimensional torus that can be understood as a nonlinear deformation of the cohomological equation (so the nondegeneracy assumption plays a major role in the corresponding analysis), and a new KAM theorem tailored for the study of divergence-free fields in three dimensions whose Poincar\'e map cannot be computed explicitly.


Introduction
The magnetohydrodynamics (MHD) equilibria in a smooth bounded domain Ω ⊂ R 3 are often described by the solutions to the system of equations Here the vector field B(x) is the steady magnetic field of a perfectly conducting plasma, the scalar function P (x) is the plasma pressure and N (x) denotes the outer unit normal at a boundary point x ∈ ∂Ω.It is also customary to assume that P is constant on the boundary ∂Ω, so we can set (1.1d) P = 0 on ∂Ω .
Equations (1.1), which we will referred to as the fixed boundary problem, model the equilibrium configurations of a plasma confined in the fixed magnetic domain Ω.Because of this connection, the most interesting case is when Ω is a toroidal domain, that is, when the boundary of the domain is diffeomorphic to a torus.It is worth noting that these equations also describe stationary solutions to the 3D Euler equations in fluid mechanics, with B playing the role of the velocity field of the fluid and P being the negative of the Bernoulli function (see e.g.[1]).
In this paper we will be concerned with the problem of finding toroidal domains that admit MHD equilibria whose pressure is constant on the boundary but not in the interior.This problem can be traced back at least to the work of H. Grad in the 1960s, who conjectured [18] that no smooth solutions fibring the domain with toroidal surfaces exist unless the domain is axially symmetric.An important somewhat related result, due to Bruno and Laurence [8] in the 1990s, is the existence of weak solutions with nonconstant stepped-pressure in nonsymmetric toroidal domains that are small perturbations of an axisymmetric solid torus.A very illuminating numerical implementation of this model suggesting the existence of stepped-pressure equilibria in toroidal domains far from axisymmetric was developed in [21,22].See also [10] for the construction of smooth solutions with nonconstant pressure, but subject to a small force, on toroidal domains that are close to be axisymmetric.However, Grad's influential conjecture remains wide open.A comprehensive recent account of the subject can be found in [19,20].
It is important to remark that in the context of plasma physics, the pressure profile P (usually described as a function of the magnetic flux) and the rotational transform profile (which is, essentially, the ratio of the components of the frequency vector of the magnetic field on each invariant toroidal surface) are prescribed a priori, as inputs of the computations.The relevance of these prescriptions lays in the fact that they yield uniqueness in the solutions of Equations (1.1), if they exist.In contrast, in this work we will be interested in the existence of stepped-pressure MHD equilibria in domains very far from symmetric, without imposing any profile, and then the pressure and the rotational transform profiles show up a posteriori.The same comment applies to the free boundary problem that is next stated.
Another related equation that appears in plasma physics, particularly concerning the design of plasma confinement devices for nuclear power generation, describes a free boundary steady state surrounded by vacuum with an external current J ext .In terms of the interior and exterior magnetic fields, B and B ext , this system reads in this case as For simplicity, we have assumed that P = 0 on ∂Ω.The jump condition (1.2f) uses that the (hydrodynamic) pressure in the vacuum is simply 1  2 |B ext | 2 .In this work we shall assume that the external current is a current sheet, i.e., (1.2h) where dS is the surface measure on the boundary ∂Ω of a certain domain Ω enclosing Ω and J : ∂Ω → R 3 is a tangent vector field on the surface.We want to stress that in applications in the context of plasma confinement (stellarators), the external current J ext is generated by a certain number of current coils, but the techniques we use in this work do not allow us to deal with this case (see Remark 1.4).We will additionally impose the tangency condition (1.2i) B • N = 0 on ∂Ω and refer to the system of equations (1.2) as the free boundary problem.
We observe that a related, although different, free boundary problem for MHD equilibria has been considered in the literature, see e.g.[29,3] and references therein.In these works, the problem consists in determining the shape of an axisymmetric plasma region, which is surrounded by vacuum and contained in a given axisymmetric shell.It is assumed that the shell is a perfect conductor, so no external current is introduced in the vacuum region; moreover the vacuum field is assumed to be tangent to the boundaries of the shell and of the plasma, and the continuity of the hydrodynamic pressure is also imposed (in order to have a weak solution of the MHD equations across the plasma boundary).In contrast, in the free boundary problem that we consider above, the vacuum vessel is not a perfect conductor, so this leads to consider MHD equilibria on a fixed plasma region Ω (the fixed boundary problem), and look for an external current J ext in the vacuum whose corresponding magnetic field B ext forces the plasma to be in equilibrium on the whole space (the reason for the boundary conditions (1.2e) and (1.2f)).
The main result of this article is the existence of piecewise smooth MHD equilibria with nonconstant stepped-pressure in a wide range of toroidal domains, which can be very different from an axisymmetric domain.The same philosophy works, with only minor modifications, for the fixed and free boundary problems.The equilibria we construct are not C 1 , like those of Bruno and Laurence for almostaxisymmetric domains, and in fact they feature singular current sheets (cf.Remark 2.4).The toroidal domains we consider can even be knotted in an arbitrarily complicated fashion.Specifically, the result applies to any toroidal domain with an analytic boundary satisfying a certain nondegeneracy assumption, which enables us to employ KAM-theoretic techniques in a certain step of the proof.
1.1.Nondegenerate toroidal domains.To define the nondegeneracy condition that our toroidal domains must satisfy, we need to introduce some notation.We recall that by toroidal domain we mean that the domain is diffeomorphic to a solid torus.Firstly, recall that Hodge theory ensures the existence of a unique (modulo a multiplicative constant) solution to the boundary problem curl h = 0 in Ω , div h = 0 in Ω , h • N = 0 on ∂Ω on a toroidal domain Ω ⊂ R 3 .We refer to h as the harmonic field of Ω.
Beltrami fields are solutions to the equation for some nonzero real constant λ.The space of Beltrami fields on a toroidal domain Ω is infinite dimensional.Specifically, as the curl defines a self-adjoint operator with dense domain in the space L 2 div,h (Ω) of square-integrable divergence-free fields that are L 2 -orthogonal to the harmonic field h, it is standard that there is an orthogonal basis of L 2 div,h (Ω) consisting of Beltrami fields with zero harmonic part, i.e., solutions {B n } ∞ n=1 to Equation (1.3) for certain nonzero constants {λ n } ∞ n=1 that satisfy the additional constraint , there is a unique Beltrami field with parameter λ and prescribed harmonic part, that is, a unique solution to the boundary problem (1.3) subject to the additional constraint All along this paper we shall use the term eigenvalue when referring to the constant proportionality factor λ of a Beltrami field, even taking into account that, strictly speaking, λ is an eigenvalue of curl only when it belongs to the aforementioned discrete set of points {λ n } ∞ n=1 .In both cases (that is, when B has a prescribed harmonic part and when B = B n for some n), the trace to the boundary of the Beltrami field is a smooth vector field tangent to the (embedded, nonsymmetric, possibly knotted) toroidal surface ∂Ω.In any case, ∂Ω is an invariant torus of B. Now recall that this invariant torus is Diophantine with frequency vector ω ∈ R 2 if there exist global coordinates ϕ : ∂Ω → T 2 , with T := R/2πZ, such that the restriction of the field B to ∂Ω reads in these coordinates as (1.4) and ω is a Diophantine vector.This means that there exist constants γ > 0 and τ > 1 such that for any vector with integral components k ∈ Z 2 \{0}.It is well known that the Diophantine property (possibly with different γ) of the frequency vector ω is independent of the choice of coordinates ϕ.We recall that {∂ ϕ1 , ∂ ϕ2 } stands for the holonomic (or coordinate) basis associated to the coordinate system ϕ : ∂Ω → T 2 ; as usual, the notation ∂ ϕj means that the coordinate basis vector field can be identified with the partial derivative operator (under the usual interpretation of vector fields as differential operators acting on functions).
In this paper, we will say that a toroidal domain Ω ⊂ R 3 is nondegenerate of type I or II if there is a Beltrami field on Ω for which the boundary ∂Ω is a Diophantine invariant torus and if the determinant of certain 2 × 2 constant matrices is not zero (type I) or not equal to a certain constant depending on the Beltrami field (type II).To streamline the Introduction, the expression of these matrices (which are the average on T 2 of matrix-valued quantities involving the specific Beltrami field, the associated Diophantine frequency vector and the linearizing coordinates ϕ) is relegated to Definitions 3.1 and 4.1 in the main text.To get some intuition about the meaning of this condition, recall that Beltrami fields appear in the plasma physics literature as force-free fields with constant factor, so the nondegeneracy condition can be heuristically understood as the existence of a generic force-free field on the domain whose rotational transform is irrational on the boundary.For concreteness, we shall refer to a Beltrami field with this property as a nondegenerate Beltrami field of type I or II .Some observations are in order.Firstly, since there are infinitely many curl eigenfields that do not necessarily vanish on the toroidal surface ∂Ω, and since the set of Diophantine vectors has full measure, it is natural to conjecture that a "generic" toroidal domain should be nondegenerate, in this sense.However, genericity questions for vector-valued eigenvalue problems are remarkably hard [30,14] and we have not been able to rigorously establish this claim.A particular case that one can study in detail is the class of thin toroidal domains, which one can understand as thickenings of an arbitrary closed curve in space.There one has a good understanding on the structure of harmonic fields, which can be used to rigorously show that thin toroidal domains are indeed nondegenerate.Details are given in Proposition 7.1.A concrete consequence of this fact is that there certainly are nondegenerate toroidal domains of any knot type, which are obviously not small perturbations of an axisymmetric domain.The boundary can be chosen to be smooth, and even analytic.
1.2.Statement of the main results.We are now ready to present our main result about the existence of MHD equilibria in nondegenerate toroidal domains.
In the context of the fixed boundary problem, the result can be stated as follows: Theorem 1.1 (Fixed boundary MHD equilibria).Let Ω 1 ⊂ R 3 be a nondegenerate toroidal domain of type I with analytic boundary, and let B 1 be a nondegenerate Beltrami field of type I with eigenvalue λ 1 on Ω 1 , in the sense of Definition 3.1.Then, for any N > 1 and almost all (λ 2 , . . ., λ N ) ∈ R N −1 with λ j = λ k , the following statements hold: (i) There exists a collection of "nested" nondegenerate toroidal domains of type I {Ω k } N k=2 with analytic boundary, all of them diffeomorphic to Ω 1 , with the property that Ω k−1 ⊂ Ω k for all 2 k N (see Figure 1).(ii) There is a piecewise smooth MHD equilibrium (B, P ) in the fixed domain Ω N , which satisfies Equations (1.1) in the weak sense.In particular, P + 1 2 |B| 2 is continuous across the surfaces ∂Ω k for 1 k N − 1. (iii) For each 1 k N , the magnetic field and the pressure satisfy Here {p k } N −1 k=1 are distinct nonzero constants, p N := 0 and we have set Remark 1.2.The full measure set of values (λ 2 , . . ., λ N ) for which the theorem holds is described as the whole R N −1 minus a finite set of hyperplanes; the values of λ j that must be avoided are given by Equation (3.5), and depend on the nondegenerate Beltrami field B 1 and the pressure jumps in a very nontrivial way.
Likewise, for the free boundary problem, our main result can be stated as follows.
In the statement, since B is a Beltrami field, it solves Equation (1.2a) in Ω for any constant P , whose value is inessential and can be set to be 0. Since the plasma pressure of the vacuum is 0 as well, this yields free boundary force-free equilibria.J dS, where dS is the surface measure on the boundary of an analytic domain Ω that is diffeomorphic to Ω and encloses Ω, and where J is an analytic tangent vector field on ∂Ω , such that (B, B ext , J ext , Ω) is a solution of the free boundary problem (1.2) with P = 0 in Ω.
Remark 1.4.The external surface current J ext can be approximated (in the sense of distributions) by a current Jext whose support consists of finitely many closed curves on ∂Ω , see e.g.[17].Then, Jext represents the current produced by a number of external coils (and hence it cannot be longer represented as J dS for some smooth tangent vector field on ∂Ω ).However, (B, Bext , Jext , Ω) is not a solution of the free boundary problem (1.2) because the boundary terms Bext • N and | Bext | 2 − |B| 2 on ∂Ω are small but not necessarily zero.
Remark 1.5.A minor modification of the proof of Theorem 1.3 using Theorem 5.1, allows us to construct free boundary MHD equilibria with constant pressure P = c = 0 in Ω, provided that |c| is small enough.In this case, the boundary condition (1.2f) is replaced by the identity |B ext | 2 − |B| 2 = 2c on ∂Ω, and Ω needs to be a nondegenerate toroidal domain of both types I and II.Since the pressure of the vacuum region is 0 and c = 0, this modified construction yields stepped-pressure free boundary equilibria in R 3 .
A consequence of these theorems and of the above discussion about nondegenerate domains is the existence of piecewise smooth MHD equilibria with nonconstant pressure and (fixed or free) toroidal boundaries of any knot type.A precise statement is given in Corollary 7.3.For the benefit of the reader, in Section 2 we recall the definition of weak solutions to the system (1.1a)-(1.1c),which is required to make sense of MHD equilibria that are only piecewise smooth.
It is worth mentioning that a minor modification of the proof of Theorem 1.1 permits to prove the existence of Lipschitz continuous force-free fields with a nonconstant proportionality factor on toroidal domains of complicated geometry.Details are provided in Theorem 8.1.This is interesting because, in a certain precise sense, smooth force-free fields with a nonconstant factor are rare, as discussed in [27,15].1.3.Strategy of the proof.The strategy of the proofs of Theorems 1.1 and 1.3 is similar, so let us focus on the former.The basic idea behind Theorem 1.1 is motivated by the work of Bruno and Laurence [8] on MHD equilibria on small perturbations of an axisymmetric toroidal domain.The perturbative construction they use in their proof, however, hinges strongly on having approximately axisymmetric solutions, where one can obtain very precise information about the solutions and their trajectories, and cannot be extended to deal with toroidal domains that are not approximately symmetric.
To explain the gist of our approach, let us stick to the simplest case, N = 2.The case of an arbitrary N 2 is obtained by repeating the process N − 1 times.Our initial observation (Lemma 2.1) is that, if we have two Beltrami fields B 1 , B 2 defined on two disjoint domains Ω 1 , Ω 2 := Ω 2 \Ω 1 , with Ω 1 ⊂ Ω 2 , one can define a piecewise smooth MHD equilibrium on the domain Ω 2 , with a certain piecewise constant pressure function P , provided that the difference We start by choosing B 1 as a nondegenerate Beltrami field in the toroidal domain Ω 1 , so that the analytic surface ∂Ω 1 is a Diophantine invariant torus of B 1 .To construct a Beltrami field B 2 in an exterior neighborhood of ∂Ω 1 , we use a version of the Cauchy-Kovalevskaya theorem for the curl operator [13] (see also Appendix B) with a Cauchy datum given by an analytic vector field tangent to ∂Ω 1 .A key aspect of this theorem is that one can only grant the existence of a local solution to the equation provided that the Cauchy datum satisfies an additional constraint.When one takes this constraint into account, showing that |B 1 | 2 − |B 2 | 2 is constant on ∂Ω 1 becomes equivalent to proving the existence of an analytic solution to a certain nonlinear Hamilton-Jacobi equation on T 2 .
The key difficulty of the problem is that, as the toroidal domains we consider are far from the axisymmetric case, we cannot extract from the equations enough information about the trajectories of the vector fields.The first manifestation of this difficulty is that we have not found a way of effectively using trajectories to analyze the aforementioned Hamilton-Jacobi equation.Instead, we have shown that one can exploit the fact that the restriction B 1 | ∂Ω1 is conjugate to a Diophantine rotation to regard the equation as a nonlinear perturbation of the cohomological equation which appears in KAM theory.With this approach, we eventually establish the existence of analytic solutions by means of a Newton quadratic scheme (Theorem 5.1).The next step is to show that the resulting field B 2 does in fact have an invariant torus enclosing a toroidal domain Ω 2 ⊃ Ω 1 , which permits to make sense of the basic geometric configuration used to construct the MHD equilibrium.To this end, we prove that ∂Ω 1 is a twist (in the KAM theoretic sense) invariant torus of B 2 , so that it is accumulated by a set of Diophantine analytic invariant tori.However, the key difficulty is that we cannot compute a good approximation for the trajectories of B 2 .This means that we do not have enough information to apply the existing KAM theorems for divergence-free fields (see e.g.[26,23,16,24]), which are based on studying the Poincaré map of the field on a transverse section.
To solve this problem, we establish a KAM theorem for divergence-free vector fields in R 3 with two key features that make it rather different from other KAM theorems in the same context [6,7,23].First, it applies to vector fields which do not need to be approximately integrable or in Birkhoff normal form.Secondly, the twist condition is written solely in terms of the vector field and of the approximate invariant torus.An additional advantage is that the formulas take a particularly simple form when the field is Beltrami.Recall that a KAM theorem for perturbations of integrable volume-preserving diffeomorphisms was obtained in [9,31,32].
1.4.Organization of the paper.After recalling the definition of weak MHD equilibria, in Section 2 we prove a lemma ensuring that one can construct piecewise smooth MHD equilibria by gluing two Beltrami fields defined on non-intersecting domains with a common boundary component, provided that the boundary traces of these Beltrami fields satisfy a certain constraint.The main arguments of the proofs of Theorems 1.1 and 1.3 are presented in Sections 3 and 4. For clarity of exposition, however, the two essential technical points of proof (which are of independent interest) are relegated to Sections 5 and 6.Specifically, in Section 5 we solve, using a cohomological equation, the Hamilton-Jacobi equation associated with the constraint that we came across in Section 2. Also, in Section 6 we prove our new KAM theorem for divergence-free vector fields in R 3 .Section 7 is devoted to rigorously proving that thin toroidal domains of any topology are generically nondegenerate (of type I and II).The existence result for Lipschitz-continuous forcefree fields with a nonconstant factor is presented in Section 8.The paper concludes with two technical appendices.In the first appendix we show that Beltrami fields are analytic up to the boundary if the domain is analytic, and in the second we record certain results for Beltrami fields that we proved in [13,16,12] and which are relevant for the problem under consideration.

Construction of weak MHD equilibria from Beltrami fields
In this section we introduce the definition of a weak MHD equilibrium.We say that a pair (B, P ) of class, say, L 2 (Ω) is a weak solution to the stationary MHD equations in Ω if for any vector field w ∈ C 1 c (Ω) and any scalar function φ ∈ C 1 (Ω).Of course, if B and P are continuously differentiable, this is equivalent to saying that they satisfy Equations (1.1a)-(1.1c) in Ω.As usual, the subscript c means that we take functions with compact support.Lemma 2.1.Let {Ω k } N k=1 be N 2 bounded domains in R 3 with smooth connected boundaries.Assume that these domains are nested in the sense that is a piecewise smooth MHD equilibrium on Ω N with piecewise constant pressure Here c 0 is any real constant (in particular, it can be chosen so that Remark 2.2.As B k is a Beltrami field defined in smooth domains and tangent to the boundary, it is standard [11] that B k is smooth up to the boundary.Therefore, the constraint (2.1) makes sense pointwise.A related analytic regularity result up to the boundary, which will be needed later on, is proved in Appendix A.
Remark 2.3.The result and the proof remain valid when λ k = 0 if the corresponding vector field B k is additionally assumed to be divergence-free in Ω k .
Proof.To keep the formulas as simple as possible, we will prove the result for N = 2; the general case is analogous.We start by noticing that, for all where we have used that div B 1 = div B 2 = 0 in their respective domains and Hence div B = 0 in the sense of distributions.
Let us now take an arbitrary vector field w ∈ C 1 c (Ω 2 ).We can write Integrating by parts, and using that B 1 • N = 0 on ∂Ω 1 , we easily obtain To pass to the last equation we have used the well known identity for Beltrami fields Analogously, using the same identity for B 2 on Ω 2 and that B 2 • N = 0 on the boundary, we can compute the term I 2 .Notice that by the connectedness of the boundaries of Ω 1 and Ω 2 we have ∂Ω 2 = ∂Ω 1 ∪ ∂Ω 2 , and that the outward pointing normal vector of Ω 2 on ∂Ω 1 is −N .We thus obtain The surface integral is taken only on ∂Ω 1 because w = 0 on ∂Ω 2 .
Putting together these computations and using the boundary condition (2.1), we finally conclude that . It then follows that (B, P ) is a weak solution of the MHD equations in Ω 2 , as claimed.
Remark 2.4.It is easy to check that the plasma current J := curl B of the solution constructed in Lemma 2.1 is the vector-valued distribution where dS k and N k are the area measure and outer unit normal on ∂Ω k .The current sheet terms appearing in this formula are a consequence of the discontinuity of the magnetic field across the surfaces ∂Ω k .

Fixed boundary equilibria: proof of Theorem 1.1
In this section we show how to implement the strategy discussed in the Introduction to prove the first main result of this paper, modulo some technical results that will be presented in later sections.
Here ω is the frequency vector of B on ∂Ω, ϕ are the linearizing coordinates in Equation (1.4), G is the metric matrix (or first fundamental form) of the surface ∂Ω in the coordinates ϕ, and R(ϕ) is the unique zero mean solution to the equation on the torus where X := B| ∂Ω and the constant κ is chosen so that T 2 (κ − |X| 2 ) dϕ = 0.By rescaling the Beltrami field, we can henceforth assume that the above constant is κ = 1 (which has already been considered when writing Equation (3.1)).
The Beltrami field and the analytic toroidal domain satisfying the nondegeneracy assumption are B 1 and Ω 1 , and the corresponding eigenvalue is λ 1 .Clearly, the nondegeneracy of B 1 is invariant under rescaling the field by a nonzero constant.

Construction of the first layer.
As ∂Ω 1 is an invariant torus of B 1 with Diophantine frequency vector ω (1) , it is standard [25] that one can parametrize the invariant torus by an embedding K 1 : T 2 → R 3 satisfying the equation is the pointwise derivative of K 1 in the direction of ω (1) .In this picture, of course, the invariant torus is the image ).Let us emphasize from the beginning that parametrizing the invariant tori by embeddings is essential for the KAM theorem that we will prove in Section 6 (Theorem 6.2), which will play a key role in this proof later on.
Since the boundary of the domain Ω 1 is analytic, Theorem A.1 implies that B 1 is analytic up to the boundary.A theorem of Herman [32,33] then ensures that the linearizing parametrization K 1 is analytic, or in other words, there are analytic coordinates ϕ : Here G K1 := DK 1 DK 1 is the matrix representation of the pullback of the Euclidean metric to the surface ∂Ω 1 obtained using the embedding K 1 , so it is a positive definite 2 × 2 symmetric (nonconstant) matrix.Theorem 5.1 implies that, for any constant c 2 that is small enough in absolute value, there is an analytic vector field X 2 on ∂Ω 1 such that: (i) The 1-form that is dual to X 2 with respect to the metric on ∂Ω 1 induced by the Euclidean metric in R 3 , is closed.(ii) The pointwise norm of X 2 , computed with the induced metric on ∂Ω 1 , satisfies the equation for some constant bounded as |b 2 | C|c 2 |. (iii) The vector field X 2 depends continuously on the parameter c 2 , in the C r -topology of vector fields for any r, and is Diophantine with frequency vector 1) .
In particular, one can write where in what follows O(c 2 ) stands for a quantity (which may vary from line to line) whose C r norm is bounded by C|c 2 |, for any fixed integer r.We remark that the fact that the frequency vectors ω (2) and ω (1) are proportional, cf.Equation (3.2), is a consequence of the proof of Theorem 5.1.It is not a necessary condition, but it is convenient to write all the estimates of the proof of Lemma 5.4 in a simple way, and to keep the same Diophantine constants.

Now we consider the Cauchy problem
curl for some nonzero constant λ 2 = λ 1 that we will fix later.Since ∂Ω 1 and X 2 are analytic, and the 1-form dual to X 2 is closed, Theorem B.1 (which is a sort of Cauchy-Kovalevskaya theorem for the curl operator) implies that there exists a unique analytic solution to this Cauchy problem in a neighborhood of ∂Ω 1 .Eventually, we will only be interested in the behavior of the solution outside Ω 1 .
By construction, ∂Ω 1 is a Diophantine invariant torus of the vector field B 2 .We claim that it is twist (in the KAM sense; see Definition 6.1) for almost all choices of the constant λ 2 .This implies that ∂Ω 1 is accumulated by a set of Diophantine invariant tori of B 2 contained in R 3 \Ω 1 .Since B 2 is analytic, these tori are analytic as well.This follows from Corollary 6.3 to the KAM theorem for divergence-free fields that we shall prove in Section 6.
Let us denote by K 2 : T 2 → R 3 an embedding that is a linearizing parametrization of the invariant torus ∂Ω 1 with frequency vector ω (2) of the vector field B 2 .Then, we can introduce coordinates (which we still denote by ϕ) such that X 2 = B 2 | ∂Ω1 becomes the Diophantine linear field 2 ∂ ϕ2 .In general, K 2 is different from the parametrization K 1 that linearizes X 1 , but it follows from the previous discussion that both parametrizations differ by a higher order correction, i.e., (3.3) Since B 2 satisfies the equation curl B 2 = λ 2 B 2 , an easy computation shows that the following identity holds: , where DB 2 is the Jacobian matrix of B 2 and × denotes the vector product, both computed in Cartesian coordinates.
Proof.The proof is straightforward.Indeed, To invoke Corollary 6.3, we must check the twist condition (cf.Definition 6.1).This involves computing a two-component vector field (or 2 × 1 matrix) appearing in Equation (6.3), which in this case takes the form Here DB 2 and B 2 are evaluated at K 2 (ϕ) and the normal vector is defined in terms of K 2 as in Definition 3.1.Observe that DK 2 is a 3 × 2 matrix.
Since the vector field B 2 × n 2 is tangent to ∂Ω 1 and perpendicular to B 2 , we infer that there is a nonvanishing vector (given by a 2 × 1 matrix) where we have used Lemma 6.4 to pass to the second equality.It is clear from this expression that the vector A 2 (ϕ) only depends on the way the torus ∂Ω 1 is embedded in R 3 , on the Diophantine vector ω (2) , on the parametrization K 2 linearizing X 2 , and on the eigenvalue λ 2 .
According to Definition 6.1, the invariant torus ∂Ω 1 of B 2 is twist if the twist constant is nonzero.Here F 2 (ϕ) is the function defined as the projection of the field α 2 in the (ω (2) ) ⊥ direction, i.e., This function is nonvanishing because the field DK 2 α 2 is perpendicular to B 2 | ∂Ω1 , so α 2 and ω (2) cannot be proportional at some point of T 2 .
Arguing in the same way, we obtain an analogous expression for the twist constant T 1 of the invariant torus ∂Ω 1 for the vector field B 1 .As ∂Ω 1 is an invariant torus for B 1 and B 2 , and as the corresponding parametrizations K 1 and K 2 and Diophantine vectors ω (1) and ω (2) differ just by an error of order c 2 by Equations (3.2)-(3.3),we conclude that where This shows that ∂Ω 1 is a twist invariant torus of B 2 for almost all choices of λ 2 .
Hence, we can take a Diophantine analytic invariant torus Σ 2 of the vector field B 2 , lying outside Ω 1 , which is η-close to ∂Ω 1 .By this we mean that, for any fixed r, there is a diffeomorphism Ψ 1 of R 3 which maps ∂Ω 1 into Σ 2 and which is close to the identity as Ψ 1 − id C r < η.The invariant torus Σ 2 is then the boundary of a toroidal domain Ω 2 ⊃ Ω 1 .It is easy to check that the matrix M 2 in Equation (3.1), associated to the vector field B 2 | ∂Ω2 is related to the matrix M 1 of B 1 | ∂Ω1 as As M 1 is invertible and ∂Ω 1 is accumulated by Diophantine invariant tori of B 2 , we can therefore take η (and |c 2 |) small enough so that M 2 is invertible too.We then conclude that Ω 2 is a nondegenerate toroidal domain of type I.

Conclusion of the proof.
As Ω 2 is another nondegenerate toroidal domain of type I, we can repeat the argument to construct a vector field B 3 in a neighborhood of ∂Ω 2 that solves the equation for some constant λ 3 = λ 2 , and the Cauchy datum X 3 satisfies with arbitrarily small constants c 3 and b 3 = O(c 3 ).Here X 2 := B 2 | ∂Ω2 .Again, one can pick c 3 and λ 3 so that ∂Ω 2 is a twist Diophantine invariant torus of B 3 .Therefore there is an analytic Diophantine invariant torus of B 3 , which is the boundary of another nondegenerate toroidal domain of type I Ω 3 ⊃ Ω 2 .
This process can be iterated N − 1 times to obtain a family of (analytic) nested tubes {Ω k } N k=1 , different constants λ k , small constants c k , b k and vector fields To construct a weak solution (B, P ) of the MHD equations in the toroidal domain Ω N , we set The constant p 1 is arbitrary, and the constants p k are defined in terms of c j , b j as Note that, generically, p k = p j if k = j.A straightforward application of Lemma 2.1 shows that (B, P ) is a piecewise smooth MHD equilibrium with all the properties stated in Theorem 1.1.
Remark 3.3.In view of potential applications of our method of proof in the context of plasma physics, we want to highlight that the toroidal surfaces ∂Ω 2 , . . ., ∂Ω N we construct in the proof of Theorem 1.1 are not too distorted from one to the other, and their separation is small.Quantifying such distortions is feasible (but not trivial) using quantitative versions of the Cauchy-Kovalevskaya and the KAM theorems (see also Remark 4.2).

Free boundary equilibria: proof of Theorem 1.3
We first introduce the class of toroidal domains that we consider in Theorem 1.3.
Definition 4.1.A toroidal domain Ω ⊂ R 3 is nondegenerate of type II if there exists a Beltrami field B on Ω (with eigenvalue λ) such that the following conditions hold: (i) ∂Ω is an invariant torus of the field with Diophantine frequency vector ω.
(ii) The twist constant T (see Definition 6.1) satisfies is the linearizing embedding of ∂Ω in coordinates ϕ (so, in particular, ∂Ω = K(T 2 )), n := ∂ ϕ1 K ×∂ ϕ2 K is a normal vector and the R 2 -valued function α(ϕ) is defined as The nondegeneracy of B is invariant under rescaling the field by a nonzero constant.
The proof of Theorem 1.3 follows the same strategy as the proof of Theorem 1.1, although it is easier because it does not make use of the Hamilton-Jacobi equation we study in Section 5. Therefore, to make the presentation easier, we will use the same notation as in the previous section without further mention.
As the analytic toroidal domain is nondegenerate in the sense of Definition 4.1, there exists a Beltrami field B in Ω, satisfying the equation curl B = λB for some nonzero constant λ and the boundary condition B • N = 0, such that ∂Ω is a Diophantine invariant torus of B with frequency vector ω and its twist constant satisfies (4.1) We focus on the special case of a globally force-free equilibrium with P = 0, whose magnetic field B is a Beltrami field.In Remark 1.5 we explain how to modify this construction to obtain more general equilibria.It is then clear that (B, B ext , Ω) are a solution to the equations for a free boundary MHD equilibrium with external current J ext if the external magnetic field and the external current satisfy curl To ensure that J ext is a current sheet, we also aim to construct a toroidal domain Ω ⊃ Ω and a tangent vector field J on ∂Ω such that J ext = J dS , (4.2e) where dS is the surface measure on ∂Ω .Note that the tangent vector field J must be divergence-free with respect to the induced metric on ∂Ω because Equation (4.2a) implies that, in the sense of distributions, 0 = div J ext = (div ∂Ω J) dS .
Thus proving Theorem 1.3 boils down to constructing a domain Ω and an analytic divergence-free tangent vector field J on ∂Ω such that the solution to the exterior div-curl problem (4.2a)-(4.2c) on R 3 \Ω satisfies the conditions (4.2d)-(4.2f).
To construct solutions to this overdetermined system, we follow the same philosophy as in Section 3. Let X := B| ∂Ω be the restriction of the Beltrami field B on the boundary of the domain Ω.Observe that X is analytic in view of Theorem A.1 and the associated 1-form X is closed on ∂Ω by Theorem B.1.Therefore, the Cauchy problem has a unique analytic solution in a small tubular neighborhood U of ∂Ω as a consequence of Theorem B.1.By construction, h • N = 0 and |h| 2 = |B| 2 on ∂Ω.Lemma 2.1 then ensures that the field with U := U \Ω, is a weak solution to the stationary MHD equations in the toroidal domain Ω 2 := U ∪ Ω.
Proceeding just as in Subsection 3.2, the twist constant T h of the invariant torus ∂Ω of the harmonic field h can be readily shown to be , where T is the twist constant of B (cf. Definition 6.1).This is simply Equation (3.4), where we have set λ 2 = 0 because the field h is harmonic and c 2 = 0 because The nondegeneracy assumption of type II, i.e., Equation (4.1), ensures that T h = 0. Thus Corollary 6.3 implies that ∂Ω is accumulated (in both components of its complement) by analytic Diophantine invariant tori of h.We can therefore choose an analytic domain Ω ⊃ Ω whose boundary is one of these invariant tori.
To conclude, let us now define the vector field for x ∈ R 3 \Ω.As h is divergence free in Ω \Ω and tangent to ∂Ω and ∂Ω , an elementary computation shows that div B ext = 0 in the sense of distributions.Furthermore, the corresponding current is readily computed using that and N is the outer unit normal on ∂Ω .Therefore, (B ext , J, Ω ) satisfies the system (4.2),so Theorem 1.3 follows.
Remark 4.2.Quantitative versions of the Cauchy-Kovalevskaya theorem (Theorem B.1) and the KAM theorem (Theorem 6.2) provide an estimate for the separation that we can obtain between the current sheet ∂Ω and the domain Ω in terms of the Diophantine constants of the frequency vector ω and of the analyticity radii and the analytic norms of the different objects that appear in the construction (namely, the Beltrami field B| ∂Ω and the linearizing embedding K of the invariant torus ∂Ω).However, this is a hard job that will be the object of future contributions.

Solving a Hamilton-Jacobi problem via the cohomological equation
Let Ω be an analytic nondegenerate toroidal domain of type I in R 3 .By definition, there exists a Beltrami field B in Ω that satisfies the equation curl B = λB for some constant λ, ∂Ω is a Diophantine invariant torus of B, and the corresponding matrix M defined in Equation (3.1) is invertible.Arguing as in the beginning of Section 3.2 we infer that B is analytic up to the boundary and there are analytic coordinates ϕ : ∂Ω → T 2 such that where ω ∈ R 2 is a Diophantine frequency vector.With a slight abuse of notation, in this section we will use the same name for a quantity on ∂Ω (a function or a vector field) and for its expression in these coordinates.
In this section, if Z 1 , Z 2 are two vector fields on T 2 , |Z 1 | and Z 1 • Z 2 denote the norm of Z 1 and the scalar product of Z 1 and Z 2 , respectively, computed with respect to the metric on ∂Ω induced by the Euclidean metric, which we write in the coordinates ϕ.As before, [f ] T 2 will denote the mean of a function f in ϕcoordinates, i.e., The main result of this section is the following theorem.This result is used in the proof of Theorem 1.1 to construct Cauchy data which satisfy the constraint equation (2.1) appearing in the Cauchy-Kovalevskaya theorem for the curl operator.
Theorem 5.1.Let Ω and Y be as above.Then, for any constant c with small enough absolute value, there is a nonnegative constant b and an analytic vector field X on ∂Ω of the form where ∇ denotes the gradient operator on ∂Ω associated with the induced metric, such that: (ii) X is analytically conjugate to a linear field with Diophantine frequency vector ω = (1 + c) 1/2 ω.(iii) For any fixed r, the scalar analytic function H on ∂Ω and the constants b, a 1 , a 2 are bounded as Moreover, X depends continuously on the parameter c in the C r -topology of vector fields, for any fixed r.
Remark 5.2.We are working in the analytic category because we need analytic solutions to apply the Cauchy-Kovalevskaya theorem in the proof of Theorem 1.1.
In this category, we can prove Theorem 5.1 using a quadratic Newton scheme.
Using instead a Nash-Moser iteration, one can prove a completely analogous result in the C r setting, for large enough r.
Remark 5.3.The expression of X guarantees that the dual 1-form of X (computed with respect to the metric induced from R 3 ), which we denote by X , is closed: dX = 0.This condition is also required to apply the Cauchy-Kovalevskaya theorem.We also observe that the introduction of nonzero constants a 1 , a 2 in the expression of X is important to prove Lemma 5.4 below, where they are crucial to ensure the solvability of certain cohomological equations that appear in the Newton's iterative scheme.
Proof.As in Section 3, we can rescale B so that In coordinates ϕ, condition (i) for the vector field X is equivalent to picking the analytic function H : T 2 → R and the constants a = (a 1 , a 2 ) so that the equation is satisfied.Here L ω ≡ ω 1 ∂ ϕ1 + ω 2 ∂ ϕ2 is the derivative in the ω direction and we have set Throughout this proof, we use the shorthand notation To study the existence of solutions to Equation (5.2) it is convenient to define the nonlinear operator The definition of the constant b ensures that [T c (H, a)] T 2 = 0.
To study the operator T c we will make use of the Banach space Ḣρ of holomorphic functions H on the complex strip (5.3) ∆(ρ) := {ϕ : Re ϕ ∈ T 2 , |Im ϕ| < ρ} that have zero mean on {Im ϕ = 0} (i.e., [H] T 2 = 0).This space is endowed with the supremum norm H ρ := sup ϕ∈∆(ρ) |H(ϕ)|.We will also denote by Ḣρ an analogous space of vector or matrix-valued functions.With some abuse of notation, we still use the notation ϕ for the complexification of the toroidal coordinates ϕ, and the same name for a function (vector or matrix) on T 2 and its complexification on ∆(ρ).Since the induced metric on ∂Ω is analytic and Y is also an analytic vector field, it then follows that there is some ρ 0 > 0 such that T c defines a map T c : Ḣρ → Ḣρ for all ρ < ρ 0 .
To solve the equation (5.4) T c (H, a) = 0 , we will crucially use the additional requirement that X is analytically conjugate to the linear field (1 + c) 1/2 ω.More precisely, let us consider the equation where Φ(ϕ) := ϕ + v(ϕ) is a diffeomorphism of T 2 and v : T 2 → R 2 .We denote the LHS of this equation by R c (H, v, a).
Our goal is to find analytic solutions (H, v, a) to Equations (5.4) and (5.5) when |c| is small enough.Notice that (5.5) automatically guarantees that condition (ii) is satisfied.

We shall now use the notation of this lemma without further notice. As
ρ |c| .It is obvious that there is c 0 > 0 small enough such that the assumption (5.8) holds for all |c| < c 0 (of course, the smallness assumption on v 0 is also satisfied because v 0 ρ = 0).It remains to check the generic condition on the matrix M (0) , cf.Equation (5.15).Since Φ 0 = id, an easy computation shows that the columns of the 2 × 2 matrix M (0) are given by the vectors In terms of the positive definite symmetric matrix G describing the metric on ∂Ω in the ϕ-coordinates, it is immediate to check that Equation (5.15) where R(ϕ) is the unique zero mean solution to the equation This condition is immediately satisfied, by Definition 3.1, if B is a nondegenerate Beltrami field of type I on the toroidal domain Ω.
For any ρ ∈ (0, ρ), we can then conclude from Lemma 5.4 that there exists a unique triple (H, v, a) ∈ Ḣρ × Ḣρ × R 2 in a neighborhood of (0, 0, 0) such that Equations (5.4) and (5.5) hold, provided that |c| is small enough.It is clear that (H, v, a) depends continuously on c.
The bound (iii) then follows from the estimate (5.9) below, the usual Cauchy estimate and the obvious bound 5.1.Existence of solutions of the Hamilton-Jacobi equation.In this section we prove the basic lemma used to establish the existence of analytic solutions to Equations (5.4)- (5.5).To this end, note that Equation (5.5) reads as Here and in what follows, when the operator L ω acts on vector or matrix-valued functions, its action is understood componentwise.
To solve the system (5.7) we shall use Newton's quadratic scheme and the Rüssmann estimates for analytic solutions to cohomological equations [28].We recall that the constants γ > 0 and τ > 1 in the proof appear in the definition of the Diophantine vector ω, and that one can assume γ 1 without any loss of generality.
Lemma 5.4.Let us take c ∈ [− 1 2 , 1  2 ] and consider a triple (H 0 , v 0 , a 0 ) ∈ ( Ḣρ , Ḣρ , R 2 ).For any ρ ∈ (0, ρ), if v 0 ρ and (5.8) are smaller than a certain constant ε 0 > 0 that depends on ρ but not on c, and if the approximate solution (H 0 , a 0 , v 0 ) satisfies the generic assumption given by Equation (5.15) below, then there exists a unique solution (H, v, a) ∈ Ḣρ × Ḣρ ×R 2 to Equations (5.4) and (5.6) (or, equivalently, (5.7)) bounded as Proof.To set a quadratic Newton iteration, we introduce corrections (ξ, η, α) to the approximate solution so that is a solution to the equations modulo a quadratic error, which is bounded by CE 2 0 (precise estimates will be shown later).We also take a constant b 0 ensuring that [T c (H 0 , a 0 )] T 2 = 0, and introduce the correction where the constant β will be fixed later.
Setting E 0 H := T c (H 0 , a 0 ) and E 0 v := R c (H 0 , v 0 , a 0 ), we then obtain (ξ, η, α) as solutions to the linearized equations In Equation (5.10), the vector field X 0 is defined as and the constant β will be chosen later to ensure the solvability of the equation (that is, so that a certain zero mean condition holds).In Equation (5.11), the symbol D is used to denote the Jacobian matrix of a vector field.
Taking the pullback of Equation (5.10) with the diffeomorphism Φ 0 := id + v 0 , defining the function ξ := ξ • Φ 0 , and using that In this equation, I denotes the 2 × 2 identity matrix.We also observe that if v 0 ρ is small enough, the matrix I + Dv 0 is invertible.The second summand, which denotes the action of the vector field −2(I + Dv 0 ) −1 E 0 v (understood as a first order differential operator) on the function ξ, is in fact a quadratic term; precise estimates will be given below.Thus, we can drop this term and consider the following equation: Following Zehnder [34], to study Equation (5.11) we define a new function η as η =: (I + Dv 0 ) η .
Computing the Jacobian matrix of the equation that defines E 0 v , one obtains the identity Plugging this expression into Equation (5.11), and dropping the term (DE 0 v ) η, which is quadratic, we can write Summarizing, we have changed the original linearized system of equations by the equivalent linear cohomological Equations (5.12) and (5.13).Choosing the constant β in Equation (5.12) so that T 2 = 0 , Equation (5.12) admits a unique zero-mean solution ξ, depending on the constant vector α, which is of the form Here and the constants β 0 , β (with i = 1, 2) guarantee that all the above functions of the form Φ * 0 • • • have zero mean.This ensures that the action of the operator L −1 ω (mapping functions of zero mean to functions of zero mean) is well defined.Note, in particular, Next, let us plug the expression for ξ in Equation (5.13) and consider the 2 × 2 matrix-valued function M (0) whose columns are the vector fields The solvability of Equation (5.13) then follows if and only if one can pick a vector α ∈ R 2 such that (5.14) This linear equation has a solution if and only if the matrix M (0) satisfies the invertibility condition (5.15) det (I + Dv 0 ) −1 M (0) • (id + v 0 ) T 2 = 0 .This is the generic assumption appearing in the statement of the lemma.Note that Equation (5.15)only depends on (H 0 , v 0 , a 0 ), on the vector field Y and on the domain Ω.
We have then proved that, fixing the constants α and β as above, there is a unique solution (ξ, η) to the linearized equations (5.12)-(5.13).Now let us estimate the analytic norms of these solutions to show that this scheme is indeed quadratic, and that it can be iterated because the norms of the corrected approximate solutions are uniformly bounded.For this we use Rüssmann estimates [28]: if A ∈ Ḣρ and ω is a Diophantine vector, for each δ > 0 the cohomological equation L ω B = A has a unique solution B ∈ Ḣρ−δ that can be bounded as

The constant C is independent of δ
In what follows, let us fix a small constant 0 < δ < 1, and δ < ρ (which will measure the loss of analytic band in the iteration) and assume that E 0 < ε 0 (cf.Equation (5.8)), where ε 0 satisfies the smallness condition (5.16) By assumption, v 0 ρ < ε 0 , and hence (5.16) implies that v 0 ρ δ, which guarantees, by the Cauchy estimate, that I + Dv 0 is close to the unit matrix.
A straightforward computation using Rüssmann estimates and the Cauchy estimate for derivatives then implies Here we have used that the condition v 0 ρ δ ensures that the diffeomorphism Φ 0 is close to the identity, and hence the constant C can be taken independent of v 0 .The constant does depend on Y and Ω, though.
Analogously, the constant β in Equation (5.12) can be bounded as In these bounds, the constant C only depends on H 0 ρ , |a 0 |, Y and Ω.
In order to show that this scheme can be iterated, let us check that the norms of the corrected approximate solution (H 1 , v 1 , a 1 ) and the constant b 1 remain uniformly bounded.Indeed, if H 0 ρ + |a 0 | + |b 0 | < κ 0 for some positive constant κ 0 , we can easily derive that where in the last inequality we have used the assumptions (5.8) and (5.16).Analogously, if E 0 is small enough (cf.Equation (5.16)), The generic assumption (5.15) is then satisfied in the iteration because the difference M (1) − M (0) is bounded as Therefore, if δ is small enough, satisfies an invertibility condition by hypothesis.
To complete the proof of the proposition, we have to check that the new errors ) are quadratic with respect to the errors E 0 H and E 0 v .This follows from the fact that Cγ −6 δ −9−6τ E 2 0 , so the scheme is indeed quadratic.In particular, Cγ −6 δ −9−6τ E 2 0 .so the new error E 1 is smaller than ε 0 because of the smallness condition (5.16).
It is now standard that the scheme can therefore be iterated to yield a unique solution (H, v, a) ∈ Ḣρ × Ḣρ × R 2 to Equations (5.4) and (5.6) that is bounded as in Equation (5.9), with ρ := ρ − 8δ.This is an easy consequence of the previous estimates and the following well known lemma: for some constant C > 0, positive reals a, b, and small constant δ 0 , with The result is then proven.

A KAM theorem adapted to divergence-free vector fields
In this section we prove a KAM theorem that applies to divergence-free vector fields in R 3 that are not necessarily close to integrable ones or in Birkhoff normal form.Another key technical advantage of this result is that the twist condition is written in terms of the (approximately) invariant torus and the vector field itself, so one does not need any fine information about the trajectories of the vector field.The proof follows the parametrization method as presented in [25] in the context of Hamiltonian systems.
It is convenient to study the reducibility properties of the invariant tori using embeddings of T 2 → R 3 , denoting by ϕ the natural coordinates on the torus.As before, in this section D denotes the Jacobian matrix of a vector-valued function and | • | stands for its norm computed with the induced metric.
6.1.Statement of results.Let B be an analytic divergence-free vector field in R 3 and let ω ∈ R 2 be a frequency vector satisfying the Diophantine condition (1.5).Let K : T 2 → R 3 be an analytic embedding, which parametrizes a torus T := K(T 2 ) ⊂ R 3 that we will eventually assume to be approximately invariant under B with frequency ω.The associated error, or defect, of invariance is measured using the function where, as before, L ω K(ϕ) := DK(ϕ)ω is the derivative of K in the direction of ω.
We still denote by the matrix representation of the pull-back to T 2 of the Euclidean metric.It is obviously a nonsingular 2 × 2 matrix because K is an embedding.As K and B are analytic, they can be analytically continued to a complex strip that we call ∆(ρ) of the form (5.3).
An important ingredient of a KAM theorem is its twist condition, which in this case we define in terms of the embedding K and the vector field B as follows: Definition 6.1.The embedding K : T 2 → R 3 is twist for the vector field B if the twist constant Here the 2 × 1 matrix A is defined as |n(ϕ)| 2 and n is the vector field normal to the torus T given by (6.4) n(ϕ) The main result of this section, which is of interest in itself, is the following KAM theorem.We recall that two tori are said to be η-close if there is a diffeomorphism Ψ of R 3 mapping one into the other and such that Ψ − id C r < η, where r 4 is any fixed integer.Theorem 6.2.Let B, ω, ρ and K be as above.Assume that the embedding K is twist with respect to the vector field B. If the invariance error E ρ is small enough, there is a constant λ * and an analytic Diophantine invariant torus T * of B with frequency vector ω The following corollary, which is a standard consequence of Theorem 6.2, is the result we employed in the proof of Theorem 1.1.Corollary 6.3.Assume that the (analytic, divergence-free) vector field B has an invariant torus T with Diophantine frequency vector ω.If T is twist (in the sense of Definition 6.1), then it is accumulated (in both connected components of the complement R 3 \T ) by a set of Diophantine analytic invariant tori of B.
Proof.Simply apply Theorem 6.2 to the triple (B, ω , K), with a Diophantine frequency vector ω = ω which is very close to ω and has the same Diophantine constants (γ, τ ).This ensures the invariance error is as small as one wishes.The accumulation property of the statement follows from the fact that the set of Diophantine numbers with fixed (γ, τ ) has positive measure in any neighborhood of ω.The invariant torus T with frequency vector (1 + λ )ω that one obtains with Theorem 6.2, lies in the exterior component of R 3 \T or in the interior one depending on whether is smaller or bigger than ω2 ω1 (because of the twist condition).
To prove Theorem 6.2 we will iteratively correct the embedding and the frequency vector by means of the Newton method.Denoting the corrected quantities by K(ϕ) := K(ϕ) + ∆ K (ϕ) , one is led to choosing ∆ K (ϕ) as a solution of the linearized equation Our next goal is to analyze this equation, which will involve developing a geometric setting in which the analytic properties of the equation are laid bare.The most efficient way to do this is by first considering the (trivial) case where T is an actual invariant torus of B. This case will be considered in the next subsection.Subsequently, in Subsection 6.3 we will refine this approach to deal with the approximately invariant case, which will enable us to prove Theorem 6.2.6.2.Geometric study of the invariant case.In this subsection we shall assume that T is an invariant torus of the (divergence-free) vector field B with frequency vector ω ∈ R 2 , which we parametrize by the map K : T 2 → R 3 .The invariance equation reads as (6.6)L ω K(ϕ) = B(K(ϕ)) .
The columns of the matrix DK are obviously a basis of the tangent space T K(ϕ) T .With n(ϕ) being the normal vector (6.4), the key geometric observation is that the frame (6.7) greatly simplifies the analysis of the operator R.More precisely, in this frame (6.7), the linear operator R has a triangular structure that reduces the study of Equation (6.5) to that of two cohomological equations with constant coefficients.
To prove this, we start by computing the action of the operator R on the frame (6.7).Firstly, by taking derivatives of both sides of Equation (6.6), we obtain that L ω DK(ϕ) = DB(K(ϕ))DK(ϕ) , which implies that R(DK) = 0. To compute R( n(ϕ) |n(ϕ)| 2 ), we use the following lemma: Lemma 6.4.If Equation (6.6) is satisfied, Then the elementary identity holds if U is a 3 × 3 traceless matrix and v 1 , v 2 are 3 × 1 vectors.If we take U as DB(K(ϕ)), and note that DB has zero trace because B is divergence-free, the lemma then follows.
We are ready to compute the action of the linear operator R on the normal vector: Here we have used Lemma 6.4 to pass to the last line.It is straightforward to check that this expression can be written in the frame (6.7) as (6.9)R n(ϕ) |n(ϕ)| 2 = DK(ϕ)A(ϕ) , where the 2 × 1 matrix A(ϕ) is (6.10) Observe A(ϕ) is therefore just as in Equation ( 6.3).
6.3.Proof of Theorem 6.2.Our goal now is to use the frame introduced in the previous section to analyze Equation (6.5), which determines the small corrections ∆ K and δ ω .The problem ultimately boils down to studying the inverse of the operator R.
Taking derivatives in Equation (6.1), we obtain (6.11) Arguing as in the proof of Lemma 6.4 but using this equation instead of (6.6), we prove the following: Lemma 6.5.The normal vector n(ϕ) satisfies the equation This allows us to compute the quantity R( n(ϕ) |n(ϕ)| 2 ) as in (6.8): Here A(ϕ) is the vector (6.10) and The following lemma shows how to get an approximate solution of Equation (6.5), modulo a quadratic error, using solutions to a pair of cohomological equations: Lemma 6.6.Suppose that the functions ξ 1 , ξ 2 on T 2 satisfy the cohomological equations solves Equation (6.5) modulo a quadratic error.More precisely, Remark 6.7.We observe that, at least formally, the quantities E 1 (ϕ) and E 2 (ϕ) are quadratic in E(ϕ).This is obvious except, perhaps, for the presence of the constant δ ω in Equation (6.15).The fact that this constant is of order [E] T 2 is a consequence of the solvability condition that the average of Equation (6.15) must be zero (see the discussion after Equation (6.21)).
To solve the cohomological equations (6.15) and (6.16), we need to ensure that the RHS of these equations have zero mean.In the case of Equation (6.16), a simple computation shows that where we have used that n(ϕ) L ω K(ϕ) = 0.The last integral can be equivalently written as where dS := |n| dϕ is the induced area form on T , N := n |n| is the unit normal and Ω is the domain in R 3 bounded by T .To obtain the last equality, we have also used that B is divergence-free.We thus conclude that [n E] T 2 = 0, so there exists a unique zero mean solution ξ2 to Equation (6.16).The general solution to this equation is Plugging the expression for ξ 2 in Equation (6.15), we obtain that the average of the RHS of this equation is We now make the additional assumption that δ ω and ω are collinear, so δ ω := Λω for some Λ ∈ R.This is crucial to preserve the Diophantine properties of the frequency vector.The twist condition in Definition 6.1 then ensures that there is a unique pair ([ξ 2 ] T 2 , Λ) ∈ R 2 for which the quantity (6.21) is 0. It is immediate to check that the solution ( We thus obtain an expression for the corrected embedding K(ϕ) := K(ϕ) + ∆ K (ϕ).The invariance error associated with the corrected embedding can be easily computed: Therefore, the new error is quadratic in E(ϕ) by Lemma 6.6.
In view of the above estimates, it is now well known that, if a certain smallness assumption is satisfied, the associated quadratic scheme starting with (K 0 , Λ 0 ) := (K, 1) ∈ Ḣρ × R + converges to some (K * , Λ * ) ∈ Ḣρ × R + .The embedded torus T * := K * (T 2 ) is then a Diophantine invariant torus of B with frequency vector ω * := Λ * ω.The details are standard (see e.g.[25]), and go essentially as in Section 5.1, so we will just sketch the argument.

Thin toroidal domains are generically nondegenerate
As discussed in the Introduction, it is reasonable to expect that a "generic" toroidal domain satisfies both nondegeneracy conditions (cf.Definitions 3.1 and 4.1), but proving generic results for vectorial problems is often extremely hard.Our objective in this section is to prove an analog of this result in the class of thin toroidal domains, where one can analyze the harmonic field (and other Beltrami fields) in detail.
Given a closed smooth curve γ : T → R 3 , let us denote by Ω(γ, ε) the toroidal domain (or tube) of thickness ε defined by this curve: The main result of this section can then be stated as follows.When we say that the result holds for almost all small enough ε, it means that there exists some ε 0 > 0 and a subset Z ⊂ [0, ε 0 ] of measure zero such that the result holds for all ε ∈ (0, ε 0 ]\Z. Proposition 7.1.For a generic curve γ (in the sense of a dense subset in the C r topology, for any r 3), the toroidal domain Ω(γ, ε) is analytic and nondegenerate of type I and II for almost all small enough ε > 0.
Let us next prove condition (ii) in Definition 3.1.In the following computations we shall use the formulas and results obtained in [16] without further mention.Consider the coordinates (α, r, θ) : Ω 1 → T × (0, 1) × T that parametrize Ω 1 (we are assuming, without any loss of generality, that the length of the core curve γ is 2π); in particular, ∂Ω 1 corresponds to the surface {r = 1}.In these coordinates, Y := B| r=1 takes the form where τ is the torsion of γ, and O(ε) stands for a quantity whose C m norm is bounded by C|ε|, for any fixed integer m.Using the expression of the metric on ∂Ω 1 induced from R 3 , in coordinates (α, θ), we can readily compute |Y | 2 , which leads to Moreover, we can compute the linearizing coordinates (ϕ 1 , ϕ 2 ) in terms of (α, θ) as and hence Y takes the following form in these coordinates: This implies that the frequency vector of Y (which is Diophantine) is given by Accordingly, we infer from Equation (7.2) that the function R in Definition 3.1 solves the cohomological equation Putting together all these computations, we obtain a matrix M , cf.Equation (3.1), of the form where G ε is the matrix of the metric (which depends on ε) in ϕ-coordinates.It is obvious that M is invertible because the (inverse) metric matrix G −1 ε is positive definite, which shows that Ω 1 is nondegenerate of type I.In fact, [16,Theorem 7.8] ensures that the twist constant T of the invariant torus ∂Ω 1 of B 1 (see Definition 6.1) is where c γ is certain explicit constant that depends on the curve γ through its curvature and torsion but not on ε nor λ.This constant is nonzero for a generic curve γ [16, Lemma 7.9].
To conclude, let us establish property (ii) in Definition 4.1.Straightforward computations using the formulas for thin tubes derived in [16] show that for some ε-independent constant C.Here α is the R 2 -valued function introduced in Definition 4.1, and not an angular coordinate.One can then use (7.3) to see that Therefore, provided that the nonzero constant λ satisfies |λ| < C|c γ |ε 4 for a certain positive constant C independent of ε.
Remark 7.2.The nondegenerate Beltrami field B in the proof of Proposition 7.1 has eigenvalue λ 1 = O(ε 3 ).In particular, the ratio between the L 2 norms of the plasma current curl B L 2 and the magnetic field B L 2 , in the toroidal domain Ω(γ, ε), is small.Since the Beltrami equation is linear, we can safely normalize the magnetic field so that B L 2 = 1, thus implying that the plasma current is small, which is a desirable feature for stellarator design [4].
As any toroidal domain is isotopic to a thin tube, we infer that there are MHD equilibria of the kind described in Theorems 1.1 and 1.3 that have arbitrary topology: Corollary 7.3.There exist piecewise smooth MHD equilibria with fixed or free toroidal boundaries of arbitrary topology.More precisely, given any toroidal domain Ω 0 ⊂ R 3 , one can assume that the domains Ω k of Theorem 1.1 and Ω, Ω of Theorem 1.3 are diffeomorphic to Ω 0 .
Proof.It is an immediate consequence of Theorems 1.1 and 1.3 and of the fact that, by Proposition 7.1, a generic thin tube is nondegenerate of type I and II.

Lipschitz continuous force-free fields with nonconstant factor
A force-free field in a domain Ω is a vector field B that satisfies the equations curl B = f B in Ω , div B = 0 in Ω , B • N = 0 on ∂Ω for some scalar function f .In the context of hydrodynamics, these fields are called Beltrami fields with nonconstant factor [27,15].
It is obvious that a force-free field satisfies the MHD equations in Ω with P = 0.This is used to define force-free fields of low regularity.Specifically, one says that a vector field B ∈ L 2 (Ω) is force-free if in Ω, where the factor Proof.We use the same construction as in the proof of Theorem 1.1, but we consider the particular case where c k := 0 for all k = 2, • • • , N .This obviously implies that b k = 0 as well.The effect of this choice of constants is that the vector field B k in a neighborhood of ∂Ω k−1 is constructed using the Cauchy-Kovalevskaya theorem with Cauchy datum given by B k−1 | ∂Ω k−1 , so ∂Ω k−1 is no longer a discontinuity surface of the magnetic field, and that we do not need to use Theorem 5.1.
For almost all choices of the constant λ k , the Diophantine invariant torus ∂Ω k−1 one obtains is twist.This ensures the existence of a family of nested toroidal domains {Ω k } N k=1 as above and of a weak solution (B, P ) to the MHD equations with constant pressure P = 0.In each set Ω k \Ω k−1 , B satisfies the equation curl B = λ k B and is analytic up to the boundary; in particular, B is Lipschitz continuous on Ω N .The plasma current curl B is given by Remark 2.4, and the singular terms supported on the toroidal surfaces vanish because B is continuous.This formula shows that B is in fact a force-free field with piecewise constant factor f as above.The theorem is then proven.
Remark 8.2.The assumption that the matrix M is invertible that appears in the definition of nondegenerate toroidal domains of type I is not used to prove Theorem 8.1 because here we do not need Theorem 5.1.Thus, the theorem holds under the slightly weaker assumption on the domain Ω 1 that it admits a Beltrami field for which the boundary is a Diophantine invariant torus.
Remark 8.3.In view of Proposition 7.1, Theorem 8.1 implies the existence of Lipschitz-continuous force-free fields with nonconstant factor on thin tubes of any topology.
Therefore, P maps each differential form ψ ∈ H 1 (Ω, Λ • ) to the projection of its boundary trace to the space of H 1/2 sections of the vector bundle E over the boundary whose section at x is E x .Equivalently, one can understand P as the action of a zeroth order pseudodifferential operator (whose principal symbol is an orthogonal projector) on the usual boundary trace ψ| ∂Ω ∈ H 1/2 (∂Ω, Λ • ).
We are interested in the boundary value problem for the Dirac operator defined as Dψ − λ ψ = G in Ω , (A.1a) P ψ = 0 on Ω .(A.1b) The key observation now is that, for each x ∈ ∂Ω and each ξ ∈ T * x ∂Ω, the operator a(x, ξ) := γ(N (x))γ(ξ) ∈ L(Λ • ) maps E x onto E x and viceversa.It is then well known [2, Corollary 3.18] that (A.1) defines an elliptic boundary value problem.For the benefit of the reader, let us mention that this is a straightforward consequence of the fact that the principal symbol of the Calderón projector associated to the boundary problem (A.1) is [|ξ|I − ia(x, ξ)] /(2|ξ|).Since (A.1) is an elliptic boundary value problem, the theory of Boutet de Monvel [5] ensures that ψ is analytic up to the boundary if G ∈ C ω (Ω, Λ • ).This implies that the vector field B is analytic up to the boundary.In order to see this, let ψ be the differential 1-form dual to B, which we can regard as an element of H 1 (Ω, Λ • ).If G ∈ C ω (Ω, Λ • ) denotes the 1-form dual to F , the equations div B = 0 and curl B − λB = F read in terms of the associated forms as As is the identity, we infer that ψ satisfies Equation (A.1a) in Ω.Furthermore, the condition B • N = 0 amounts to saying that the differential form ψ (which is of

3. 1 .
Nondegenerate toroidal domains of type I. Let us begin by defining the class of toroidal domains that we consider in Theorem 1.1.Definition 3.1.A toroidal domain Ω ⊂ R 3 is nondegenerate of type I if there exists a Beltrami field B on Ω such that the following conditions hold: (i) ∂Ω is a Diophantine invariant torus of the field.(ii) The 2 × 2 constant matrix (3.1)

Ω(B 2 |B| 2
⊗ B) • ∇w − 1 div w dx = 0 and Ω B • ∇φ dx = 0 for any vector field w ∈ C 1 c (Ω) and any scalar function φ ∈ C 1 (Ω).The strategy of the proof of Theorem 1.1 can be readily adapted to show the existence of Lipschitz-continuous force-free fields with nonconstant factor on toroidal domains of any topology.More precisely, one can prove the following: Theorem 8.1.Let B 1 be a nondegenerate Beltrami field of type I with eigenvalue λ 1 on an analytic toroidal domain Ω 1 .For any N 2 and almost all distinct constants {λ k } N k=2 , there exists a family of nested analytic toroidal domains {Ω k } N k=1 as in Theorem 1.1 and a Lipschitz continuous vector field B satisfying the equation curl B = f B , div B = 0