{'fraction_non_alphanumeric': 0.04656544759633928, 'fraction_numerical': 0.011799151442897718, 'mean_word_length': 3.9508723201111016, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 25, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimal shape of stellarators for magnetic confinement fusion

We are interested in the design of stellarators, devices for the production of controlled nuclear fusion reactions alternative to tokamaks. The confinement of the plasma is entirely achieved by a helical magnetic field created by the complex arrangement of coils fed by high currents around a toroidal domain. Such coils describe a surface called"coil winding surface"(CWS). In this paper, we model the design of the CWS as a shape optimization problem, so that the cost functional reflects both optimal plasma confinement properties, through a least square discrepancy, and also manufacturability, thanks to geometrical terms involving the lateral surface or the curvature of the CWS. We completely analyze the resulting problem: on the one hand, we establish the existence of an optimal shape, prove the shape differentiability of the criterion, and provide the expression of the differential in a workable form. On the other hand, we propose a numerical method and perform simulations of optimal stellarator shapes. We discuss the efficiency of our approach with respect to the literature in this area.


Motivations: towards a shape optimization problem
Nuclear fusion is a nuclear reaction involving the use of light nuclei. In order to produce energy by nuclear fusion, high temperature plasmas 1 must be produced and confined. For these reactions to occur, the nuclei must get close to each other at very small distances. They must therefore overcome the Coulomb repulsion. This happens naturally in a plasma during collisions if the energy of the nuclei is sufficient. This is the objective of devices called tokamaks, steel magnetic confinement chambers that allow a plasma to be controlled in order to study and experiment with energy production by nuclear fusion. The magnetic confinement technique allows to maintain a sufficient temperature and density of the plasma, in an intense magnetic field. The simplest configuration for the magnetic field is the toroidal solenoid; this is the configuration found in most current experiments.
Unfortunately, the magnetic field is not uniform in general, which causes a vertical drift of the particles, in opposite directions for the ions and for the electrons. This charge separation creates a vertical electric field which, in turn, causes the particles to drift out of the torus. This phenomenon dramatically reduces the confinement. To get around this obstacle, the effect of such drifts is canceled by giving a poloidal component 2 to the magnetic field: the field lines are wound on nested toroids. Thus, the particles, following the magnetic field lines, have their vertical drift cancelled at each turn. In a tokamak, the poloidal magnetic field is created by a toroidal electric current circulating in the plasma. This current is called plasma current.
A possible alternative to correct the problems of drift of magnetically confined plasma particles in a torus is to modify the toroidal shape of the device, by breaking the axisymmetry, yielding to the concept of stellarator. A stellarator is analogous to a tokamak except that it does not use a toroidal current flowing inside the plasma to confine it. The poloidal magnetic field is generated by external coils, or by a deformation of the coils responsible for the toroidal magnetic field. This system has the advantage of not requiring plasma current and therefore of being able to operate continuously; but it comes at the cost of more complex coils (non-planar coils) and of a more important neoclassical transport [15].
The confinement of the plasma is then entirely achieved by a helical magnetic field created by the complex arrangement of coils around the torus, supplied with strong currents and called poloidal coils.
Despite the promise of very stable steady-state fusion plasmas, stellarator technology also presents significant challenges related to the complex arrangement of magnetic field coils. These magnetic field coils are particularly expensive and especially difficult to design and fabricate due to the complexity of their spatial arrangement.
In this paper, we are interested in the search for the optimal shape of stellarators, i.e., the best coil arrangement (provided that it exists) to confine the plasma. In general, two steps are considered: first, the shape of the plasma boundary is determined in order to optimize the physical properties, among which the neoclassical transport and the magnetohydrodynamic (MHD) stability. In a second step, we search for the coil shapes producing approximately the "target" plasma shape resulting from the previous step.
In this article, we focus entirely on the second step, assuming that the target magnetic field B T is known. It is then convenient to define a coil winding surface (CWS) on which the coils will be located (see Figure 1). The optimal arrangement of stellarator coils corresponds then to the determination of a closed surface (the CWS) chosen to guarantee that the magnetic field created by the coils is as close as possible to the target magnetic field B T . Of course, it is necessary to consider feasibility and manufacturability constraints. We will propose and study several relevant choices of such constraints in what follows.

State of the art and main contributions of this article
The question of determining the best location of coils around a stellarator, reformulated as an optimal surface problem, is a major issue for the construction of stellarators with efficient confinement properties. The physical and mathematical literature dedicated to plasmas is rich of references on this issue. We mention hereafter a non-exhaustive list of various important contributions around this problem. Let us first mention [19], where all the basic theoretical elements to understand the modeling of stellarator magnetic fields are gathered.
Regarding optimal design issues, let us distinguish between several optimization/optimal control approaches and modeling choices. Each discrete stellarator coil can be represented as a closed one-dimensional curve embedded in R 3 [39,40,38]. In these references, several optimization methods are tested among which the steepest descent and Newton like methods.
Another common choice consists in using the aforementioned CWS, in other words to define a closed toroidal winding surface enclosing the plasma surface on which all coils lie. Two kinds of issues related to the optimal design of stellarators can then be addressed. The simplest is to assume the CWS to be given, and to look for currents on this surface generating the desired magnetic field for confining the plasma. Indeed, in the limit of a large number of coils, a set of discrete coils can be described by a continuous current density on the CWS. Let us mention NESCOIL [23,29], where the current potential representing a surface current distribution is sought such that the normal component of the magnetic field vanishes in a least-squares sense at the plasma boundary. In the same vein, REGCOIL [21] improves the NESCOIL approach by adding a Tikhonov regularization term in the minimization functional whereas COILOPT [34] uses an explicit representation of modular coils on a toroidal winding surface. A review of such approaches can be found in [12]. Recently, a similar problem where an extra Laplace forces penalization term is taken into account has been investigated in [31].
A much more difficult problem is to determine the CWS and the density current distribution at the same time. This is expected to improve the performances of the resulting device. On the other hand, this approach requires solving a dual optimization problem, including a rather challenging surface optimization problem. This is the main purpose of this article. In the following we mention some of the many contributions on this topic and position our contribution through this literature. In [27], this problem is modeled using a cost functional written as the weighted sum of four terms: the first one is the surface-integrated-squared normal magnetic field on the desired plasma surface. The second is the opposite of the total volume enclosed by the coil-winding surface, acting to enforce the coil-plasma separation. The third one is a measure of the spectral width of the Fourier series describing the coil-winding surface. This allows to overcome the non-uniqueness of the Fourier series representation of the coil-winding surface. The last one is the L 2 norm of the current density, allowing to obtain coils with good manufacturing properties. It is important to note here, and this is related to the motivation for this paper, that the approach developed in [27] rests upon a (truncated) Fourier series parameterization of the surface equation. The authors thus compute derivatives of their cost with respect to these Fourier coefficients.
In [28], a more complex model involving a drift kinetic equation is considered and similar shape optimization issues are investigated.
In what follows, we propose a continuous approach, which does not rely on any parameterization of the surfaces involved. We use the notion of Hadamard variation and shape derivative. We rigorously analyze, in a continuous framework, the sensitivity with respect to the domain of a REGCOIL-type cost. We thus obtain intrinsic expressions with respect to any parametrization. This makes our approach flexible and the formulas obtained by using developments of the parametric equation of the surfaces in Fourier series can be adapted without any difficulty to other choices of parametrization. We also propose several choices of manufacturability terms in the cost functional and discuss their relevance.
The issues addressed in the following as well as our main contributions are summed-up hereafter: • Modeling of the problem (Section 1.4). Using the CWS concept, we propose a continuous formulation of the question of the best coil arrangement as a shape optimization problem, regardless of any surface parametrization. In particular, several choices of manufacturing constraints are proposed. They are integrated to the cost functional using a penalization/regularization term. From the mathematical point of view, the main issue comes to minimize a functional involving the trace of the solution of an elliptic partial differential equation (PDE) on a manifold, under geometrical constraints involving the distance to this manifold.
• Analysis of the shape optimization problem (Sections 2 and 3). Having in mind the determination of an efficient algorithm for finding an optimal form for the above problem, we focus mainly on two questions. The first one is dedicated to the existence of an optimal shape (Section 2.2). In this context, the developed approach is not completely standard and requires to carefully establish semicontinuity properties of the trace of the solution of the PDE on manifolds satisying a uniform regularity property. The second one concerns the establishment of optimality conditions using the notion of form derivative (Section 3.2). Here again, due to the particular nature of the PDE at stake, the classical approach cannot be used in a direct way and many adaptations are necessary. We establish a workable expression of this derivative, which is the basis of the numerical approaches developed in the next section.
• Numerical implementation (Section 4). A relevant aspect of this paper is that the study of the sensitivity of the studied criterion to a variation of the shape of the stellarator is carried out without using any parameterization of the surface to be designed. As a result, the sensitivity relations obtained at the end of the previous step are totally intrinsic with respect to any parameterization of the surface. As a consequence, we can apply the more robust "optimize then discretize" approach, instead of a "discretize then optimize" procedure as in most of the methods implemented for this application. The shape derivatives constitute the basis of a quasi-Newton optimization method that we implement by using a parametric representation of the surface in terms of Fourier series.

Notations
In what follows, the notation S is used to denote a C 1,1 toroidal surface 3 in R 3 , equipped with the Riemannian metric induced by the canonical embedding i S : S → R 3 , i.e., the scalar product between two vectors v and w tangent to S at a common point is equal to v, w , where ·, · denotes the Euclidean scalar product in R 3 . We also denote by µ S the associated Riemannian volume form, which coincides with the two-dimensional Hausdorff measure on S ([11, Theorem 2.10.10] and [10, Theorem 2 in Section 2.2]). We write V to denote the bounded domain of R 3 such that S = ∂V . Throughout this article, we use the following notation: • For n ∈ N * = {1, 2, . . . } and any integer m ≥ n, H n denotes the n-dimensional Hausdorff measure in R m ; • P denotes a smooth toroidal domain 4 of R 3 standing for the plasma domain; • X(S) denotes the set of smooth tangent vector fields on S; • F S = L 2 (Γ(T S)) denotes the completion of X(S) for the inner product • × denotes the cross product in R 3 ; • given a function F : R n1 → R n2 and x ∈ R n1 , DF (x) denotes the n 1 × n 2 Jacobian matrix of F . In the case where n 1 = n 2 , |DF | stands for the absolute value of the determinant of F . The symbol D x is used to denote the Jacobian operator with respect to the (vector) variable x; • ∇ S denotes the tangential gradient to S in R 3 , defined for every differentiable function f : R 3 → R by ∇ S f = ∇f − ∇f, ν ν on S, where ν stands for the outward normal vector to V . Similarly, the notation div S stands for the tangential divergence given by div S θ = div θ − Dθν, ν on S, where θ is a vector field on R 3 ; • The norm on X(S) induced by the inner product ·, · F S is denoted · F S ; • F 0 S is the closure under the norm · F S of the (tangential) divergence-free vectors of X(S); • The flat two-dimensional torus is denoted by T = (R/Z) 2 . X(T ), F T and F 0 T are defined similarly to what has been done above; • The Hausdorff distance d V and the signed distance b V from V are defined as: 3 By toroidal surface, we mean here the range of the toroidal solenoid by a homeomorphism. In what follows, we will rather consider smooth toroidal surfaces, where the wording "smooth" refers to at least C 1,1 regularity. 4 toroidal domain stands for any three-dimensional domain whose boundary is a toroidal surface • The reach of V [11] is given by More explanations are provided in Appendix B. For more exhaustive informations about this notion, we refer to [11] and [7, Sect. 6.6]; • Given a differentiable vector field θ : R 3 → R 3 , we denote by e(θ) the symmetric part of the Jacobian matrix Dθ, that is, • For two Banach spaces E and F , we denote by L(E, F ) the Banach space of continuous linear maps from E to F and by L(E) the Banach space of continuous endomorphisms; • the adjoint of a linear operator L is denoted by L † ; • If A and B denote two matrices in M 3 (R), we define their doubly contracted product as • I 3 denotes the identity matrix in R 3 .

Modeling: towards a shape optimization problem
Since we are interested in solving a shape optimization problem whose unknown is the coil winding surface S, we are led to make some assumptions on S motivated by the application under consideration. In particular, we assume in what follows that the distance d(S, P ) between S and the plasma domain P is uniformly bounded from below, namely, we fix δ > 0 and we require that We now introduce the main operator we will deal with, which plays a crucial role in electromagnetism: the so-called Biot and Savart operator. This operator associates with each current distribution on S the corresponding magnetic field in P . It can be considered as a kind of inverse of the curl operator.
Definition 1 (The Biot and Savart operator BS S [9]). Let S be a smooth two-dimensional manifold and X belong to F S . Let δ S denote the single layer distribution supported on S defined by Let u denote the unique distributional solution of the PDE that falls off at infinity, i.e., Then the Biot and Savart operator is defined as the map BS S : F S −→ L 2 (P, R 3 ) associating with X the restriction of u to the plasma domain P . By introducing the kernel K given by Remark 1. According to (H dist,P,δ ), the restriction of K to S × P is uniformly bounded. By standard regularity results for parameterized integrals, the mapping P y → BS S (X)(y) is smooth and the operator BS S , seen as going from F S to C k (P, R 3 ) with k ∈ N ∪ {+∞}, is continuous. As a consequence, the operator In what follows, we will use several times that for every x, y, h ∈ R 3 with x = y we have Computation of the optimal current j. In view of modeling the optimal design problem we will deal with, let us now introduce a target magnetic field B T ∈ L 2 (P, R 3 ). The target magnetic field B T being given, we model the optimal design of a stellarator problem as a kind of regularized least square problem, where one aims at determining both the current j and the manifold shape S leading to the magnetic field closest to B T on S. To this aim, and according to the REGCOIL procedure [21], we introduce the shape functional C defined, for every closed smooth two-dimensional manifold S, as where λ > 0 denotes a regularization parameter.
The shape optimization problem. To state the shape optimization problem that we will consider, let us first define the set of admissible manifolds. We gather hereafter several conditions evoked previously that we will take into account in the search of the CWS.
• Topology and uniform boundedness. To preserve the topology of the device (see Footnote 3), we will only consider CWSs that are two-dimensional closed toroidal manifolds. Moreover, we will fix a compact set D of R 3 and we require the CWS to be contained in D.
• Uniform distance constraint of the coils to the plasma. To build the vacuum vessel around the plasma, we will assume that the CWS satisfies assumption (H dist,P,δ ).
• Manufacturing cost. In order to avoid irregular shapes that are too difficult to build, we will assume that the CWS has a minimal regularity, say C 1,1 , and a minimal reach condition. More precisely, we will assume that the reach of the CWS is uniformly bounded from below by some r min > 0. We recall that this condition imposes that the curvature radii (where they can be defined) are larger than r min and that there is no bottleneck of distance smaller than 2r min (see, e.g., [1, Figure 3]). To sum-up, S is a C 1,1 closed toroidal surface such that Reach(S) r min > 0.
(H reach,rmin ) As it will be emphasized in what follows, the regularity assumption is actually a consequence of the reach constraint: indeed, the class of sets satisfying a "Reach" constraint is closed in a sense to be specified later and all elements are of class C 1,1 .
Other constraints such as a bound on the two-dimensional Hausdorff measure H 2 (S) of S (in other words the perimeter of the stellarator in R 3 ) will also be considered: To sum-up, let us introduce the admissible set of shapes we will deal with in what follows: Note that the reach condition has been imposed on the surface and not only on the volume. The resulting shape optimization problem we will consider reads In the two following sections, we investigate two important aspects of the shape optimization problem (P shape ). The first one concerns the existence of optimal shapes and is investigated in Section 2. The second one is related to the derivation of first order optimality conditions, at the heart of the algorithms implemented in the last section of this article. To this aim, we apply in Section 3.2 the so-called Hadamard boundary variation method recalled in Section 3.1.

Existence of an optimal current for a given shape (Solving of Problem (P S ))
We first establish that the infimum defining (P S ) is in fact a minimum. Moreover, the minimizer is unique.
The optimization problem (P S ) has a unique minimizer j S . Moreover, one has Proof. First observe that F 0 S is a Hilbert space and that the mapping F 0 is lower semicontinuous for the strong topology of F 0 S . Indeed, let (j n ) n∈N denote a sequence of F 0 S converging to j ∈ F 0 S . According to (H dist,P,δ ), by using the dominated convergence theorem and since K is uniformly bounded in S × P , one has It follows that the functional to minimize is lower semicontinuous (and even continuous) in F 0 S , whence the existence of a unique minimizer j S for Problem (P S ).
Since BS S is continuous, its adjoint BS † S is well defined on F 0 S . It is hence standard that the first order optimality condition for this problem reads The operator BS † S BS S is compact and symmetric. Besides its spectrum is positive and we can therefore consider its resolvent for negative real numbers −λ with λ > 0, so that The expression of C(S) given in (4) follows from a straightforward computation.

Remark 2.
When confronted with the numerical implementation of the shape optimization, motivated by the structure of F 0 S and the properties of the in vacuo Maxwell-equations, we will find it useful to: • optimize on a closed affine subset j a S +F 0 S ⊂ F 0 S instead of the entire set F 0 S . We refer to Section 4.1.3 for further details; • replace the target magnetic field in L 2 (P, R 3 ) by its normal component on the plasma surface (thus, by an object in L 2 (∂P, R)). Indeed, a divergence-free vector field on a 3D domain (in absence of electric currents in the plasma) is nearly entirely characterized by its normal component on the boundary. Further details are given in Section 4.1.2 and Appendix A.3.
Nevertheless, such changes have a minor impact on the theoretical discussion on the shape optimization process and we believe that, for the sake of clarity, it is better to postpone the details about such modifications to Section 4.

Existence of an optimal shape
Theorem 1. The shape optimization problem (P shape ) has at least one solution.
The proof follows the direct method of the calculus of variation. Most of the compactness on our set of admissible shapes comes from the bounded reach assumption (H reach,rmin ). In particular, the following Lipschitz estimate is crucial. Lemma 2 (Theorem 2.8 of [4]). Let V ⊂ R n be a nonempty set such that Reach(∂V ) r min and H n (∂V ) = 0. Then, for every h ∈ (0, r min ), the gradient ∇b V of the signed distance function is 2 rmin−h -Lipschitz on the tubular neighborhood U h (∂V ).
With this estimate, we can state the following compactness result.
Proof. Compactness properties among Hausdorff distances from sets of uniformly positive reach are well known and remain valid for the signed distance (see, e.g., [7,Chapter 6]). Besides, as stated in [4], the convergence property holds true for the strong topology of C 1,α (for α < 1) and for the weak topology of W 2,∞ in a tubular neighborhood of S ∞ . As a consequence, d(S ∞ , P ) δ and Reach(∂V ∞ ) r min . In particular, thanks to Lemma 2, b V∞ is C 1,1 on U r (S ∞ ). Finally, the convergence of H 2 (S n ) to H 2 (S ∞ ) follows from standard results on the continuity of S → H 2 (S) (see [4] or [13]).
The end of this section is devoted to the proof of Theorem 1. Let (S n ) n∈N = (∂V n ) n∈N be a minimizing sequence for Problem (P shape ). Denote by S ∞ a closure point of this sequence in the sense of Lemma 3. In what follows, we will still denote by (S n ) n∈N the converging subsequence introduced in Lemma 3.
We will proceed by showing a semicontinuity property of the criterion, namely that Let j n denote the minimizer of Problem (P S ) for the surface S n , whose existence is provided by Lemma 1.
The idea is to consider a volume integral as approximation of the surface integral in the same spirit as in [6]. For this purpose, we need to extend locally j n to a volume around S n . Notice that, without loss of generality, S n is contained in U rmin (S ∞ ) for every n, which implies, in particular, that ∇b V∞ is everywhere defined and Lipschitz continuous on S n . Let h > 0 be a small constant to be fixed later and define the map where A h (S n ) denote the image of T n . Notice that T n is a bijection between (−h, h) × S n and A h (S n ) if the latter is contained in U rmin (S ∞ ) (see Figure 2).
where ∇ Sn (∇b V∞ )(x 0 ) is a 3 × 3 matrix, according to the notation introduced in Section 1.3, and T x0 S n is identified with a linear subspace of R 3 . We can identify DT n (t 0 , x 0 ) with a 3 × 3 matrix by choosing an orthogonal basis on T x0 S n and its determinant, denoted |DT n (t 0 , x 0 )| in what follows, is independent of such a choice.
Using the regularity of ∇b V∞ near the surface S ∞ , one can prove the following crucial estimate.
Proof. By Lemma 3, b V∞ is C 1,1 in a neighborhood of S ∞ , thus ∇b V∞ | Sn is Lipschitz continuous for n large enough with a Lipschitz constant independent of n and, in particular, there exists C > 0 such that for all n large enough |∇ Sn (∇b V∞ )| L ∞ (Sn) C.
) delimitated by the blue dashed (resp., purple dotted) curves. The black arrows represent the field ∇b V∞ and the red ones represent ∇b Vn . Note that both V n and V ∞ are on the right of the figure. Besides, the linear mapping R x0 : is the unit normal outward vector. Hence, its determinant is equal to 1. Since the determinant is where the reminder term is uniformly bounded with respect to n for n large enough. This concludes the proof.
In what follows ε > 0 is a small parameter to be fixed and h is as in the statement of Lemma 4. We shall also assume that A h (S n ) ⊂ U rmin (S ∞ ) for every n. Notice that, as soon as ε < 1, for n large enough T n is a diffeomorphism. Thus we can define on A h (S n ) the projector p Sn onto S n along the field ∇b V∞ by requiring that p Sn coincides with the S n -component of the inverse of T n (see Figure 3). This allows us to introduce n , defined on A h (S n ) bỹ Using [7, Chap. 7, theorem 8.5], for any m ∈ N ∪ {∞} large enough, we get where the limit is a consequence of Lemma 3.
Using the change of variable formula (also known as area formula for Lipschitz functions), one gets for n ∈ N ∪ {∞} large enough, every f ∈ L 1 (S n ), and every κ ∈ (0, h), which also rewrites Let κ be in (0, h) and η > 0 be small enough so that κ − η > 0 and κ Using that with Equation (6), we obtain which ensure that n belongs to L 2 (S ∞ , R 3 ). (Notice however that n is not necessarily in F 0 S∞ , as it is neither, in general, a tangent vector field nor a divergence free one.) Equation (9) The remaining two steps of the proof consist first in showing the semicontinuity property and then in checking that j ∞ belongs to F 0 S∞ . Notice that, even if we have defined the operator BS S∞ only among the vector fields tangent to S ∞ , by a slight abuse of notation it still makes sense to consider BS S∞ j ∞ , defined using formula (2), even without having checked that j ∞ is in F S∞ .
Both steps rely on the following lemma.

Lemma 5.
Given C > 0 and ε > 0, there exists N ∈ N such that for every r > 0 and every f ∈ C 1 (U r (S ∞ )) such that f L ∞ (Ur(S∞)) ≤ C and f is C-Lipschitz continuous on U r (S ∞ ), we have Proof. Let r > 0 and f ∈ C 1 (U r (S ∞ )) be such that f L ∞ (Ur(S∞)) ≤ C and f is C-Lipschitz continuous on U r (S ∞ ). Up to taking h small enough, we can assume that for n large enough.
As above, consider κ ∈ (0, h), 0 < η < min(κ, h − κ), and n large enough so that (8) holds true. By (7), we have where, for notational simplicity, we write |DT n | for |DT n (T −1 n (x))| and |DT ∞ | for |DT ∞ (T −1 ∞ (x))|. Hence, where we added and subtracted 1 2κ Aκ−η(Sn) f (p S∞ (x)) n (x)|DT ∞ | dx to get We are going to show that A 1 and A 2 can be made arbitrarily small by suitably choosing κ and η (depending only on C and not on the specific function f ) and letting n be large enough.
The term A 2 can be estimated using the inequality f L ∞ (Ur(S∞)) ≤ C, as follows: where the factor 4 comes from the fact that the measure of (−κ − η, −κ + η) ∪ (κ − η, κ + η) is equal to 4η. Notice that j n L 1 (S n ) ≤ j n F Sn H 2 (S n ) is bounded uniformly with respect to n, so that A 2 can be made arbitrarily small by choosing η small enough (depending only on C).
Let us now focus on the term A 1 . Since and because of (11), it follows that for n large enough. Hence, we have the estimates Hence A 1 can be made arbitrarily small choosing ε and then η small enough, and letting n large enough.
Let us start the proof of (10) by comparing BS Sn and BS S∞ . Given y ∈ P , one has Notice that |K(·, y)| is bounded in a neighborhood of S ∞ , uniformly with respect to y ∈ P , since sup (x,y)∈Sn×P |K(x, y)| 1 δ 2 . Moreover, for every y ∈ P and every ρ > 0, the map x → D x K(x, y) is upper bounded by 4 ρ 3 outside U ρ (P ) according to (3). Assume that h < δ, so that A h (S n ) is at distance at least δ − h from P for every n. Consider ρ < δ − h and a Lipschitz neighborhood N of R 3 \ U δ−h (P ) not intersecting U ρ (P ). Since the geodesic distance in N is equivalent to the restriction to N of the standard Euclidean distance, we deduce that there existsC > 0 independent of y such that We deduce from Lemma 5 that for every ε > 0 there exists N > 0 such that for any integer n > N , and, in particular, Using the compactness of BS S∞ , we have This concludes the proof of (10), since ε is arbitrary.
To conclude the proof, it remains to check that j ∞ belongs to F 0 S∞ . By weak convergence of n to j ∞ and according to Lemma 5, where we used that j n is orthogonal to ∇b Vn everywhere on S n . According to Lemma 2, moreover, lim n→∞ ∇b V∞ − ∇b Vn L ∞ (S∞,R 3 ) = 0, and we conclude that j ∞ , ∇b V∞ L 2 (S∞,R 3 ) = 0 since the sequence ( j n L 2 (S∞,R 3 ) ) n∈N is bounded. This proves that j ∞ is a vector field tangent to S ∞ .
To prove that j ∞ is divergence free (in distributional sense), we have to check that j ∞ is orthogonal to Indeed, this characterization of divergence-free vector fields follows from the Hodge decomposition (see Appendix A.1). For g ∈ C ∞ (R 3 ), since div j n = 0 on S n , one has 0 = Sn j n , ∇ Sn g dµ Sn = Sn j n , (∇g − ∇g, ∇b Vn ∇b Vn ) dµ Sn .
Set G n := ∇g − ∇g, ∇b Vn ∇b Vn for n ∈ N ∪ {∞}. Notice that G n converges uniformly to G ∞ in a neighborhood of S ∞ . Hence, again using Lemma 5, This concludes the proof of Theorem 1.
3 Shape differentiation for Problem (P shape )

Reminders on the Hadamard boundary variation method
Let us recall hereafter some notions of topology on sets of regular domains defined in terms of particular perturbations called identity perturbations. The latter are of the form τ = Id +θ, where θ is small enough in a suitable sense. More precisely, according to the approach developed in [26,25], one defines Let us recall that the space W k,∞ (R 3 , R 3 ) endowed with the norm is a Banach space. Choosing τ in T 2,∞ allows to preserve the topological and regularity properties of sets we are interested in, as highlighted in the next result.
The first statement of this lemma follows from standard arguments. For instance, the property "τ (∂Ω 0 ) = ∂(τ (Ω 0 ))" directly results from the fact that τ defines a homeomorphism. The second one comes from a direct application of the Banach fixed-point theorem. We refer to [7, Chapter 4] for more explanations.
As a consequence of the lemma, T 2,∞ induces a topology on the set O 1 of open sets of R 3 whose boundary belongs to C 1,1 (according for instance to [26, Assertion 2.52]). Given Ω 0 ∈ O 1 , let V (Ω 0 , ε) be the set of domains of the type (Id +θ)Ω 0 with θ W 2,∞ ε and ε > 0 small enough so that Id +θ is a diffeomorphism (Lemma 6). Then, one defines a topology on the space O 1 with the help of the neighborhood basis given by the sets of V (Ω 0 , ε). Furthermore, it is shown in [24] that every neighborhood of Ω 0 in V (Ω 0 , ε) is metrizable with a Courant-type distance (induced by that associated with · W 2,∞ in T 2,∞ ) and has the structure of complete separable manifold.
Let us conclude this section by recalling the notion of shape differentiability.

Definition 2. A shape functional Ω → J(Ω) is said to be shape differentiable at Ω (in the sense of Hadamard)
in the class of domains with C 1,1 boundary whenever the underlying mapping with Ω θ = (Id + θ)(Ω), is differentiable in the sense of Fréchet at θ = 0. The corresponding differential dJ(Ω), · is the so-called shape derivative of J at Ω and, by definition of Fréchet differential, the following expansion holds: In the next section we study the shape differentiability of the cost C. In order to fit Definition 2, C is implicitly identified with a functional V → C(∂V ) = C(S) on the set of C 1,1 toroidal domains.

Shape derivative of the cost functional C
This section and the next one are devoted to the computation of the shape derivative of the functional C. Theorem 2. Let S = ∂V ∈ O ad . Let Z P ∈ L(L 2 (P, R 3 ), F S ) and Z P , a bilinear mapping from L 2 (P, R 3 ) × F 0 S into F S , defined by Z P (k) = P K(·, y) × k(y) dµ P (y), where for i ∈ {1, 2, 3}, (X 2 ) i: denotes the i-th line of X 2 seen as a column vector, and ν denotes the outward normal vector to S = ∂V .

Remark 3.
The proof of this result relies crucially on the expression of the magnetic field provided through the Biot and Savart operator BS S (see Definition 1). In general, in many shape optimization problems involving PDEs on bounded domains, PDEs are interpreted as implicit equations on the deformation variable θ and on the state variable. They are in general taken into account by applying the implicit function theorem which also provides an expression for the material (or Lagrangian) derivative of the state with respect to the deformation (see, e.g., [16,Chapter 5]). In the present case, dealing with the Biot and Savart operator comes to consider a PDE on an unbounded domain. The approach we have chosen here, instead, is based on the integral representation of the state variable (the magnetic field here). To establish the above result, we use suitable changes of variables that allows us to rewrite the criterion as an integral over a fixed domain and derive it as a parameterized integral with respect to θ. Although the principle of this calculation is simple, its implementation is not straightforward.

Proof of Theorem 2
For the sake of notational simplicity, the inverse of a group element ϕ ε will be denoted with a slight abuse of notation by ϕ −ε := (ϕ ε ) −1 .
Let S and θ be as in the statement of the theorem. Assume for now that the criterion C is shape differentiable at S. We will comment on this assumption at the end of the proof. In what follows, we concentrate on the computation of the shape derivative in the direction θ.
Since C is shape differentiable at S, we infer that with S ε = (Id +εθ)S.
Step 1: a change of variable. In order to compute C(S ε ), we need to compute some kind of derivative of BS S ε and its adjoint. Nevertheless, we aim to overcome the fact that the domain of BS S ε depends on ε.
Notice that, according to the discussion in Section 3.1, the mapping ϕ ε = Id +εθ induces a bijection between X(S) and X(S ε ). Nevertheless ϕ ε does not map F 0 S into F 0 S ε . This leads us to introduce the linear mapping where J(µ S , µ ε S )ϕ ε denotes the Jacobian determinant 5 of ϕ ε (see Appendix C for further details and the explicit expression of J(µ S , µ ε S )ϕ ε ). The following result will be crucial in what follows since it confirms that Φ ε is indeed a diffeomorphism preserving divergence-free vector fields.

Lemma 7.
For every ε small enough, Φ ε is a diffeomorphism from F 0 S to F 0 S ε . Proof. Since ϕ ε is an orientation preserving diffeomorphism, one has J(µ S , µ ε S )ϕ ε > 0. Besides, As a consequence, Φ ε defines a diffeomorphism from F S to F S ε . We are left to prove that it preserves divergence-free vector fields. According to the Hodge decomposition (see Appendix A.1), it is enough to check Using the change of variables formula (cf. (28)), one has, for every X ∈ X(S), where the notation ϕ ε, * stands for the conormal derivative of ϕ ε . Then X is divergence-free if and only if Φ ε (X) is. The lemma is thus proved.
Step 2: computation of the variation of j. Since we prefer to avoid dealing with operators defined on S ε , we will use Φ ε to relate F 0 S and F 0 S ε . Let us first compute (Φ ε ) † . Let j ∈ F S and g ∈ F S ε . One has We thus infer that (Φ ε ) † is given by Let j ε := Φ −ε (j S ε ). According to Lemmas 1 and 7, j ε is well defined and belongs to F 0 S . To compute the differential of j ε , it is convenient to introduce the operators According to the optimality condition (5) on j S ε , j ε is uniquely characterized by the identity Let us now compute the first order variation of j ε . To this aim, we use the expansion obtained in [16,Lemma 5.4.15]. Recall that the notation div S θ stands for the tangential divergence of θ on S.

Lemma 8. Let S and θ be chosen as above. Then, one has
where, for every j ∈ F 0 S and y ∈ P , and e(θ) is defined as in (1).
Proof of Lemma 8. Let us start with L ε . Given j ∈ F 0 S and y ∈ P , we have Moreover, it can be easily checked that the reminder term of this expansion grows at most linearly with respect to j F S . Regarding Q ε , a similar reasoning using (14) yields concluding the proof.
Combining all the results above, we now compute the sensitivity of j ε with respect to ε. The following result is an immediate consequence of Lemma 8 and (13).
Step 3: computation of the cost functional derivative. Recall that By differentiating this expression and according to Proposition 1, we get Note that Using the fact that ∂C ∂j (j S ) = 0 since j S is the minimizer of j →C(S, j), we get In what follows, we will use the identity stated in the following lemma.

Lemma 9.
Let j ∈ F 0 S , k ∈ L 2 (P, R 3 ), and θ be as in the statement of Theorem 2. Then Proof. The proof follows from straightforward computations, by combining the Fubini theorem with standard properties of the scalar triple product 6 .
By combining (15), (18), and Lemma 9, one computes To conclude this computation, observe that for all vectors u and v in R 3 , We thus obtain where X 1 and X 2 have been introduced in the statement of the theorem. Now, according to [16,Prop. 5.4.9], the shape differential dC(S), θ above can be recast as where, for i ∈ {1, 2, 3}, (X 2 ) i: denotes the i-th line of X 2 seen as a column vector, [(X 2 ) i: ] S is the tangential part of (X 2 ) i: defined as [(X 2 ) i: ] S = (X 2 ) i: − [(X 2 ) i: ] S , ν , and κ denotes the mean curvature 7 on S. The expected formula is obtained by noting that each line of X 2 is tangential (in other words, normal to ν). Indeed, this follows from the definitions of the mapping Z P , the function j S , and the fact that I 3 − νν T corresponds to the matrix of the orthogonal projection onto S. To conclude this proof, it remains to investigate the shape differentiability of S → C(S). Let us introduce τ θ = Id +θ where θ is chosen as in the statement of the theorem. It is straightforward to show that the real number C(τ θ (S)) can be written as a smooth function of integrals written on the fixed domain S, for which the integrand depends regularly on θ. Indeed, this can be straightforwardly obtained by replacing Id +εθ by Id +θ in the reasoning above, and mimicking the associated computations leading to (13), (16) and (17). This yields to the expansion with S ε = (Id +εθ)S and the shape differentiability of C hence follows.

Numerical implementation
The results obtained above are intrinsic in the sense that they do not depend on the parametrization of the objects (surfaces, magnetic field, electric current, . . . ). There are several ways to represent them numerically. We have chosen to use what is, to the best of the authors' knowledge, the classical approach in the stellarator community. In particular: • Surfaces, vector fields and magnetic fields are represented by Fourier coefficients. We detail the parametrization in Section 4.1.
• Stellarator symmetry is imposed on all the objects. We refer to [19,Section 12.3] and [8] for details and justifications of the stellarator symmetry.
• As mentioned in Remark 2, BS S is slightly modified. Not only the optimization space F 0 S is replaced by a suitable affine subspace of it, but also we restrict the image of BS S to the plasma boundary. We provide further details in Sections 4.1.2 and 4.1.3 and Appendix A.3.

Surface representation
We represent a toroidal surfaces S as the image of the two-dimensional flat torus T = (R/Z) 2 by an embedding Stellarators often exhibits a discrete symmetry by rotation. For example W7X is invariant by the rotation of angle 2π/5 along the vertical axis and NCSX has an invariance by the rotation of angle 2π/3. To reduce the complexity, we only represent one module of the surface and we denote by N p the number of modules needed to generate the entire surface (using rotations of angle 2π/N p ). We introduce the cylindrical coordinates (R, ϕ, Z). We will make the assumption of no toroidal folding, i.e., that the intersection of each half plane {ϕ = constant} with S is a single loop. We express ψ in cylindrical coordinates (R(u, v), 2πv Np , Z(u, v)) as Then we develop R and Z in Fourier components and we impose the stellarator-symmetry Note the absence of sin terms for R and cos terms for Z. For the numerical simulation, we truncate the number of Fourier components in (19) and (20).

Remark 5.
The cost considered in this paper only depends on the surface (and is independent of its parametrization ψ). On the toroidal direction, we have already imposed that ϕ = 2πv/N p . On the other hand, we can compose ψ with any diffeomorphism f v : R/Z → R/Z on the poloidal direction u. Namely, ψ(f v (u), v) and ψ(u, v) have the same image for a fixed v. Thus our problem is invariant under the action of a smooth family of diffeomorphisms. This extra degree of freedom has two consequences: • If we use a regular discretization for the surface (ψ( i nu , j nv )) i,j of size n u × n v , we need |∂ u ψ| and |∂ v ψ| to be as regular as possible.
• As we take a finite number of harmonics, we would like to "compress" as much as possible the information on the shape by using low harmonics.
This problem has been study in the plasma community in [17] and gave rise to the notion of spectrally optimized Fourier series. Nevertheless, we would like to highlight that this approach is extrinsic (it depends on the parametrization) and should not be used for other purposes than fixing the gauge invariance. We have not implemented it since our numerical results empirically already provided a reasonable regularity on the poloidal parametrization.

Magnetic field representation
In the previous sections, we represented the target magnetic field as a three-dimensional vector field in the plasma domain. Nonetheless, thanks to the structure of Maxwell's equations, the magnetic field inside P is nearly entirely determined by its normal component along ∂P . It is thus possible to work with a scalar quantity on a surface instead of a three-dimensional vector field on a volume. Indeed, let us introduce the line integral (also called circulation) of the magnetic field along one toroidal turn. By Stoke's theorem (also known as Ampere circuital law in electromagnetism), this is equal to the total flux of the electric current across any surface enclosed by the above-mentioned toroidal loop. This quantity is called the total poloidal current and is denoted by I p .
As proved in Appendix A.3, I p and the normal component of the magnetic field across ∂P characterize completely the magnetic field inside the plasma.
As a consequence, it is reasonable to minimize with the total poloidal current of j fixed, where ν denotes the outward normal unit vector to ∂P . This idea has been used by physicists for a long time, for example [23,21]. Besides, if we consider two currents distribution j and on two toroidal surfaces S andS outside of P with the same total poloidal currents, the induced magnetic field in P satisfies We provide mathematical proofs of these facts in Appendix A.3.
We also use a normal target magnetic field that respects the stellarator symmetry, that is, B m,n sin(2π(mu + nv)).
As before, we truncate the Fourier series to obtain a numerically tractable expression.

Current-sheet representation
As mentioned in the previous section, we need to parameterized all divergence-free vector field on S with a fixed total poloidal current I p . In Appendix A.2 we prove that The following lemma describes how embeddings induce isomorphisms between F 0 T and F 0 S . Lemma 10. Let ψ : T → R 3 be an embedding with S = ψ(T ) and consider Then Ψ induces an isomorphism between F 0 T and F 0 S . The proof is completely similar to that of Lemma 7. Let us suppose now that (u, v) are poloidal and toroidal coordinates for the parameterization ψ, that is, is a loop doing exactly one poloidal turn (and 0 toroidal ones); • Γ t : R/Z t → ψ(0, t) ∈ S is a loop doing exactly one toroidal turn (and 0 poloidal ones).
Besides, as is it in general the convention in the dedicated literature, we assume that ψ is orientation reversing, meaning that with ν the outward normal vector field.
Proof. Remark that div Ψ(X) = 0 ensure that the flux across any loop depends only on the isotopic class of the loop considered. Recall that the flux of Ψ(X) across some loop Γ is given by where the choice of the sign in −ν is due to to the convention (22). Thus, the flux across Γ t (the poloidal flux) is The computation of the flux across Γ p is analogous.
Thanks to this lemma, in order to minimize on F 0 S we fix I p and I t and minimize with respect to Φ. Indeed, I p is fixed by the toroidal circulation of the target magnetic field (see A.3), whereas I t is usually set to 0. This second condition is necessary to ensure the existence of "poloidal coils". Otherwise, no closed field lines would realize one poloidal turn and zero toroidal ones. Thus, the set of admissible currents is described by We say that Φ is the scalar current potential. By stellarator symmetry, its expansion in Fourier series is Φ(u, v) = m 0 n∈Z Φ m,n sin(2π(mu + nv)).
Let us denote j a,S = Ψ( It is straightforward that the affine decomposition J adm (S) = j a,S +F 0 S is compatible with Φ ε (cf. (12)), meaning that for any shape deformation θ, Thus, we can consider the restriction of BS S toF 0 S that we will denote BS S . Let BS † S be its adjoint (inF 0 S ) andπ S the orthogonal projector defined in F 0 S ontoF 0 S . Then Lemma 1 holds and the expression of the unique minimizer is given bŷ
The full code is available on our gitlab 8 under MPL 2 license. The constraints on the perimeter, the reach and the plasma-CWS distance are implemented as a nonlinear penalization cost which blows up rapidly once the values exceed (or subceed) a given threshold. We refer to the code documentation for further details 9 .

Numerical results
In what follows, the data used for the simulations come from the NCSX stellarator equilibrium known as LI383 [37]. We will also use as reference CWS the one used in the original REGCOIL paper [21].
We present here four simulations. We used either λ = 2.5e −16 or λ = 5.1e −19 as regularization parameter in the expression of the cost C. We mesh the CWS and the plasma surface with 64 × 64 grids. The scalar current potential Φ is developed in Fourier series up to order 12 in both directions. The optimization is performed with up to 2000 steps of the BFGS algorithm. In every simulation we implemented a penalization on the perimeter of the CWS (penalization above 56m 2 ) and on plasma-CWS distance (penalization under 20cm). We also implemented a reach penalization for two simulations (penalization under 7.69cm). Let us call Ref the initial CWS. We use DP to refer to the simulations with distance and perimeter penalization and DPR for those with additional reach penalization. The numerical results are summarized in Tables 1 and 2. Figure 7 illustrates the convergence history of the implemented optimization algorithm.  Remark 6. Without penalization on the reach, one naturally obtains better results (as less constraints are applied on the set of admissible shapes). Nevertheless, such an approach seems a very bad idea: • theoretically, since the existence of an optimal shape is guaranteed only for bounded reach, • numerically, because sharper and sharper "spikes" appear, as shown in Figure 4. Those spikes can be arbitrary long while still keeping a finite perimeter (and encapsulated volume).  To obtain the complete CWS, it is enough to make two rotations of angle 2π/3 around the principal axis of the stellarator.

A Some differential geometry
In this section, we recall some basics fact about differential geometry and vector fields on toroidal surfaces and domains.

A.1 Hodge decomposition
We recall in this part some notions of differential geometry and in particular of Hodge theory. We refer to [20, Chapter 3] and [22] for details and precise definitions in the smooth setting. Although we are only interested in C 1,1 manifolds in this article, note for the sake of completeness that details on Hodge theory for Lipschitz manifolds can be found for instance in [35]. The Hodge decomposition is a powerful tool which gives an orthogonal decomposition of the space of square integrable p-forms on a Riemannian closed manifold M as where B p is the L 2 -closure of {dα | α ∈ Ω p−1 (M )}, B * p is the L 2 -closure of {d * β | β ∈ Ω p+1 (M )} (d * is the coderivative), and H p is the set {ω ∈ Ω p (M ) | ∆ H ω = 0} of harmonic p-forms with ∆ H the Hodge Laplacian.
We apply this result to the simple case of 1-forms on a two-dimensional closed Riemannian manifold S. We recall a few basics facts: Figure 7: History of convergence for the implemented optimization algorithm. From left to right, evolution of the costs (left), distance and perimeter constraints (middle) and the curvature constraint (right) along the optimization process. From top to bottom: configurations DPR and DP, respectively for Tables 1 (λ = 2.5e −16 ) and 2 (λ = 5e −19 ).
• 1-forms and vector fields can be identified thanks to the Riemannian metric. This isomorphism is called the musical isomorphism and we denote by X b the 1-form defined as the image of a vector field X thought the musical isomorphism. Conversely w # denotes the vector field which is the image of the 1-form ω.
• The divergence of a vector field X is −d * X b .
We want to show that the space of "divergence-free" 1-forms (i.e., ker d * ) coincides with B * 1 ⊕ H p . It is clear that the latter space is contained in ker d * | Ω 1 (M ) . Conversely, for every exact form ω, i.e., such that ω = df with f ∈ C ∞ (M ), one has d * ω = d * df = ∆ H f . We recall that the Hodge Laplacian coincides with the Laplace-Beltrami operator on 0-forms. But ∆ H f = 0 implies that f is constant on each connected component, thus d * ω = 0 implies that ω = 0. As a result the space of divergence-free 1-forms is B * 1 ⊕ H p = (B 1 ) ⊥ . Equivalently, the space of divergence-free vector fields coincides with the orthogonal to {∇f | f ∈ Ω 0 (M )}. In Appendix A.2 we give an explicit description of (B * 1 ⊕ H p ) # for the two-dimensional flat torus.

A.2 Divergence-free vector field on a flat torus
Let T = (R/Z) 2 be the flat torus with Cartesian parametrization (u, v). We want to characterize the set of divergence-free vector fields on T . As explained in A.1, we only need to characterizes B * 1 (T ) and H 1 (T ).
is a two-dimensional vector space as the first Betti number for a 2D torus satisfies b 1 = 2. We easily compute Using the musical isomorphism, we deduce that all divergence-free vector fields in L 2 have the form given in Equation (21).

A.3 Poisson equation on a toroidal 3D domain
Given a toroidal 3D domain P , we want to study the Maxwell equations in vacuum inside P . We introduce a toroidal loop Γ inside P and denote by I p the electric current-flux across any surface enclosed by Γ. By the conservation of charges (div j = 0), this quantity is well defined. By smoothness of the Biot and Savart operator and of the plasma boundary ∂P , all functions considered in this appendix may be assumed to be C ∞ .

Lemma 12.
Let g be the normal magnetic field on ∂P (i.e., the normal component of B| ∂P ). Then g and I p determine completely the magnetic field B in P . Besides, there exists a constant C > 0 such that for every other magnetic fieldB with the same total poloidal current, |B −B| L 2 (P,R 3 ) C|g −g| L 2 (∂P ) whereg is the normal component ofB| ∂P .
Before going to the proof of this statement, we emphasize that the structure of the space L 2 (P, R 3 ) is well understood and admits a generalized Hodge decomposition (P is not a closed manifold, thus part A.1 does not apply) from which the lemma follows easily. Such a decomposition is proved, for example, in [3]. For completeness, we provide the following proof.
For simply connected domains of R 3 , the complex is an exact sequence, meaning that Im(∇) = ker(curl) and Im(curl) = ker(div). For a 3D toroidal domain, the dimension of the quotient space ker(curl) Im(∇) is always one. This is a consequence of the De Rham cohomology of P . We refer to [22,Diagram 16.15] for further details.
Thus ker(curl) ⊂ Im(∇), i.e., there exists X ∈ C ∞ (P, R 3 ) such that X ∈ Im(∇) and curl X = 0. Without loss of generality, we can suppose that div X = 0. Indeed, it is enough to consider X = X − ∇ζ with ζ solution of the Poisson equation ∆ζ = div X in P, ζ = 0 on ∂P.
To have an intuition, the reader can think of the vector field X = e θ R in R 3 \ {R = 0} in cylindrical coordinates (R, θ, z). This vector field is divergence and curl free but is not in the image of a gradient.
We recall Maxwell's equations for a the static magnetic field in vacuum: div B = 0 in P.
Using Stoke's theorem, the line integral of B along Γ is given by the total flux I p of electric currents across any surface enclosed by Γ. In particular the contribution of the term ∇ξ to I p is zero, yielding The quantity Γ X, Γ |Γ | dµ Γ is nonzero, since otherwise, by the De Rham isomorphism, X would be in Im ∇. Thus, α is uniquely determined by I p , since X does not depend on B.
Equation (24) Thus ξ is determined by g and α as the unique solution of a Laplace equation with Neumann boundary conditions. Finally, let B = ∇ξ + αX andB = ∇ξ + αX with ξ andξ the solutions of equation (25) corresponding to g andg, respectively. The difference δ = ξ −ξ is solution of ∆δ = 0 in P, ∂ n δ = g −g on ∂P.

B Reach constraint and sets of positive reach
In this section, we gather some reminders about the notion of reach. We refer to [7, Chapter 6, Section 6] for more exhaustive explanations around this notion.
Recall first that, if V is a nonempty subset of R n , its skeleton, denoted by Sk(V ), is the set of all points in R n whose projection onto V is not unique. The set V is said to have a positive reach whenever there exists h > 0 such that every point v of the tubular neighborhood U h (V ) has a unique projection point on V .
Recall that the definition of U h (V ) is provided in Section 1.3. One thus defines the reach of V as Reach(V ) = sup{h > 0 | (26) is satisfied}.
An equivalent definition of the reach writes where B h (v) denotes the Euclidean open ball centered at v with radius h.
The notion of reach is actually closely related to the so-called uniform ball condition. The next result make this relationship precise. Theorem 3 (Theorems 2.6 and 2.7 in [4]). Let Ω be an open subset of R n with a nonempty boundary.
• If there exists h > 0 such that Ω satisfies a uniform ball condition, namely then ∂Ω has a positive reach which is larger than h and the Lebesgue measure of ∂Ω in R n is equal to 0. Furthermore, ∂Ω is a C 1,1 hypersurface of R n .
• If ∂Ω has a positive reach and if its Lebesgue measure in R n is equal to 0, then it satisfies the ball condition (27) for every h ∈ (0, Reach(∂Ω)) and in particular, ∂Ω is a C 1,1 hypersurface of R n .

C Jacobian determinant and changes of variables on manifolds
We recall here some basic results about integration on manifolds which can be found in [2] or [33] for example. Let M and N be two compact Riemannian n-dimensional manifolds with volume forms µ M and µ N . Let ϕ : M → N be an orientation preserving diffeomorphism. Then, for any v ∈ C 1 (N ), In the particular, when M and N are closed 2-dimensional submanifolds of R 3 of class C 1,1 , and ϕ is of the type ϕ = Id +θ with θ ∈ W 2,∞ (R 3 , R 3 ) and θ W 2,∞ (R 3 ,R 3 ) < 1 (so that ϕ defines a diffeomorphism in W 2,∞ (R 3 , R 3 )), one has J(µ M , µ N )ϕ = det(Id +Dθ)|((Id +Dθ) ) −1 ν| with ν the outward normal to M . We refer for instance to [16,Section 5.4.5] for a shape optimization oriented proof or [20, Chapter 5] for a more differential geometry oriented presentation.