Ali Elarif

Blaise Faugeras

Francesca Rapetti

Tokamak free-boundary plasma equilibrium computation using nite elements of class C 0 and C 1 within a mortar element approach

Keywords: Tokamak, equilibrium, reduced Hsieh-Clough-Tocher nite element, Newton method, geometrical coecients tokamak, équilibre, elements nis de Hsieh-Clough-Tocher réduits, méthode de Newton, coecients geométriques

The numerical simulation of the equilibrium of the plasma in a tokamak as well as its self-consistent coupling with resistive diusion should benet from higher regularity of the approximation of the magnetic ux map. In this work, we propose a nite element approach on a triangular mesh of the poloidal section, that couples piece-wise linear nite elements in a region that does not contain the plasma and reduced Hsieh-Clough-Tocher nite elements elsewhere. This approach gives the exibility to achieve easily and at low cost higher order regularity for the approximation of the ux function in the domain covered by the plasma, while preserving accurate meshing of the geometric details in the rest of the computational domain. The continuity of the numerical solution at the coupling interface is weakly enforced by mortar projection. A new technique for the computation of the geometrical coecients is also presented.

Introduction

The way a magnetic eld inuences the transport properties of charged particles is of high interest for a wide spectrum of physical systems and areas. It is a complex challenging problem that goes beyond the purpose of the present paper. However, its suggests some key ingredients that need to be correctly treated from the mathematical and numerical points of view.

Among the ingredients we rst have the computation of the equilibrium of a plasma in a Tokamak (see for example [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF] and the references therein). It is a free boundary problem described by the Grad-Shafranov equation in axisymmetric conguration (see more details in [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF]). The right-hand side of this equation is a nonlinear source, which represents the toroidal component of the plasma current density. On the numerical side, the use of nite elements (FEs) enables an accurate resolution of complicated geometric features of realistic tokamak devices, while Newton methods are involved to solve the resulting nonlinear nite-dimensional systems thus ensuring fast convergence. The Newton schemes are highly non-trivial, since the domain covered by the plasma is unknown and depends non-linearly on the poloidal magnetic ux, ψ, the primal unknown. Boundary conditions at innity are consistently incorporated through boundary integral equations [START_REF] Faugeras | FEM-BEM coupling methods for tokamak plasma axisymmetric free-boundary equilibrium computations in unbounded domains[END_REF]. Another ingredient is the correct simulation of the transport in the plasma which is rather complex due to the presence of turbulence at the plasma edge and of thermal exchanges with the chamber wall (see [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF], [START_REF] Falchetto | The European Integrated Tokamak Modelling (ITM) eort: achievements and rst physics results[END_REF]).

An aspect linking the two previous ones is the precise computation of the magnetic conguration and of 1D averaged quantities, which is the object of the present work. In tokamaks the energy and particle transport is anisotropically distributed, that is, along magnetic surfaces it is much greater than across the magnetic surfaces. We can thus assume that densities and temperatures are constant on each magnetic surface. If we label by ρ the magnetic surface S, we dene the average over S of a given quantity u by u = ∂ V ( S u dV ) = 1 V S u dS |∇ρ| where V = ∂ ρ V and V is the volume enclosed inside the surface S. The derivative ∂ ρ V can be computed once suitable 1D proles (such as, for example, ∂ ρ ψ or ∂ ψ ρ) are known. These 1D proles, among which the so-called geometric coecients have to be computed with care. They are important functionals of the solution of such equilibrium problems that are essential to incorporate resistive diusion eects into plasma evolution modeling [START_REF] Grad | Classical diusion in a tokomak[END_REF][START_REF] Hinton | Theory of plasma transport in toroidal connement systems[END_REF][START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF][START_REF] Jardin | Computational methods in plasma physics[END_REF]. Also many plasma characteristics (e.g. the so-called safety factor or the average current density prole), important to quantify stability or for monitoring during the experiment, are dened as integrals involving the gradient ∇ψ of the poloidal ux.

In order to improve the description of ∇ψ we wish to introduce in certain parts of the computational domain, FE functions that are not only continuous, but have also rst order continuous derivatives. For this purpose, instead of relying on composite meshes [START_REF] Heumann | A nite element method with overlapping meshes for freeboundary axisymmetric plasma equilibria in realistic geometries[END_REF], we consider a non-overlapping domain decomposition formulation of the physical problem and couple, by means of mortar projection [START_REF] Bernardi | A new nonconforming approach to domain decomposition: The mortar element method[END_REF][START_REF] Bernardi | Basics and some applications of the mortar element method[END_REF], reduced Hsieh-Clough-Tocher (rHCT) FEs [START_REF] Clough | Finite element stiness matrices for analysis of plates in bending[END_REF][START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] in the plasma domain to piece-wise linear Lagrange FEs in the exterior domain. These rHCT FEs have the advantage to be dened on triangles and to be less expensive from a computational point of view than quintic FEs used in [START_REF] Jardin | A triangular nite element with rst-derivative continuity applied to fusion mhd applications[END_REF]. Concerning the computation of geometric coecients, we propose and compare two strategies. The rst one, rather usual in equilibrium codes, relies on the explicit computation of iso-contours. The second one which proves to be very ecient is based on the coarea formula and motivated by results from [START_REF] Drescher | A high order method for the approximation of integrals over implicitly dened hypersurfaces[END_REF].

The paper is organised as follows. Section 2 briey introduces the equations modelizing the equilibirum of the plasma in a tokamak. In Section 3 we introduce the associated variational formulation in a non-overlapping domain decomposition framework. Section 4 exposes the nite elements used and the proposed coupling method. In Section 5 we move to the fully discretized problem and the implemented Newton method. Section 6 presents the numerical methods pro-Inria posed for the computation of geometric coecents, and in Section 7 we show numerical results and as an example apply the method to the computation of an equilibrium in the JT60-SA tokamak.

2 Free-boundary Grad-Shafranov equation The equations which govern the equilibrium of a plasma in presence of a magnetic eld in a tokamak are the solenoidal condition and Ampère's law in the whole space (including the plasma) and the force balance in the plasma itself, which read, respectively,

div B = 0, curl 1 µ B = J, grad p = J × B, (1)
where p is the plasma kinetic pressure, B is the magnetic induction, J is the current density and µ the magnetic permeability. These equations are sucient for the modeling of the plasma static equilibrium. To simulate the plasma quasi-static evolution, the set of equations ( 1) has to be completed with also Faraday's law in all the conducting structures, and Ohm's laws in the plasma, coils and passive structures. Assuming axial symmetry in the geometry of the tokamak, we may introduce a cylindrical coordinate system (r, ϕ, z), such that r = 0 is the major axis of the tokamak torus. We recall that the transformation x = r cos ϕ and y = r sin ϕ allows to pass from cylindrical coordinates (r, ϕ, z) to Cartesian ones (x, y, z). Equations ( 1) are reformulated in a (poloidal) section ϕ = constant of the tokamak, making the hypothesis that the scalar eld p does not depend on the angle ϕ, thus grad p belongs to the poloidal (r, z)-plane.

The classical primal unknowns for toroidal plasma equilibria described by (1) are the poloidal magnetic ux ψ = ψ(r, z), the pressure p and the diamagnetic function f . The poloidal magnetic ux ψ := rA•e ϕ is the scaled toroidal component (ϕ-component) of the magnetic vector potential A, such that B = curl A, where e ϕ denotes the unit vector for the ϕ coordinate. Note that A is divergence-free (Coulomb gauge) by construction. The diamagnetic function f = rB • e ϕ is the scaled toroidal component of the magnetic induction B. The magnetic induction B can thus be represented as

B = B pol + B tor , B pol = 1 r ∇ ψ × e ϕ , B tor = B ϕ e ϕ = f r e ϕ
as it is assumed to be independent of the angle ϕ, where ∇ is the gradient operator in the r, z coordinates. The previous relations show that the magnetic surfaces are generated by the rotation of the iso-ux lines around the axis of the torus. We refer to standard text books, e.g. [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF], [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF], [START_REF] Wesson | Tokamaks. The International Series of Monographs in Physics[END_REF], [START_REF] Peter | Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF], [START_REF] Goedbloed | Advanced magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF] and [START_REF] Jardin | Computational methods in plasma physics[END_REF] for the details and state in the following only the nal equations. We introduce

Ω ∞ = [0, ∞] × [-∞, ∞],
the positive half plane, to denote the meridian plane that contains the poloidal section of the tokamak, centered at the origin. In axisymmetric coordinates, force balance, the solenoidal condition and Ampère's law in (1) yield the following equation for the ux

ψ(r, z) in Ω ∞ -∆ * ψ = J ϕ , (2)
where J ϕ e ϕ is the toroidal component of J, and the second order elliptic dierential operator -∆ * is dened by

-∂ r 1 µ(ψ)r ∂ r ψ -∂ z 1 µ(ψ)r ∂ z ψ := -∆ * ψ (3)
where µ(ψ) is the magnetic permeability, which is equal to µ 0 , the constant permeability of the vacuum, everywhere except in the possibly existing iron parts of the tokamak (see Figure 1). The geometry of the tokamak determines various subdomains (see Figure 1) which are then used to specify J ϕ accordingly: -Ω Fe ⊂ Ω ∞ denotes those parts of Ω ∞ made of iron where the permeability µ is not constant and given as a (non-linear) function of ψ, namely µ(ψ

RR n°9364 r Ω Fe z Ω c i passive structures Ω p Ω L ∂Ω L 0
) = µ F e (|∇ψ| 2 r -2 ); if Ω Fe = ∅, then µ = µ 0 everywhere; -Ω ci ⊂ Ω ∞ , 1 ≤ i ≤ N c
, denotes the intersection of the ith coil with the poloidal plane. We suppose that ith coil has n i wire turns and cross section area |Ω ci |;

-Ω L ⊂ Ω ∞ , denotes the domain bounded by the limiter, thus the domain accessible by the plasma;

-Ω p ⊂ Ω L , denotes the domain covered by the plasma and the boundary ∂Ω p is the outermost closed ψ-isocontour contained within the limiter region Ω L .

The toroidal component of the current density J ϕ is zero everywhere outside the plasma domain and the poloidal eld coils (and possibly the passive structures). In the coils Ω ci , we set J ϕ = Ii |Ωc i | , where I i is the total current (in At, ampère turns) in the ith coil. In the static modeling, I i is constant whereas in the quasi-static case I i is related to assigned tensions v i (t) in the supplies and to mutual and self inductance via electric circuit equations, as explained in [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF]. In the passive structures, we set J ϕ = 0 for static equilibrium computations whereas for the quasi-static evolution of the equilibrium we need to set J ϕ = -(σ/r)∂ t ψ where σ is the electric conductivity of the passive structures. In the plasma domain Ω p , the equations (1) imply that both the pressure p and the diamagnetic function f are constant on each ψ-isoline, i.e., p = p(ψ) and f = f (ψ). One then deduces the so-called Grad-Shafranov-Schlüter equilibrium equation in the plasma [START_REF] Grad | Hydromagnetic equilibria and force-free elds[END_REF], [START_REF] Shafranov | On magnetohydrodynamical equilibrium congurations[END_REF], [START_REF] Lüst | Axialsymmetrische magnetohydrodynamische Gleichgewichtskongurationen[END_REF] -

∆ * ψ = rp (ψ) + 1 µ 0 r f f (ψ) (4)
where the right-hand side of equation ( 4) is the toroidal component J ϕ of the plasma current density. Functions p and f f , that are zero outside Ω p , are non-linear with respect to ψ. They Inria can be either reconstructed starting from eld measures (this is the goal of the inverse equilibrium reconstruction problem) or supplied as data (in the direct equilibrium problem). In the latter case, we need to introduce few additional notations. The plasma domain Ω p (ψ) is unknown and depends non-linearly on the poloidal ux ψ, in other words, Ω p = Ω p (ψ) is a functional of the poloidal ux ψ (and we have a free-boundary problem). The dierent characteristic shapes of Ω p (ψ) are illustrated in Figure 1: the boundary of Ω p (ψ) either touches ∂Ω L the boundary of Ω L (limiter conguration) or contains one or more saddle points of ψ (divertor conguration). In the latter case, the boundary of the plasma domain is named the magnetic separatrix. The saddle points of ψ, denoted by (r X , z X )=(r X (ψ), z X (ψ)), are called X-points of ψ. The plasma domain Ω p (ψ) is the largest subdomain of Ω L bounded by a closed ψ-isoline in Ω L and containing the magnetic axis (r a , z a ). The magnetic axis is the point (r a , z a ) = (r a (ψ), z a (ψ)), where ψ has its global maximum (or minimum, depending on axis positive direction) in Ω L . For convenience, we introduce also the coordinates (r b , z b ) = (r b (ψ), z b (ψ)) of the point that determines the plasma boundary. Note that (r b , z b ) is either an X-point of ψ or the contact point with ∂Ω L .

The domain of p and f f is the interval [ψ a , ψ b ], with the scalar values ψ a and ψ b being the ux values at the magnetic axis and at the boundary of the plasma (supposing ψ a < ψ b ):

ψ a (ψ) := ψ(r a (ψ), z a (ψ)) , ψ b (ψ) := ψ(r b (ψ), z b (ψ)) . (5)
Since the domain of p and f f depends on the poloidal ux itself, it is more practical to supply these proles as functions of the normalized poloidal ux ψ N (r, z):

ψ N (r, z) = ψ(r, z) -ψ a (ψ) ψ b (ψ) -ψ a (ψ) . (6)
These two functions, subsequently termed S p and S f f , have, independently of ψ, a xed domain [0, 1]. They are usually given as piecewise polynomial functions. Another frequent a priori model is

S p (ψ N ) = λ β r 0 (1 -ψ α N ) γ , S f f (ψ N ) = λ(1 -β)µ 0 r 0 (1 -ψ α N ) γ , (7)
with r 0 the characteristic radius (in meters) of the tokamak vacuum chamber and α, β, γ ∈ R given parameters. The parameter β is related to the poloidal beta [6, p. 15], whereas α and γ describe the peakage of the current prole, λ is a scaling parameter related to the total plasma current. In the following we use the formulation

J ϕ = λ( r r 0 A(ψ N ) + r 0 r B(ψ N )) with A = r 0 S p and B = 1 r 0 µ 0 S f f given functions on [0, 1].
Equilibrium equation ( 2) can be either considered in the whole poloidal plane Ω ∞ , with lim ||(r,z)||→+∞ ψ(r, z) = 0 as boundary condition at innity and ψ(0, z) = 0 at the axis r = 0, or in a restricted bounded domain Ω re , as in Figure 2 (right). In this second case, the boundary of Ω re can be viewed as a measurement contour Γ including Ω L , and possibly some of the coils and passive structures, with ψ(r, z) assigned for all points (r, z) on the boundary of Ω re . Suitable initial conditions on ψ should be added in case we were interested by the quasi-static evolution problem (see [START_REF] Faugeras | An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code[END_REF]). In this paper, we restrict ourselves to static equilibrium computations for an iron-free tokamak.

RR n°9364

As we are going to present later a discretization scheme for the problem (2) that employs dierent approximation spaces on the tokamak poloidal section Ω ∞ , we formulate the variational problem in a bounded spatial domain Ω ⊂ Ω ∞ directly in a non-overlapping domain decomposition framework.

Variational formulation

There are two possible ways to dene a computational domain Ω ⊂ Ω ∞ to be triangulated. It can either be a restricted domain Ω re enclosed in the measurement contour Γ with Ω L ⊂ Ω re , as in Figure 2 (right), or a suciently large semi-circle centered at the origin and containing the whole poloidal section of the tokamak, as in Figure 2 (left). The rst case is called the bounded domain case and the second the ABB domain case, that bears its name from Albanase, Blum, de Barbieri, who rst introduced in [START_REF] Albanese | On the solution of the magnetic ux equation in an innite domain[END_REF] the boundary integral method on the semi-circle used to take into account conditions at innity. In the following, we work with the ABB domain. To dene the ABB domain, we choose a semi-circle γ of radius ρ γ surrounding all the coil domains Ω ci . The ABB domain Ω ⊂ Ω ∞ , that is used for computations, has boundary ∂Ω = γ ∪ γ 0 , where γ 0 := {(0, z), -ρ γ ≤ z ≤ ρ γ }. We have now to select the functional space where the unknown ψ lives. We recall that the coordinates (r, z, ϕ) vary in Ω = Ω × [0, 2 π] ⊂ R 3 and ψ(r, z) = rA ϕ (r, z) where A = (0, 0, A ϕ ) is the magnetic vector potential. The vector potential A has to be selected in H(curl, Ω), the space of vector elds in L 2 ( Ω) 3 with curl in L 2 ( Ω) 3 . We thus need to have: (ii) A and curl A in L 2 ( Ω) 3 , therefore

(i) [A × n] I = 0,
Ω |A ϕ | 2 r dr dz dϕ = 2 π Ω |A ϕ | 2 r dr dz = 2 π Ω |ψ| 2 1 r dr dz < ∞ and Ω |curl A| 2 r dr dz dϕ < ∞ that is 2 π Ω 1 r 2 |∇ (r A ϕ )| 2 r dr dz = 2 π Ω |∇ ψ| 2 1 r dr dz < ∞.
We thus remark that A ∈ H(curl, Ω) if and only if ψ ∈ H 1 (Ω) where

L 2 * (Ω) = {g : Ω → R, g 2 * ,Ω := Ω g 2 1 r dr dz < ∞} and H 1 (Ω) = {u ∈ L 2 * (Ω), ∇ u ∈ L 2 * (Ω) 2 } is the Hilbert space endowed with the norm u 2 H 1 (Ω) = u 2 * ,Ω + |u| 2 H 1 (Ω) being |u| 2 H 1 (Ω) = ∂ r u 2 * ,Ω + ∂ z u 2 * ,Ω . For ψ ∈ H 1 (Ω)
the trace on γ 0 vanishes in the following sense [START_REF] Jiang | Articial boundary conditions for axisymmetric eddy current probe problems[END_REF]

lim r→0 + {r}×[-ργ ,ργ ]∩Ω ψ(r, z) 2 1 r 2 dz = 0 .
To formulate (2) (see Figure 2 (left)) as a variational problem in a non-overlapping domain decomposition framework, we set Ω = Ω in ∪ Ω ex where Ω in is a bounded domain containing Ω L (see Figure 2 (right)) and Ω ex is the complement of Ω in in Ω (see Figure 2 (center)). The boundary of Ω in is denoted I, to recall that it is an interface between the two subdomains Ω in , Ω ex , on which we will impose, at the discrete level, the continuity of ψ, in a weak sense, that is through a mortar-like L 2 projection [START_REF] Bernardi | A new nonconforming approach to domain decomposition: The mortar element method[END_REF]. Note that I = Ω ex ∩ Ω in . Let us now introduce the

functional space V = {(v, w) ∈ H 1 (Ω ex ) × H 1 (Ω in ), v |γ0 = 0, v |I = w |I } .
Continuity is required for ψ in Ω in in order to have meaningful ψ a and ψ b that appear in the denition of Ω p and ψ N [7, Remark I.5, page 18]. The weak formulation of ( 2) is:

Find ψ = (ψ ex , ψ in ) ∈ V such that a(ψ, s) := a ex (ψ ex , v) + a in (ψ in , w) = (I, s) ∀s = (v, w) ∈ V 0,I (8)
where

V 0,I = H 1 0 (Ω ex ) × H 1 0,I (Ω in ).
In [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF], we have set

a ex (ψ, v) := Ω ex 1 µ r ∇ ψ • ∇ v drdz + c(ψ, v), a in (ψ, w) := Ω in 1 µ 0 r ∇ ψ • ∇ w drdz -J p (ψ, w) ,
where

J p (ψ, w) := Ωp(ψ) λ r r 0 A(ψ N ) + r 0 r B(ψ N ) w drdz , (I, s) := Nc i=1 I i |Ω ci | Ωc i (χ Ω ex v + χ Ω in w) drdz , (9)
with (I, s) containing the expression χ Ω ex v + χ Ω in w to take into account the presence of coils in Ω in and Ω ex (here, χ D is the characteristic function of a set D).

RR n°9364

The bilinear form c(•, •) is dened as

c(ψ, ξ) := 1 µ 0 Γ ψ(x)N (x)ξ(x)dS(x)+ 1 2µ 0 Γ Γ (ψ(x) -ψ(y))M (x, y)(ξ(x) -ξ(y)) dS(x) dS(y) , (10)
and accounts for the boundary conditions at innity [START_REF] Albanese | On the solution of the magnetic ux equation in an innite domain[END_REF], with x = (x r , x z ), y = (y r , y z ) and

M (x, y) = k x,y 2π(x r y r ) 3 2 2 -k 2 x,y 2 -2k 2 x,y E(k x,y ) -K(k x,y ) N (x) = 1 x r 1 δ + + 1 δ - - 1 ρ Γ and δ ± = x 2 r + (ρ Γ ± x z ) 2 .
Here, K and E are the complete elliptic integrals of rst and second kind, respectively, and

k x,y = 4x r y r (x r + y r ) 2 + (x z -y z ) 2 .
We refer to [START_REF] Grandgirard | Modélisation de l'équilibre d'un plasma de tokamak[END_REF]Chapter 2.4] for the details of the derivation. Rigorous existence and uniqueness assertion for the general case are still an open problem. See [START_REF] Temam | Remarks on a free boundary value problem arising in plasma physics[END_REF][START_REF] Berestycki | On a free boundary problem arising in plasma physics[END_REF][START_REF] Blum | Existence and control of plasma equilibrium in a tokamak[END_REF][START_REF] Masa | On an equilibrium of the plasma in a tokamak with a limiter[END_REF] for some theoretical work related to such results.

Coupling dierent nite elements on non-overlapping meshes

In the domain containing the plasma, we wish to have a nite element approximation ψ h for the poloidal ux ψ that is not only continuous but has also continuous gradient ∇ ψ h . This is possible if we use the reduced or minimal Hsieh-Clough-Tocher (rHCT) nite element space on a triangular mesh over Ωin . The rHCT triangular element is one of the simplest elements which provide continuous dierentiability of the approximated solution ψ h . For the denition of a rHCT nite element on a triangle we refer to [START_REF] Clough | Finite element stiness matrices for analysis of plates in bending[END_REF], [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] and to [START_REF] Bernardou | Basis functions for general Hsieh-Clough-Tocher triangles, complete or reduced[END_REF] for more details on the basis functions' denition. This regularity is not necessary in the external domain therefore we couple rHCT nite elements in Ω in with continuous piece-wise linear nite elements on a mesh of triangles of Ω. We thus introduce τ ex (resp. τ in ) a mesh of triangles that covers Ωex (resp. Ωin ). The two meshes τ ex , τ in are shape regular and quasi-uniform. They are supposed to match at the common interface I, that is (τ ex ) |I = (τ in ) |I . We assume that I is a polygonal with nodes and edges in τ ex (this is adopted to write the discrete equivalent of the matching condition at interface I). We denote by h ex (resp. h in ) the maximum element diameter in τ ex (resp. τ in ) and set h = max(h ex , h in ).

The nite element spaces

Now, let us introduce the nite element spaces, that we use in the simulations. Locally, on one triangle, they are triples (T, P loc (T ), Σ(T )) where T denotes a triangle of the mesh, P loc (T ) the local space of functions dened on that triangle and Σ(T ) a set of unisolvent degrees of freedom for the functions in the local space (see [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF]). The indices i, i + 1, i + 2, in Denition 2 below, take values 1, 2, 3. When i + 1 > 3 (resp. i + 2 > 3), we replace it by [(i + 1) mod 3] + 1 (resp.

[(i + 2) mod 3] + 1).

Denition 1 Let us denote by

T = [V 1 , V 2 , V 3 ] the triangle of vertices V 1 , V 2 , V 3 . The P 1 La- grange nite element associated with T is the triple (T, P 1 (T ), Σ 0 (T )) where Σ 0 (T ) = {σ i : v → v(V i )} i=1,2,3

Inria

Tokamak free-boundary plasma equilibrium computation using nite elements of class C 0 and C 1 11

B B B V V V 1 2 3 N N N 1 2 1 3 3 2 B B B V V V 1 2 3 N 1 2 1 3 V 4 N 2 N 4 N 3
= [V 1 , V 2 , V 3 ] is cut into three triangles B i : each B i = [G, V m , V ]
having vertices in V m , V with m, ∈ {1, 2, 3} \ {i} and at the barycenter G (denoted by the small circle at the interior) of the triangle T . We can thus reconstruct the height N i = ψ h (V i ) of ψ at the three vertices V i of T and the tangent plane to the surface ψ h at the vertices V i , as generated by

∂ r ψ h (V i ), ∂ z ψ h (V i ).
The coupling with a piece-wise linear element (right).

Denition 2 Let us denote by

T = [V 1 , V 2 , V 3 ] the triangle of vertices V 1 , V 2 , V 3 . The triangle T is divided into three subtriangles B i = [G, V i+1 , V i+2 ]
where G is the barycenter of T (see Figure 3). The reduced Hsieh-Clough-Tocher (rHCT) nite element associated with T is the triple (T, P loc (T ), Σ(T )), where P loc (T ) is the polynomial space of functions w ∈ C 1 (T ) such that,

w |Bi ∈ P 3 (B i ) and (∂ n w) |bi ∈ P 1 (b i ) for any edge b i ∈ ∂B i ∩ ∂T.
Here above, n is the outward normal vector to ∂T , b i the edge on ∂T that does not insist in the vertex V i and (∂ n w) |bi the normal derivative of w at b i , and Σ(T ) is the set dened by the following functionals:

ζ i,00 : w → w(V i ), (11)
ζ i,10 : w → (grad w)(V i ) • (V i+1 -V i ), (12)
ζ i,01 : w → (grad w)(V i ) • (V i+2 -V i ), (13)
where (gradw

)(V i ) • (V i+1 -V i ) is the directional derivative of w on [V i , V i+1 ]. Note that (grad w)(V i ) • (V i+1 -V i ) = (∂ n w) |bi+1 (V i ) -(∂ n w) |bi (V i ).
Globally, over the meshes τ ex and τ in , the nite element spaces are, respectively,

V ex = {v ∈ C 0 (Ω ex ), v |γ0 = 0, v |T ∈ P 1 (T ), ∀ T ∈ τ ex } V in = {w ∈ C 1 (Ω in ), w |T ∈ P loc (T ), ∀ T ∈ τ in }
with the space P loc (T ) dened in Denition 2 . We then denote by V ex ∂ (resp. V in ∂ ) the trace space of V ex (resp. V in ) on the closed polygonal line I. For functions in V ex , the degrees of freedom are given in Denition 1 . One degree of freedom is associated with each node V k ∈ τ ex , so the total number of degrees of freedom in V ex is equal to N ex , the total number of mesh nodes in τ ex , divided into N ex ∂ nodes on the interface I and the remaining N ex o nodes on Ωex \ I. For RR n°9364 functions V in , the degrees of freedom are given in Denition 2. Three degrees of freedom are associated with each node V k ∈ τ in , so the total number of degrees of freedom in V in is equal to 3N in , being N in = N in o + N in ∂ the total number of mesh nodes in τ in divided into N in ∂ nodes on the interface I and the remaining N in o nodes on Ωin \ I. By the denition of V ex and V in there are nite element spaces

V ex • and V in • such that V ex = V ex • ⊕ EV ex ∂ and V in = V in • ⊕ EV in ∂
, where E denotes the trivial extension operators. The elements of V ex

• and V in • have vanishing Dirichlet trace on I.

The discrete coupling condition

Let us denote by {v ex i } i=1,N ex the basis of V ex in duality with the degrees of freedom of Denition 1 associated with vertices V i ∈ τ ex and {w in j } j=1,3N in that of V in in duality with the degrees of freedom of Denition 2 associated with vertices V j ∈ τ in . Then, if u ex and u in represent the vectors gathering the values of degrees of freedom of ψ ex h ∈ V ex and ψ in h ∈ V in we have the decomposition

u ex = (u ex • , u ex ∂ ) and u in = (u in • , u in ∂ ) where u ex • (resp. u in • ) and u ex ∂ (resp. u in ∂ ) are the degrees of freedom in V ex • (resp. V in • ) and V ex ∂ (resp. V in ∂ ).
Since the dierential operator -∆ * is only second order, there is no transmission condition on (∂ n ψ h ) I to be included in the discrete space [START_REF] Bernardi | A new nonconforming approach to domain decomposition: The mortar element method[END_REF]. Examples of MEM coupling between rHCT FEs on both subdomains can be found to solve problems involving four-order operators (we refer to [START_REF] Marcinkowski | Mortar methods for some second and fourth order elliptic equations[END_REF] for convergence studies with application to elasticity). In plasma physics, the requirement of C 1 continuity all over Ω could be necessary when discretizing, for example, reduced MHD models which are characterised by stream function formulations and thus require a proper discretization of four-order operators (see [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF]).

The interpolation case

The imposition of the transmission condition [ψ h ] I = 0 by interpolation means to ask for the equality

(ψ ex h ) |I (V ) = N ex ∂ j=1 (u ex ∂ ) j v ex j (V ) = N in ∂ k=1 (u in ∂ ) k w in k (V ) = (ψ in h ) |I (V ) , ∀ V ∈ I. (14)
If the meshes τ ex , τ in are coincident at I, then [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF] becomes

u ex ∂ = u in ∂ .
If the meshes τ ex , τ in were not coincident at I, relation [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF]

for V ≡ V i ∈ (τ ex ) |I would read P u ex ∂ = D u in ∂ where (P) ij = v ex j (V i ) = δ ij , i, j = 1, N ex ∂ (D) i,k = w in k (V i ), i = 1, N ex ∂ , k = 1, N in ∂ .
By looking at Figure 4, we remark that relation [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF] allows to pass the information on the values of ψ h at the interface I (left-hand side of the gure) without inuencing the value of ∂ n ψ h (right-hand side of the gure). Indeed, with matching grids and a point-wise coupling, the block contained in u in ∂ and associated with the degrees of freedom ( 12), (13) of Denition 2 at V i ∈ I is not involved to dene u ex ∂ .

The L 2 -projection case

The imposition of the transmission condition [ψ h ] I = 0 in a weak sense, by mortar L 2 -projection, means to ask that

I h u ex h z h dI = I h u in h z h dI, ∀ z h ∈ M h (15)
Inria where I h denotes the description of the interface in terms of edges of the slave side mesh (this will become clear with the denition of M h ) and M h is a space of Lagrange multipliers, the choice of which is of crucial importance in the error analysis. For simplicity, we omit the foot-index h on I (also because we consider matching grids at the common interface). The mathematical rationale behind the imposition of the transmission condition by L 2 -projection (rather than the more natural condition of pointwise continuity ( 14) at a chosen set of grid nodes on I) becomes clear from the convergence analysis, shortly presented in Remark 2. This type of condition was rstly studied by Bernardi, Maday and Patera [START_REF] Bernardi | A new nonconforming approach to domain decomposition: The mortar element method[END_REF], who introduced the mortar element method (MEM) to generalize the spectral element method to geometrically nonconforming partitions, to subdomains with dierent resolutions (polynomial degrees) on subdomain interfaces, and also to allow the coupling of spectral element methods with other methods, such as, e.g., the h-version of the nite element method. Its generality and exibility go far beyond these two specic examples (see [START_REF] Bernardi | Basics and some applications of the mortar element method[END_REF]).

Here, we dene the mortar multiplier space as M h = {ξ h ∈ C 0 (I) : ξ h|e ∈ P 1 (e) , ∀ e ∈ (τ ex ) |I }. It coincides with the whole trace space V ex ∂ , due to the fact that I is a closed curve. In the mortar terminology, the choice of M h confers to Ω in the master role (resp. to Ω ex the slave role) with respect to passing the information (value of ψ h ) across I. The left-hand-side of [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF] raises no diculty to be computed since the two discrete functions u ex h , z h live on the same mesh inherited from Ωex on I. On the contrary, the right-hand-side involves discrete functions that live on dierent meshes. Whenever the meshes τ ex , τ in coincide or not at I, to compute eciently the integrals in the coupling condition we rely on numerical quadrature. Relation (15) links the block u ex ∂ to the block u in ∂ by the matrix relation P u ex ∂ = D u in ∂ where this time

(P) i,j = I v ex ∂,i v ex ∂,j dr dz i, j = 1, N ex ∂ , (D) i,k = I v ex ∂,i w in ∂,k dr dz i = 1, N ex ∂ , k = 1, N in ∂ , with I represented in τ ex .
Locally, on each edge e ∈ I, the 2 entries of the block u ex ∂ are dened in terms of the 6 entries of the block u in ∂ . By looking at Figure 5, we can see that relation [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF] allows to pass the information on the values of ψ h at the interface I (left-hand side of the gure) improving the behavior of ∂ n ψ h (right-hand side of the gure). Indeed, with the relation [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF], the block u in ∂ is now completely involved in dening u ex ∂ . We recall that there is no condition imposed on [∂ n ψ h ] |I , however a good behavior of ∂ n ψ h passing through the interface I is important to keep under control the consistency error (as explained later in Remark 2). The behavior of ∂ n ψ h at the interface I cannot be improved further by considering for example quadratic multipliers. According to RR n°9364

The discrete problem

We set

X h = {u h ∈ L 2 (Ω) : u h|Ω in = u in h ∈ V in , u h|Ω ex = u ex h ∈ V ex }
, where h = max(h ex , h in ) and the discrete space V h = {u h ∈ X h : u ex h|γ0 = 0 and (15) }.

The discrete problem to solve reads: Find

ψ h ∈ V h such that a(ψ h , s h ) = ( I, s h ) ∀s h = (v h , w h ) ∈ V ex • × V in • . (17)
The bilinear and linear forms a(., .), ( I, .) are dened as for the problem (8) and evaluated in [START_REF] Peter | Principles of magnetohydrodynamics: with applications to laboratory and astrophysical plasmas[END_REF] for the functions in the discrete space.

Remark 2 The presence of the weak coupling condition (15) prevents V h from being a subspace of V, i.e., we are using a non-conforming method to approximate the solution of the problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF].

The second Strang lemma allows to derive the following error bound for such an approximation

||u -u h || 1, * ≤ c inf z h ∈V h ||u -z h || 1, * + sup z h ∈V h I ∂ n u [z h ] I dI ||z h || 1, * (18)
where

||.|| 1, * is the broken norm ||z h || 2 1, * = ||z ex h || 2 1,Ω ex + ||z in h || 2 1,Ω in .
In the right-hand side of (18), the rst term represents the best approximation error of ψ (that is, the distance between the exact solution ψ and the nite-dimensional space V h ) and the extra error, the second, involving interface jumps, is known as the consistency error. The consistency error is related to the variational crime on the conformity property, due to the fact of dealing with a discrete space V h ⊂ V. The error estimate [START_REF] Grad | Classical diusion in a tokomak[END_REF] is optimal if each term on the right-hand side can be bounded by the norm of local errors arising from the approximation of ψ in Ω ex and Ω in , in an additive fashion. In this way, we can take advantage of the local regularity of the exact solution as well as

Inria

Tokamak free-boundary plasma equilibrium computation using nite elements of class C 0 and C 1 15 the approximation properties enjoyed by the local subspaces. The second term is optimal owing to the special choice of the mortar. In fact, due to the orthogonality property [START_REF] Freidberg | Ideal Magnetohydrodynamics[END_REF], we can subtract from ∂ n ψ an arbitrary function z h ∈ M h . It is therefore important that such functions allow for an optimal approximation of the normal derivative of ψ at the interface. In the present cas, with the choice we made, we have both ∂ n ψ in h ∈ P 1 (I) and M h = P 1 (I), with I described in the same way by τ ex and τ in . [START_REF] Bernardou | Basis functions for general Hsieh-Clough-Tocher triangles, complete or reduced[END_REF] The discrete problem in matrix form Let us rst focus on the discrete problem in Ω in , that reads : Find ψ h ∈ V in , with ψ h = g D on I, such that

Ω in 1 µ 0 r ∇ψ h •∇w h dr dz- Ωp(ψ h ) J ϕ (ψ N,h , r) w h dr dz = Ωc i ⊂Ω in I i |Ω ci | Ωc i w h drdz , ∀w h ∈ V in • , ( 19
) with J ϕ (ψ N , r) = λ( r r 0 A(ψ N ) + r 0 r B(ψ N )) and g D is a given function representing a Dirichlet data. With N in = N in • + N in ∂ , the space V in (resp. V in • ) has dimension n = 2N in + N in • (resp. 3N in = n + N in ∂ ). At a point (r, z) ∈ Ω in , we have ψ h (r, z) = n j=1 ψ j w j (r, z) + N in ∂ =1 g D, v D, (r, z),
where {w j } j (resp., {ψ j } j ) is the set of n basis functions (resp., n degrees of freedom) for the rHCT nite element space V in • and {v D, } the set of N in ∂ basis functions for rHCT nite elements associated with the values of the discrete function at the mesh boundary nodes of I. The reals {g D, } are the values of the Dirichlet boundary condition at the mesh nodes on I.

The normalized discrete ux is ψ N,h (r, z) = (ψ N (ψ h (r, z)), ψ a (ψ h ), ψ b (ψ h )) with ψ a (ψ h ) = ψ h (r a , z a ), ψ b (ψ h ) = ψ h (r b , z b ). ( 20
)
The critical points (r a , z a ) and (r b , z b ) are not necessarily located at nodes of the mesh, as is the case with piece-wise linear nite elements, and in this work we use a Newton method to nd them by solving ∇ψ h (r, z) = 0. We use a six-point quadrature rule of order 4 (see [START_REF] Lyness | Moderate degree symmetric quadrature rules for the triangle[END_REF]) on each sub-triangle for the approximation of integrals in [START_REF] Grad | Hydromagnetic equilibria and force-free elds[END_REF]. A diculty which arises with rHCT nite elements is the computation of the iso-contours and therefore of the plasma boundary Γ p (ψ h ). Contrary to what has been done in [START_REF] Heumann | Quasi-static free-boundary equilibrium of toroidal plasma with CEDRES++: Computational methods and applications[END_REF][START_REF] Faugeras | An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code[END_REF] with piece-wise linear nite elements, we do not compute explicitely the plasma boundary and do not consider a specic quadrature rule for elements T such that the intersection T ∩ Ω p (ψ h ) is neither empty nor the whole triangle. We do not look for such elements T , but rather set A(ψ N,h (r, z)) = 0 and B(ψ N,h (r, z)) = 0 for points (r, z) outside the plasma domain. Once (r b , z b ) and (r a , z a ) are known, the plasma domain can be dened as {(r, z) ∈ Ω, 0 ≤ ψ N,h (r, z) ≤ 1 and ∇ψ h (r, z) • (r a -r, z a -z) ≥ 0}. The integral of the current density term over the plasma domain is then approximated with the same quadrature rule as for the linear term.

Let us denote by ψ the vector gathering the n degrees of freedom {ψ i } i and g D the vector of size N in ∂ collecting the values of the Dirichlet boundary conditions. The rst term in [START_REF] Grad | Hydromagnetic equilibria and force-free elds[END_REF] leads to the linear expression Aψ + A D g D , where A is the n × n stines matrix and A D is the n × N in ∂ RR n°9364 matrix associated with the boundary conditions. The second term leads to the denition of a vector J(ψ) of size n with components

(J(ψ)) i = Ωp(ψ h ) J ϕ (ψ N,h , r) w i dr dz, i = 1, n, (21)
that depend non-linearly on ψ. Using the quadrature method explained above and denoting by T i the subset of elements T of τ in to which the node of the mesh associated with the i th degree of freedom belongs, and by x q,T = (r q,T , z q,T ), ω q,T , the quadrature points and weights in the element T , the integral is computed as

J(ψ) i = T ∈Ti q J ϕ (ψ N,h (x q,T ), r q,T ) w i (x q,T ) ω q,T . (22)
In order to derive a Newton's method, we need to compute the Jacobian matrix Jac ψ (ψ) which follows from the chain rule now that all notations are set:

[Jac ψ (ψ)] ij = T ∈Ti q ω q,T w i (x q,T ) ∂J ϕ ∂ψ N (ψ N,h (x q,T ), r q,T ) [ ∂ψ N ∂ψ ∂ψ h ∂ψ j + ∂ψ N ∂ψ a ∂ψ a ∂ψ j (ψ h ) + ∂ψ N ∂ψ b ∂ψ b ∂ψ j (ψ h )] = T ∈Ti q ω q,T w i (x q,T ) ∂J ϕ ∂ψ N (ψ N,h (x q,T ), r q,T ) 1 (ψ b -ψ a ) [w j (x q,T ) + ψ h (x q,T ) -ψ b (ψ b -ψ a ) w j (x a ) - ψ h (x q,T ) -ψ a (ψ b -ψ a ) w j (x b )] (23)
where w j (x a ) (resp. w j (x b )) are null for indices j not corresponding to the element T a (resp., T b ) where x a = (r a , z a ) (resp., x b = (r b , z b )) lies. Nevertheless Jac ψ (ψ) has entries at indices (i, j) not local to an element. Each element T is linked to T a (resp. T b ) through the w i (x q,T )w j (x a ) (resp. w i (x q,T )w j (x b )) terms. This is not standard but the matrix can still be assembled element by element as it is usually done with nite elements. For the right-hand-side of equation ( 19), we introduce the matrix L in of size n × N in c with entries

L in i,j = 1 |Ω ci | Ωc i w in j dr dz, Ω ci ⊂ Ω in ,
and the vector U in I of size N in c holding the currents I i of the coils Ω ci ⊂ Ω in . Newton's iterations for equation [START_REF] Grad | Hydromagnetic equilibria and force-free elds[END_REF] in its fully discretized form, say e(ψ) = 0 with e(ψ) := Aψ + A D g D -J(ψ) -L in U in I , can be written as

ψ k+1 = ψ k -[e ψ (ψ k )] -1 e(ψ k ), [e ψ (ψ)] = A -Jac ψ (ψ) . (24)
In order to write the matrix form of the discrete problem in the whole domain, instead of a Dirichlet type datum g D on the boundary I, we have to impose the mortar coupling condition. We denote by X the reduced variable such that ψ = Q X with Q the coupling matrix. In detail, we have

ψ =     u ex • u ex ∂ u in • u in ∂     =     I 0 0 0 0 P -1 D 0 I 0 0 0 I       u ex • u in • u in ∂   = Q X,
with the matrices P and D dened in the previous section. So, J(ψ) = J(Q X) = H(X) and when we derive J(ψ) with respect to X we have D X H(X)dX = Jac ψ (ψ)QdX. Newton's Inria iterations for problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] in its fully discretized form, say e(X) = 0 with now e(X)

:= Q [(A + C) Q X -J(ψ) -L U I ]
, can be written as

X k+1 = X k -[e X (X k )] -1 e(X k ), [e X (X)] = Q [(A + C) -Jac ψ (ψ k )] Q . (25)
Here, L (resp. U I ) now takes into account the contributions associated with all the coils (resp. all the currents I i of the coils) contained in Ω in and in Ω ex . See the section on numerical results for a convergence history of the proposed Newton's method.

Computation of geometric coecients

Important outputs of an equilibrium computation are 1D proles computed from integrals of the form

A α,β (y) = Cy r α |∇ψ h | β ds (26)
dened on iso-contours C y = {x ∈ Ω p , ψ N (x) = y}, for y ∈ [0, 1], in the plasma domain. Note that the iso-contour C y is a priori not given by an explicit parametrization, but implicitly as level set of a scalar function. The 1D proles, such as [START_REF] Jardin | Computational methods in plasma physics[END_REF], are essential to incorporate resistive diusion eects into plasma evolution modeling. Also many plasma characteristics (e.g., the safety factor or the averaged current density prole), that are important to quantify stability or for monitoring during the experiment, are dened as integrals over iso-contours [START_REF] Blum | Numerical Simulation and Optimal Control in Plasma Physics with Applications to Tokamaks[END_REF]. These 1D proles are used for example in the denition of averaged quantities

< r α |∇ψ h | β >= A α+1,β-1 A 1,-1 .
With P 1 Lagrange FE, high accuracy in computing these proles can be achieved only by relying on a very ne mesh. The use of rHCT FEs enables to compute a smooth poloidal magnetic eld and to locate precisely the position of the magnetic axis and of the X-point independently of the mesh. At a xed number of degrees of freedom, it improves the precision of the computation of these 1D proles with respect to that with P 1 Lagrange FEs.

The iso-contour method

A drawback of using rHCT FE is that the computation of precise iso-contours for ψ is much more involved than when using P 1 Lagrange FEs. Nevertheless, it is possible and we can proceed as follows: Given the magnetic axis and the point dening the plasma boundary, iterate, until convergence (up to a xed threshold), the following steps.

RR n°9364

1. Find a rst point x 0 y on the levelset Cy, where ψy := ψa + (ψ b -ψa)y thanks to a Newton's method along the ray (xa, x b ). 2. Construct isocontour points by rst following the direction t orthogonal to ∇ψ h with a step a, namely, set x k+1,0 y = x k y + at.

3. Correct x k+1,0 y to obtain a point located exactly on the levelset, with iterations

x k+1, +1 y = x k+1, y + δs ∇ψ h (x k+1, y )
where δs is the Newton's increment for ψ h (x k+1, y + s ∇ψ h (x k+1, y

)) -ψy = 0.

If ||x k+1, +1 y -x k+1, y
|| 2 inferior to a xed threshold ε, stop and proceed with next point, otherwise go back to step 2 and reduce step size a.

Typically in our numerical tests (see Section 7) the step a is set to a fraction of h in , a = hin 4 , the threshold is ε = 10 -12 and the Newton loop in step 3 above converges in 2 iterations.

Once the levelset points are found, the integrals (26) are computed using the trapezoidal quadrature rule. This method works perfectly well and gives precise results but the computation of iso-contours can be time consuming. On the other hand these computations are completely independent one from the other and can thus be easily parallelized.

The weak formulation method

In [START_REF] Drescher | A high order method for the approximation of integrals over implicitly dened hypersurfaces[END_REF] it is proposed to use a weak formulation based on the coarea formula (see Theorem 3.2.12 in [START_REF] Federer | Geometric measure theory. Die Grundlehren der mathematischen Wissenschaften[END_REF]) to compute integrals such as [START_REF] Jardin | Computational methods in plasma physics[END_REF]. The coarea formula underlines the relationship between integrals on iso-contours and integrals on Ω ⊂ R d , here d = 2. Given a function f : Ω → R integrable and ψ N : Ω → [0, 1] Lipschitz continuous, the iso-contour integrals are dened as g f,ψ N (y) := x∈Ω, ψ N (x)=y f (x) dx. The coarea formula states that g f,ψ N ∈ L 1 (0, 1) and

1 0 g f,ψ N (y) dy = Ω f (x) |∇ψ N (x)| dx. ( 27
)
As done in [START_REF] Drescher | A high order method for the approximation of integrals over implicitly dened hypersurfaces[END_REF], from the identity g f,ψ N (y) λ(y) = g f λ(ψ N ),ψ N (y) and ( 27), we have

1 0 g f,ψ N (y) λ(y) dy = 1 0 g f λ(ψ N ),ψ N (y) dy = Ω f (x) λ(ψ N (x)) |∇ψ N (x)| dx. (28)
By combining [START_REF] Lyness | Moderate degree symmetric quadrature rules for the triangle[END_REF] and ( 26), we have :

1 0 A α,β (y)λ(y)dy = Ωp(ψ h ) r α |∇ψ h | β λ(ψ N,h )|∇ψ N,h |drdz, ∀λ ∈ L 2 (0, 1) (29)
In [START_REF] Drescher | A high order method for the approximation of integrals over implicitly dened hypersurfaces[END_REF], to approximate L 2 (0, 1), the space of all polynomials of degree less than or equal to P in [0, 1] is used as discrete space, with a basis of Legendre polynomials. This enables to have a diagonal (P + 1) × (P + 1) mass matrix corresponding to the left-hand side of (29) but on the other hand requires the use of an expensive high order quadrature formula to compute terms corresponding to the right-hand side of (29) since Legendre polynomials up to degree 30 or 40 have to be used. In this paper we have succesfully tested the use of a cubic spline basis which does not require such high order quadrature formulas.

We propose an equivalent formulation of (26) which does not involve anymore the explicit computation of iso-contours. Numerically the method consists in writing A α,β (y) = i a i α,β λ i (y),

Inria

Table 1: Error values in H m (Ω in )-norm, m = 0, 1, 2, for the linear problem (30) using rHCT FEs on rened meshes.

h ||u -u h || L 2 ||u -u h || H 1 ||u -u h || H 2 0.
44194174e-1 0.99885478e-5 0.15105871e-2 0.33870552e-0 0.22097087e-1 0.12031353e-5 0.37243609e-3 0.16757605e-0 0.11048543e-1 0.14814052e-6 0.92608080e-4 0.83400905e-1 0.55242717e-2 0.18396474e-7 0.23099315e-4 0.41611326e-1 where λ i are basis Legendre polynomials or spline functions on [0, 1], and nding the vector of coecients a by inverting the linear system M a = b where

M i,j = 1 0 λ i (y)λ j (y)dy and b i = Ωp(ψ h ) r α |∇ψ h | β λ i (ψ N,h )|∇ψ N,h |drdz
For this last surface integral, as it is done for the equilibrium equation, we do not resolve the plasma domain Ω p (ψ h ) but only check whether the quadrature points are inside or not.

Numerical results

We present some numerical results which highlight the features of the proposed method. We start by checking the correct implementation of rHCT FEs, then show an example of the use of the non-overlapping MEM for the JT60-SA tokamak. The section ends with the computation of several geometric coecients and 1D output proles. These last experiments are conducted with the code NICE, a C++ implementation of the methods for axisymmetric free boundary plasma equilibria described in [START_REF] Faugeras | An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code[END_REF].

Validation of rHCT FE implementation

Since our aim is to combine P 1 Lagrange FEs with rHCT FEs, we can in general only expect that the global error reduction at each mesh renement step is comparable to the error reduction of a solution without MEM, relying on P 1 Lagrange FEs everywhere. To keep the high-order precision everywhere, we would have had to combine rHCT FEs on the exterior with high order rHCT FEs on the interior, with a large increase in the size of the system to solve. On the other hand, the use of rHCT FEs in Ω in really improves the quality of the computation of interesting physical paramaters, as we will see in the following. To validate the precision of the rHCT in Ω in , we consider the linear problem

-∇ • (∇ u) = f in Ω in = [0, 1] 2 u = u D on I = ∂Ω in . ( 30
)
The mesh τ in is a triangular one, satisfying the regularity requirements stated in [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] for the error estimate ||u -u h || m = O(h 3-m ), m = 0, 1, 2, with u h the rHCT FE approximated solution of problem [START_REF] Masa | On an equilibrium of the plasma in a tokamak with a limiter[END_REF]. The data f and u D are consistent with u(r, z) = x 4 (y -1) 2 + y 4 (x -1) 2 being the solution of [START_REF] Masa | On an equilibrium of the plasma in a tokamak with a limiter[END_REF]. Convergence rates (see Figure 6 and the detailed values in Table 1) are in agreement with the theoretical ones.

RR n°9364

1.e-8

1.e-7

1.e-6

1.e-5

1.e-4

1.e-3

1.e-2

JT60-SA tokamak test case

Numerical solution Numerical results for problem [START_REF] Ciarlet | The Finite Element Method for Elliptic Problems[END_REF] are obtained for the conguration presented in Figure 2. The tokamak JT60SA, where JT60 stands for Japan Torus-60 and SA indicates the upgraded model with Superconducting coils, is the largest machine before ITER and is intended to be able to run with the same D-shaped plasma as ITER. The construction of the JT60SA ocially began in 2013 and its assembly was completed between spring and summer 2020. The rst plasma in this machine is planned for the end of 2020, it thus constitutes a perfect benchmark conguration for which numerical simulations can compute and analyse a plasma equilibrium before this equilibrium is physically generated in the machine.

Starting from the JT60-SA tokamak machine description (coils, limiter, vaccum vessel and passive structure) and a coupling interface contour , a mesh is generated using the mesh generator Triangle [START_REF] Shewchuk | Triangle: Engineering a 2D Quality Mesh Generator and Delaunay Triangulator[END_REF]. It is made of 6123 triangles and 3198 vertices for the external domain and 7275 triangles and 3721 vertices for the internal domain. The coupling interface contour is made of 165 vertices. Computations are initialized with an inverse static computation in the full ABB domain using P 1 FE as implemented in the code NICE [START_REF] Faugeras | An overview of the numerical methods for tokamak plasma equilibrium computation implemented in the NICE code[END_REF]. It consists in nding the currents in PF coils which enable to have a desired plasma shape. Functions p and f f are given through parameterization [START_REF] Blum | Existence and control of plasma equilibrium in a tokamak[END_REF] where α = 2, γ = 0.8 and β = 0.5. The scaling parameter λ is computed such that the plasma current is I p = 5.5 × 10 6 A. This provides an initial P 1 FE equilibrium from which the coupled rHCT-P 1 solution is initialized.

It is essential to verify numerically the implementation of the computation of the Jacobian matrix given by [START_REF] Heumann | A nite element method with overlapping meshes for freeboundary axisymmetric plasma equilibria in realistic geometries[END_REF]. This is done by comparison with nite dierences and we check that

e(ε) := || J(ψ + εh) -J(ψ) ε -Jac ψ (ψ)h|| = O(ε) (31)
where the perturbation vector h is chosen randomly. Table 2 gives a typical example of such a numerical test which enables to make sure that the computation of the Jacobian is correct.

The error e starts by decreasing as ε decreases. As it is typical for very small ε values, the error increases due to the accumulation of round-o errors. This is conrmed by the good convergence of the Newton iterations (25) (see Table 3) for the resolution of the non-linear coupled problem.

The non-overlapping MEM for the P 1 -rHCT FEs' coupling works ne to reconstruct the plasma equilibrium in the tokamak (see Figure 7). Inria ε e 10 0

1.13214 × 10 +1 10 -1 4.94353 × 10 -2 10 -2 1.32271 × 10 -2 10 -3 3.15033 × 10 -3 10 -4 6.87336 × 10 -5 10 -5 6.61486 × 10 -7 10 -6 7.23428 × 10 -6 10 -7 7.80983 × 10 -5 10 -8 1.21922 × 10 -3 10 -9 5.06769 × 10 -3

n -X n-1 ||/||X n-1 ||.
With the non-overlapping MEM we are able to introduce FE functions in Ω in that are continuous with continuous derivatives, hence the location of critical points is no more restricted to a nite number of points. Indeed, we have that the plasma axis position (r a , z a ) and the X-point position (r b , z b ) are at points other than mesh nodes of τ in , as shown with a zoom of the solution ψ h in Figure 8, left and right, respectively.

To highlight the inuence of the continuous derivatives we compare the results with the MEM that uses piece-wise linear FEs (P 1 ) instead of the rHCT FEs. In Figure 9 the behavior of the poloidal component 1 r ∇ψ × e ϕ of the magnetic induction is shown as a function of r along the chord {(r, z a ), r min ≤ r ≤ r max }, with r min (resp., r max ) the minimum radius (resp., the maximum radius) of the tokamak vacuum chamber. The staircase eect associated with the discontinuity of the derivatives for piece-wise linear FEs is replaced by a smooth prole when rHCT FEs are adopted in Ω in . Geometric coecients and 1D output proles After an equilibrium is computed one can proceed with the evaluation of several output quantities which are used either to characterize the plasma (e.g the safety factor which plays a role in stating the stability of the plasma) or geometric coecients which are use in transport models which are used to compute the quasistatic evolution of the plasma. Here we concentrate on a few 1D proles for the computation of which we show the benet of using rHCT FE instead of P 1 FE.

For

ψ N ∈ [0, 1], f (ψ) = S f (ψ N ) is computed by integration of B. S f (ψ N ) = [(B 0 r 0 ) 2 -2(ψ b -ψ a )λµ 0 r 0 1 ψN B(x)dx] 1/2 (32)
where B 0 is the vacuum toroidal eld at r = r 0 . For simplicity in what follows we will note f (ψ N ) for S f (ψ N ).

RR n°9364

Let us dene a discretization of the unit interval [0, 1] by S + 1 values ψ 0 N = 0, . . . , ψ S N = 1. These points are taken as abscissa for all computed 1D proles. We also note Ω Again in this last surface integral and as is done for the equilibrium equation we do not resolve the domain Ω ψ s N but only check if the quadrature points are inside or not. Now that the f (ψ s N ), φ(ψ s N ) and ρ(ψ s N ) proles are known we can proceed with other quantities. The safety factor is computed as

q(ψ s N ) = A -1,-1 (ψ s N ) f (ψ s N ) 2π
and the ψ derivative of the plasma volume as Among the many other quantities that can be computed out from an equilibrium conguration we are only going to examine one more, the magnetic shear prole dened as sh := ρ q ∂q ∂ρ which is dicult to compute smoothly because of the q derivative. This derivative is approximated using centered nite dierences.

Figure 10 shows these computed proles for 4 dierent cases: the P 1 FE case and three rHCT cases, using the isocontour method or the weak formulation method with either Legendre polynomials or splines in order to compute integrals A α,β . As expected proles computed from the rHCT solution are smoother than those computed from the P 1 solution. This is particularly clear for example for gm 2 close to the magnetic axis (ψ N = 0) or for the magnetic shear. Proles computed from the rHCT solution with three dierent numerical methods dier very little. Concerning the computation time it can be noted that these runs were done on a laptop with two dual cores using OpenMP to parallelize the isocontour computations. Despite this, the weak formulation method is still faster than the isocontour method.

Conclusion

In this paper we propose a numerical method which enables the computation of a high order approximation of the equilibrium of the plasma in a tokamak. The method relies on the decomposition of the computational domain in two subdomains, the external one, not containing the plasma and in which C 0 piece-wise linear Lagrange FEs are used, and the internal one, containing the plasma in which C 1 rHCT FEs are used. The coupling between the two subdomains is done thanks to a mortar element method for which we showed that the L 2 -projection method enables to ensure an excellent transmission of information from one domain to the other. At the discrete level we propose a Newton method to solve the coupled nonlinear problem and give numerical evidence of its convergence. We also propose two dierent approaches to compute the geometric coecients which are essential outputs of a plasma equilibrium computation. The second one, based on a variational formulation over the plasma domain is fast and precise. It does not require the explicit computation of iso-contours as the rst method does.

RR n°9364

Figure 1 :

 1 Figure 1: Left: Geometric description of the tokamak in the poloidal plane. Middle and right: Sketch for characteristic plasma shapes. The plasma boundary touches the limiter (middle) or the plasma is enclosed by the separatrix, a ux line that goes through an X-point (right).

Figure 2 :

 2 Figure 2: The global ABB conguration containing the poloidal section of the JT60SA tokamak (left), as union of two subdomains, namely the external part not containing Ω L (center) and the internal part containing Ω L (right).

  that is the continuity of the tangential trace of A across any interface I ⊂ Ω. This condition is associated with [B•n] I = 0, the continuity of the normal component of B across I. Since A = A ϕ e ϕ whereas n, the normal vector to I, is in the poloidal section Ω, it is sucient to have A ϕ continuous across I. The transmission condition [ 1 µ curl A × n] I = 0, associated with Inria [H × n] I = 0, the continuity of the tangential trace of the magnetic eld H = (1/µ)B across I, is embedded in the variational formulation and can be recovered by integration by parts.

Figure 3 :

 3 Figure 3: A visualization of the rHCT reconstruction of the function ψ in a triangle T ∈ τ in (left). The mesh triangle T= [V 1 , V 2 , V 3 ] is cut into three triangles B i : each B i = [G, V m , V ]having vertices in V m , V with m, ∈ {1, 2, 3} \ {i} and at the barycenter G (denoted by the small circle at the interior) of the triangle T . We can thus reconstruct the height N i = ψ h (V i ) of ψ at the three vertices V i of T and the tangent plane to the surface ψ h at the vertices V i , as generated by ∂ r ψ h (V i ), ∂ z ψ h (V i ). The coupling with a piece-wise linear element (right).

Figure 4 :

 4 Figure 4: The prole of ψ h (left) and ∂ n ψ h (right) at the interface edge mid-points, with pointwise interpolation on matching grids.

Figure 5 :

 5 Figure 5: The prole of ψ h (left) and ∂ n ψ h (right) at the interface edge mid-points, with mortar L 2 -projection by linear mutipliers.

Figure 6 :

 6 Figure 6: Error decay in H m (Ω in )-norm, m = 0, 1, 2, for the linear problem (30) using rHCT FEs on rened meshes.

  ψ s N = {x ∈ Ω p , ψ N (x) ≤ ψ s N } the domain bounded by C ψ s N .The toroidal ux coordinate is dened and computed as ρ(ψ s N ) = φ(ψ s N )

  ψ s N ). Once this is known we easily obtain the ρ derivative of the plasma volume ∂V ∂ρ Two important geometric coecients are gm 1 :=< 1 r 2 > and gm 2 :=< |∇ρ| 2 r 2 > computed asgm 1 (ψ s N ) = A -1,-1 (ψ s N ) A 1,-1 (ψ s N ) , gm 2 (ψ s N ) = ( ∂ρ ∂ψ (ψ s N )) 2 A -1,1 (ψ s N ) A 1,-1 (ψ s N ).

Figure 7 :Figure 8 :

 78 Figure 7: Iso-contours of ψ h in the whole domain Ω (left) and in Ω L (right), computed by the proposed non-overlapping MEM to have P 1 FEs in Ω ex and rHCT FEs in Ω in .

Figure 9 :

 9 Figure 9: Poloidal magnetic eld ||∇ψ(r, z a )||/r as a function of r in the vacuum vessel. In black the P 1 solution and in red the rHCT solution

Figure 10 :

 10 Figure 10: The prole of ∂ψ ∂ρ (ψ N ) (top, left), q(ψ N ) (top, right), gm 1 (ψ N ) (center, left), gm 2 (ψ N ) (center, right), sh(ψ N ) (bottom, left) and of ∂V ∂ρ (ψ N ) (bottom, right). Black: 101 iso-contours method from the P 1 solution. Red: 101 iso-contours method from the rHCT solution. Cyan: weak formulation method with Legendre polynomials up to degree 40 and a quadrature rule of order 19. Blue: weak formulation method with 43 splines and a quadrature of order 10.

Table 2 :

 2 Convergence of the error e (Eqn (31)) for the computation of the Jacobian Jac ψ (ψ)

	n residual relative error
	1 3.47561 × 10 -2
	2 2.91101 × 10 -4
	3 3.32183 × 10 -7
	4 1.96881 × 10 -13

Table 3 :

 3 Convergence history of Newton iterations, iteration number n and residual relative error ||X

RR n°9364

route des Lucioles -BP 93 06902 Sophia Antipolis Cedex Publisher Inria Domaine de Voluceau -Rocquencourt BP 105 -78153 Le Chesnay Cedex inria.fr ISSN 0249-6399

Acknowledgements

This work is supported by the French National Research Agency grant SISTEM (ANR-19-CE46-0005-03). The rst and second authors have worked within the framework of the EURO Fusion Consortium and have received funding from Euratom research and training programs 2014-2018, 2019-2020, under the grant agreement N.633053. The views and opinions expressed herein do not necessarily reect those of the European Commission. The third author thanks the project team CASTOR at INRIA for the delegation in 2020, during which this work was completed. Inria


Keyan Bennaceur

Patrice Jacques

Fabien Portier

Patrice Roche

D C Glattli

Unveiling quantum Hall transport by Efros-Shklovskii to Mott variable range hopping transition in Graphene

Keywords: numbers: 73.43.-f, 72.20.My, 71.30.+h

come

Since its recent discovery [START_REF] Novoselov | [END_REF]2], the anomalous Quantum Hall Effect (QHE) displayed by relativistic like electrons in a Graphene monolayer has been mostly investigated to search for quantum Hall ferromagnetism [3] or Fractional Quantum Hall effect [4,5] in very low disorder samples to favor interactions effects. Here we address the opposite regime where disorder is strong enough to hide electron interactions, as it is the case for standard exfoliated Graphene monolayers deposited on the oxide layer of a silicon substrate. Graphene offers a new set of parameters to revisit the quantum phase transition of localization in the Quantum Hall regime. In particular, we show that the screening of interactions by the 300 nm close back-gate provided by highly doped silicon allows to observe for the first time the transition from Efros-Shklovskii (E-S) to Mott Variable Range Hopping (VRH) in QHE for large localization length. The universal scaling exponents of quantum localization at Hall plateau transitions deduced from our whole set of temperature and bias current data, together with the E-S to Mott VRH transition, definitely validate the Polyakov-Shklovskii (P-S) suggestion that VRH transport is sufficient to describe the quantum Hall electrical transport [7].

The integer Quantum Hall effect occurs whenever Landau Levels (LLs) which form due to quantization of carrier cyclotron orbits are fully filled. In a conventional 2D-electron gas (2DEG) positive energy LLs with E n = eB/m * (n + 1/2), n integer, lead to LLs filling factors ν = k, and quantized Hall resistance R H = h/ke 2 with k a positive integer and spin degeneracy lifted assumed. In a Graphene monolayer, carriers obey a ultrarelativistic Dirac equation. LLs have both positive and negative energies E n = ±hV F /l c (2n) where V F is the Fermi velocity and l c = ( /eB) the magnetic length. Because of the existence of a zero kinetic energy LL and the non lifted four fold valley and spin degeneracy, the filling factors series leading to Hall plateaus be-comes ν = ±4(k + 1/2). The finite width of the quantized Hall plateaus occurs because of localized states in the bulk provided by potential disorder whose energies spread around the unperturbed LL energies. For filling factors yielding the Hall plateau series mentioned above, the Fermi level E F lies between two unperturbed LLs, i.e. on localized state energies. The longitudinal conductance σ xx vanishes at zero temperature, preventing electron backscattering through the bulk and ensuring perfect quantization of the Hall current circulating on the sample edges. The transition between two Hall plateaus occurs when the Fermi level lies in the middle of a disorder energy broadened LL. For an infinite size sample, the localized state size ξ diverges at a single energy E c ≈ E n resulting in backscattering and a longitudinal conductance peak. According to the quantum localization theory [8]

ξ ∼ |E -E c | -γ ∼ |ν -ν c | -γ
with γ ≃ 7/3, a value confirmed in many experiments on conventional 2DEGs [10,13].

Deducing the scaling exponent from transport measurement requires a model linking transport quantities to ξ. In general one can write σ xx = f (ξ/L(T )) where f is a universal scaling function and L(T ) a characteristic length depending on the conduction mechanisms at finite temperature. Let us first discuss bulk conduction on the Hall plateaus. It is generally accepted that transport occurs via phonon assisted inelastic transitions between localized states, the so-called variable range hopping mechanism. For non-interacting electrons, the VRH Mott's law gives

σ xx ∝ T M T exp(-(T M /T ) 1/3 ) (1)
or equivalently σ xx ∝ ( ξ LM (T ) ) 2 exp(-(L M (T )/ξ) 2/3 ), which defines the characteristic length L(T ) labeled as L M (T ) = 1/πg(ε F )k B T , where g(ε) is the energy independent density of states at the Fermi energy. However, in the QHE regime screening is poor and Coulomb repulsion must be included. One thus enters the Efros-Shklovskii (E-S) VRH regime, where the density of states [14] , with both the length L E-S (T ) = 4πε 0 εk B T /Ce 2 and the energy k B T 0 = Ce 2 /4πεε 0 ξ given by the Coulomb energy. C ≃ 6.2 is a numerical constant [15]. Measuring T 0 thus allows to determine ξ and to probe the scaling law far from the conductance peaks. Still in the same regime, passing a current I trough the sample while keeping a fixed low temperature gives an E-S VRH like law for σ xx where the current plays the role of the temperature. This is the P-S model [16] which uses the effective electronic temperature k B T → eE H ξ/2 where the local Hall electric field E H is proportional to the current I. This leads to

g(E) ∝ |E -E F | yields σ xx ∝ T 0 T exp(-(T 0 /T ) 1/2 ) (2) or σ xx ∝ ( ξ LE-S ) exp(-(L E-S /ξ) 1/2 )
σ xx ∝ exp(-((E 0 /E H ) 1/2 ) (3)
and

σ xx ∝ exp(-((E 1 /E H ) 1/3 ) (4)
for E-S and Mott's VRH respectively. Probing the localization scaling law for large localization length requires understanding the conduction mechanism close to the conductance peaks, in the plateau transition region where the Fermi level approaches a unperturbed LL. Historically, the first conduction mechanism proposed for conductance peaks was a metallic regime. The Pruisken model [17] sets the characteristic length L(T ) as the phase coherence length L φ (T ) = (Dτ φ ) 1/2 . Here D is a diffusion constant and the phase coherence time τ φ ∝ T -p follows a non-universal power law with T (p = 2 accounts for most observations). The localization scaling exponent γ can be indirectly accessed by the temperature dependence of the Full Width at Half Maximum (FWHM) ∆ν of the conductance peaks. The latter is obtained when ξ(∆ν/2) ≃ L φ (T ) giving ∆ν = (T /T 1 ) κ with the non universal exponent κ = p/2γ. However Polyakov and Shklovskii proposed that the VRH regime should last in the plateau transition region and the FWHM obtained from ξ(∆ν/2) ≃ L E-S (T ) (or T ≃ T 0 ) ) giving ∆ν = (T /T 1 ) κ and the now universal κ = 1/γ. Here k B T 1 = Ae 2 /4πεε 0 ξ with A a numerical constant. Similarly, the dependence of the FWHM with bias current using P-S model is ∆ν = (I/I 1 ) µ with µ = 1/2γ = κ/2 while using the phase coherence length approach µ = p/4γ. From the latter discussion we see that, as p = 2 is a reasonable exponent for τ φ , the scaling law of the FWHM with temperature can not discriminate between the two scenariios nor that with bias current.

A further criterion is thus needed to discriminate the two scenarios, and that is addressed experimentally in this Letter. The idea originates from the Aleiner Shklovskii (A-S) [18] prediction that a cross-over from E-S to Mott VRH occurs when interactions are screened, for example by a gate parallel to the 2DEG at a distance d. This requires ξ > 2d which is likely to occur on conductance peaks for sample size ≫ 2d. In the Mott VRH regime the FWHM conductance peak now becomes ∆ν ∝ (T /T 2 ) κ with κ = 1/2γ and ∆ν ∝ (I/I 1 ) µ with µ = 1/3γ. On the contrary in the Pruisken scenario screening is not expected to impact the temperature dependence of the FWHM. This yet never observed E-S to Mott cross-over in the QH regime would definitely establish the P-S scenario, that VRH describes transport almost everywhere in the QHE regime even close to the maximum of the conductance peaks and that the phase coherence length approach is not appropriate.

Previous measurements performed in conventional 2DEGs, including Si-MOSFETs, and InAs/InGAAs or GaAs/AlGAAs heterojunctions have been able to probe the scaling exponent of ξ. Experiments using direct determination of ξ from the E-S VRH [13] and even more directly by geometrical comparison with sample width [11] have given γ ≃ 2.3. Probing the scaling law using the conductance peak width is less direct and showed a dispersion in the extracted values of κ. Works combining temperature and bias current have shown excellent agreement with the P-S model [11,13,19,21]. Recently the scaling law has been studied in Graphene using temperature, but no current bias study was done [20]. Except on the n = 0 LL level the results were found compatible with γ = 2.3. So far no experiments in the QHE have shown the E-S to Mott cross-over. It has been only observed in highly disordered 2D electron systems in zero field [6]. Here, the cross-over occurs not because of screening but for energies well above the Coulomb gap, restoring a constant density of states.

In this paper we present a complete set of data in temperature and bias current on the QHE regime performed on Graphene monolayers. The silicon back-gate with d = 300nm gives for the first time access to the cross-over from Mott to ES VRH regime on conductivity peaks for the highest filling factors where ξ is large. All the scaling exponents γ, κ and µ agree with the P-S and A-S prediction in the screened and unscreened regime. This provides a definitive confirmation of these models. The scaling exponent found to be γ ≃ 2.3 is the same for first two n = ±1, ±2 LLs including the n = 0 LL where no anomalous behavior is observed as found in Ref. [20].

Four samples (S1 to S4) have been fabricated using exfoliation of natural graphite flakes [22] . All samples have been deposited on the 300nm thick oxide layer of highly doped silicon wafer which serves as a back gate. Contacts where made using e-beam lithography and evaporating 5/70 nm Ti/Au in high vacuum. After processing, S1 and S2 were covered with PMMA while S3 and S4 were not. Samples S3 and S4 were heated up to 450K in cryogenic vacuum during several hours until the four points resistance reached a steady high value signaling low dopant concentration. Mobilities at 1.10 12 cm -2 are around 3000 cm 2 V -1 s -1 for S1 and S2, 6000 cm 2 V -1 s -1 for S3 and 10000 cm 2 V -1 s -1 for S4. Sample S4 showed the quantum Hall plateaus series and a Raman spectrum of a monolayer although the gate efficiency was two times smaller than expected. It is likely a twisted bilayer.

The longitudinal resistivity ρ xx was recorded while varying density for about 50 temperature values ranging from 1.6K (S4) or 4K (S3) to 300K and similar runs are repeated for different fixed high magnetic fields from 6 to 17 Tesla. These extensive measurements have been done on two samples (S3 and S4). Measurements as a function of bias current were also performed ranging from 10nA to 100µA at fixed temperature (4.2K or 1.6K). These measurements have been done on the four samples. Four point measurements were done on S3 (R xx ) and S4 (R xx and R hall ) whereas 3 point measurements were performed on the other two samples. In order to extract the resistivity from the resistance the aspect ratio was estimated using a numerical electrostatic calculation. For all measurements the density is slowly swept (less than 2 V/mn ) to avoid hysteretic effects due to trapped charges in silicon oxide. To get rid of the weak remaining hysteresis the data are taken always in the same sweeping direction. Most of the data presented in this Letter are from sample S3 and some of the data from sample S4. Other samples were used to confirm the reproducibility of the physical properties.

We now present the experimental results. An example of ρ xx measurements is shown in Figure 1 first focus on the Hall plateau regions. The variation of σ xx with T was shown to follow accurately the E-S VRH law [7], see supporting material [23]. This is well obeyed from 1.6K to ≃ 80K throughout the ν = ±2 plateaus for all magnetic fields. We emphasize that both the simple activated law and the 2D Mott's VRH law yields poor fits. Above 100K departure from the E-S VRH law signals thermal activation to the next LLs. From the fit combining activated and E-S VRH law we can extract both the VRH temperature T 0 as shown in Fig. 2 and the activation energy ∆. The observed ∆ are smaller than those given by the Dirac equation in magnetic field

∆ = E n+1 -E n = 1 2 ( √ n + 1- √ n) 2e v 2
F B (solid lines) because of LL disorder broadening. At 17T, 200K and 70K LL broadening are found respectively for sample S3 and S4. The values for S3 are comparable with those obtained in [24], see [23] for a complete set of data.

The dependence with current bias shows accurate agreement with the E-S VRH like law, Eq.3, following the P-S prediction. Comparing Eqs. 2 and 3 we can write ( εH0 εH ) α1 = ( T0 T ) α2 + k, where k is a constant. If the E-S law is obeyed for bias currents one should have α 1 = α 2 = 1/2, where α 2 = 1/2 was already established. In order to check this, we define the effective temperature T ef f (I) such that ρ xx (I) = ρ xx (T ef f ). On Fig. 1, right inset, T ef f is plotted as a function of the bias current in the logarithmic scale for ν = ±2, ±6. It is clear that below 100K (no thermal activation) T ef f ∝ I showing that α 1 = α 2 and a VRH like law for current is well obeyed by σ xx . Fig. 2 shows T σ xx (T ) using Eq.3 for a continuous series of filling factors except on the very maximum of the conductance peaks. From these measurements we can extract the localization length via T 0 using: ξ(ν) = Ce 2 /4πǫ r ǫ 0 k B T 0 (ν). In the inset of Fig. 2 the smallest localization length ξ min found in the middle of the plateaus is plotted as a function of the magnetic field. The solid lines show the magnetic length l c for comparison. In sample S3, ξ min ∼ 40 nm for ν = ±2 at B = 16.5 Tesla is around seven times larger than l c whereas in sample S4 ξ min is of the order of l c . Smaller localized states can be expected due to higher mobility of S4 or, if it is a twisted bilayer, due to screening of the silicon oxide charge impurities by the lower layer. ξ values in S3 are consistent with measurements of Ref. [25] but quite below those found in Refs. [20,26].

In Fig. 3, lower graph, we show ξ(ν) extracted from a E-S VRH analysis of the data for continuous values of the filling factor and for different magnetic fields. The line at ξ = 2d = 600nm signals the limit of validity for ξ extracted from the E-S VRH law. Indeed for larger ξ we expect screening of the interactions and a cross over from E-S to Mott VRH law. This is what is observed as shown by a log plot of T σ xx as a function of 1/ √ T on the upper part of 3. For ξ < 600nm a linear variation is found while for ξ > 600nm the variation is no longer linear and is well fitted by the Mott's law. This yet never observed cross-over from E-S to Mott's law is one of the two main results of our experiment.

Before going further in the study of the Mott's regime, it is important to determine the universal scaling exponent γ in the E-S VRH regime for which ξ can be reliably known. ν -ν C for ν C ≃ ±2 for samples S3 and S4. The linear variation indicates that γ = 2.3 is a reasonable exponent. This result is in agreement with Ref. [20].

The second important result of this work comes from the study of the scaling exponent of the FWHM ∆ν of resistivity peaks between Hall plateaus with both temperature and bias current. Here, as ξ > 2d in this regime, Mott's VRH law is obeyed and we expect the exponent values will allow to discriminate between the E-S VRH and the phase coherence length scenario. The FWHM of resistance peaks are plotted for both S3 and S4 on Fig. 4 on the lower left and lower right part respectively for ν = [±2, ±6], [±10, ±6], [-2, +2]. The figures clearly show a universal behavior of ∆ν at temperature below 100K and bias current below 10µA. κ is found to be equal to 0.23 ± 0.02 and µ to 0.13 ± 0.01 which are both in good agreement with Mott's VRH confirming Mott's law at the edge of the quantum Hall plateaus. If the phase coherence length approach was relevant, we would have found κ = 0.42 and µ = 0.21. ] and ν = [±2, 0] at different magnetic field, the black curve is for S4 while the others are for S3. Lower part: FWHM of the ρxx peaks for temperature measurement (left) and for bias current measurement (right). Blue, green and red curves are data for the peaks between ν = ±10 and ±6, ν = ±6 and ±2, and ν = -2 and 2 respectively. Squares are for S3 and triangle for S4. The solid red and dashed black lines correspond to the expected exponent for E-S scenario and for the A-S scenario in the Mott regime respectively. The Pruisken scenario is also represented by the dashed black lines.

To conclude we have studied the quantum localization in the quantum Hall regime using Graphene monolayers. The standard localization length scaling exponents was found for all studied Landau Levels. More important, our first observation of an Efros-Shklovskii to Mott VRH cross-over in the quantum Hall regime found on conductance peaks allows to discriminate between the Polyakov-Shklovskii and the Pruisken phase coherence length scenarios describing the conduction on the plateau transition.

FIG. 1 :

 1 FIG. 1: Longitudinal resistance ρxx versus filling factor at 16.5T. On the left the temperature is varied from 4.2K to 300K and on the right the bias current from 10nA to 100µA at 4.2K. Inset: effective temperature T ef f (I) (see definition in text) versus bias current I for different filling factors on plateaus.

FIG. 2 :

 2 FIG.2: Square root of E-S VRH temperature T0 versus filling factor for sample S3 at different magnetic fields ranging from 6.6T to 16.5T. The inset shows ξmin as a function of the magnetic field. Green, red, and blue data are for ν = ±2,±6, and ν = ±10 respectively. Squares are for sample S3 and triangles for S4. The solid curve corresponds to lc.

Fig. 4 ,FIG. 3 :

 43 FIG.3: Lower graph: localization length ξ versus filling factor for S3 for different magnetic fields. The horizontal solid line ξ = 600nm shows the threshold above which screening of interactions by the gate plays a role. Upper figure: plots of T ρxx versus 1/ √ T at 14.8T for various filling factors near ν = -6 and ξ above and below 600nm. The red solid curve is a E-S law's fit of the data while the green solid curve uses Mott's law.

3 FIG. 4 : 3 0

 343 FIG. 4: Upper figure: T -2.3 0as a function of ν -νC for ν = [±2, ±6] and ν = [±2, 0] at different magnetic field, the black curve is for S4 while the others are for S3. Lower part: FWHM of the ρxx peaks for temperature measurement (left) and for bias current measurement (right). Blue, green and red curves are data for the peaks between ν = ±10 and ±6, ν = ±6 and ±2, and ν = -2 and 2 respectively. Squares are for S3 and triangle for S4. The solid red and dashed black lines correspond to the expected exponent for E-S scenario and for the A-S scenario in the Mott regime respectively. The Pruisken scenario is also represented by the dashed black lines.

We thank J-N Fuchs and M.Goerbig for fruitful discussions as well as K. Novoselov for help in the early stage of the project. We thank the Quantronics group and especially D. Vion for help in nanofabrication. The CNano IdF grant JoseqhGraphn is acknowledged.


Sandrine Chaignon

Laetitia Virlouvet

Marie Javelle

Valérie Méchin

Marie-Pierre Jacquemot

Thierry Balliau

Charlotte Henriet

Llorenç Cabrera-Bosquet

Sandrine Balzergue

Etienne Delannoy

Johann Joets

Michel Zivy

Marie-Laure Martin-Magniette

Claude Welcker

Christophe Sallaud

François Tardieu

Sylvie Coursol

MYB31 triggers water stress memory responses in maize

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Water deficit directly impacts the ability of plants to intercept and convert light into biomass. Because leaf growth is one of the first processes affected by water deficit, many physiological studies concentrated in short-term responses and associated mechanisms. They demonstrated the roles of cellular and metabolic processes such as changes in cell turgor, hydraulic conductance and cell wall plasticity. However, our understanding of how water deficit impacts cell wall biosynthesis is still fragmentary. Here, we report that ZmMYB31, a R2R3-MYB transcription factor implicated in the regulation of lignin biosynthesis [1][2], is induced by water deficit in the growing zone of maize leaves. Consistent with this finding, we also showed that ZmMYB31 colocates with quantitative trait loci for growth responses to water deficit and anthesis-silking interval (ASI) under water deficit conditions. To assess the potential contribution of ZmMYB31 to maize responses to water deficit, we generated transgenic maize overexpressing ZmMYB31. We showed that increasing ZmMYB31 expression in maize maintained ASI under water deficit conditions. On the basis of these data, we generated transcriptomic and proteomic data to detect transcripts and proteins that show substantial changes in abundance in the leaf growing zone of three representative independent transgenepositive homozygous maize lines compared with that of transgene-negative sibling plants grown in PhenoArch under well-watered and water deficit conditions. We identified 352 genes that were transcriptionally regulated by ZmMYB31 under well-watered conditions and displayed opposite transcription levels under water deficit conditions, suggesting that ZmMYB31 may function in a water stress memory response pathway, which likely helps the transgenic maize enduring the water deficit stress. Our studies further revealed co-expression relationships between these genes and others from the proteome and a cell wall data sets, suggesting potential regulatory influence.

Funding acknowledgment: French National Research Agency (ANR-10-BTBR-01-Amaizing), FranceAgrimer and LabEx Saclay Plant Sciences-SPS (ANR-10-LABX-0040-SPS).

Références bibliographiques