Some inverse problems around the tokamak Tore Supra

We consider two inverse problems related to the tokamak \textsl{Tore Supra} through the study of the magnetostatic equation for the poloidal flux. The first one deals with the Cauchy issue of recovering in a two dimensional annular domain boundary magnetic values on the inner boundary, namely the limiter, from available overdetermined data on the outer boundary. Using tools from complex analysis and properties of genereralized Hardy spaces, we establish stability and existence properties. Secondly the inverse problem of recovering the shape of the plasma is addressed thank tools of shape optimization. Again results about existence and optimality are provided. They give rise to a fast algorithm of identification which is applied to several numerical simulations computing good results either for the classical harmonic case or for the data coming from \textsl{Tore Supra}.


Introduction
A very challenging potential application of inverse boundary problems for elliptic equations in doubly connected domains is related to plasma confinement for thermonuclear fusion in a tokamak [12].
First of all the tokamak is a concept invented in the Soviet Union and introduced in the 1950s. The meaning of the word is, after translation from the Russian words, "toroidal chamber with magnetic coils" [50] and simply refers to a magnetic confinement device with toroidal geometry (see Figure 1). Today, tokamaks are the most used devices for the fusion experiments and needless to say represent the most suitable approach for their control. In the present paper we focus on the case of the tokamak Tore Supra built at the CEA/IRFM Cadarache (France).
In order to control the plasma position by applying feedback control, it is necessary to know its position in a very short time, or in other words in a time which has to be smaller than the sampling frequency of the plasma shape controller (some microseconds). As a consequence, the poloidal flux function has to be fast identified and computed in a very effective way since its knowledge in the domain included between the exterior wall of the tokamak and the boundary of the plasma is sufficient to determine completely the plasma boundary. Insofar as the only data available are the ones obtained by measurements uniformly distributed on the entire exterior wall of the tokamak, namely the poloidal flux function and its normal derivative, the identification of the plasma boundary may be viewed as the solution of a free boundary problem.
The computation of the poloidal flux function in the vacuum and the recovery of the plasma shape have already been extensively studied (see for example [12]). Consequently, several numerical codes have been developed and we mention here the code Apollo actually used for the tokamak Tore Supra [46] based on an expansion of Taylor and Fourier types for the flux. Others expansions may be found in the literature such as the one making use of toroidal harmonics involving Legendre functions [3]. Naturally, those series expansions are truncated for computations and the coefficients are determined so that they fit to the measurements.
Needless to say, the inverse problem consisting in identifying Dirichlet as well as Neumann data inside a domain from such overdetermined data on part of the external boundary has already generated a lot of methods. We refer for instance to Tikhonov regularization [49], iterativ method [37], conformal mappings [27], integral equations [45], quasi-reversibility [13,36,40] or level-set [6,14,15].
The purpose of this paper is to formulate alternative and original methods for the resolution of the free boundary problem. Two different approaches are considered: a first one making use of Complex Analysis tools and a second one based on Shape Optimization. Up to the knowledge of the authors, such approaches have never been carried out before in the case of the tokamak (although many methods have been developed to solve some close inverse problems, see for instance [6,13,14,27,45]).
Firstly, we are interested in the following inverse problem : assume that the domain under study, denoted by Ω l , is the vacuum located between the outer boundary Γ e of the tokamak and the limiter Γ l (see Section 2.2) in such a way that it refers to a fixed annular domain, or more precisely, to a conformally equivalent doubly-connected domain. The poloidal flux and its normal derivative being given on Γ e , we want to recover those magnetic data on Γ l inside the device. This amounts to solve a Cauchy problem from overdetermined data on part of the boundary. This problem is known to be ill-posed since the work of Hadamard. However sufficient conditions on the available data on Γ e , together with a priori hypotheses on the missing data, may be provided for continuity and stability properties to hold. In order to deal with these constraints, a link is established between the equation satisfied by the poloidal flux in the vacuum and the conjugate Beltrami equation [9]. As a result of this part, we will be able to provide the flux and its normal derivative everywhere in the domain Ω l then more particularly on the limiter Γ l itself.
In a second time, we investigate the inverse problem consisting in recovering the shape of the plasma inside the tokamak, from the knowledge of Cauchy data on the external boundary Γ e (or on Γ l if the resolution of the inverse problem of last paragraph has been performed first). Indeed, if the plasma domain were known, the poloidal flux u inside the tokamak would be determined until the boundary of the plasma Γ p by an elliptic partial differential equation, namely the homogeneous magnetostatic equation for the poloidal flux, with Cauchy data on Γ e . The question of the global existence of a solution for such a system is open and probably quite difficult. The point of view developed in this paper does not need such a result. Indeed, given that the shape of the plasma is a level set of the poloidal flux u, the question of the determination of the shape of the plasma inside the tokamak comes down intrinsically to solve an overdetermined partial differential equation. Noticing that, in general, an overdetermined condition can be interpreted as the first necessary optimality conditions of a Shape Optimization problem (see for instance [2,30]), we chose to see the shape of the plasma as a minimizer of a shape functional in a given class of admissible domains. Furthermore, we use a standard shape gradient algorithm to compute numerically the boundary Γ p .
We point out that the Shape Optimization part has led to numerical developpements that have been used in the present paper. However an effective numerical method based on the Complex Analysis part is under progress. This latter would extend some recent results obtained for the Laplace operator [27] to the case of more general elliptic operators and could probably help to initialize the part treating with the shape of the plasma (see discussion at end of Section 3.3).
The overview of the article is as follows. Section 2 is devoted to general notations and a description of a model of plasma equilibrium. The equations governing this latter allow to characterize the boundary Γ p of the plasma as a level line of the flux in the domain bounded by the limiter. A brief description of the geometry of Tore Supra is provided too. Our main existence and stability results for the Cauchy problem explained above are stated in Section 3. We introduce some generalized harmonic functions associated with the problem, which in fact belong to generalized Hardy spaces of an annulus. A useful density result in such classes is given which enables the resolution of the Cauchy problem formulated as a bounded extremal one. Finally, in Section 4, we introduce a Shape Optimization problem whose solution is supposed to be the shape of the plasma. We first prove the existence of minimizers for such a problem and state the associated necessary first order optimality conditions. In particular, it is shown that, under a regularity assumption for the optimum, the solution of this problem verifies the system characterizing the shape of the plasma. Some numerical simulations, based on an optimization algorithm and the use of the shape derivative, are included at the end of the paper and permit to recover in a satisfying way  We denote by (r, ϕ, z) the three-dimensional cylindrical coordinates system where r is the radial coordinate, ϕ is the toroidal angle and z is the height; e r , e ϕ and e z will denote the axis unit vectors. Thus, given a generic vector A, its component along the unit vectors will be denoted by A r , A ϕ and A z , respectively, so as to have Since the tokamak is an axisymmetric toroidal device, we may assume that all magnetic quantities do not depend on the toroidal angle ϕ. That involves that a plasma equilibrium may be studied in any cross section (r, z), named poloidal section (see Figure 1). Now writing B for the magnetic induction and u for the poloidal magnetic flux, and denoting by ∇ · B = ∂Br ∂r + ∂Bz ∂z , ∇ × B = ∂Bz ∂r − ∂Br ∂z and ∇u = ( ∂u ∂r , ∂u ∂z ), it is easy to derive from Maxwell's equations [12], especially from Gauss's law ∇ · B = 0, the following equations Combining those relations with Ampere's law in the vacuum region located between the outer boundary of the tokamak and the limiter (see Section 2.2), i.e we obtain the following equation for the poloidal flux u in the vacuum It is a linear homogeneous second order elliptic equation in divergence form in two dimensions. Observe that this latter does not stand in the whole space located inside the tokamak owing to the presence of a current density j in the plasma. Here, denoting by j T its toroidal component, equation (3) becomes the so-called Grad-Shafranov [47] equation which is a non-linear elliptic equation in two dimensions. However observe that in the framework if this paper, we only consider equation (3).

The geometry of Tore Supra
We present here the main features of the tokamak Tore Supra built at CEA/IRFM Cadarache in France in a poloidal cross section (r, z). All numeric quantities for the distances are given in meters.
Here Γ e stands for the outer boundary of the tokamak. It corresponds to a circle with center (R, 0) = (2.42, 0) and radius ρ e = 0.92. Inside the tokamak, the first material piece encountered is the limiter Γ l . It corresponds to a closed curve and its goal is to prevent any interaction between the plasma and Γ e . In a nutshell, the plasma cannot pass through Γ l . Inside the limiter stands the plasma, referring to the domain Ω p with a boundary ∂Ω p = Γ p .
Furthermore we add two notations. In the following we will denote by Ω v the domain occupied by the vacuum, i.e the one included between Γ p and Γ e . Moreover Ω l will denote the domain included between Γ l and Γ e . Then refering to (3), the equation satisfied by the poloidal magnetic flux u in Ω v , whence a fortiori in Ω l , is ∇ · ( 1 r ∇u) = 0. It should be noticed that in Tore Supra, the domain Ω l is not exactly an annulus as it can be observed on Figure 6. But as explained in [9, Section 6], all properties of functions solutions to (6) below are preserved by composition with conformal maps. Accordingly, the results of the next section will always refer to a annular domain like in Figure 2.

Definition of the plasma boundary Γ p
A commonly accepted definition of the plasma boundary Γ p may be given in these terms [12]: it is the outermost closed magnetic surface entirely contained in the limiter Γ l . Moreover it corresponds to an equipotential of the function u. In fact, two situations can happen for the plasma equilibrium in tokamak devices: • the limiter case. Here Γ p and Γ l have at least one common point. The equipotential's value of u defining the plasma boundary Γ p is given by the maximum of u on the limiter Γ l . At all contact points, the limiter and the plasma are tangent.
• the divertor case. Γ p and Γ l have no point in common. In this case, the plasma is strictly included inside the limiter and its boundary has an X-point (hyperbolic point).
Let us specify that we made no attempt to study the divertor case. In fact, the paper focuses only on the limiter case, more particularly when the plasma and the limiter have only one point in common. However it should be noticed that this study could easily be applied in the case of several common points. We shall make use of this observation without further notice.
3 How to recover magnetic data in the tokamak Tore Supra ?
Here the differential operator ∇· and ∇ are understood in the sense of distributions (with respect to the real variables r, z in R 2 ). In this section we are dealing with the Cauchy problem described in the introduction. This problem is, as defined by Hadamard, ill-posed for magnetic data u and ∂ n u prescribed on Γ e not only in Sobolev spaces W 1−1/p,p (Γ e ), 1 < p < +∞ but also in L p (Γ e ), where they still make sense [9]. However, well-posedness may be ensured if L 2 -norm constraints is added on Γ l . It follows that the extrapolation issue from boundary data turned out as a best approximation one on Γ e , still named bounded extremal problem (BEP). Note that from now on we restrict ourselves to the particular Hilbertian framework, i.e for p = 2, in order to deal with any possible boundary data on Γ e having finite energy. Remark 1. When σ is constant, this problem reduces to solving a Cauchy problem for the Laplace operator, that is to say to recover the values of a holomorphic function in a domain of analyticity from part of its boundary values. It is well worth noting that this problem has been extensively studied in Hardy spaces on simply and doubly connected domains (see [7,8,25]).

From Grad-Shafranov equation to conjugate Beltrami equation
To study the Cauchy problem for (6), our approach proceeds via a complex elliptic equation of the first order, namely the conjugate Beltrami equation: where ν ∈ W 1,∞ (Ω l ; R) satisfies and where we use the standard notations From this point, it is easy to verify that (7) decomposes into a system of two real elliptic equations of the second order in divergence form. Indeed, assume first that f = u + iv is a solution to (7) with real valued u and v. Putting this in (7) yields that u satisfies (6) while with σ = 1−ν 1+ν . Moreover, from that definition of σ, we obtain that (8) implies (5). Conversely, assume that u is real valued and satisfies (6). By putting ν = 1−σ 1+σ , a simple computation shows that (7) is equivalent to the system This system admits a solution when Ω l is simply connected. Indeed the differential form Consequently it is closed and applying Poincaré's Lemma yields the existence of a real valued function v, unique up to an additive constant, such that f = u + iv satisfies (7) with the definition of ν given above. Again, conditions (5) implies (8).
Note that when Ω l is doubly connected, such like in the situation of Tore Supra, the existence of v as a real single-valued function may only be local. Indeed, from Green's formula applied to u and to any constant function in Ω l , we get But, since (10) may be rewritten on ∂Ω l as ∂v ∂t = σ ∂u ∂n where ∂ ∂t stands for the tangential partial derivative on ∂Ω l , it holds that Hence v is clearly multiplied-valued if Γe σ ∂u ∂n = 0 and f = u + iv a fortiori too. However it is always possible to define u as the real part of a single valued function f in Ω l . See for example [25] where the holomorphic case is processed by the aid of the logarithmic function (well-known to be multiplied-valued in the unit disk). We only mention, without giving more details, that our situation may be similarly solved by making use of toroidal harmonics [4].

Generalized Hardy spaces
Assume that ν ∈ W 1,∞ (Ω l ; R) verifies (8). When handling boundary data in the fractional Sobolev space W 1/2,2 (∂Ω l ), it is a result from [16] that the Dirichlet problem for (7) admits a unique solution in W 1,2 (Ω l ). But as mentioned before, the assumptions on the boundary data are relaxed in a manner that they now belong to the Lebesgue space L 2 (∂Ω l ). In this case, it is obvious that the solution of the Dirichlet problem for (7) has no more reason to belong to W 1,2 (Ω l ) in general. On the other hand the Dirichlet problem will rather have a solution in some generalized Hardy spaces whose definition will follow.
This focuses the attention on the fact that those generalized Hardy spaces are the natural spaces when trying to solve inverse problem linked to (7) for L 2 boundary data (see [23] for the simply connected case). Moreover, in Section 3.3, we will show that the Cauchy problem may be reformulated as a bounded extremal one admitting a unique solution in those appropriate Hardy classes.
Let us denote by D R,ρ and T R,ρ the disk and the circle centered at (R, 0) of radius ρ. Remember that from now Ω l stands for the annular domain shown in Figure 2. Thus we have T R,ρe = Γ e . Definition 3.1. If ν ∈ W 1,∞ (Ω l ; R), we denote by H 2 ν (Ω l ) the generalized Hardy space which consists in Lebesgue measurable functions f on Ω l , solving (7) in the sense of distributions in Ω l and satisfying where Equipped with the norm defined by (11), We remark immediately that when ν = 0 and Ω = D = D 0,1 , H 2 ν (Ω l ) is nothing but the classical H 2 (Ω l ) space of holomorphic functions on the unit disk bounded in norm L 2 on T = T 0,1 (see [22,26]). Recall just that it consists in the functions f ∈ L 2 (T) which Fourier coefficients of negative order vanish. Most of the properties of generalized Hardy spaces on simply connected domains derive from those of the classical H 2 spaces and still hold for multiply connected domains. Basically, the main idea relies on the connection between functions f ∈ H 2 ν (Ω l ) and functions ω satisfying in the distributional sense, for and such that the condition (11) still holds. That connection has been introduced in [11] and leads to .
in other words ∂ω = αω. From the definition of ω, it is obvious that (11) still holds for this latter. By the same way, and noticing that f = ω + ν ω √ 1 − ν 2 , the converse is valid.
Furthermore solutions of type ω, which satisfy condition (11), can be represented as ω = e s g where s is continuous on Ω l and g ∈ H 2 (Ω l ) [9,11]. This remark allows to establish that on simply connected domains, classical and generalized Hardy spaces share similar properties and those latter are naturally exportable to multiply connected domains. We recall in the following result most of them 1. any function f ∈ H 2 ν (Ω l ) has a non-tangential limit almost everywhere on ∂Ω l = Γ e ∪ Γ l , called the trace of f and denoted by trf , which belongs to L 2 (∂Ω l ); 2. trf L 2 (∂Ω l ) defines an equivalent norm on H 2 ν (Ω l ); , trf cannot vanish on a subset of ∂Ω l with positive measure unless f ≡ 0; 5. each f ∈ H 2 ν (Ω l ) satisfies the maximum principle, i.e |f | cannot assume a relative maximum in Ω l unless it is constant.
We refer to [9] for the proofs in simply connected domains.
These properties are necessary to prove the results of next section, i.e. the density of traces of function in H 2 ν (Ω l ) in L 2 (I) whenever I is a subset of non-full measure of the boundary ∂Ω l = Γ e ∪ Γ l , which are the key point to solve extremal problems with incomplete boundary data.

Density results and bounded extremal problem
This section is devoted to the resolution of the Cauchy problem formulated in the introduction.The first step which needs to be established is a density result of fundamental importance which asserts that if I ⊂ ∂Ω l has positive measure as well as J = ∂Ω l \ I, then every L 2 -complex function on I can be approximated by the trace of a function in H 2 ν (Ω l ). Remember that in the situation under study in this paper, the subset I of the boundary ∂Ω l corresponds a priori to the whole boundary Γ e . But in order to ensure density results as well as better algorithm's convergence (see the end of Section 3.3), we restrict from now on the set of available magnetic data to a strict subset, still named I of the outer boundary Γ e . We mention that the case where I is effectively the whole boundary Γ e is already under study.
Thus the following result is a natural extension of [9, Theorem 4.5.2.1] for an annular domain.
Theorem 3.2. Let I ⊂ ∂Ω l = Γ e ∪ Γ l be a measurable subset such that J = ∂Ω l \ I has positive Lebesgue measure. The restrictions to I of traces of H 2 ν -functions are dense in L 2 (I).
From the magnetic data on I, one can compose the function F d = u + iv where v = Γe σ ∂u ∂n according to the formulation of (10) on ∂Ω l . As a consequence of Theorem 3.2, it exists a sequence of functions (f n ) n 0 ∈ H 2 ν (Ω l ) which trace on I converges to F d in L 2 (I). Suppose now that lim n→+∞ trf n L 2 (J) = +∞. Up to extracting a subsequence, we may assume that (trf n ) n 0 is bounded in L 2 (J). But (trf n ) n 0 is already bounded in L 2 (I) by assumption. Consequently it remains bounded in L 2 (∂Ω l ) and then in H 2 ν (Ω l ) by point 2 of Proposition 2. As H 2 ν (Ω l ) is a Hilbert space, there exists a subsequence (trf np ) p 0 converging weakly toward some f ∈ H 2 ν (Ω l ). Now, by restriction, it is obvious that (trf np |I ) p 0 converges weakly to f |I in L 2 (I). Since F d is the strong limit of (trf n ) n 0 in L 2 (I), we conclude that F d = f a.e on I. For this reason F d can either already be the trace on I of a H 2 ν -function, or trf n L 2 (J) → +∞ as n → +∞. This behaviour of the approximant on J shows that the Cauchy solution is unstable, which in fact explains that the inverse problem is ill-posed. Note that in general the function F d does not coincide with the trace of a function belonging to H 2 ν (Ω l ) insofar as magnetic data proceeds from sensor subject, like for all mechanical engineering, to roundoff errors.
Therefore, in order to prevent such an unstable behaviour, an upper bound for the L 2 norm of the approximation on J will be added. Thus the Cauchy problem may be expressed as a bounded extremal one in appropriate Hardy classes. In practical terms, define, for M > 0 and φ ∈ L 2 (J), In light of this definition, it is now possible to find a unique solution for the minimization problem formulated on the class B Ω l .
For every function F d ∈ L 2 (I), there exists a unique solution g * ∈ B Ω l such that Proof. Since B Ω l is clearly convex and the norm . L 2 (∂Ω l ) lower semi-continuous, it is enough to show that B Ω l is closed in L 2 (I) to ensure the existence and the uniqueness of g * . Indeed, there is a best approximation on any closed convex subset of a uniformly convex Banach space which is in this case L 2 (I) [10, Part 3, Chapter II, 1, Propositions 5 and 8].
Let us now prove that, if By making use of the triangle inequality, for all λ such that 0 < λ < 1. Now, taking advantage of the assumption, we have that for λ > 0 sufficiently small Re (g * + λh) − φ L 2 (J) ≤ M . To recap, g * + λh ∈ B Ω l and F d − (g * + λh) L 2 (I 1 ) < F d − g * L 2 (I 1 ) . This contradicts the optimality of g * .
An explicit formulation of g * has been established in [1] in the case ν = 0 which is still valid in the present situation. Indeed, denote by P ν the orthogonal projection from L 2 (∂Ω l ) onto trH 2 ν (Ω l ) (the operator P ν is the natural extension of the classical Riesz projection from L 2 (T) onto trH 2 (D) [22]). The unique solution of the bounded extremal problem, when F d / ∈ B Ω l , is given by where λ > 0 is the unique (Lagrange-type) parameter such that the constraint on J is saturated, in other words Re g * − φ L 2 (J) = M and χ J is the characteristic function of J.
Otherwise, the case where F d ∈ B Ω l corresponds to λ = −1. The behaviour with respect to λ of the error e(λ) = F d − g * 2 L 2 (I 1 ) and the constraint M is discussed in [1]. Observe that g * may be computed constructively if a complete family of trH 2 ν (Ω l ) is known. Indeed, the operator P ν can be expressed thanks to the conjugation operator which, for every function in H 2 ν (Ω l ), associates the trace of Re f to the trace of Im f . A formula for this latter may be found in [5]. Thus it is sufficient to expand any function in L 2 (∂Ω l ) on a basis for which the conjugation operator is easily computed. This method is described in [23] for simply connected domains.
Finally we briefly describe a algorithm for solving numerically the bounded extremal problem. It has been developed in [17] and consists in taking advantage of the fact that initially the magnetic data are available on the entire outer boundary Γ e . As a first step compute the function F d thank the magnetic data on Γ e (as it is explained just after Theorem 3.2) and split the measurement set into two disjoint parts Γ e = I 1 ∪ I 2 . Secondly, solve the bounded extremal problem with respect to χ I 1 F d and a constraint M 1 > 0. One obtain a solution g * 1 which depends of M 1 . From this latter compute the real number M 2 = Argmin . And lastly solve the bounded extremal problem with respect to χ I 1 F d and the constraint M 2 > 0. This cross validation procedure provides good results in the harmonic case, i.e. for the Laplace's equation.
We conclude this section by claiming that the solution g * of the bounded extremal problem could be a starting point of the next section which focuses on the recovery of the plasma boundary by technics of Shape Optimization. Indeed g * provides magnetic data, i.e. u and ∂ n u, at every point of Ω l . It is then possible to initiate the free boundary problem with an outer boundary Γ E located in Ω l , which is then closer to the potential boundary Γ p of the plasma. Moreover this new boundary Γ E can still be Γ l itself.
An other feature, under the assumption that the point of contact M 0 between the plasma and the limiter is known in advance, is to provide the constant c = g * (M 0 ) corresponding to the value of the flux at this point.

Recovering the shape of the plasma inside the tokamak
Tore Supra

A Shape Optimization frame
In this section we make use of the notations introduced on Figure 2. In particular, Γ e denotes the external boundary of the disk centered at (r, z) = (R, 0) with given radius ρ e > 0, representing the external wall of the tokamak Tore Supra. We are interested in the determination of the shape of the plasma inside the tokamak. In the following analysis, we will assume that the plasma is composed of one simply connected component, whose boundary is denoted Γ p . Let us introduce Ω, the domain located between the two closed curves Γ e and Γ p . Physically, as it is underlined in Section 2.3, the only interesting case for our study is the limiter one. It means that there exists a point of the plasma boundary denoted by M 0 , where the plasma meets the toroidal belt limiter Γ p and the poloidal magnetic flux is maximal. We have seen in Section 2.3 that the boundary Γ p may be seen as a level set of the solution. In other words, we look for a domain Ω and its boundary Γ p = ∂Ω\Γ e such that where u still denotes the poloidal magnetic flux inside the tokamak, ∇ denotes the nabla operator with respect to the variables (r, z), and σ is a given analytic function such that 0 < C 1 < σ < C 2 in Ω. We recall that M 0 is such that u(M 0 ) = max x∈Ω u(x) (M 0 is not unique a priori). Indeed, let us recall that in the case of the tokamak Tore Supra, one has (see Section 3.2) σ(r, z) = 1 r and Ω ⊂ D R,ρ .
Although a refined study has been led in Section 3 to obtain explicit weak solutions to (7), whence to (6), with boundary data u 0 and u 1 belonging to a larger space as usually (the trace of a generalised Hardy space), and since we want to use standard technics of Shape Optimization, we make the classical assumptions on the boundary data, that is u 0 ∈ W 1/2,2 (Γ e ) and u 1 ∈ W −1/2,2 (Γ e ).
In other words, Ω can be seen as the solution, if it exists, of the "free boundary problem" where c > 0 is a constant chosen so that M 0 ∈ Γ p . Written under this form, this problem looks like a free boundary problem and we could have a temptation to apply Holmgren's uniqueness or a Cauchy-Kowalewski's type theorem. Nevertheless and to our knowledge (see for instance [34]), until yet, that is only possible to deduce from these theorems, the existence of a neighborhood of the external boundary Γ e where is defined a local solution of the three first equations of (15), but the existence of a closed level set Γ p inside the domain delimited by Γ e seems to be a difficult question. Let us say one word on the uniqueness of solutions for this "free boundary problem", that is the identifiability of the interior boundary curve Γ p from one pair of Cauchy data on the exterior boundary curve. This question has been already studied, for instance in [38,39] and the proof fits immediately into our frame: u 0 and u 1 being fixed, if two domains Ω 1 and Ω 2 respectively associated with the interior plasma boundary Γ 1 p and Γ 2 p satisfy the overdetermined System (15), and if Γ 1 p and Γ 2 p have both a Lipschitz regularity, then Ω 1 = Ω 2 .
Notice that a close family of "free boundary problems" have been studied in the past, referring to overdetermined like problems with additional Bernoulli type boundary condition. However, in these problems, the Bernoulli condition applies to the free boundary and permits in general to ensure some kind of regularity of this latter and its uniqueness. The case presented in this section is quite unusual in the context of Shape Optimization since the overdetermined condition applies on the fixed boundary Γ e . Let us notice nevertheless that this kind of inverse problem has been strongly studied with different approaches and tools (conformal mapping, integral equations, quasi-reversibility, level-sets and so on) and we mention for instance [6,13,14,27,45]. About the literature on free boundary problem with a Bernoulli condition, we refer to [28,29,31,32,33,35,41,42].
Furthermore and with respect to the physical experiments, we will make the following assumptions: coming from the physical frame. Let us notice that, thanks to these additional informations, if we assume moreover that the domain Ω satisfies the interior sphere property condition, then the application of Hopf's maximum principle yields min Ω u = min Γe u 0 and ∂u ∂n (argmin u) = u 1 (argmin u) < 0.
In the same way, one has also max Ω u = max Γp u = c, and ∂u ∂n | Γp > 0.
For the reasons mentioned before, we decided to see this "free boundary problem" as a Shape Optimization one. Indeed, let us introduce, Ω and c > 0 being fixed, the solution u Ω of the following mixed Dirichlet-Neumann elliptic problem Let us denote by D the closure of D R,ρ , i.e. the compact set which boundary is Γ e . If the capacity of the compact set D\Ω is strictly positive and u 1 belongs to W −1/2,2 (Γ e ), then this problem has obviously a unique solution, by virtue of Lax-Milgram's Theorem. Now, a natural idea to solve the overdetermined system (17) consists in introducing a least square type shape functional J through the Shape Optimization problem where O ad denotes the set of admissible shapes that has to be clarified. The natural questions arising here are: • How to choose O ad to ensure the existence of solutions for Problem (18)?
In the following section, we will bring some partial answers to these questions.

An existence result
This section is devoted first to the statement of an existence result for the Shape Optimization problem (18) and secondly to the writing of the first order necessary optimality conditions of this problem. We need to define some convergence notions for the elements of O ad , providing a topology within this class. where d H (K 1 , K 2 ) = max(ρ(K 1 , K 2 ), ρ(K 2 , K 1 )), for any (i, j) ∈ {1, 2} 2 , ρ(K i , K j ) = sup x∈K i d(x, K j ), and ∀x ∈ D, d(x, K i ) = inf y∈K i d(x, y).
• converging to Ω in the sense of characteristic functions if for all p ∈ [1, +∞), • converging to Ω in the sense of compacts if 1. ∀K compact subset of D, K ⊂ Ω ⇒ ∃n 0 ∈ N * , ∀n n 0 , K ⊂ Ω n .
It is in general a hard task to get existence results in Shape Optimization. Not only, many such problems are often ill-posed (see for instance [2,30]), but if the class of admissible domains is too large, a minimizing sequence of domains may converge to a very irregular domain, for which the solution of the partial differential equation may exist in a very weak sense. For these reasons, a partial solution consists in restricting the class of admissible domains, by assuming some kind of regularity. For that purpose, let us define the notion of ε-cone property, introduced in [18].
Let us make precise the class of admissible domains. For a given η > 0, we call D η , the annular compact set D\B(X 0 , η), where X 0 = (0, R) is the center of the disk D.  To prove this proposition, we follow the direct method of calculus of variations. Although this technic is standard, we hope that this result, in the particular frame of plasma Physics, is new.
Proof. Let (Ω n ) n 0 be a minimizing sequence for Problem (18) and let Γ p,n be the part of the boundary contained in Ω n , that is ∂Ω n \Γ e . Since for any n ∈ N, Ω n is included in the compact set D η , there exists Ω ⋆ such that (Ω n ) n 0 converges to Ω ⋆ for the Hausdorff metric and in the sense of characteristic functions (see [30,Chapter 2]). Let us denote by Γ ⋆ p its internal boundary. Furthermore, since any Ω n satisfies an uniform ε cone property, the sequence (Ω n ) n 0 converges to Ω ⋆ in the sense of compacts and necessarily, Ω ⋆ satisfies the uniform ε-cone property [30,Theorem 2.4.10].
To prove the existence result, we need to prove at least the lower semi-continuity of the criterion J in the class O ε,η ad . Let us first prove that the sequence (u Ωn ) n 0 is bounded in W 1,2 (D η ). Notice that due to the homogeneous Dirichlet boundary condition on the lateral boundary Γ p,n , we can extend u Ωn by c outside Γ p,n . So we can consider that the functions are all defined on the box D η (or D) and the integrals over Γ p,n and over D η (or D) will be the same. For that purpose, let us multiply Equation (17) by u Ωn − c and then integrate by part. From Green's formula, we get C 1 Ωn |∇u Ωn | 2 Ωn σ|∇u Ωn | 2 = Γe σu 1 (u Ωn − c), and the right hand side is uniformly bounded with respect to n, due to the fact that the sequence (J(Ω n )) n 0 is convergent and hence bounded. Now, it remains to prove that the sequence (u Ωn ) n 0 is bounded in L 2 (Ω n ) (in the sense that Ωn u Ωn (x) 2 dx is uniformly bounded with respect to n). This is actually a consequence of the non positivity of u 1 . Indeed, the use of Hopf's Theorem yields immediately that u Ωn attains its maximum value at x m ∈ Γ e ∪ Γ p,n and that ∂u Ωn ∂n (x m ) > 0 and hence, max Ωn u Ωn = u Ωn | Γp,n = c. Notice that the regularity condition on Ω n required to apply Hopf's lemma is satisfied by choosing every Ω n having a regular enough boundary, which is clearly not restrictive. Let us introduce u η/2 , the solution of (17) set in the interior of D η/2 . A comparison between u Ωn and u η/2 , using Hopf's maximum principle ensures that for any n ∈ N, u η/2 u Ωn c a.e. in D. Thus, β = max( u η/2 L 2 (D) , πρ 2 c) is an uniform bound of the L 2 (D) norm of u Ωn , which proves the W 1,2 (D) boundedness of u Ωn , and shows at the same time that (u Ωn ) n 0 is bounded in W 1,2 (D η ).
As a result, using Banach Alaoglu's and Rellich-Kondrachov's Theorems, there exists a function u ⋆ ∈ W 1,2 (D η ) such that u Ωn It remains to prove that u ⋆ is the solution of Equation (17) on Ω ⋆ . Let us notice that, once this assertion will be proved, the continuity of the shape functional J will immediately follow, because of the above convergence result. Let us write the variational formulation of (17). For any function ϕ satisfying ϕ ∈ W 1,2 (D η/2 ), ϕ = 0 on Γ e ∪ Γ ⋆ p , and for all n ∈ N, the function u Ωn verifies D σ∇u Ωn · ∇ϕ = Γe σu 1 ϕ.
It suffices to use the weak H 1 -convergence of u Ωn to u ⋆ to show that u ⋆ satisfies also (19).
To conclude, it remains to prove that u ⋆ = c on Γ ⋆ p . This is actually a consequence of the convergence of (Ω n ) n 0 into Ω ⋆ and the fact that Ω ⋆ has a Lipschitz boundary. We refer to Theorem 2.4.10 and Theorem 3.4.7 in [30]. Finally, let us notice that the inclusion Γ e ⊂ ∂Ω ⋆ is verified. Indeed, it follows from the stability of the inclusion for the Hausdorff convergence. To be convinced, one way consist in considering that Γ e is the boundary of an open set Σ that does not contain D, and to replace the constraint Γ e ⊂ ∂Ω by Σ ⊂ Ω, where Ω = Ω ∪ Γ e ∪ Σ. Hence, the stability of the inclusion of open sets for the Hausdorff convergence applies and the conclusion follows.
Notice that the existence result do not guarantee that the solution u Ω ⋆ of (17) set in the optimal domain Ω ⋆ satisfies the Cauchy conditions on Γ e , that is the nonhomogeneous Dirichlet and Neumann boundary conditions on Γ e .

First order optimality conditions
Let us now write the first order optimality conditions for Problem (18). For that purpose, we will first compute the shape derivative of the criterion J. The major difficulty when dealing with sets of shapes is that they do not have a vector space structure. In order to be able to define shape derivatives and study the sensitivity of shape functionals, we need to construct such a structure for the shape spaces. In the literature, this is done by considering perturbations of an initial domain; see [2,30,43,48].
Practically speaking, it is preferable to enlarge the class of admissible domains to write the first order optimality conditions, even if there is a risk of loosing the existence result. Indeed, it would be else strongly difficult for instance to take into account the fact that any admissible domain must satisfy an ε-cone property. For this reason, let us define Let Ω ∈ O ad given. Let us consider a regular vector field V : R 2 → R 2 with compact support, which does not meet the fixed boundary Γ e . For small t > 0, we define Ω t = (I + tV )Ω, the image of Ω by a perturbation of Identity and f (t) := J(Ω t ). We recall that the shape derivative of J at Ω with respect to V is We will denote it by dJ(Ω; V ). Let us now compute this shape derivative. The classical formulae yield whereu denotes the shape derivative of u, and is solution of the system The proof of such an assertion may be found for instance in [24, Section 3.1]. It is more convenient to work with another expression of the shape derivative and to write it as a distribution with support Γ p . For that purpose, let us introduce the adjoint state p defined formally as the solution of the system    ∇ · (σ∇p) = 0, Proposition 4. Assume that ∂Ω is C 2 (in other words, Γ p is C 2 ). Then, the criterion J is differentiable with respect to the shape at Ω and one has for all admissible perturbation As a consequence, if Problem (18) has a solution Ω ⋆ having a internal boundary Γ ⋆ i which is a C 2 closed curve, then Ω ⋆ solves the overdetermined problem (15).
Proof. The shape differentiability is standard (see [20,30]) under the assumption that the boundary of Ω is regular enough. Let us now prove the expression of the shape derivative. First multiply equation (21)  The conclusion follows easily. Now, let us investigate the first order optimality conditions. Since Γ i,⋆ is Lipschitz and surrounds a domain whose interior is nonempty, we deduce that if V belongs to the cone of admissible perturbations, then so does −V . Hence, ∂p ∂n ∂u Ω ∂n = 0 on Γ ⋆ p , and it follows that there exists a subset of Γ ⋆ p of nonzero surface measure on which u Ω ⋆ = c and ∂p ∂n ∂u Ω ⋆ ∂n = 0. Furthermore, by virtue of Holmgren's theorem, necessarily ∂u Ω ⋆ ∂n = 0 on such a subset. Hence ∂p ∂n = 0 on this subset, and since one has also p = 0 on Γ ⋆ p , a second application of Holmgren's theorem yields p ≡ 0. In other words, u Ω ⋆ = u 0 , which is the desired result.

Numerical computations
The method presented here is a gradient method, based on the use of the shape derivative dJ(Ω; V ) computed in Proposition 4. It consists basically in finding a deformation of a given domain Ω guaranteeing the decrease of the criterion J(Ω). To sum up, the algorithm writes: • Initial guess Ω 0 ∈ O ad is given; • Iteration k: Ω k+1 = (Id + t k V k )Ω k where V k is a vector field shrewdly chosen and t k the time step; Practically speaking, we are looking for a vector field V k such that and p k is the adjoint state, solution of Problem (22) set on the domain Ω k . One possible choice for V k is For a precise description of this method or others like the level set method and examples of applications in the domain of Shape Optimization, we refer for instance to [2,14,21,19,44].
To ensure (23), let us write V k as the solution of the problem written under variational form and Ω k (∇V k : ∇ϕ where the test functions ϕ describe the space of functions H 1 (Ω k ) 2 verifying ϕ = 0 on Γ e , and the symbol ":" denotes the tensorial notation used to represent the doubly contracted tensorial product, i.e. for A and B, two smooth vector fields of R 2 , We observed that both approach gave similar results. Before presenting our numerical results on the tokamak Tore Supra, let us begin with some numerical tests, in order to verify the efficiency of our method. A first test case has been done, taking σ = 1, and using the fact that the solution of the Laplacian operator is know on an annular domain. Indeed, let us assume that Ω = D 2 \D 1 where D 1 (respectively D 2 ) denotes the disk centered at the origin of radius R 1 > 0 (respectively R 2 > 0). Then, the solution of is given by Note that the normal derivative of u on ∂D 2 is clearly .
A first test can then be made, looking for the internal radius R 1 , assuming that the normal derivative of u on ∂D 2 is constant equal to c 1 > 0 and the expected solution is obviously The results for this case are plotted on Figure 3.
A second test case has been built, by choosing an arbitrary domain Ω c and a function u 0 , computing the solution u of (17), and choosing for u 1 , the quantity ∂u ∂n . Running the program with these definitions of u 0 and u 1 permits to compare the computed domain with the theoretical domain Ω c . The results for this case are plotted on Figure 4.
Let us now present the simulations of the plasma shape in the tokamak Tore Supra. Measures of u 0 and u 1 on the external boundary of the tokamak, made by Physicists from CEA/IRFM, are the starting point of our algorithm. We interpolated these data thanks to cubic splines to have a complete definition of these functions (see Figure 5). All the simulations presented here have been realized using the software Matlab. With these functions, the algorithm provided the following results plotted on Figure 6. Notice that we did not impose in our algorithm for the plasma to stay inside the limiter, but we fortunately observed it. The point marked by a cross is a control point and the free boundary has to contain this point (it is on the limiter). These conditions being fulfilled, validate our method.   (15) with an arbitrary constant c 1 ≫ c yields a solution u and the associated optimal domain is Ω 1 = {u 0 u c 1 }. In fact, the true solution is provided by Ω ⋆ = {u 0 u u(M 0 )}, where M 0 is a point of the limiter where u reaches its maximal value (that is c !). It can be seen as a continuation method up to a flux level set touching the limiter. Indeed, it is very easy to show that, as a consequence of the maximum principle, Ω ⋆ satisfies (14). On Figure 7 is an example of the solution obtained by running the program with c 1 = 0.2, the expected constant being closed to 0.153.