Nicolas Crouseilles
email: crouseil@math.unistra.fr

Ahmed Ratnani
email: ratnani@math.unistra.fr

Eric Sonnendrücker

An Isogeometric Analysis Approach for the study of the gyrokinetic quasi-neutrality equation

Keywords: isogeometric analysis, NURBS, gyrokinetic quasi-neutrality equation, numerical simulations

In this work, a new discretization scheme of the gyrokinetic quasi-neutrality equation is proposed. It is based on Isogeometric Analysis; the IGA which relies on NURBS functions, seems to accommodate arbitrary coordinates and the use of complicated computation domains. Moreover, arbitrary high order degree of basis functions can be used. Here, this approach is successfully tested on elliptic problems like the quasi-neutrality equation.

Introduction

Nowadays, the modeling of magnetized plasmas is a key issue for controlled thermonuclear fusion. In practice, the study of such plasmas requires solving the Maxwell equations coupled to the computation of the plasma response. Different ways are possible to compute this response: the fluid or the kinetic description. Obviously solving the full Vlasov equation involves the discretization of the six-dimensional phase space, which is a challenging problem. On the other hand, the fluid approach seems to be insufficient when one wants to study the behavior of zonal flow, or the interaction between waves and particles for example (see [START_REF] Dimits | Comparisons and physics basis of tokamak transport models and turbulence simulations[END_REF][START_REF] Hammett | Fluid models for Landau damping with application to the ion-temperature-gradient instability[END_REF]).

In the context of strongly magnetized plasmas however, the motion of the particles is particular since it is confined around the magnetic field lines; the frequency of this cyclotron motion is faster than the frequencies of interest. Hence, averaging the Vlasov equation over the cyclotron motion reduces the dimensionality, and numerical simulations, even if they remain very costly, become possible (see [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Hatzky | Energy conservation in a nonlinear gyrokinetic particlein-cell code for ion-temperature-gradient-driven modes in θ-pinch geometry[END_REF]). These simulations are performed using particles methods or a phase space grid (Eulerian methods) but the computation of the electric potential is always performed on a physical 3D grid. Moreover, the configuration of a tokamak is such that the physics is highly anisotropic and structures along the magnetic field lines are quite larger than across the magnetic field lines. In order to reduce the numerical effort which is huge anyway, it is quite important to use this information in order to define the grid resolution in each direction. Indeed, if the grid is aligned or almost aligned with the magnetic field lines, it is possible to use a lot fewer points in the parallel direction than in the transverse direction.

It is of great importance to develop a quasi-neutrality solver that is flexible with respect to geometry and that can provide high order accuracy. In this work, we are investigating an approach which can accommodate arbitrary coordinates and a complicated geometry of the computation domain. In [START_REF] Czarny | Bézier surfaces and finite elements for MHD simulations[END_REF] Czarny and Huysmans used Bézier elements for MHD simulations. Bézier surfaces are the most basic tool in (CAD) computer aided design. However, it cannot preserve the exact geometry. The Isogeometric Analysis (IGA), which has been introduced recently by Hughes et al. [START_REF] Hughes | Analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement[END_REF], seem to provide all these features. The IGA relies on NURBS functions, which are a generalization of Spline functions and provides an exact modeling of large classes of computational domains including conics and all spline surfaces. Moreover they rely on a cartesian grid of the parameter space and are fairly easy to use even using spline basis functions of arbitrary degree. Moreover for domains that can be represented using a periodic angular variable as is the case for the poloidal plane of the tokamak, we were able to develop a fast solver that is comparable in computation time, for any spline degree, with the spectral and Finite Difference solver used in Gysela [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF].

The results presented in this work use a bean shaped domain corresponding to a poloidal cut of the tokamak plasma. The next step will be to couple this field solver with a gyrokinetic solver. Moreover, being based on the projection of the approximated function into a finite dimensional space, it can not only deal with differential operators but also integral operators such as those involved in an exact computation of the double gyroaverage, which is up to now approximated by the transverse Laplacian.

Finally, let us mention a point which is particularly useful for parallel computations. Due to the adiabatic assumption for electrons, a nonlocal term intervenes in the equation which is very penalizing for massive parallelization of full gyrokinetic codes. In this work, we propose an algorithm which enables an interesting decoupling of the quasi-neutrality equation. Indeed, by decomposing the electric potential between its average on a magnetic surface and the difference between itself and this average, it is possible to solve the quasi-neutrality equation by: first, solving a 1D radial ordinary differential equation, and second, solving N ϕ poloidal 2D equations (where N ϕ is the number of poloidal planes). This latter equations are solved using the NURBS approach. The two equations are local and seems suitable for massively parallel computations.

The rest of the paper is organized as follows: After some recalls about the quasi-neutral approximation of Poisson's equation in gyrokinetics, the IGA approach is presented and applied to simple examples and the fast solver is described in the case of polar coordinates. Then, the last part is devoted to the numerical solution of the quasi-neutrality equation satisfied by the electric potential.

Quasi-neutrality equation

The Poisson equation enables to determine the electric potential φ as a function of the distribution function

∇ 2 φ = -4π|e|(n i -n e ).
where n i (resp. n e ) stands for the ion density (resp. the electron density). In classical tokamak plasmas, the Debye length is one order smaller than the Larmor radius so that (see [START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF]), the quasi-neutrality equation is given by

n i (x) = n e (x), (2.1)
The ion density n i is evaluated at particles position, and can be computed from the solution of the gyrokinetic equation (posed on guiding-center position). At the first gyrokinetic order, n i can be written

n i (x) = J (f + g)(x, v)dv, J (f )(x) = 1 2π 2π 0 f (x + ρ)dϕ, ρ = |ρ|(cos ϕ, sin ϕ), (2.2)
where f stands for the gyrocenter distribution and g is the first order correction, which can be approximated (as in [START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF][START_REF] Lee | On high-order corrections to gyrokinetic Vlasov-Poisson equations in the long wavelength limit[END_REF][START_REF] Hahm | Nonlinear gyrokinetic equations for tokamak microturbulence[END_REF]) by

g = ∂ µ F M (φ -J (φ)), with F M = n 0 √ 2πT exp(-mµ/(2T )),
where n 0 = n 0 (r) is an equilibrium density and T = T (r) is an equilibrium ion temperature. These two profiles are generally given. Hence, injecting this last expression into (2.2) leads to

n i (x) = ni (x) + n 0 T (φ(x) -φ(x)), (2.3)
where µ)dµ is the second gyroaverage transformation of φ. In Fourier variables, this term has a compact form using the function Γ 0 (b) = +∞ 0 exp(-x)J 2 0 (bx)dx, where J 0 is the Bessel function. Then (2.3) becomes in Fourier variables

ni (x) = J (f )(x, v)dv and φ(x) = ∞ 0 J 2 (φ)(x) exp(-
n i = ni + n 0 φ T i (1 -Γ 0 (b)), with b = k 2 ⊥ ρ 2 i , ρ 2 i = T /B 2 .
By expanding the Gamma function Γ 0 using a Padé approximation (see [START_REF] Dubin | Nonlinear gyrokinetic equations[END_REF][START_REF] Hahm | Nonlinear gyrokinetic equations for tokamak microturbulence[END_REF][START_REF] Lee | On high-order corrections to gyrokinetic Vlasov-Poisson equations in the long wavelength limit[END_REF]), we obtain

n i = ni + 1 B 2 ∇ ⊥ • (n i ∇ ⊥ φ) .
For the electron density, an adiabatic assumption is often performed so that we can write the following equality

n e (x) = n 0 + n 0 T e (φ -φ ),
where the operator • is an integration over constant magnetic surfaces. The quasi-neutrality equation then reads,

- 1 B 2 ∇ ⊥ • (n i ∇ ⊥ φ) + n 0 T e (φ -φ ) = ni -n 0 .
Linearization around the equilibrium density n 0 together with the approximation B ≈ B 0 are usually performed (see [START_REF] Hahm | Nonlinear gyrokinetic equations for tokamak microturbulence[END_REF][START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]); these simplifications provides only a radial dependence for the anisotropic factor

- 1 B 2 0 ∇ ⊥ • (n 0 ∇ ⊥ φ) + n 0 T e (φ -φ ) = ni -n 0 . (2.4)
In the rest of the paper, the equation (2.4) is intended to be solved.

NURBS

Let us first recall some important properties on B-splines and their generalization to NURBS. These will be used to define the computational domain and also as basis functions for our Finite Element formulation.

Splines

Let T = (t i ) 1 i N +k be a non-decreasing sequence of knots.

Definition 1 (B-Spline) The i-th B-Spline of order k is defined by the recurrence relation:

N k j = w k j N k-1 j + (1 -w k j+1 )N k-1 j+1
where,

w k j (x) = x -t j t j+k-1 -t j N 1 j (x) = χ [tj ,tj+1[ (x)
We note some important properties of a B-splines basis:

• B-splines are piecewise polynomial of degree p = k -1

• Positivity • Compact support; the support of N k j is contained in [t j , .., t j+k ]
• Partition of unity :

N i=1 N k i (x) = 1, ∀x ∈ R • Local linear independence • If a knot t has a multiplicity m then the B-spline is C (p-m) at t Let (P i ) 1 i N ∈ IR d be
a sequence of control points, forming a control polygon.

Definition 2 (B-Spline curve) The B-spline curve in IR d associated to T = (t i ) 1 i N +k and (P i ) 1 i N is defined by :

M(t) = N i=1 N k i (t)P i

Fundamental geometric operations

Once a B-spline curve is defined by a number of control points N , a degree k and a set of knots T . It is possible to define the same curve using more knots or a higher degree. Mesh refinement starting from a given spline or NURBS curve for the boundary can be performed in this way. Let us briefly sketch the algorithms for performing these operations. After modification, we denote by N , k, T the new parameters. (Q i ) are the new control points.

Knot insertion

One can insert a new knot t, where t j t < t j+1 . For this purpose we use the DeBoor algorithm:

N = N + 1 k = k T = {t 1 , .., t j , t, t j+1 , .., t N +k } α i =    1 1 i j -k + 1 t-ti t i+k-1 -ti j -k + 2 i j 0 j + 1 i Q i = α i P i + (1 -α i )P i-1

Order elevation

We can elevate the order of the basis, without changing the curve. Several algorithms exist for this purpose. We used the one by Huang et al. [START_REF] Huang | Fast degree elevation and knot insertion for B-spline curves[END_REF].

k = k + m m i = m i + m N = N + ms
Differential coefficients are defined as P l i :

P l i =      P i l = 0 1 ti+p-t i+l ( P l-1 i+1 -P l-1 i ) l > 0, t i+k-1 > t i+l 0 l > 0, t i+k-1 = t i+l β i = i l=1 m l , 1 i s -1, et α i = i l=1 k-1-l k-1+m-l , 1 i k -2
We present the algorithm by [START_REF] Huang | Fast degree elevation and knot insertion for B-spline curves[END_REF]:

1. Compute P j 0 , 0 j k -1 et P i β l ,1 l s -1 , k -m l i k -1 2. Compute Q j 0 = j l=1 ( k-l k+m-l ) P j 0 , 0 j k -1 3. Compute Q j β l +ml = j l=1 ( k-l k+m-l ) P j β l , 1 l s -1 , k -m l i k -1 4. Compute Q k-1 β l +ml+i = Q (k-1)
β l +ml , 1 l s -1 , 1 i m 5. Compute Q 0 i
Note that there exist other algorithms such those given by (see [START_REF] Prautzsch | A fast algorithm to raise the degree of B-spline curves[END_REF][START_REF] Piegl | The NURBS Book[END_REF] and others). The one given in [START_REF] Huang | Fast degree elevation and knot insertion for B-spline curves[END_REF] is more efficient and much more simple to implement. We can also use a more sophisticated version of this algorithm to do the insertion of new knots while elevating the degree.

Refinement strategies

Refining the grid can be done in 3 different ways. This is the most interesting aspects of B-splines basis.

• using the patch parameter h, by inserting new knots. This is the h-refinement, it is the equivalent of mesh refinement of the classical finite element method.

• using the degree p, by elevating the B-spline degree. This is the p-refinement, it is the equivalent of using higher finite element order in the classical FEM.

• using the regularity of B-splines, by increasing / decreasing the multiplicity of inserted knots. This is the k-refinement. This new strategy does not have an equivalent in the classical FEM.

J.A. Evans et al. [START_REF] Evans | n-Widths, supinfs, and optimality ratios for the k-version of the isogeometric finite element method[END_REF] studied the k-refinement using the theory of Kolmogorov n-widths. As we will see in this article, the use of this strategy can be more efficient than the classical p-refinement, as it reduces the dimension of the basis.

NURBS

Let ω = (ω i ) 1 i N be a sequence of non-negative reals. The NURBS functions are defined by a projective transformation:

Definition 3 (NURBS) The i-th NURBS of order k associated to the knot vector T and the weights ω, is defined by

R k i = ω i N k i N j=1 ω j N k j .
Notice that when the weights are equal to 1 the NURBS are B-splines.

Definition 4 (NURBS curve)

The NURBS curve of order k associated to the knot vector T , the control points (P i ) 1 i N and the weights ω, is defined by

M(t) = N i=1 R k i (t)P i Definition 5 (NURBS surface)
The NURBS surface of order k associated to the knot vectors {T (1) , T (2) }, the control points (P i,j ) 1 i N1,1 j N2 and the weights {(ω (1) , ω (2) }, is defined by

M(t (1) , t (2) ) = N1 i=1 N2 j=1
R i,j (t (1) , t (2) )P i,j

with R i,j (t (1) , t (2)

) = R (1)
i (t (1) )R

(2) j (t (2) ) Remarks NURBS functions inherit all B-splines properties. Remark that in the interior of a knot span, all derivatives exist, and are rational functions with non vanishing denominator. We present here the definition of the perspective mapping. We construct the weighted control points

P ω i = (ω i x i , ω i y i , ω i z i , ω i ).
Then we define the B-spline curve in four-dimensional space as

M ω (t) = N i=1 N k i (t)P ω i .
For fundamental geometric operations on NURBS curves, we use the latest transformation and algorithms on B-spline curves. NURBS functions allow us to model, exactly, much more domains than B-splines. In fact, all conics can be exactly represented with NURBS. For more details, see [START_REF] Piegl | The NURBS Book[END_REF].

Examples

Circle

One can draw a circle using only 9 control points, and the parameters:

N = 9, p = 2 , T = {000, 1 4 1 4 , 1 2
i P i ω i 1 (1, 0) 1 2
(1, 1)

1 √ 2 3 (0, 1) 1 4 (-1, 1) 1 √ 2 5 (-1, 0) 1 6 (-1, -1) 1 √ 2 7 (0, -1) 1 8 (1, -1) 1 √ 2 9
(1, 0) 1

2 circles

To draw a circle inside another circle we need the following parameters:

N 1 = 9, p 1 = 2, T (1) = {000, 1 4 1 4 , 1 2
1 2 , 3 4
3 4 , 111}
and

N 2 = 3, p 2 = 1, T (2) = {00, 1 2 , 11}.

Ellipse

To draw an ellipse we can change the scaling of control points coordinates.

Modeling a Tokamak

Starting with an ellipse, one can modify the control points and the weights to have a simplified model of tokamak. We can approximate a domain, using data points. We can refer to the works of [START_REF] Ma | NURBS Curve and Surface Fitting for Reverse Engineering[END_REF][START_REF] Heidrich | Fitting Uncertain Data with NURBS[END_REF].

The IGA approach for elliptic problems

The first step to use NURBS as finite element basis functions, is to define the patch, i.e the parametric domain, and therefore the mapping F , that maps the patch into the physical domain.

Patch

When the edge can be described as a NURBS curve, we can construct an initial patch, by adding some control points that rise from the definition of the knot vectors used. From this coarse mesh, we can then, use h/p/k refinements to create the grid. We can also use multiple patches to describe more complex domains [START_REF] Cottrell | Isogeometric Analysis, toward Integration of CAD and FEA[END_REF]. NURBS models suffers from the fact that control points must lie in a rectangular grid. Hence, we will have a lot of superfluous control points, that might exist only to satifsy this constraint. Sederberg et al. [START_REF] Sederberg | T-spline Simplification and Local Refinement[END_REF] defined the notion of T-splines that allow us to reduce the number of those control points. In [START_REF] Dörfel | Adaptive isogeometric analysis by local h-refinement with T-splines[END_REF] Dörfel et al. used T-splines for local h-refinement in isogeometric analysis.

Grid generation

For this purpose, we use alternatively h and p-refinement. The minimal degree of the basis functions is imposed by the domain design. When inserting knots, we can use uniformly-spaced knots or non uniformlyspaced ones.

Variational formulation of an elliptic problem

A general elliptic problem writes

   -∇ • (A∇u) + cu = f, Ω u = g, Γ D ∇u • n = k, Γ N (4.5)
where

Γ D Γ N = ∂Ω and Γ D Γ N = ∅.
A is some given tensor and c a given nonnegative function.

Let define the spaces

X = H 1 (Ω), S = {v/v ∈ H 1 (Ω), v| Γ D = g}, and V = {v/v ∈ H 1 (Ω), v| Γ D = 0}.
The variational formulation of (4.5) is:

Find u ∈ S such that: a(w, u) = l(w), ∀w ∈ V (4.6) where a(w, v) = Ω A∇v • ∇w + Ω c v w, and
l(w) = Ω f w + Γ N kw.

The discrete variational formulation

Let X h , V h be finite linear sub-spaces of X , V, where h is intended to tend to zero. S h is a finite dimensional approximation of S.

The discrete variational formulation is

Find

u h = v h + g h ∈ S h such that: a(w h , u h ) = l(w h ), ∀w h ∈ V h , (4.7)
where g h , h are given, and g h ∈ S h . By re-ordering the basis functions (NURBS) to have : R i | Γ D = 0, ∀i ∈ {1, .., nn D }, we then obtain a construction of a basis of V h i.e. :

∀w h ∈ V h , w h = n-n D i=1 [w h ] i R i g h ∈ S h such that [g h ] i = 0, ∀i ∈ {1, .., n -n D } and so g h = n i=n-n D +1 [g h ] i R i . Then, u h ∈ S h is of the form u h = n-n D i=1 [u h ] i R i + n i=n-n D +1 [g h ] i R i .
We obtain classically the linear system:

Σ [u h ] = L,
where the following notations have been used

Σ = (Σ i,j ) 1 i,j n-n D , L = (L i ) T 1 i n-n D , [u h ] = ([u h ] i ) T 1 i n-n D .

Computing element integrals

Let Q be a cell in the physical domain. Q is the parametric associated cell such that Q = F ( Q). let J F be the jacobian of the transformation F , that maps any parametric domain point (ξ, η) into physical domain point (x, y). For any function v of (x, y) we associate its representation in the parametric domain

v((ξ, η)) := v • F ((ξ, η)) = v((x, y)).
The basis functions R i will not be affected by these changes, the reader can always know if the we are working in the physical or parametric domain thanks to (x, y) = F (ξ, η), x = α(ξ, η) and y = β(ξ, η).

We then have, for v, w defined on

F ( Q) = Q, with J F -1 the jacobian of F -1 , Q (A∇v) • ∇w dx dy = Q=F -1 (Q) (AJ F -1 ∇ v) • (J F -1 ∇ w) 1 det J F -1 dξ dη = Q ∇ v T J T F -1 A T J F -1 ∇ w det(J F ) dξ dη. (4.8)
Introducing the notations

α 1 = ∂α ∂ξ , α 2 = ∂α ∂η , β 1 = ∂β ∂ξ , β 2 = ∂β ∂η ,
we have for the determinant of the Jacobian det(J

F ) = α 1 β 2 -α 2 β 1 and for J F -1 J F -1 = 1 det(J F ) β 2 -β 1 -α 2 α 1 . Now, let us define Θ = J T F -1 A T J F -1 , so that (4.8) becomes Q (A∇v) • ∇w dx dy = Q ∇ v T Θ ∇ w det(J F ) dξ dη.
When A is the identity matrix the explicit form of the Θ matrix writes:

Θ = 1 det(J F ) 2 α 2 2 + β 2 2 -α 1 α 2 -β 1 β 2 -α 1 α 2 -β 1 β 2 α 2 1 + β 2 1 .
For the term of order zero, in the left hand side, we have :

Q cvw dx dy = Q c v w det(J F ) dξ dη.
For the right hand side part, the computations are easier and lead to

Q f v dx dy = Q f v det(J F ) dξ dη,
whereas for the Neumann part, we have

∂Q∩Γ N k v dΓ = ∂ Q∩ Γ N k v det(J F |Γ )d Γ.
In practice, we use Gauss-Legendre quadrature to compute those integrals. We give here the number of used quadrature points, per direction, with respect to the spline degree. Notice that, Hughes et al. [START_REF] Hughes | Efficient quadrature for NURBS-based isogeometric analysis[END_REF], began the study of efficient numerical quadrature formulae for NURBS. Hence, we get,

Σ i,j = Ω A∇R i • ∇R j dx dy + Ω c R i R j dx dy = Q Q ∇ R T i Θ ∇ R j det(J F ) dr dθ + Q c R i R j det(J F ) dr dθ , (4.9)
and

L i = Ω f R i + Γ N hR i = Q Q f R i det(J F ) dξ dη + ∂ Q∩ Γ N k R i det(J F |Γ )d Γ . (4.10)
As one can see, the matrix Σ is the sum of the stiffness matrix (terms Ω A∇R i • ∇R j ), and the mass matrix (terms Ω c R i R j ).

A fast solver for polar coordinates

Let us now see how the specific structure of our problem allows us to derive a specific solver which is a lot faster than using a generic sparse matrix solver. For our quasi-neutrality equation, the mapping F is the mapping defining polar coordinates F (r, θ) = (r cos θ, r sin θ). Then we have det(J F ) = r. Moreover the matrix A is of the form a(r)I and c ≡ c(r) is only a function of r. The matrix Θ becomes

Θ = a(r) 0 0 a(r) r 2 .
On the other hand the basis functions can be written as products of functions of r only and of θ only. Let us also introduce a numbering using the two indices of our logical patch Q in (r, θ). Thus the degree of freedom with index i will be associated with the grid index (i r , i θ ) and the basis function R i is such that R i (r, θ) = R ir (r) R i θ (θ). Plugging this into the element integrals defining the matrix Σ in (4.9) we get

Q ∇ R T i Θ ∇ R j det(J F ) dξ dη + Q c R i R j det(J F ) dξ dη = Q a(r) ∂ R i ∂r ∂ R j ∂r + a(r) r 2 ∂ R i ∂θ ∂ R j ∂θ r dr dθ + Q c(r) R i R j r dr dθ = a(r) R ′ ir (r) R ′ jr (r) r dr R i θ (θ) R j θ (θ) dθ + a(r) r 2 R ir (r) R jr (r) r dr R ′ i θ (θ) R ′ j θ (θ) dθ + c(r) R ir (r) R jr (r) r dr R i θ (θ) R j θ (θ) dθ. (4.11)
These formulas lead us to two observations. First the elementary matrices decouple into products of integrals in r and integrals in θ so that the final matrices can be written in a Kronecker product structure (explanations and applications of this structure can be found in the article by Van Loan [START_REF] Van Loan | The uniquitous Kronecker product[END_REF] and references therein). The exploitation of this structure for developing fast solvers was developed in [START_REF] Lynch | Direct solution of partial difference equations by tensor product methods[END_REF]. Second the matrices in θ are simple mass and stiffness matrices with no varying parameter inside, so that if the mesh is uniform in θ the matrix will be circulant. One way to express the Kronecker product structure is to write the unknown degrees of freedom and the right-hand-side as matrices where the terms correspond to the indices (i r , i θ ). We denote those respectively by U and F . Then, in our case, the linear system can be written

K ar U M θ + M ar U K θ + M cr U M θ = F, (4.12)
where K ar is the weighted stiffness matrix in r corresponding to the terms a(r) R ′ ir (r) R ′ jr (r) r dr, M ar is the weighted mass matrix in r corresponding to the terms a(r) r 2 R ir (r) R jr (r) r dr, M cr is the weighted mass matrix in r corresponding to the terms c(r) R ir (r) R jr (r) r dr, K θ is the stiffness matrix in θ corresponding to the terms R ′ i θ (θ) R ′ j θ (θ) dθ and M θ is the mass matrix in θ corresponding to the terms R i θ (θ) R j θ (θ) dθ. The latter two matrices K θ and M θ are circulant, which means that they can be both diagonalized in the same orthonormal basis corresponding to the Fourier modes. This can be expressed by

M θ = P Λ M P * , K θ = P Λ K P * ,
where Λ M and Λ K are the diagonal matrices of the eigenvalues and a multiplication by P corresponds to the normalized Fast Fourier Transform and a multiplication by P * to its inverse. This can be exploited for the fast solution of (4.12) bringing the solution of a linear system arising from a 2D problem to sets of smaller systems corresponding to 1D problems. The procedure can be performed with the following algorithm: 1. Multiply system (4.12) on the right by P (amounts to a 1D FFT on each line of F ). Then we get

K ar U P Λ M + M ar U P Λ K + M cr U P Λ M = F P. (4.13)
2. Note that a multiplication on the right by the diagonal matrix of eigenvalues corresponds to multiplying each column of the matrix by the corresponding eigenvalue, which implies that (4.13) corresponds to uncoupled problems on each of the columns of Û = U P . So denoting by Û1 , . . . , Ûn the columns of the matrix Û and by F1 , . . . , Fn the columns of the matrix F = F P , (4.13) becomes for each column j,

(λ Mj (K ar + M cr ) + λ Kj M ar ) Ûj = Fj (4.14)
This is a set of banded systems of size the number of points in the r direction, that can be solved very efficiently using the LAPACK routines DPBTRF for the Cholesky factorization that is only called once at the beginning and then DPBTRS for the solution at each time step.

3. Compute U = Û P * by inverse FFT of the lines of Û .

Let us now compute the cost at each time step for a N r × N θ mesh, disregarding the cost of the Cholesky factorization for the systems in r which needs only to be performed once for a many time steps computation. The algorithm consists of three steps: 1) N r FFTs which need O(N θ log 2 N θ ) operations, 2) N θ up and down sweeps of a Cholesky decomposed banded system which cost O(N r ) each, 3) N r inverse FFTs which cost O(N θ log 2 N θ ) operations. So all together the cost is O(N r N θ log 2 N θ ) operations, which is almost optimal. This algorithm uses the structure of the system in an optimal manner and only works on dense matrices. A generic sparse systems solvers could not do this.

Numerical results

We have tested this new approach to solve an elliptic partial differential equation, using the analytic solution u(r, θ) = sin(2πr) sin(2πθ), which solves :

-∇ 2 u(r, θ) + u(r, θ) = F (r, θ), Ω (4.15)
where Ω = [0, 1] × [0, 1], ∇ 2 is the Laplacian cartesian. The boundary conditions are :

u(r = 0, •) = u(r = 1, •) = 0 (4.16)
and periodic boundary condition on θ.

In Figures 8 and9, we compare the time CPU spent to solve the linear system, using this approach, spsolve, from scipy based on SUPERLU solver, for the classical formulation (detailed in the sections before). The test was done on grids 128 × 128 and 256 × 256, using different spline order. As one can see, the new method does not really depend on the spline degree, the reason is that, for 1D problem, the bandwidth is equal the spline order, whereas it is a quadratic function of the order, in the classical approach. On a grid of 512 × 512, we spent 0.3 sec to solve the linear system, using cubic splines. Using splines of order 8, it took only 0.321 sec.

Numerical validation

In all tests we have treated, we consider an elliptic problem under Dirichlet boundary condition, i.e u(x) = 0, ∀x ∈ ∂Ω.

Spline degree

Test case 1: Order of convergence for Poisson in polar coordinates

For the validation, we solved the polar coordinate elliptic problem, which fits into our generic elliptic problem (4.5) with A the identity tensor and c = 0,

-∇ 2 r,θ φ(r, θ) = n(r, θ),
with the following analytical solution:

φ(r, θ) = sin π r 2 max -r 2 min (r 2 -r 2 min ) .
Injecting the solution into the elliptic problem enables to compute the right hand side n which is given in input of the code. This test enables to check the L 2 norm of the error (in log scale) between the numerical and analytical solution, with respect to the parameter h = max{diam( Q)}, for different orders of the basis functions (from order 2 to order 6). This is shown in Figure 10. We verify that the slopes of the different curves correspond to the order of the basis functions. In particular, for high orders, the machine precision is achieved. We also plot in Figure 10 the CPU time as a function of the parameter h = max{diam( Q)}, for different orders of the basis functions (from order 2 to order 6). These two last figures gives information for choosing a priori the best compromise between precision (which order of the basis function should be chosen) and the CPU time.

Test case 2: Chaotic solution

Following [START_REF] Nishimura | A finite element Poisson solver for gyrokinetic particle simulations ina global field aligned mesh[END_REF], we first test our solver on a chaotic solution on a polar coordinate Laplacian -∇ 2 φ = n. This analytic solution takes the form

φ math (r, θ) = sin(2πξ) + ǫ M B M sin(2πM ξ) l A l cos(lθ + Θ l ), (5.17)
where 0 ≤ ǫ ≤ 1, ξ = (rr min )/(r maxr min ), A l and B M are random numbers which range between 0 and 1, with |M |, |l|(≤ 20 for the first simulation and ≤ 40 for the second one). Finally, the phase Θ l is also given by a random number in 0, 2π. The right hand side of -∇ 2 φ = n is given by

n(r, θ) = - 1 r ∂ r (r∂φ math ) - 1 r 2 ∂ 2 θ φ math = - 1 r σ r (r, θ) l A l cos(lθ + Θ l ) + sin(2πξ) + ǫ M B M sin(2πM ξ) l l 2 r 2 A l cos(lθ + Θ l ), with σ r (r, θ) = 2π ∆r cos(2πξ) - 4π 2 r ∆r 2 sin(2πξ) + ǫ M B M 2πM ∆r cos(2πM ξ) - 4π 2 M 2 r ∆r 2 sin(2πM ξ)
In our numerical experiments, ǫ is taken equal to 0.4. Note that the solution satisfies the homogeneous Dirichlet condition at r = r min = 0.2 and r = r max = 0.4 and periodic boundary conditions in the θ direction. We show on Figure 11, the L 2 norms (in log scale) of the difference between the analytical and numerical solutions as a function of the parameter h = max{diam( Q)}. The same observations as before are also available for this test; the slopes of the curves corresponds to the order of the basis, even when a very large number of modes is considered in the solution. In this kind of test high order basis functions enable to better capture the very fine scales of the solution. Comparing with Nishimura results [START_REF] Nishimura | A finite element Poisson solver for gyrokinetic particle simulations ina global field aligned mesh[END_REF], we see in Table 1 that the new method allows us to reach same error order, with much fewer degrees of freedom, by increasing the spline order. Let us notice, that to reach an error of 10 -6 , Nishimura used a grid of N = 154, 560, with our method, we only need 8 times less (N = 20064) degrees of freedom by using sixth or seventh order for spline basis.

Spline degree

Numerical solution of the quasi-neutrality equation

This section is devoted to the numerical approximation of the quasi-neutrality equation (2.4) by applying the method introduced above. To that purpose, we will consider a given right hand side (that is to say the distribution function of the ions is given and the coupling with the gyrokinetic equation is not considered), and we focus on the computation of the solution of (2.4) in polar coordinates:

Given n 0 (r), T e (r) two radial profiles and F (r, θ, ϕ), solve for φ(r, θ, ϕ) satisfying

-∇ ⊥ • (n 0 ∇ ⊥ φ) + n 0 T e (φ -φ ) = F (r, θ, ϕ). (6.18)
The • operator refers to the magnetic flux average

φ (t, r) = 1 J(r, θ)dθdϕ φ(t, r, θ, ϕ)J(r, θ)dθdϕ, (6.19)
with J(r, θ) a jacobian defining the poloidal geometry, whereas the given right hand side F (r, θ, ϕ) reads

F (r, θ, ϕ) = (n i (r, θ, ϕ) -n 0 (r)) ,
where ni is the ion density, T e the electronic temperature and n 0 is a correction term taking into account the presence of the electrons. The first term on the left hand side is known as the polarization term which corresponds to the difference between the guiding-center density and that of particles. The boundary conditions are supposed to be 2π-periodic in the θ, ϕ variables, whereas the radial boundary conditions are imposed by the Dirichlet condition and writes φ(r = r min , θ, ϕ) = 0 and φ(r = r max , θ, ϕ) = 0 .

The main goal of this section consists in the numerical solution of (6.18). More precisely, we want to derive an efficient method for (6.18), that is to say a method which is as local as possible. Indeed, the principal difficulty is the average term φ which is totally nonlocal since it couples every values of θ and ϕ. The nonlocality is an important obstacle for an efficient parallelization, but also when one wants to use Fourier transforms in the θ, ϕ variables. We propose a decoupling approach that enables to decompose the solution of (6.18) into two local problems, one of them reduces to a two-dimensional elliptic type problem of the form (4.5).

The rest of this section presents different ways to solve the quasi-neutrality equation, pointing out the advantages and disagrements of the solvers. Then, some numerical results are shown to compare the performance of the solvers.

The decoupling approach

In this subsection, we propose a new method to allow us the derivation of a local algorithm for (6.18) in the ϕ variable. An efficient algorithm can then be developed for its resolution. To this purpose, we first restrict in this section to a Jacobian such that J(r, θ) = but the computations for arbitrary Jacobians J(r, θ) are performed in Appendix A. The main ingredient consists in the decoupling of (6.18) into two local equations: one equation on φ which is given by (6. [START_REF] Ma | NURBS Curve and Surface Fitting for Reverse Engineering[END_REF]), and one equation on Φ = φφ .

The first equation is simply derived by averaging the left hand side of the quasi-neutrality equation (6.18) with respect to θ and ϕ variables

1 (2π) 2 ∇ ⊥ • [n 0 (r)∇ ⊥ φ] dθdϕ = n 0 (r) r + n ′ 0 (r) ∂ r φ + n 0 (r)∂ 2 r φ . (6.20)
Finally, introducing F (r) = 1/(2π) 2 F (r, θ, ϕ)dθdϕ, we have

∇ ⊥ • [n 0 (r)∇ ⊥ φ] = ∇ ⊥ • [n 0 (r)∇ ⊥ φ ] = F , (6.21)
from which we can deduce a 1D equation for φ = φ (r)

-n 0 (r)∂ 2 r φ + n 0 (r) r + n ′ 0 (r) ∂ r φ = F (r), (6.22)
Then, we want to derive an equation satisfied by Φ = φφ . This can be done easily by adding and removing the term

∇ ⊥ • [n 0 (r)∇ ⊥ φ ] in (6.18) -∇ ⊥ • [n 0 (r)∇ ⊥ (φ -φ )] + n 0 T e [φ -φ ] -∇ ⊥ • [n 0 (r)∇ ⊥ φ ] = F (r, θ, ϕ).
By introducing the new unknown Φ = φφ in this last equation, we get

-∇ ⊥ • [n 0 (r)∇ ⊥ Φ] + n 0 T e Φ -∇ ⊥ • [n 0 (r)∇ ⊥ φ ] = F (r, θ, ϕ). (6.23)
Thanks to (6.21), we can write the equation which is satisfied by Φ(r, θ, φ) = φφ

-∇ ⊥ • [n 0 (r)∇ ⊥ Φ] + n 0 T e Φ = F (r, θ, ϕ) -F (r). (6.24)
More precisely, the equation to solve is

-n 0 ∂ 2 r Φ - n 0 r + n ′ 0 ∂ r Φ - n 0 r 2 ∂ 2 θ Φ + n 0 T e Φ = F (r, θ, ϕ) -F (r). (6.25)
Thanks to these straightforward computations, we totally decoupled (6.18) by introducing the new unknowns Φ and φ . Moreover, equation (6.18) is strictly equivalent to the equations (6.22)- (6.24). This new formulation provides a system of equation which does not include nonlocal term any more. The algorithm is then the following Algorithm

• solve the 1D equation (6.22) to get φ (r)

• solve the 2D equation (6.24) ∀ϕ to get Φ(r, θ, ϕ)

• compute φ = Φ + φ
Now, an important point consists in a good choice of the numerical approximation to solve the 2 equations (6.22) and (6.24). A two-dimensional solver will provide Φ whereas a one-dimensional solver for (6.22) leads to φ . By the summation of these two solutions, we can obtain the desired solution φ(r, θ, ϕ).

First approach: spectral + finite differences

The first approach to solve (6.24) consists in the use of a Fourier transform in the θ direction (Φ is periodic in the θ variable since φ is). The projection of Φ(r, θ, ϕ) onto the Fourier space writes Φ(r, θ, ϕ) = N m=1 Φ m (r, ϕ)e imθ . so that equation (6.25) is rewritten

-n 0 ∂ 2 r Φ m + n 0 r + n ′ 0 (r) ∂ r Φ m + n 0 1 T e + m 2 r 2 Φ m = S m (r, φ), (6.26)
where

S m (r, φ) = F (r, ϕ) if m = 0 F 0 (r, ϕ) -1 (2π) 2 F 0 (r, ϕ)dϕ else,
is the Fourier coefficient of order m of the right hand side of (6.24). We are faced to an ODE for each Fourier mode m; one classical way to solve this problem is the use of a finite difference scheme in the r direction to solve (6.26); this leads to the constitution of a tridiagonal linear system, the resolution of which allows to get the Fourier modes Φ m (r, ϕ). An inverse Fourier transform leads to Φ(r, θ, ϕ).

The resolution of the ODE (6.22) can be performed in the same way as (6.26), using finite differences of order two. Remark 6.1 From the CPU time point of view, the present approach enables us to solve the total quasineutrality equation (i.e. to compute φ = φ + Φ) with a cost of order O(N 2 ln(N )) for a N × N grid in (r, θ). Our NURBS based solver has the same complexity with only the bandwith of the systems in r that increase with the degree of the splines.

Second approach: FEM

The second option discretizes directly equation (6.24), that considers two dimensions (r, θ) and a parameter ϕ. The finite element method described in section 3 is performed here.

For the ODE (6.22), a high order method is required not to penalize the high order achieved in the two-dimensional solution. Here, we propose a collocation method introduced in [START_REF] De Boor | Collocation at Gaussian points[END_REF].

Recalling, that the equation 6.24 can be written in the form:

-∇ ⊥ • (A ∇ ⊥ Φ) + c Φ = g (6.27)
where,

A = n 0 0 0 n 0 , c = n 0 T e
, and, g = F (r, θ, ϕ) -F (r). (6.28)

Numerical results

In the present section, we compare the different approximation of the quasi-neutrality equation (6.18) we proposed in the previous section. The different methods are studied on analytic tests, since we impose a right hand side which is consistent with a solution of (6.18).

Test case 1: In this test, we consider the full quasi-neutrality equation (6.18) in which the profile n 0 and T e are supposed constant equal to 1, with the following solution

φ(r, θ, ϕ) = sin 2π L (r -r min ) u(θ, ϕ), θ, ϕ ∈ [0, 2π], r ∈ [0.2, 0.8],
with L = r maxr min and u(θ, ϕ) = 1/π 2 (cos 2 θ cos 2 ϕ). This form of solution imposes the average to be φ (r) = sin 2π L (rr min ) . Injecting this solution in (6.18) leads to an analytical right hand side which is used to recover numerically the true solution. This analytical way enables comparison with the analytical solution. Let us detail the computations leading to the right hand side. First, we compute the three component of

∇ 2 ⊥ φ -∂ 2 r φ = 4π 2 L 2 sin 2π L (r -r min ) u(θ, ϕ), - 1 r ∂ r φ = 2π L cos 2π L (r -r min ) u(θ, ϕ), - 1 r 2 ∂ 2 θ φ = 2 r 2 π 2 sin 2π L (r -r min ) cos(2θ) cos 2 ϕ,
Then, we can compute F which will be used in the right hand side of equations (6.22) and (6.24)

F (r) = 4π 2 L 2 sin 2π L (r -r min ) - 2π rL cos 2π L (r -r min ) .
We also compute F (r, θ, ϕ) for the resolution of (6.24)

F (r, θ, ϕ) = sin 2π L (r -r min ) u(θ, ϕ) 4π 2 L 2 + 1 + 2π L cos 2π L (r -r min ) u(θ, ϕ) + 2 r 2 π 2 sin 2π L (r -r min ) cos(2θ) cos 2 ϕ -sin 2π L (r -r min ) = sin 2π L (r -r min ) u(θ, ϕ) 4π 2 L 2 + 1 + u(r, θ) sin 2π L (r -r min ) - 1 r 2π L cos 2π L (r -r min ) -F (r) -sin 2π L (r -r min )
We then test the decoupling method for which the spectral approach, the IGA-FEM solution is compared to the analytic solution. For the first approach, we fix the number of points in the angular directions θ, ϕ to 64 whereas the number of points in the radial direction is modified to recover the order of the finite difference method. We show in the following figures comparisons between the different approaches we propose. On Figures 12, we plot the L 2 norm errors between the numerical and analytical solution. We first plot the error for the elliptic part (6.24) for which the behaviour is similar to the curves obtained in section 4. Then, we focus on the total error of the solvers for (6.18). We observe that the curves are nearly the same; this is due to the collocation method which is very precise. The error made in this part is then negligeable compared to the error made in the elliptic part. We can also remark that the error for the classical approach spectral + finite differences is of order 2 (due to the second order finite differences). This approach needs to consider a lot of point in the radial direction (2048) to reach a competitive accuracy, whereas the IGA-FEM approach presents very precise results even considering cubic or quartic order.

Test case 2: In this test, some anisotropy is introduced in the Laplacian term by initializing n 0 . Moreover, a profile is imposed to T e . These two functions are solutions of a differential equation. The radial profiles of the ion temperature T e (r) as well as the radial density profile n 0 (r) are deduced by numerical integration of their gradient profile given by 1 T e (r)

dT e (r) dr = -cosh -2 (3(r -r p )), (6.29)
with r p = 0.5. These parameters are employed in [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. Hence, imposing analytical form of the solution, φ(r, θ, ϕ) = C sin ϕ cos(4θ)(rr min ) 6 (r maxr) 6 , (6.30)

with C a constant chosen such that max[φ] = 1, we deduce an expression of the right hand side of (6.18) which is given to the code to compute back an approximation of the solution. The right hand side F is given by

F (r, θ, ϕ) = -∂ 2 r φ -∂ r φ 1 r + n ′ 0 n 0 - 1 r 2 ∂ 2 θ φ + φ T e ,
where the different term are

∂ r φ = 6C sin(4θ) sin ϕ (r -r min ) 5 (r max -r) 6 -(r -r min ) 6 (r max -r) 5 , ∂ 2 r φ = 30C sin(4θ) sin ϕ(r -r min ) 4 (r max -r) 4 (r max -r) 2 - 12 5 (r -r min )(r max -r) + (r -r min ) , ∂ 2 θ φ = -16φ.
As in the previous tests, we are interested in the L 2 norm of the error as a function of the size of the mesh. Results are presented in Table 2 in which the L 2 norm of the error is given as a function of the number of points and the degree of the splines. We also give the results for the standard approach where finite difference and FFT are used. In this test also, the IGA-FEM approach has a very good behaviour since with 32 points per direction; to achieve the same precision the finite difference+FFT approach needs 1024 points per direction. Note a kind of saturation of the error: this is due to the computation of the solution of the ODE 6.29. Indeed, T e and n 0 are numerical solution of an ODE and are not analytical. We verify that when the ODE is solved precisely (by refinement), the correct orders are recovered.

On Figure 13, the solution together with the error between the analytical solution and the numerical one are plotted for the IGA-FEM approach, using degree 2 and 32 points per direction.

Test case 3: non-circular cross section Following the Nishimura idea [START_REF] Nishimura | A finite element Poisson solver for gyrokinetic particle simulations ina global field aligned mesh[END_REF], we have generated a turbulence test over our non-circular tokamak model. We give the result of such test, with l = m = 10 and NURBS degree p = 4, in Figure 14.

For such a domain, the use of the standard finite elements method, would need a great number of degrees of freedom, moreover we couldn't be able to represent exactly the tokamak. Table 2: Test 2 QN : Number of degree of freedom and L 2 norm of the error for some spline degrees and for the method finite-difference+FFT (FD+FFT). In the third column, the ODE (6.29) is solved using 10 3 points whereas in the last column, 10 4 are used. A Appendix: Decoupling approach for J(r, θ) 21

Figure 1 :

 1 Figure 1: (left) A B-spline curve and its control points, (right) B-splines functions used to draw the curve. N = 9, k = 2 , T = {000, 1 4

Figure 2 : 2 ,

 22 Figure 2: Illustration of h-refinement with p = 2, T = {000, 111}, T = {000, 1 2 , 111} and T = {000, 1 2 , 3 4

1 2 , 3 4 3 4 ,

 144 111}. Control points and weights are given in the following table:

Figure 3 :

 3 Figure 3: Domain plot and the exterior control points.

Figure 4 :

 4 Figure 4: Poloidal plane designing of a tokamak. (left) a tokamak model, (right) tokamak after h-refinement.

Figure 5 :

 5 Figure 5: Grid generation. 1 st line: (left) after h-refinement, N 1 = 17, N 2 = 5, (right) after p-raffinement, p 1 = p 2 = 3. 2 nd line: (left) after h-refinement p 1 = p 2 = 3, (right) using unstructured mesh, p 1 = p 2 = 2.

Figure 6 :

 6 Figure 6: Mapping from the patch to the physical domain: (left) initial patch, (right) patch after h-refinement in the η direction

Figure 7 :

 7 Figure 7: The number of used Gauss-Legendre points with respect to the spline degree.

Figure 10 :

 10 Figure 10: Validation test : (left) L 2 error norm, (right) CPU time.

Figure 11 :

 11 Figure 11: Nishimura test: L 2 error norm as a function of h; (left) for 20 modes, (right) 40 modes.

Figure 12 :

 12 Figure 12: Test case 1: L 2 error norm, (left) for the elliptic part, (right) for the global problem ODE + FEM.

  Spline degree Degrees of freedom L 2 error norm L 2 error norm "refined"

Figure 13 :

 13 Figure 13: Test 2 QN: numerical solution with NURBS (left) and the difference with the analytical solution (right). 32 points are used per direction, order 3.

Figure 14 : 3 . 2 3 . 3 7 4 4 . 5 13 6

 14323374513 Figure 14: Turbulence test: (left) the numerical solution, with l = m = 10 and p = 4, (right) the difference uu h for a turbulence test in the tokamak model, with l = m = 2 and p = 4.

Table 1 :

 1 Nishimura test: Number of degree of freedom and L 2 norm of the error for each spline degree.

		Degrees of freedom L 2 error norm
	2	17408	7.20 10 -3
	3	18060	1.07 10 -3
	4	18720	1.71 10 -4
	5	19388	2.84 10 -5
	6	20064	4.81 10 -6
	7	20748	1.67 10 -6

Conclusion

This work presented a new approach for the resolution of elliptic type operator, with an application to the quasi-neutrality equation. The IGA-FEM approach seems to be very efficient to deal with this kind of problem, but generalization to a large class of operators can also be dealt with following the FEM approach.

A Appendix: Decoupling approach for J(r, θ)

The starting equation writes

with F = (n in e )/n 0 and φ (r) = φ(r, θ, ϕ)J(r, θ)dθdϕ J(r, θ)dθdϕ .

The appendix is devoted to the extension of the decoupling approach to (r, θ) depending jacobians J. As in the case where J(r, θ) = r, we integrate (A.31) with respect to θ, ϕ to get a one-dimensional equation on φ

Let us remark that this operation is transparent for the operator • . Moreover, we are looking for an equation satisfied by Φ = φ -φ from (A.31):

Let us come back to (A.32) in order to derive an equation satisfied by φ . From (A.32), we have

By introducing the notation h(r) = φφ , we can derive an equation for φ (r)

To conclude, we present the decoupling algorithm which enables to solve (A.31)

• (c) compute (and store) tmp = (Φ -Φ )


Julien Clement
email: jclement@liafa.jussieu.fr

Carole Delporte-Gallet

Hugues Fauconnier

Mihaela Sighireanu
email: sighirea@liafa.jussieu.fr

Mode d'emploi pour la v érification des protocoles de population †

Keywords: protocole de population, abstraction par comptage, model-checking, Spin, Prism

Cet article est consacré à la vérification basée sur le modèle (model-checking) du modèle des protocoles de population introduit par Angluin et al. [AAD + 06]. Ces deux dernières années, la vérification des protocoles de population par model-checking a fait l'objet de nombreuses études et de nouveaux outils ont été proposés ([PLD08, BFPvdP08, LPSZ09, SLDP09, CMS10]). Nous montrons dans cet article que, dans une certaine mesure, les outils classiques de model-checking tels que Spin et Prism peuvent être utilisés pour effectuer les vérifications. Pour cela, nous appliquons la technique d'abstraction par comptage [GS92] pour obtenir des modèles abstraits de protocoles de population qui peuvent être vérifiés efficacement par les outils de model-checking existants. Le problème essentiel pour la complexité de la vérification des protocoles de population concerne la condition d'équité forte. Cette dernière ne peut être utilisée directement avec Spin même pour des exemples de relativement petite taille. Nous montrons qu'on peut cependant remplacer dans de nombreux cas cette équité par une équité faible efficacement vérifiable en Spin. Plus notable encore, nous montrons que la vérification avec la condition d'équité forte est équivalente à un problème de vérification d'un modèle probabiliste. Ainsi, le model-checker probabiliste Prism s'avère être un outil de vérification adapté aux protocoles de population.

Introduction

Population protocols [AAD + 06] are a formal model for sensor networks. Since such networks should become a part of the every day life, it is important to ensure their correctness. This task is now possible in a fully automatic way due to the important progress of verification techniques. Proposed in the early eighties, the model-checking technique [START_REF] Emerson | Characterizing correctness properties of parallel programs using fixpoints[END_REF][START_REF] Queille | Iterative methods for the analysis of petri nets[END_REF] allows to check complex correctness properties for models with finite number of states. The system to be verified is formalized in some high level formal language, e.g., CSP or Petri Nets, to obtain a model. The requirements of the system are specified using some logical language, e.g. LTL [START_REF] Pnueli | The temporal logic of programs[END_REF], to obtain a property. Model-checking techniques explore exhaustively the model in order to check that each state of this model satisfies the given property (a reference for model-checking techniques is [START_REF] Clarke | Model Checking[END_REF]). High performance model-checkers, e.g., Spin [START_REF] Holzmann | Spin model checker : the primer and reference manual[END_REF], are now available. They are able to deal with systems of more than 10 20 states. Moreover, the model-checking techniques have been extended to deal with more complex finite models like the probabilistic ones, e.g., in the Prism [START_REF] Kwiatkowska | Probabilistic symbolic model checking with PRISM : A hybrid approach[END_REF] model-checker.

Unfortunately, some models of realistic systems have a number of states which exceeds the capacity of finite model-checkers. Thus, the research have focus on designing abstraction techniques to reduce the model-checking problem of huge systems to the model-checking of smaller systems. One of the first abstraction techniques proposed is the counting abstraction [START_REF] German | Reasoning about systems with many processes[END_REF]. It reduces systems with many identical processes running in parallel to a system which keeps track only of the number of processes which are in some particular state. This abstraction has shown its efficiency for the model-checking of safety properties in cache coherence protocols.

In this work, we apply the counting abstraction technique and finite model-checking tools to verify population protocols. Due to lack of space, we refer the reader to the long version [START_REF] Clement | Guidelines for verification of population protocols[END_REF] for formal definitions, proofs and experimental results.

The model of population protocols (PP) [AAD + 06] involves individual agents with a very simple behavior, which can be seen as a finite state machine. An important property is the uniformity of the protocol, i.e., the fact that the protocol description is independent of the number of agents (called also the population size) or their identity. When two agents come into range of each other ("meet"), they can exchange information. The agents are anonymous and move in an asynchronous way. At the beginning of the protocol, each agent receives a piece of input. The goal of the protocol is to stabilize each agent into a state in which it outputs the value to be computed by the protocol. Thus, the specification of a population protocol includes the "correct stabilization" (CS) property which requires that each agent stabilizes its computation to an output value which corresponds to a function on the input assigned to agents.

Our contribution is to highlight the use of the existing model-checking techniques to check the CS property for protocols with finite fixed sizes. For this, we use the counting abstraction to reduce the PP model to a vector addition system model, which is a model computationally equivalent to classical place/transition Petri Nets. We show that this abstraction allows the verification of the CS property in a more efficient manner. We also give sufficient conditions to check the CS property under a fairness constraint weaker than the one required in [AAD + 06]. These theoretical results allow us to verify, using the Spin tool, a benchmark of PP models for population size greater than 10 3 . In a second approach, we highlight that the fairness of the PP model can be exactly captured for the CS property in a probabilistic model where probabilities strictly greater than some ε > 0 are assigned to transitions. Then, we show that Prism can deal very efficiently with the model-checking of this property on the abstract model of PP.

Population Protocol Model and its Abstraction

In this section, we briefly introduce the population protocol model. More details are available in [AAD + 06]. For a population of k agents, Π = {π 1 , . . . , π k } denotes the set of agents.

Definition 1 A population protocol (PP) is specified as a six-tuple P = (Q, X,Y, ι, ω, δ) which contains a finite set Q of possible agent states, an input assignment ι : X → Q mapping the finite set X of inputs to the agents' state, an output assignment ω : Q → Y mapping the agents' states to a finite set of outputs Y , and a transition relation δ ⊆ Q 4 .

In this definition, the initial state of an agent is fixed (using ι) by an input value in X received by the agent. In each state q, an agent outputs a value in Y given by the mapping ω. The interaction between agents follows the rules described by the relation δ : if two agents in states q 1 and q 2 meet each other, they change into states q 3 resp. q 4 , where (q 1 , q 2 , q 3 , q 4 ) ∈ δ. We also use notation q 1 q 2 -→ q 3 q 4 for the elements of δ. The semantics of a PP model is given by a labelled transition system as follows. A configuration s of the protocol is a mapping Π → Q specifying the state of each agent. A transition between configurations

s ℓ -→ s ′ takes place if ∃π i , π j ∈ Π s.t. t = s(π i ) s(π j ) - → s ′ (π i ) s ′ (π j ) ∈ δ, ℓ = ({i, j},t) and ∀π ∈ Π\{π i , π j }. s(π) = s ′ (π).
From an initial configuration s 0 , a computation C from s 0 is s 0 , ℓ 0 , s 1 , ℓ 1 , . . . with s i ℓ i -→ s i+1 . One of the characteristics of the PP model is that the order in which pairs of agents interact is unpredictable. In order to model this aspect, one may imagine the presence of a scheduler which is scheduling the interactions. This scheduler may force two agents to never interact. In the presence of such scheduler, the PP has no chance to compute its goal. Thus, Angluin et al. [AAD + 06] require that the scheduler allows only strong globally fair computations.

Definition 2 A computation C is strong globally fair (GF) iff for every s and s ′ such that s -→ s ′ , if s = s i for infinitely many i in C, then s j = s ′ for infinitely many j in C i.e., s -→ s ′ ∧ ✷✸s ⇒ ✷✸s ′ .

Under the GF requirement for the scheduler, a protocol of size k stably computes a function f : X k → Y if, for every input x ∈ X k , every GF computation starting in ι( x) stably outputs f ( x), i.e., ι( x) ∧ GF ⇒ ✸✷ ω(s) = f ( x) . The important result of Angluin et al. [START_REF] Angluin | The computational power of population protocols[END_REF] is that a predicate (i.e., function with co-domain Y = {0, 1}) is stably computable by the PP model iff it can be defined as a first-order formula in Presburger logic.

To scale up the verification task, we use as abstract model for the PP model the vector addition system model [START_REF] Karp | Parallel program schemata[END_REF].

Definition 3 A vector addition system (VAS) is a pair A = (c, D) where c = c 1 , . . . , c n is a finite vector of integer variables called counters, and D ⊆ N n × Z n is a finite set of guarded translations (µ, τ) s.t. µ + τ ≥ 0.

To obtain a VAS model from a PP model, we apply the counting abstraction [START_REF] German | Reasoning about systems with many processes[END_REF]. Intuitively, a configuration s is abstracted by a vector c indexed by Q such that c[q] is the number of agents in state q in s. This abstraction is possible due to the uniformity of the protocol : the behavior of an agent depends only on its state and not on its index. We prove (see [START_REF] Clement | Guidelines for verification of population protocols[END_REF]) that the counting abstraction preserves the CS property and thus the verification of the CS property on the PP model can be transferred into a verification of the same property in the abstract PP model.

Verification under Weak Fairness

The verification of the CS property on the abstract model of PP for finite population sizes can be done using existing finite state model-checkers. We consider the Spin tool [START_REF] Holzmann | Spin model checker : the primer and reference manual[END_REF] which is one of the most efficient model-checkers for finite-state systems. To check a property on a model, Spin translates the property into a Büchi automaton which size is exponential in the size of the formula. This fact prevents us to encode the GF constraint in the formula in order to select only GF computations. Fortunately, Spin provides an algorithm which verifies a property on the set of "weakly fair" computations of a model. The selection of these computations is done on-the-fly and very efficiently.

Definition 4 A computation C is weakly fair (WF) iff for every ℓ, if it exists j and s ′ i with s i ℓ -→ s ′ i ∈ δ for any i ≥ j, then s k ℓ -→ s k+1 for infinitely many k > j in C, i.e., ✸✷ℓ is enabled ⇒ ✷✸ℓ is fired.

Thus, we study the conditions on which the WF constraint can be used for the verification of the abstract PP models. The problem with the use of the WF constrain is that we can obtain false negatives since the set of WF computations includes the set of GF computations. We provide two sufficient conditions on the transition system of the abstract PP model to eliminate the false negatives. These conditions constrain the strongly connected components (SCC) of the induced transition system.

Theorem 1 (Correction of verification under WF) Let P be a PP model whose abstract model has an induced transition system A which satisfies one of the following constraints :

1. A has only sink strongly connected components. 2. All the SCC of A have the same output value. Then P stably computes f under the GF constraint if A stably computes f under the WF constraint

Probabilistic Verification

An alternative way to introduce fairness constraints in the PP model is to consider probabilistic schedulers. This alternative is considered by Angluin et al. [START_REF] Angluin | A simple population protocol for fast robust approximate majority[END_REF] which define the notion of probabilistic computation for a PP model. The simplest and more natural probabilistic scheduler proposed in [AAD + 06] schedules ordered pairs of agents.

Definition 5 A random pairing (RP) scheduler chooses, in each state of a computation, the ordered pair of agents which interact in a random independent and uniform manner from all ordered pairs of agents.

The RP scheduler does not combine well with the counter abstraction because this abstraction collapses states whose sets of agents is the same. Therefore, we have to consider a scheduler whose decision is based on the transitions of the PP which are also transitions of the abstract VAS model.

Definition 6

The random ruling (RR) scheduler chooses, in each state of the computation, one enabled transition t ∈ δ of the PP in a random uniform and independent manner. We show that, as in the non-probabilistic case, to verify that a predicate f is computed almost surely when probabilities are strictly positive consists in searching the closed strongly connected components in the transition system. Thus, we obtain the equivalence of the RR and RP schedulers w.r.t. the stably computation property as a consequence of the ergodic theorem on finite discrete Markov chains.

Proposition 1 A PP P stably computes a predicate f ⇔ P computes f almost surely with an RP scheduler.

This result is the theoretical basis for the use of probabilistic model-checkers to verify that P stably computes a predicate f even for a non randomized model of population protocols. Indeed, we have to verify that the abstract model of P computes f almost surely with some scheduler assigning strictly positive probabilities to transitions.

We use the Prism tool [START_REF] Kwiatkowska | Probabilistic symbolic model checking with PRISM : A hybrid approach[END_REF] to perform probabilistic model-checking. For each protocol, we encode its abstract VAS model in a discrete Markov chain model of Prism. The probabilities on transitions are not fixed, thus Prism associates to each transition enabled at some state a uniform distribution of probabilities. This model is checked against the property ✸✷ f with probability at least one, i.e., almost surely.

Conclusion

We demonstrate that standard, finite state model-checking tools may help to verify that PP stably computes some predicate. In addition, these tools may be used to verify other properties of population protocols, i.e., properties on computations that can be expressed by LTL formulas. Such properties may be (1) all agents are infinitely often in some state, or (2) all computations are private in some PP with privacy [START_REF] Delporte-Gallet | Secretive birds : Privacy in population protocols[END_REF].

† This work has been accepted for publication at ICDCS 2011[START_REF] Clement | Guidelines for verification of population protocols[END_REF].


Olivier De

La Barrière
email: barriere@satie.ens-cachan.fr

Carlo Appino

Fausto Fiorillo

Ragusa Carlo

Hamid Ben Ahmed

Mohamed Gabsi

Frédéric Mazaleyrat

Martino Lobue

O De La Barrière

C Ragusa

Loss separation in soft magnetic composites

Keywords: numbers: 75, 50, Bb

à 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.

I. INTRODUCTION

Three-dimensional flux paths and high frequencies are ubiquitous in modern electrical machines [1] [2] and laminated cores may often prove inadequate. In fact, while inhibited from out-of-plane magnetization, soft magnetic laminations show a rapid increase of the power loss and fall into skin effect phenomena in the presence of harmonics and/or high working frequencies [3] [4]. Soft magnetic composites (SMCs) are therefore often proposed as convenient substitutes, their use appearing especially suited for devices like permanent magnet, claw-pole, and axial flux machines.

While many literature efforts in these materials have been devoted to preparation methods and microstructural studies, with the aim of improving their manufacturing process, their mechanical properties, and their magnetic behavior [5] [6], the analysis of the magnetic properties has mostly been performed on qualitative grounds [7] [8], aiming at rapid exploitation in SMC machine core modeling [9]. For example, the important problem of loss calculation in such materials has been approached by use of polynomial interpolation methods [10] [11], which are not accurate enough to predict the loss under highly distorted induction waveforms. De Wulf, et al. [12] have resorted to the concept of loss separation and attempted to apply Bertotti's statistical theory [13] to the frequency dependence of the energy loss in different types of SMCs. Their study was, however, limited to the maximum frequency of 100 Hz and their calculations disregard the actual heterogeneous structure of the material.

In this paper we present significant new results on the energy losses and their frequency dependence (DC -10 kHz) in commercial soft Fe bonded composites. The experiments show that the sample cross-sectional area negligibly affects the loss at all frequencies. The classical loss component is therefore calculated for eddy current paths confined within the grains and the loss decomposition, leading to the determination of the hysteresis and excess losses, is accordingly performed . The excess loss is observed to depend on frequency in a quasi-linear fashion, in contrast with the square-root law usually observed in magnetic laminations. This provides evidence for correspondingly weakened dynamic homogenization of the magnetization process, ensuing from combination of strong internal demagnetizing effects and lack of long-range eddy currents.

I. EXPERIMENTAL RESULTS AND DISCUSSION

Two commercial SMCs (Somaloy Prototyping and Somaloy 110i1P) have been investigated. The first one (Material A) has the following physical properties: density δ = 7450 kg/m 3 , electrical resistivity ρ = 2.80⋅10 -4 Ωm, average particle size <s> = 114 µm. The second material (B) has quite smaller particle size (<s> = 29.5 µm ) and lower filling factor (density δ = 7260 kg/m 3 , electrical resistivity ρ = 76⋅10 -4 Ωm). Ring samples of mean diameter 50 mm and 25 mm were prepared for magnetic testing. The larger rings had cross-sectional area S 1 = 5 × 5 mm 2 , the smaller ones S 2 = 2.5 × 2.5 mm 2 . Hysteresis loop and loss measurements under sinusoidal induction were carried out, starting from quasi-static conditions (5 Hz for A and 20 Hz for B), up to f = 10 kHz for peak polarization values J p = 0.50, 1.0, 1.25 T. A broadband hysteresisgraph-wattmeter with digital control of the induction waveform endowed with a 5kVA DC-20 kHz CROWN power amplifier was employed. Sample heating was prevented by keeping the sample immersed in a water bath.

A remarkable outcome of these experiments is that the energy loss W(f) does not depend, at all frequencies, on the cross-sectional area of the sample (Fig. 1, SMC type A). The same occurs in the SMC type B. This finding implies that grain-to-grain conductivity phenomena, either due to random metallic contacts through the binder [START_REF] Nord | Recent Developments of Electrical Drives[END_REF] or capacitive effects at high frequencies can be disregarded from the viewpoint of losses. With the eddy current paths confined within the particles, a good simplification is introduced in the loss analysis, because one can concentrate on the behavior of the individual particles and their statistics. To start with, we calculate the classical loss component W cl , which in this case is the contribution that would originate if all particles were magnetized at the same uniform rate imposed to the sample as a whole. Scanning electron microscope images permit one to appreciate size and morphology of the particles and their distribution. The cross-sectional area S part of each particle is determined from a number of micrographs using an image processing software. Their shape is irregular (Fig. 2), but, for the sake of calculation, they are approximated as rectangles. The width and length of these equivalent rectangles are identified by performing extended micrographic analysis and by taking the mean values of particle area and aspect ratio. The Poisson equation for the magnetic field is then solved by a finite element procedure, [12] in order to compute the eddy current loss in the rectangles. The results for materials A and B are shown in Figs. 3 and4, respectively, where the decomposition of the measured loss W(f) into the hysteresis Whyst and excess Wexc(f) loss components is put in evidence. Whyst is an experimental quantity, obtained by extrapolating the measured W(f) to f=0.

The strong reduction of Wcl(f) ensuing from the reduction of the particle size is apparent. By associating it to a concurring reduction of Wexc(f), one gets that, in spite of higher Whyst, the total loss W(f) in the fine-grained SMC type B becomes lower than in the coarse-grained SMC type A beyond about 2 kHz.

As known, the excess loss arises because the dynamic magnetization process, like the quasi-static one, is inhomogeneous in nature. However, since the eddy-current counterfields become stronger with increasing frequency and they have to be compensated by the applied field in order to keep a given magnetization rate ܬ p ሶ , progressive homogenization of this process is concurrently expected. This concept is made quantitative in terms of active magnetic objects (MOs) and their frequency dependence by the statistical theory of losses [13], which provides the following general expression for the excess loss

2 2 16 ( , ) 2 1 1 o exc p o o p p o o GSV W J f n V J fJ n V σ   = + -       (1)
where G = 0.1356, S is the cross-sectional area of the sample, V o is a parameter related to the distribution of the local pinning fields for the MOs, and n 0 is the number of MOs that are simultaneously active within the sample cross-section under quasi-static magnetization. It is generally found in magnetic laminations that n 0 is at most a few units [13]. A basic tenet of the theory is that the number of active MOs varies with frequency according to the linear law n = n 0 + H exc /V 0 , where H exc = W exc /4J p is the excess field to be dynamically applied in order to counteract the mesoscopic (i.e. circulating close to the moving domain walls) eddy current fields. By taking the experimental W exc (f), as obtained by subtracting W hyst + W cl (f) from the measured W(f) values (see Figs. 3 and4), the n versus H exc behavior is obtained, as shown for the A and B type SMCs in Fig. 5 for J p = 1T. The experimental n 0 and V o values found in this way are introduced in Eq. ( 1) and the so calculated W exc (f), added to W hyst and the previously obtained W cl (f), results into the fitting loss curve for W(f) shown in Figs. 3 and4. A striking result of this analysis is that n0 is orders of magnitude larger than in soft magnetic laminations. n 0 is a quantity that refers to the whole crosssectional area of the sample [13] and it appears here that the magnetization process in SMCs tends to simultaneously invest a large number of particles at the same time already at low frequencies. It can therefore be assumed that there should be at least one active magnetic object per particle at zero excess field (i.e. under quasi-static excitation). If we take the type A SMC sample with S =25 mm2 and average grain-size area <s>2 = 0.0130 mm 2, one gets the mean number of particles <npart> in the sample cross section <npart> = 1924, close to the number n0 of MOs (Fig. 5). For this material, the observed dependence of n on the excess field (i.e. frequency), shown for Jp = 1.0 T in Fig. 5, suggests that each particle should contain a few MOs, which will be dynamically activated on increasing the frequency. The type B SMC shows a remarkably higher no value, as expected from a correspondingly lower grain size. In both cases we observe that the condition no >> Hexc/Vo is approximately satisfied in a good part of the investigated frequency range. A first order expansion of Eq. ( 3) correspondingly provides Wexc∼∝ f. This contrasts with usual Hexc/Vo >> no condition already found at low frequencies in magnetic laminations [13], which leads to the typical Wexc∼∝ f1/2 law.

III. CONCLUSIONS

The measurement of magnetic losses from DC to 10 kHz in commercial soft magnetic composite materials and the analysis of the results -centered on the concept of loss decompositionhighlights peculiar features of these materials. It is verified that the eddy currents are by and large hindered from grain-to-grain circulation at all frequencies, with ensuing benefits from the viewpoint of energy dissipation. This fact permits one to simplify the problem of loss calculation, whereby the classical loss component can be analytically derived from knowledge of the distribution of the particle size. A general formulation for the excess loss component, as provided by the statistical theory of losses, can then be applied and shown to describe the experimental results and their dependence on frequency. It is induced that the magnetization process in SMCs is characterized by a good degree of homogenization already at very low magnetizing frequencies. The lack of macroscopic eddy current patterns and the related counterfields eventually results, however, into a faster than usual increase of the excess loss component with frequency and into its correspondingly good contribution to the total loss up to the maximum investigated frequency of 10 kHz.

Figure captions

Fig. 1: Energy loss versus frequency measured at peak polarization values J p = 0.5 T, 1.0 T, 1.25 T in the type A soft magnetic composite (δ proto =7450 kg/m 3 ). The results are obtained on two ring samples of crosssectional areas S 1 = 5 × 5 mm 2 and S 2 = 2.5 × 2.5 mm 2 , respectively.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Micrograph for the type B SMC. The classical eddy current loss component Wcl(f) is calculated assimilating this structure with an array of rectangular grains.

Fig. 4 :

 4 Fig.4: Same as Fig.3for the SMC type B.

Fig. 5 :H

 5 Fig. 5: Number of active magnetic objects n versus the excess field H exc = W exc /4J p at J p = 1.0 T in the two investigated SMCs. The straight lines provides best fitting with the predicted relationship n = n 0 + H exc /V 0 .


A Sinatra

E Witkowska

J.-C Dornstetter

Li Yun

Y Castin

Limit of Spin Squeezing in Finite Temperature Bose-Einstein Condensates

Keywords: numbers: 03.75.Gg, 42.50.Dv, 03.75.Kk, 03.75.Pp, 03.75.Mn

We show that, at finite temperature, the maximum spin squeezing achievable using interactions in Bose-Einstein condensates has a finite limit when the atom number N → ∞ at fixed density and interaction strength. We calculate the limit of the squeezing parameter for a spatially homogeneous system and show that it is bounded from above by the initial non-condensed fraction.

Atomic clocks based on cold alkali atoms in two hyperfine states a and b are widely used as frequency standards. When atoms in uncorrelated quantum states are used, the clock precision is limited by the so-called projection noise, resulting from the quantum nature of the collective spin S, i.e. the sum of the effective spin 1/2 of each atom. This limit is actually already reached in most precise clocks [1]. Spin squeezing [2] amounts to creating quantum correlations among the atoms so as to increase the precision of the atomic clock beyond this standard quantum limit. The relative improvement on the variance of the measured frequency ∆ω 2 ab defines the spin squeezing parameter ξ 2 [3]. Spin squeezing in atomic ensembles was first obtained by quantum non-demolition measurements [4,5]. Recently a significant amount of spin squeezing (e.g. 6 or 8 dB) has been achieved using atoms in a resonant optical cavity [6] or exploiting atomic interactions in bimodal Bose-Einstein condensates [7][8][9]. The ultimate limits of the different paths to spin squeezing are still an open question. We determine here the influence of the non-condensed fraction for spin squeezing schemes using Bose-Einstein condensates [8][9][10].

A central issue is the scaling of the squeezing for large atom numbers. Most studies are based on a two-mode description [2]. In this case the squeezing parameter optimized over time ξ 2 best tends to zero (infinite metrology gain) for N → ∞ as ξ 2 best ∼ N -2/3 . The first analysis of squeezing at finite temperature [11] used a large N expansion in a Bogoliubov-like approach and could not predict any deviation of spin squeezing from the twomode model. Here, using fully non-perturbative semiclassical field simulations and a powerful formulation of Bogoliubov theory in terms of the time dependent condensate phase operator [12], we find on the contrary a dramatic effect of the multimode nature of the field: For a spatially homogeneous system in the thermodynamic limit, the two-mode scaling ξ 2 best ∼ N -2/3 turns out to be completely irrelevant, and the spin squeezing has a finite optimal value that we determine analytically.

The physical problem that we face is the dynamical evolution of a finite temperature Bose condensed gas after a pulse π/2 that puts each atom in a coherent superposition of two internal states a and b. This produces a non-equilibrium state that has a non-trivial evolution due to the atomic interactions inside each internal state. For simplicity, it is assumed that there is no cross-interaction between a and b atoms.

Semi-classical field simulations -In Fig. 1 we compare the two-mode theory with semi-classical field simulations at finite temperature in a trap. The gas is initially in state a at thermal equilibrium. In that state, we assume that thermal fluctuations dominate over quantum fluctuations and we use a classical field description [13][14][15][16] with an energy cut-off at k B T . The initial field ψ (0) a then randomly samples the thermal equilibrium classical field distribution for the canonical ensemble at temperature T . For the initially empty state b, inspired by the truncated Wigner approach [17,18] we represent the vacuum by a classical field ψ (0) b having in each mode independent Gaussian complex fluctuations of zero mean and variance 1/2. A sudden π/2 pulse mixes the initial fields ψ (0) a and ψ (0) b so that, at time t = 0 + , i.e. just after the pulse:

ψ a,b (0 + ) = 1 √ 2 [ψ (0) a,b ∓ ψ (0) b,a ] . (1)
At later times, each field evolves independently according to the non-linear Schrödinger equation

i ∂ t ψ a,b = - 2 ∆ 2m + 1 2 mω 2 r 2 + g|ψ a,b (r, t)| 2 ψ a,b .
(2) This corresponds to a harmonically trapped gas with same oscillation frequency ω and same coupling constant g = 4π 2 a/m for the two internal states, where a is the s-wave scattering length. As shown in Fig. 1, the squeezing is created dynamically by the interactions. However, even for a moderate non-condensed fraction N nc /N = 0.09, the best ξ 2 in the multimode theory is larger by more than one order of magnitude than in the two-mode theory.

In order to isolate the effect of the non-condensed fraction from other dynamical effects taking place in the trapped system, as for example the spatial dynamics of the condensate wave function [19,20], and to develop an analytical theory, we consider from now-on the homogeneous case. We first use the semi-classical field model that has the advantage that it can be simulated exactly, and we generalize the results to the case of a quantum field in the end. The real space is discretized on a lattice with unit cell of volume dV , within a volume V with periodic boundary conditions [12]. The Hamiltonian after the pulse for component a (and similarly for b) reads:

H = k 2 k 2 2m a * k a k + g 2 dV r |ψ a (r)| 4 . (3)
The fields have Poisson brackets i {ψ µ (r), ψ * ν (r ′ )} = δ rr ′ δ µν /dV with µ, ν = a or b, and a k (b k ) is the amplitude of ψ a(b) the over the plane wave of momentum k. In terms of the fields, the collective spin components are

S x + iS y = d 3 r ψ * a (r )ψ b (r ), (4)
S z = 1 2 d 3 r [ψ * a (r )ψ a (r ) -ψ * b (r)ψ b (r)] . (5)
The spin squeezing parameter ξ 2 is equivalent to the minimal variance of the spin orthogonally to its mean direction, divided by the mean spin length squared and suitably normalized. Here the mean spin is along x so that

ξ 2 (t) = ∆S 2 ⊥,min (t) S x (t) 2 × S x (0 + ) 2 ∆S 2 ⊥,min (0 + ) (6) ∆S 2 ⊥,min = 1 2 S 2 y + S 2 z -| (S y + iS z ) 2 | . ( 7
)
As a first step, we performed semi-classical field simulations for different temperatures and increasing system sizes [21]. The result (not shown) is that ξ 2 best converges to a finite value at the thermodynamic limit: N →∞, V →∞, ρ, g, T =constant, where ρ = N/V is the total density. Five independent physical parameters are in the model, /m, g or a, k B T, N and V . From dimensional analysis, ξ 2 best is a function of the three independent dimensionless quantities that one can form, N, ρa 3 , and k B T /ρg. The existence of a thermodynamic limit then implies

ξ 2 best = f ρa 3 , k B T ρg . (8)
As a second step, we performed simulations increasing the density in the weakly interacting limit [22], ρ → ∞, g → 0 with T, ρg=constant. We find that, for a given k B T /ρg, ξ 2 best then scales as 1/ρ ∝ ρa 3 . This implies:

ξ 2 best / ρa 3 = F (k B T /ρg) . (9)
In Fig. 2 we show the universal behavior (9). The circles and the squares correspond to two different values of ρa 3 in simulations.

Semi-classical field analytics -We now develop an analytical theory to explain these results. We split the fields after the pulse as ψ a = a0 √ V + ψ a⊥ and similarly for ψ b . We introduce the modulus and phase conjugate variables for the condensate modes

a 0 = e iθa N a0 , b 0 = e iθ b N b0 , (10)
and we introduce number conserving non-condensed fields Λ a and Λ b [22] that we expand over Bogoliubov modes with amplitudes c ak and c bk respectively [23]:

Λ a = e -iθa ψ a⊥ = k =0 U k c ak + V k c * a-k e ik•r √ V (11)
U k + V k = E k E k + ρg 1/4 ; E k = 2 k 2 2m . (12)
The spin raising component S + = S x + iS y is given by

S + = e -i(θa-θ b ) N a0 N b0 + d 3 r Λ * a Λ b . ( 13
)
Our strategy is to perform a double expansion of S 2 + . We will need terms up to ∼ N in the thermodynamic limit and up to order one in the non-condensed fraction N nc /N . In this framework, we can approximate S x (t) in the denominator of (7) by its value at t = 0 + so that

ξ 2 ≃ 4 N ∆S 2 ⊥,min . (14)
We sketch the main steps. In the Bogoliubov limit, the condensate phases at t > 0 obey [12] θ a - FIG. 2: Best squeezing ξ 2 best divided by ρa 3 as a function of kBT /ρg. Symbols: semi-classical field simulations with ρa 3 = 1.32 × 10 -2 (filled squares) and 1.94 × 10 -3 (disks). The thermodynamic limit is reached already for N = 3 × 10 4 except for the lowest value of kBT /ρg. Lower solid line: analytical semi-classical field result (21). Upper solid line: quantum result (25). Inset: quantum ξ 2 best (solid line) and non-condensed fraction Nnc /N (dashed line), both divided by ρa 3 , as functions of kBT /ρg. S is the multimode part of the relative phase derivative that is absent in the two-mode theory. In thermodynamic limit θ a -θ b ∼ 1/ √ N and it is sufficient to expand the exponential in (13) to second order. In the modulus of S + we expand:

θ b = (θ a -θ b )(0 + ) - ρg V t [(N a -N b ) + S] (15) (θ a -θ b )(0 + ) = -2 Im b (0) 0 N (0) a0 + O(N -1 ) (16) S = k =0 (U k + V k ) 2 |c ak | 2 -|c bk | 2 . ( 17
)
N a0 N b0 ≃ N tot 2 - 1 2 d 3 r |Λ a | 2 + |Λ b | 2 (18)
with

N tot = N + k |b (0) k | 2
is the total atom number in the semi-classical field picture .

Best squeezing -For the calculation of the best squeezing, one looks at the asymptotic behavior of ( 14) for t → ∞. One finds:

ξ 2 (t) = ξ 2 best + ρgt 2 1 + O N nc N + O 1 t 4 (19) with the best squeezing ξ 2 best = S 2 /N . ( 20
)
which remarkably only involves the multimode part (17) of the phase difference. An explicit calculation gives

ξ 2 best = 1 2ρ d 3 k (2π) 3 s 4 k n (0) k   s (0) k s 2 k 2 + s 2 k s (0) k 2   (21)
(solid line in Fig. 2). In (21), (12), and s (0) k is the equivalent quantity before the pulse obtained by replacing ρg with 2ρg in (12); n

s k = U k + V k given in
(0) k = k B T /ǫ (0) k
are the equilibrium occupation numbers of Bogoliubov modes before the pulse with ǫ

(0) k = [E k (E k + 2ρg)] 1/2 .
Squeezing time -From (19), the best squeezing is reached in an infinite time, which is a limitation of the / / h FIG. 3: "Close to best" squeezing time tη for η = 0.2. Filled squares and disks: simulations as in Fig. 2. Line: analytical prediction (22). Squares and circles: thermalization times in the simulations extracted from the decay of the contrast Sx . analytical approach. However, the numerical squeezing curve as a function of time is indeed quite flat around its minimum, so that it suffices in practice to determine the "close to best" squeezing time t η defined as

ξ 2 (t η ) = (1 + η)ξ 2
best , where η > 0. Then, according to (19), t η is finite and given by

ρg t η = 1 ηξ 2 best . (22)
The "close to best" squeezing time t η (22) for η = 0.2 is shown in Fig. 3 and compared to simulations.

A last important issue is that of thermalization, neglected in Bogoliubov theory and in our analytical treatment, but fully included in the semi-classical field simulations. Indeed it is possible to reach ξ 2 = (1 + η)ξ 2 best with ξ 2 best given by ( 21) only if t η given by ( 22) is shorter than the thermalization time

t η < t therm . (23)
In Fig. 4 we show the squeezing parameter ξ 2 and contrast S x the thermalization process that brings the system back to equilibrium after the pulse. For the squeezing, we compare the simulation with (i) the full Bogoliubov theory (without the analytic expansions) that we implement numerically for a finite size system and (ii) a Bogoliubov ergodic model [12] where the amplitudes c ak , c bk in (15) sample microcanonical distributions with number of particles and the energy {N a , E a } ({N b , E b }) fixed to a random value set by the pulse. Note that the simulation agrees with the Bogoliubov model at short times (included the "close to best" squeezing time) and then converges towards the ergodic model. We extract a thermalization time t therm from the contrast. As thermalization occurs the excited modes dephase and The longer time scale for the decay of b * 0 a 0 is set by phase spreading due to partition noise [24] plus thermal corrections.

S x = Re k b * k a k ∼ t→∞ Re(b * 0 a 0 ) . (24)
Quantum field -All our analytic calculations for the semi-classical field can be generalized to the quantum field. In particular, ( 20) and ( 22) are unchanged and (26) In practice ρa 3 < 10 -6 in present squeezing experiments so that (26) predicts ξ 2 T =0 best < ∼ 2.10 -5 . This value is very low, in particular below the limit given by particle losses [25]. Asymptotically for k B T ≫ ρg, ξ 2 best identifies with the initial non-condensed fraction. An interesting result is that at any temperature the initial non-condensed fraction is larger than ξ 2 best , see these two quantities in the inset of Fig. 2. Already for k B T /ρg = 2, ξ 2 best and N nc /N are within a factor three. A similar conclusion seems to hold in a trap, see Fig. 1.

ξ 2 best = d 3 k (2π) 3 s 4 k 2ρ n (0) k + 1 2 (s (0) k ) 2 s 4 k + s 4 k (s (0) k ) 2 -1 . (25
In conclusion we have shown that a realistic description of the limits of spin squeezing in interacting Bose-Einstein condensates has to be multimode: The best achievable spin squeezing ξ 2 best admits a finite limit for N → ∞ at fixed density and interaction strength, contrarily to the vanishing prediction of the two-mode model. We find that ξ 2 best is the product of ρa 3 and of a universal function of k B T /ρg that we calculated analytically, and is bounded from above by the initial non-condensed fraction. Our analytical treatment is restricted to evolution times smaller than the thermalization time, but this is enough to access ξ 2 best as we showed by semi-classical field simulations (that include thermalization) over a wide range of parameters.

FIG. 1 :

 1 FIG. 1: (Color online) Spin squeezing parameter ξ 2 as a function of time after the pulse mixing the states a and b, for N = 10 5 Rb atoms (s-wave scattering length a = 5.3 nm), in a harmonic trap with oscillation frequency ω/2π = 50Hz. The Thomas-Fermi chemical potential is µ = 15.36 ω. Finite temperature semi-classical field simulations with initial noncondensed fractions: (a) Nnc /N = 0.34 (black solid line), (b) Nnc /N = 0.20 (red line), (c) Nnc /N = 0.09 (blue line), corresponding to kBT /µ = 2.08, 1.53, 1.17 respectively. Dashed line: Two-mode theory for comparison.

FIG. 4 :

 4 FIG. 4: (Color online) Top: Spin squeezing ξ 2 as a function of time. Black solid line: simulation. Blue dashed line: Bogoliubov theory. Red dash-dotted line: ergodic model. Bottom: contrast Sx (solid line) and condensate contrast Re b * 0 a0 (dashed line) as functions of time. N = 3 × 10 4 , kBT /ρg = 3.16, ρa 3 = 1.32 × 10 -2 .

1 )

 1 -π ≃ 0.02344 .

EW acknowledges support from Polish Government Research Funds: N N202 104136, 2009-2011.