Nicolas Crouseilles

Sever Pierre Glanc

Adrian Hirstoaga

Eric Madaule

Michel Mehrenberger

Sever A Hirstoaga

Jérôme Pétri

A new

Keywords: -key 52, 30, GzGyrokinetics --key 52, 65, FfPlasma simulation: Fokker-Planck and Vlasov equation

come

Introduction

Kinetic models are known to be a precise tool able to describe plasma arising in a tokamak or a beam. However, this accurate description requires an important effort in terms of numerics, since the solution of a kinetic model depends not only on the space variable, but also on the velocity variable, which gives a six dimensional phase space. In addition to the filamentation phenomena and presence of small scales, the numerical approximation of the distribution function solution to the kinetic model is a huge challenge despite the current advancement in computers power.

In recent decades, several different numerical methods have been developed to this aim. Among them, the Particle In Cell (PIC) exploits the Lagrangian formulation of the kinetic equation by solving the characteristics equations. In spite of its relative low computational cost, these kinds of method suffer from inaccuracy in low density regions and low convergence rate when the number of macro-particles increases (see [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF]). On the other hand, Eulerian methods have been introduced. They are based on methods like finite volumes or finite difference methods which have been mostly developed in fluids mechanics; the main particularity is that they require a grid of the full phase space which makes them accurate indepen-dently of the value of the distribution function, but they are very costly in terms of memory and computational power. Moreover, a restriction on the time step has to be satisfied for stability reasons. An interesting compromise is the semi-Lagrangian method which tries to take benefit from the Lagrangian approach by tracking the time evolution of the characteristic equations of the kinetic model, and from the Eulerian framework by using a grid of the phase space. In particular, the so-obtained method is free from standard CFL conditions on the time step.

Different semi-Lagrangian methods have been developed over this last 40 years [START_REF] Gagné | A splitting scheme for the numerical solution of a one-dimensional Vlasov equation[END_REF], [START_REF] Cheng | The integration of the Vlasov equation in configuration space[END_REF], [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation[END_REF], to first compute the solution of the 1D × 1D Vlasov-Poisson equation and then of the guiding-center model and more recently of the 5D gyrokinetic model. In addition to the high dimensionality, the spatial geometry of the latter case is toroidal, which makes it more complicated to deal with.

While developing a new semi-Lagrangian solver, the gap between a linear Landau run obtained with the 1D × 1D Vlasov-Poisson model and a 5D gyrokinetic simulation in toroidal geometry is quite huge. Then, intermediate test cases are welcome for checking the performance and precision of the codes.

Tests for which the analytical solution is known at all times do not include usually sufficient complexity to tackle the difficulties of the final problem. Moreover, although interesting, some test cases are not published and often different parameter settings are chosen which makes cross validations with different schemes difficult. Therefore, one purpose of this work is to provide precise parameter settings for two test cases which have intermediate complexity between the analytic test and the realistic physical configuration. This enables us to validate the new conservative semi-Lagrangian method.

We consider here as a building block, a 2D guidingcenter type equation on an annulus satisfied by f (t, r, θ):

∂ t f - ∂ θ Φ r ∂ r f + ∂ r Φ r ∂ θ f = 0, t ∈ [0, T ], (1)
with (r, θ)

∈ Ω = [r min , r max ]×[0, 2π],
where the potential Φ solves a Poisson type equation. This model constitutes the basis of the two tests. We first revisit a 2D test case previously done with a PIC approach [START_REF] Pétri | Non-linear evolution of the diocotron instability in a pulsar electrosphere: 2D PIC simulations[END_REF] and detail the boundary conditions. We then consider a 4D drift-kinetic slab simulation [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF], [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF] which can be viewed as a polar guiding-center equation of type (1) coupled with a Vlasov type equation in the parallel (z, v) direction.

In these two configurations, some first results of a new conservative semi-Lagrangian method are given and compared to the standard one [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. Indeed, it is well known that this standard method, namely backward semi-Lagrangian (BSL), does not preserve the total mass exactly in the 2D general configuration. Then, in the spirit of [START_REF] Crouseilles | Conservative semi-Lagrangian method for the Vlasov equation[END_REF], we propose a new fully two-dimensional conservative semi-Lagrangian (CSL) method. More precisely, a generalization to the 2D case of [START_REF] Crouseilles | Conservative semi-Lagrangian method for the Vlasov equation[END_REF] is presented in this work to deal with model of the form (1) without considering directional splitting. Indeed, even if time splitting techniques are very efficient in the Vlasov-Poisson case, for the drift-kinetic or guidingcenter models considered here, the characteristic curves are more complicated and need to be tackled in the full 2D case. Hence, to overcome the lack of conservation of the standard BSL method, it seems interesting to develop, validate and compare a 2D conservative semi-Lagrangian method. Let us remark that similar motivations are considered in [START_REF] Yang | Conservative and non-conservative methods based on Hermite weighted essentially-nonoscillatory reconstruction for Vlasov equations[END_REF].

The rest of this work is organized as follows: first, the new conservative semi-Lagrangian method is presented. Then, we present in detail two test cases enabling validation of our new method. In section 3, the diocotron instability test case [START_REF] Pétri | Non-linear evolution of the diocotron instability in a pulsar electrosphere: 2D PIC simulations[END_REF] is revisited by detailing the radial boundary conditions and the associated instability rates. Finally section 4 is devoted to a 4D drift-kinetic test case.

A new 2D conservative method

In this section, we present the new conservative semi-Lagrangian method (CSL) in the cartesian case. We also briefly recall the standard BSL method and present how the new method can be extended to polar geometries.

We then consider a 2D advection equation of the form

∂ t f + U (x, y) • ∇ x,y f = 0, (x, y) ∈ [0, L x ] × [0, L y ], (2)
where the divergence free vector field is

U (x, y) = (U x (x, y), U y (x, y)),
and we consider an initial condition f (t = 0, x, y) = f 0 (x, y).

We first define a spatial uniform mesh x i = i∆x, ∆x = L x /N x , y j = j∆y, ∆y = L y /N y , and a time discretization t n = n∆t with n ∈ N, ∆t > 0. Then, for all (x i , y j ), we define the characteristics equation with X(t) := (x, y)(t)

Ẋ(t) = U (X(t)), X(t n+1 ) = (x i , y j ), t ∈ [t n , t n+1 ]. (3)
Let us denote by X(s; t n+1 , (x i , y j )) the solution of this differential system, for s ∈ [t n , t n+1 ].

In this subsection, a new conservative semi-Lagrangian is proposed. As in the finite volumes context, the unknowns in this case are not pointwise values but cell average values. Let us introduce a control volume A i,j = [x i-1/2 , x i+1/2 ]×[y j-1/2 , y j+1/2 ], with x i±1/2 = x i ±∆x/2 and y j±1/2 = y j ± ∆y/2 we define the unknown as

f (t) i,j = 1 ∆x∆y Ai,j f (t, x, y)dxdy. (4)
We use the fact that the volume is preserved with time.

In particular, between times t n and t n+1 , we get

( f n+1 i,j =) 1 ∆x∆y Ai,j f n+1 (x, y)dxdy = 1 ∆x∆y X(t n ;t n+1 ,Ai,j )
f n (x, y)dxdy, [START_REF] Crouseilles | 2D advective and conservative semi-Lagrangian schemes upwind Hermite type reconstruction[END_REF] where X(t n ; t n+1 , A i,j ) = X(t n ; t n+1 , (x, y)), (x, y) ∈ A i,j , X(t n ; t n+1 , (x, y)) is the solution at time t n of the following characteristics equation

Ẋ(t) = U (X(t)), X(t n+1 ) = (x, y) ∈ A i,j , t ∈ [t n , t n+1 ].
Since f n i,j is known ∀i, j of the mesh, we can reconstruct a piecewise polynomial function P i,j ; an example can be the Lagrange reconstruction which satisfies, ∀i, j

1 ∆x∆y A i+ ,j+k P i,j (x, y)dxdy = f n i+ ,j+k , ∀ , k = -d, . . . , d,
where d determines the order of the reconstruction (P i,j is a piecewise polynomial function of degree 2d in this case).

Another example of reconstruction we propose in this work is based on Hermite reconstruction. Let us introduce some notation. First, we introduce edge average values

f n i ± ,j
1 ∆y

y j+1/2 y j-1/2 f n (x i±1/2 , y)dy, f n i,j ± 1 ∆x x i+1/2 x i-1/2
f n (x, y j±1/2 )dx and pointwise values

f n i ± ,j ± f n (x i±1/2 , y j±1/2 ).
We need to express this quantity as a function of the known values f n i,j . For this, we first compute the edge values of the form

f n i ± ,j = s ± =r ± a ± f n i+ ,j , f i,j ± = s ± =r ± a ± f n i,j+ ,
and then the pointwise values

f n i ± ,j ± = s ± =r ± a ± f n i+ ,j ± = s ± =r ± a ± f n i ± ,j+ ,
with suitable coefficients a ± , = r ± , . . . , s ± . The integers r ± , s ± ∈ Z are generally linked to the degree of approximation. Now, for given values of r -≤ 0, s -≥ 0, we can compute coefficients a -, = r -, . . . , s -in the following way. We take

a -= - k=r - s -+1 j=r -, j ∈{0,k} (-j) s -+1 j=r -, j =k (k -j) , = r -, . . . , -1,
and

a -= s -+1 k= +1 s -+1 j=r -, j ∈{0,k} (-j) s -+1 j=r -, j =k (k -j) , = 0, . . . , s -.
In a similar way, for given values of r + ≤ 0, s + ≥ 0, the coefficients a + are given by (for = r + , . . . , s + )

a + = - -1 k=r + -1 s + j=r + -1, j ∈{0,k} (-j) s + j=r + -1, j =k (k -j) , = r + , . . . , 0,
and

a + = s + k= s + j=r + -1, j ∈{0,k} (-j) s + j=r + -1, j =k (k -j) , = 1, . . . , s + .
For p ∈ N * , we will take

r -= - p 2 , r + = - p -1 2 , s -= p -1 2 , s + = p 2 .
Using all these expressions, we then can have an approximation

f n i,j (u, v) f n (x i-1/2 +u∆x, y j-1/2 +v∆y), for 0 ≤ u, v ≤ 1.
For this, we can compute intermediate values

f n i (u, j ± ) f n (x i-1/2 + u∆x, y j±1/2 ), f n i (u, j) 1 ∆y y j+1/2 y j-1/2 f n (x i-1/2 + u∆x, y)dy, which reads, for k ∈ {j ± , j} f n i (u, k) = f n i -,k +u(f n i + ,k -f n i -,k )+u(1 -u)(6f n i,k -3(f n i + ,k + f n i -,k )).
We then obtain the 2D reconstruction for 0 ≤ u, v ≤ 1

f n i,j (u, v) = f n i (u, j -) + v(f n i (u, j + ) -f n i (u, j -)) +v(1 -v)(6f n i (u, j) -3(f n i (u, j + ) + f n i (u, j -))
). ( 6) Let us remark that we start to reconstruct in the x direction but it is possible to first reconstruct in the y direction by defining

f n j (k, v), k ∈ {i ± , i} (which is the counterpart of f n i (u, k), k ∈ {j ± , j}).
For each choice of p, a specific Hermite method is deduced. Among them, we distinguish odd and even values of p. In ( 6), we can observe that if p is odd, r ± = -s ± so that the reconstruction is continuous: this corresponds to PPM type reconstruction (see [START_REF] Colella | The Piecewise Parabolic Method (PPM) for gas-dynamical simulations[END_REF], [START_REF] Crouseilles | Conservative semi-Lagrangian method for the Vlasov equation[END_REF]). If p is even, the reconstruction is no longer continuous since r ± = -s ± : this corresponds to LAGH type reconstruction. We refer to [START_REF] Glanc | Approximation numérique de l'équation de Vlasov par des méthodes de type remapping conservatif[END_REF], [START_REF] Crouseilles | 2D advective and conservative semi-Lagrangian schemes upwind Hermite type reconstruction[END_REF] for more details.

The computation of X(t n ; t n+1 , A i,j ) needs some details since it is not possible to compute exactly this volume. In practice, we only compute the feet of the characteristics ending at the 4 corners of the cell A i,j : (x i±1/2 , y j-1/2 ) and (x i±1/2 , y j+1/2 ). This is done by computing X(t n ; t n+1 , (x i , y j )) and X(t n ; t n+1 , (x i±1/2 , y j±1/2 )) is deduced by interpolation. Hence, the computation of the feet of the characteristics is performed as in the standard semi-Lagrangian method. This determines a quadrilateral volume a i,j ≈ X(t n ; t n+1 , A i,j ) which will be supposed to be convex. This can be ensured providing that the time step ∆t is small enough.

The main step consists in the mesh intersection between the Lagrangian and Eulerian meshes. This step is not detailed here but we refer to [START_REF] Glanc | Approximation numérique de l'équation de Vlasov par des méthodes de type remapping conservatif[END_REF].

The last step is the reconstruction step. Considering a i,j = ∪ k, =1 a k, i,j where a k, i,j = a i,j ∩ A k, , we rewrite (5) as

f n+1 i,j = 1 ∆x∆y X(t n ;t n+1 ,Ai,j ) f n (x, y)dxdy, = 1 ∆x∆y ai,j f n (x, y)dxdy, = 1 ∆x∆y k, a k, i,j f n (x, y)dxdy. (7)
From this last expression, we apply the Green theorem which reformulates an integral over a volume into an integral over segments. Then, it comes

1 ∆x∆y k, a k, i,j f n (x, y)dxdy = ∂a k, i,j [P (x, y)dx + Q(x, y)dy],
where ∂a k, i,j denotes the boundary (which is composed of edges) of a k, i,j and P, Q satisfy -∂ y P (x, y)

+ ∂ x Q(x, y) = f n (x, y).
Finally, once the piecewise polynomial function has been determined, one can write

f n k, (x, y) = p,q c (p,q) k, x p y q , ∀(x, y) ∈ A i,j , so that a k, i,j f n k, (x, y)dxdy = p,q c (p,q) k, a k, i,j
x p y q dxdy.

We then have to compute the integral on the right hand side. To do that, we need some notations. For each subcell a k, i,j , we introduce its corners (x h , y h ), h = 0, . . . , S -1 and we suppose that x 0 = x S , y 0 = y S . Let us compute a k, i,j

x p y q dxdy for some given (p, q) using the Green theorem. The case p = q = 1 leads to

a k, i,j xydxdy = S h=1 (x h ,y h ) (x h-1 ,y h-1 )
x 2 y 2 dxdy

= 1 2 S h=1 1 0 [x h-1 + t(x h -x h-1 )] 2 ×[y h-1 + t(y h -y h-1 )] dt (y h -y h-1 ) = (y h -y h-1 ) 2 S h=1 [y h-1 (3x 2 h-1 + 2x h-1 x h + x 2 h ) +y h (x 2 h-1 + 2x h-1 x h + 3x 2 h )],
where the changes of variables x = x h-1 + t(x h -x h-1 ) and y = y h-1 + t(y h -y h-1 ) have been done. The algorithm of the conservative method to solve (2) then reads compute the 4 corners of a i,j by solving the characteristics ending at the 4 corners of A i,j backward in time, compute the intersection between the Lagrangian mesh and the cartesian axis, compute the corners of a k, i,j ((x h , y h ), h = 1, . . . , S), compute the piecewise polynomial function using a Lagrange or Hermite reconstruction, compute a k, i,j

x p y q dxdy with the Green theorem, compute f n+1 i,j with [START_REF] Davidson | Physics of non neutral plasmas[END_REF].

Standard numerical method

We recall the main steps of the 2D semi-Lagrangian method proposed in [START_REF] Sonnendrücker | The Semi-Lagrangian Method for the Numerical Resolution of the Vlasov Equation[END_REF].

To go from a known function f n i,j on the 2D mesh (x i , y j ) to f n+1 i,j , the semi-Lagrangian is based on the property that the solution is constant along the characteristics, and can be written as follows compute the feet of the characteristics X(t n ; t n+1 , (x i , y j )), solution of (3) compute

f n+1 i,j
= f n (X(t n ; t n+1 , (x i , y j ))),

I(f n )(X(t n ; t n+1 , (x i , y j )))
where I is an interpolation operator (cubic splines for example) since X(t n ; t n+1 , (x i , y j )) does not coincide with (x i , y j ) in general.

Extension to polar coordinates The two methods have been validated in classical test cases in cartesian configurations (using Vlasov-Poisson or guiding-center models). To deal with polar coordinates and non homogeneous Dirichlet boundary conditions, some adaptations are needed, similar to techniques proposed in [START_REF] Latu | Accuracy of unperturbed motion of particles in a gyrokinetic semi-Lagrangian code[END_REF] for the FSL method.

Since the new conservative method solves the conservative form of the transport equation, one has to consider the conservative solution g(t, r, θ) = rf (t, r, θ) which is obtained by multiplying the solution by the Jacobian of the transformation. Then g solves the following conservative form

∂ t g -∂ r ∂ θ Φ r g + ∂ θ ∂ r Φ r g = 0. ( 8
)
The conservative method is then applied to this equation. Since the BSL method solves the advective form, it does not require this adaptation. Another point concerns boundary conditions of the unknown. To overcome this delicate issue, instead of working with g, we work with the unknown (g -g eq ) = r(ff eq ) which satisfies homogeneous Dirichlet boundary conditions in our applications since the solution does not depart far from equilibrium. This is also applied to the standard BSL method.

Diocotron instability test case

This section is devoted to the presentation of the diocotron instability test case [START_REF] Madaule | Semi-Lagrangian simulations of the diocotron instability[END_REF] which enables us to validate our new conservative method on a non trivial 2D polar geometry and compare it to analytical results and to the BSL method. Let us mention the works [START_REF] Krane | Suppression of coherent structures in sheared plasma flows[END_REF] and [START_REF] Marchand | Numerical simulation of the excitation and evolution of high-azimuthal-mode-number coherent vortices in hollow magnetized electron columns[END_REF], in which the guiding-center model is solved numerically in its polar version, using a slightly different initial condition.

We consider ρ(t, r, θ), a solution of (1), with Φ(t, r, θ) solving the Poisson equation:

-∂ 2 r Φ - 1 r ∂ r Φ - 1 r 2 ∂ 2 θ Φ = γρ, (9)
where γ ∈ {-1, 1}, according to the charge of the considered particles. The guiding-center model ( 1)-( 9) satisfied by ρ(t, r, θ) is associated to an initial condition ρ(t = 0, r, θ) (see [START_REF] Davidson | Physics of non neutral plasmas[END_REF])

ρ(t = 0, r, θ) =    0, r min ≤ r < r -, 1 + cos( θ), r -≤ r < r + , 0, r + ≤ r ≤ r max , (10)
where = 10 -6 (linear case) or = 0.5 (nonlinear case).

We also need to define boundary conditions. In the θ ∈ [0, 2π] direction, periodic boundary conditions are used so that a Fourier transform is applied to (9) whereas different conditions will be considered in the radial r ∈ [r min , r max ] direction. We recall the notation Ω = [r min , r max ]×[0, 2π].

Radial boundary conditions for Φ

In this subsection, different radial boundary conditions are proposed to solve [START_REF] Glanc | Approximation numérique de l'équation de Vlasov par des méthodes de type remapping conservatif[END_REF]. In particular, their influence on the conservation of electric energy

E(t) = Ω r|∂ r Φ| 2 + 1 r |∂ θ Φ| 2 drdθ, (11)
and on total mass

M(t) = Ω rρ drdθ, (12)
is discussed. First, homogeneous Dirichlet boundary conditions are imposed at r max :

Φ(t, r max , θ) = 0, ∀t ≥ 0, θ ∈ [0, 2π].
Second, at r min , we consider two kinds of boundary conditions. We take either the homogeneous Dirichlet boundary condition

Φ(t, r min , θ) = 0, ∀t ≥ 0, θ ∈ [0, 2π],
or the following condition, named Neumann mode 0:

the homogeneous Neumann boundary condition at r min for the Fourier mode 0 in θ (∂ r Φ 0 (t, r min ) = 0) For these two choices of radial boundary conditions, one can prove (see [START_REF] Madaule | Semi-Lagrangian simulations of the diocotron instability[END_REF]) that the electric energy and the total mass are preserved in time

∂ t E(t) = ∂ t M(t) = 0,
which provides a good indicator to compare and validate the numerical methods. A last configuration has also been explored, in which homogeneous Neumann boundary conditions are imposed for all the modes ∂ r Φ k (t, r min ) = 0, ∀k.

Let us remark that however in this case, the electric energy and the total mass are not preserved with time at the continuous level (see [START_REF] Madaule | Semi-Lagrangian simulations of the diocotron instability[END_REF]).

For the solution ρ, Dirichlet type boundary conditions are considered: ρ(t, r ≤ r min , θ) = ρ(t, r min , θ), ρ(t, r ≥ r max , θ) = ρ(t, r max , θ).

Instability rates

In this subsection, the instability rates are computed from the linearized problem, considering the different radial boundary conditions.

Considering a radial equilibrium (n 0 (r), Φ 0 (r)) of ( 1)-( 9) linked by the Poisson equation

-∂ r (r∂ r Φ 0 ) = rn 0 ,
we then search for an approximate solution to(1)-( 9) of the form ρ(t, r, θ) n 0 (r) + n 1, (r) exp(i θ) exp(-iωt/γ) Φ(t, r, θ) γΦ 0 (r) + γ Φ 1, (r) exp(i θ) exp(-iωt/γ).

Inserting this particular solution into (1)-( 9) leads to the following system

         -iω n 1, -i ∂ r n 0 r Φ 1, + i ∂ r Φ 0 r n 1, = 0, -∂ 2 r Φ 1, - ∂ r Φ 1, r + 2 r 2 Φ 1, = n 1, . (13)
We express n 1, in terms of Φ 1, to get an equation for Φ 1, :

-iω + i ∂ r Φ 0 r (-∂ 2 r Φ 1, - ∂ r Φ 1, r + 2 r 2 Φ 1, ) = i ∂ r n 0 r Φ 1, . (14)
We can proceed in the numerical solution of the problem by making a discretization in the radial direction. With a uniform radial mesh r i = r min + i∆r, ∆r = (r maxr min )/N r and N r ∈ N * , we introduce the following notations and approximations

Φ 1, (r i ) = φ i , ∂ r Φ 1, (r i ) = φ i+1 -φ i-1 2∆r , ∂ 2 r Φ 1, (r i ) = φ i+1 -2φ i + φ i-1 ∆r 2 .
Then, ( 14) is approximated by Aφ = cBφ, with c = ω/ and where φ = (φ 0 , φ 1 , . . . , φ N ) T ∈ R Nr+1 and A, B are two matrices of size N r + 1. The problem is then reduced to find the eigenvalues c of the problem B -1 Aφ = cφ. For each , we look for the eigenvalue c which has the greatest (strictly positive) imaginary part. We refer to [START_REF] Davidson | Physics of non neutral plasmas[END_REF] and [START_REF] Madaule | Semi-Lagrangian simulations of the diocotron instability[END_REF] for more details.

Linear case.

Analytical growth rates | (ω)| of the linear phase are obtained for all the discussed boundary conditions (the initial condition being ( 10) with = 10 -6 ) and comparison is made with the BSL method on a 128×128 grid, ∆t = 0.05, r min = 1, r max = 10, γ = 1. In the following tables, comparisons between analytical and numerical growth rate for different modes and different r ± are performed -Table 1 refers to the first case (homogeneous Dirichlet boundary condition). -Table 2 refers to the second case (Neumann mode 0).

-Table 3 refers to the third case (homogeneous Neumann boundary condition).

One can observe that the numerical instability rates are very close to the analytical ones (relative error of about 0.1%).

Nonlinear case.

Let us now consider the initial condition ( 10) with = 0.5 and study the influence of radial boundary conditions.

To that purpose, we compare the new CSL method (the LAGH-9 method is used for the reconstruction) with the standard BSL method (cubic splines interpolation is used). Let us recall that the unknowns for CSL are centered at the vertices r i , θ j , following (4). We first focus on the Neumann mode 0 case: shown in Fig. 2 and3, we plot the time history of the electric energy E(t) and of the relative total mass (M(t) -M ex )/M ex (with M ex = rmax rmin 2π 0 ρ(t = 0, r, θ)rdrdθ = 28.27433389). A standard rectangle quadrature is used to approximate E and M given by ( 11) and [START_REF] Latu | Accuracy of unperturbed motion of particles in a gyrokinetic semi-Lagrangian code[END_REF]. We compare the BSL method with different meshes and the new CSL method with N r = 261, N θ = 256 and ∆t = 0.1 (N r = 261 is chosen to ensure that r ± coincide with a mesh point). First, we can see that the new CSL method preserves the total mass better than the standard BSL method. Indeed, even with a refined mesh N r = 1026, N θ = 1024 (purple curve in Fig. 3), the error in the BSL on this diagnostic is more important than the error in the CSL. The CSL method is not exactly conservative because of the boundary effect on ρ: some loss or gain of mass can occur at r min and r max . It is noted that the first order convergence of the trapezoidal quadrature rule in the r direction: refining in this direction makes the relative total mass converge to zero at the initial time. Secondly, looking at the time evolution of the energy (see 2), the same conclusions are realised since CSL is very well behaved compared to BSL: the variation of total energy is about 3×10 -2 for CSL whereas for BSL, it is about 5 × 10 -2 using 16 for both.

Then, we focus on the Neumann case. We plot the same diagnostics in Fig. 4 and5. In this Neumann configuration, the total mass and energy are not preserved at the continuous level, but it is interesting to observe how the numerical methods behave. Clearly, this is different since the CSL method does not preserve these two quantities as well as in the previous case (variations of about 0.4 for the total energy). However, it is difficult to define a good solution since the BSL method displays quite different behavior. A deeper analysis is needed and is not within the limits of this article.

Drift-kinetic Vlasov equation

This last section is devoted to the validation of the new conservative method in the context of 4D the drift-kinetic equation. Indeed, we look for f = f (t, r, θ, z, v) satisfying the following 4D slab drift-kinetic equation (see [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]) for (r, θ, z, v) ∈ Ω ×[0, L]×R. The self-consistent potential φ = φ(r, θ, z) solves the quasi neutrality equation with φ = 1 L L 0 φ(r, θ, z)dz. Periodic boundary conditions are considered in the θ and z directions whereas in the radial direction, Neumann mode 0 is imposed (see [START_REF] Latu | Improving conservation properties in a 5D gyrokinetic semi-Lagrangian code[END_REF]). For the solution f , periodic boundary conditions are imposed in the θ, z, v directions, whereas in the radial r directions, Dirichlet boundary conditions are considered.

∂ t f - ∂ θ φ r ∂ r f + ∂ r φ r ∂ θ f + v∂ z f -∂ z φ∂ v f = 0, (15)
-∂ 2 r φ + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ + 1 r 2 ∂ 2 θ φ + 1 T e (r) (φ -φ ) = 1 n 0 (r) R f dv -1, (16)
The initial function is given by

f (t = 0, r, θ, z, v) = f eq (r, v) 1 + exp - (r -r p ) 2 δr cos 2πn L z + mθ ,
where the equilibrium function is

f eq (r, v) = n 0 (r) exp(-v 2 2Ti(r) ) (2πT i (r)) 1/2 .
The radial profiles {T i , T e , n 0 } have the analytical expressions P(r) = C P exp -κ P δr P tanh( r -r p δr P ) , P ∈ {T i , T e , n 0 }, where the constants are

C Ti = C Te = 1, C n0 = r max -r min rmax rmin exp(-κ n0 δr n0 tanh( r-rp δrn 0
))dr .

Finally, we consider the parameters of [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF] (MEDIUM case)

Global algorithm

To deal with the 4D drift-kinetic model, a directional splitting is performed as in [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF]. As said before, the standard semi-Lagrangian method and the new conservative semi-Lagrangian method have common steps, in particular because in the 1D z and v configuration, we deal with the constant advection coefficient for which the BSL and CSL methods are equivalent (see [START_REF] Crouseilles | Conservative semi-Lagrangian method for the Vlasov equation[END_REF]). The only difference between the two methods lies in the 2D advection in (r, θ). We then call BSL2D the method proposed in [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF] and CSL2D the new conservative method.

The main steps of the algorithms BSL2D and CSL2D are solve the 1D advection ∂ t f + v∂ z f = 0 on ∆t/2 with the standard semi-Lagrangian method (periodic cubic splines interpolation), solve 1D advection ∂ t f -∂ z Φ∂ v f = 0 on ∆t/2 with the standard semi-Lagrangian method (periodic cubic splines interpolation), solve the quasi neutrality equation ( 16) to compute the potential Φ, with FFT in θ, z and second order finite difference in r, compute the derivatives (∂ r Φ, ∂ θ Φ, ∂ z Φ) of Φ with cubic splines in r, θ and second order finite differences in z, -2D advection in (r, θ) (1), with BSL or CSL, solve 1D advection ∂ t f -∂ z Φ∂ v f = 0 on ∆t/2 with the standard semi-Lagrangian method (periodic cubic splines interpolation), solve the 1D advection ∂ t f + v∂ z f = 0 on ∆t/2 with the standard semi-Lagrangian method (periodic cubic splines interpolation).

Numerical results

In this section, we compare the standard BSL method to the new conservative method. To do that, we consider the 4D drift-kinetic model in polar geometry ( 15)-( 16). The practical implementation is done in a modular semi-Lagrangian library [19]. We are interested in the time history of

Φ 2 (r p , θ, z)dθdz, (17)
where r p = (r max + r min )/2, of the total mass

M(t) = rmax rmin 2π 0 L 0 R f (t, r, θ, z, v)rdvdzdθdr,
and of the total energy

E(t) = rmax rmin 2π 0 L 0 R v 2 2 f (t, r, θ, z, v)rdvdzdθdr + rmax rmin 2π 0 L 0 R f (t, r, θ, z, v)Φ(t, r, θ, z)rdvdzdθdr.
The first quantity ( 17) is known to exponentially increase with time. The growth rate can be obtained by a generalization of the approach performed in [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF]. Hence, growth rates obtained in [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF] can be compared to numerical results. The total mass M(t) is preserved by the model and it is an important point for the methods to preserve this quantity as much as possible.

In Fig. 6, we plot the time history of ( 17) obtained with BSL2D and CSL2D. Different variants of the CSL2D method are shown, corresponding to different reconstructions presented in section 2: p = 3 which corresponds to third order Lagrange reconstruction, p = 4 (resp. 6) which corresponds to PPM (resp. PPM2), p = 17 which corresponds to LAGH17, and Lauritzen which corresponds to the reconstruction proposed in [START_REF] Lauritzen | A conservative semi-Lagrangian multi-tracer transport scheme (CSLAM) on the cubed-sphere grid[END_REF]. One can observe that the linear part is well reproduced by all the method since the instability rate of all the methods is very close to the analytical one proposed by [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF].

In Fig. 7, the time history of the total mass M(t) (which is theoretically preserved with time) is plotted for the different methods (BSL2D and CSL2D for different reconstructions). The final time is 8000 chosen so that the linear theory is not available any more. We can first observe that up to time 3000, the total mass is well preserved by all methods. This time interval corresponds to the linear phase in which few phenomena occur. After that time, the methods behave differently. In particular, BSL2D looses the conservation of total mass (even with a relative error of about 0.1%). On the other side, CSL2D presents better behavior (relative error of about 0.01%). Among the different reconstructions, CSL2D p = 6 (which corresponds to PPM2) behaves in the best way. Fig. 8 shows the time history of E(t) (which is theoretically preserved with time) obtained with the different methods. This invariant is quite difficult to preserve numerically and we can observe that, as for the total mass, the CSL methods display very good behavior compared to the BSL method since the total energy is better conserved (about 0.1% for CSL and about 1% for BSL). One interesting point is the very similar behavior of the total mass and total energy. It seems that the loss of the mass conservation has a strong influence on the total energy conservation. Obviously we have to go further into the comparisons but these first results are quite promising.

Conclusion

In this work, we proposed the validation of a new conservative semi-Lagrangian method CSL, on two non trivial test cases in polar geometry. The first test concerns a guiding-center model for which boundary conditions in the radial directions are discussed. It appears that some good choices lead to conservation of electric energy and total mass. This is confirmed by the numerical results in which we observe a very good behavior of the new method compared to BSL. The second test concerns a 4D drift-kinetic model which can be viewed as a first step to describe tokamak plasmas. For these two models, a linear analysis can be performed and instability rates can be obtained to validate the numerical methods in the linear phase.

The new CSL method seems to be a good alternative to the standard pointwise semi-Lagrangian method which is not conservative in these configurations. More comparisons are performed in [START_REF] Glanc | Approximation numérique de l'équation de Vlasov par des méthodes de type remapping conservatif[END_REF] and [START_REF] Crouseilles | 2D advective and conservative semi-Lagrangian schemes upwind Hermite type reconstruction[END_REF] in these two intermediate configurations (guiding-center and 4D drift-kinetic) to confirm the good behavior of the conservative semi-Lagrangian method.

Fig. 1 .

 1 Fig. 1. Illustration of the intersection between the Lagrangian mesh a k, and the Eulerian mesh Ai,j.

2π 0 ∂

 0 r Φ(t, r min , θ)dθ = 0, the homogeneous Dirichlet boundary condition at r min for the other Fourier modes ( Φ k (t, r min ) = 0, ∀k = 0) Φ(t, r min , θ) = 1 2π 2π 0 Φ(t, r min , θ )dθ .

Fig. 2 .Fig. 3 .Fig. 4 .

 234 Fig. 2. Time evolution of E(t) for different grid sizes and ∆t: Neumann mode 0 case.

Fig. 5 .

 5 Fig. 5. Time evolution of the relative total mass (M(t) -Mex)/Mex for different grid sizes and ∆t: Neumann case.

r

  min = 0.1, r max = 14.5, κ n0 = 0.055, κ Ti = κ Te = 0.27586, δr Ti = δr Te = δr n0 2 = 1.45, = 10 -6 , n = 1, m = 5, L = 1506.759067, r p = r min + r max 2 , δr = 4δr n0 δr Ti .

Fig. 6 .

 6 Fig. 6. Time evolution of Φ 2 (rp, θ, z)dθdz on 32 × 32 × 32 × 64 grid, with ∆t = 8 and vmax = 7.32.

Fig. 7 .

 7 Fig. 7. Time evolution of the total mass on 32 × 32 × 32 × 64 grid, with ∆t = 8 and vmax = 7.32.

Fig. 8 .

 8 Fig. 8. Time evolution of the total energy on 32 × 32 × 32 × 64 grid, with ∆t = 8 and vmax = 7.32.

Table 1 .

 1 Growth rate obtained by the semi-Lagrangian method and by linearized analysis with homogeneous Dirichlet boundary conditions.

		r-r+	(ω)	analytical
	2	4	5	0.2875, t ∈ [26, 55]	0.2887
	3	4	5	0.3667, t ∈ [26, 72]	0.3673
	4	4	5	0.3852, t ∈ [26, 51]	0.3840
	7	6	7	0.3424, t ∈ [26, 61]	0.3375
		r-r+	(ω)	analytical
	2	4	5	0.0719, t ∈ [93, 138]	0.0717
	3	4	5	0.2264, t ∈ [33, 87]	0.2267
	4	4	5	0.2991, t ∈ [34, 69]	0.2988
	7	6	7	0.3338, t ∈ [35, 53]	0.3307

Table 2 .

 2 Growth rate obtained by the semi-Lagrangian method and by linearized analysis with Neumann for mode 0 boundary conditions.

		r-r+	(ω)	analytical
	2	4	5	0.0672, t ∈ [89, 151]	0.0673
	3	4	5	0.2263, t ∈ [32, 89]	0.2266
	4	4	5	0.2991, t ∈ [34, 69]	0.2988
	7	6	7	0.3338, t ∈ [35, 53]	0.3307

Table 3 .

 3 Growth rate obtained by the semi-Lagrangian method and by linearized analysis with Neumann boundary conditions.


Hamida Dridi

Adel Moadhen
email: adel_moadhen@yahoo.fr

Mohamed Lazhar

Haji

L Haji

Comparative SERS study

Keywords: oxidized porous silicon, small size gold nanoparticles, SERS, Bovine Serum Albumin, Rhodamine 6G

Introduction

Optical biosensor research based on Surface Enhanced Raman Scattering (SERS) effect has attracted a great attention. It has been demonstrated that SERS can provide scattering cross sections with several orders of magnitude higher than normal Raman spectroscopy [1]. As a result, ultra low concentration even at the scale of a single molecule was detected [2,3]. Raman amplification can be explained in terms of electromagnetic [START_REF] Tourrel | Raman Microscopy: Developments and Applications[END_REF] and chemical enhancement [START_REF] Otto | [END_REF]. The first effect involves the excitation of surface plasmon (or localized surface plasmon) at or close to resonance. While, the second effect was explained by a charge transfer between the metal surface and the adsorbed probe. Based on the excitation of the Localized Surface Plasmon (LSP), SERS has become an effective tool to analyze molecules. The metallic nanoparticles LSP frequency depends on the particle type (Au, Ag, Cu...), size, shape [6,7] and as well as the refractive index of their surrounding medium, which in turn could have an influence on the SERS results [8].

Conventional SERS substrates have been often fabricated by depositing metallic nanoparticles on a planar surface [9][10][11][12]. Porous material based substrates could be interesting regarding their peculiar properties such as the large specific surface and the opened structure with pore size selective molecules filtering. Furthermore, Porous Silicon (PSi) based structures are promising for sensor application. This material was used as an efficient SERS template for various molecules detection like peptide, benzenethiol, organic molecules for different metals type [13][14][15][16][17][18][19][20][21][22]. Several parameters were studied such as porosity and thickness of the porous layer, metal nature, dipping time in the probe solution and the wavelength excitation to get the maximum of the SERS intensity [13][14][15][16][17][18]. It was reported that for PSi formed from p + type Si, the best SERS result was obtained for high porosity, [19] while for n-type Si the opposite result was obtained [13]. This difference could be attributed to the difference in PSi morphology.

It was also reported that the metal deposition rate is related to the PSi formed from n or p type Si, which influences the detection limit value [23][24]. However and up to now, the number of publications interested on bio-molecules detection based on SERS PSi is still limited. In fact, PSi was used as a SERS substrate for some bio-molecules like horseradish peroxidase (HRP) and peptide. Giorgis et al [14] have showed a detectable concentration of HRP as low as 10 -8 M. Specific complex peptide-antibody detection was performed by Virga et al., [25] which have also demonstrated a better SERS result using silver nanoparticles than silver thin film.

Furthermore, we have to mention, that surface functionalization process is interesting for some target molecules and for the metal deposition step, specifically in the case of colloidal metallic nanoparticles [14,18,26]. However, it was reported that SERS of Bovine Serum Albumin (BSA) on porous or planar substrates does not show an obvious amplification of the vibrational modes intensities [9,27]. Indeed, David et al [9] have demonstrated that optimizing the size (170 nm) and the shape of the lithographed gold nanocylinders lead to a SERS spectrum of BSA at low concentration (10 -4 M).

In this work, we report a fabrication of a SERS-active substrate obtained by immersion of oxidized porous silicon (OPSi) in a small size (5 nm) colloidal gold nanoparticles (GNP) solution for the detection of BSA proteins and R6G molecules. Firstly, we have investigated the direct interaction between GNP and BSA to verify the bio-conjugation and to adjust the Raman excitation wavelength to the LSP band. Then, SERS results of BSA proteins from two types of OPSi surfaces (silanized and unsilanized) were compared. In this case, the effect of the surface functionalization on the SERS efficiency was discussed. We have also studied the BSA solution pH effect on the SERS response. The SERS substrate activity for BSA and R6G was obtained for concentrations as low as 10 -8 M.

Experimental

1-Preparation of functionalized oxidized porous silicon samples

PSi samples were first elaborated by electrochemical anodization at room temperature and in the darkness of a highly P doped (100)-oriented silicon wafer having a resistivity of 3-7 mΩ.cm. The PSi was prepared by using a 20 % aqueous-alcoholic hydrofluoric acid solution (1:2:2). The current density and the anodization time were 80 mA/cm 2 and 90 s respectively.

After that, the PSi layers were oxidized for one hour in wet oxygen at 950°C to stabilize the surface [28][29][30] and also as a prerequisite step for the chemical functionalization [31,32]. By fitting the optical reflectivity data, we have obtained an estimated porosity of about 35 %. The OPSi layer thickness was equal to 4.8 µm and the pore diameter was of around 30 nm (determined by Scanning Electron Microscopy).

OPSi surface could be activated by different silanization molecules (based on silanol groups)

showing functional terminated groups (thiol, carboxyl, amines...). In our case, we have utilized 3-aminopropyltriethoxysilane (APTES) molecules. APTES molecules presented several functional amine terminals, which can enhance the GNP adsorption on the surface by electrostatic forces. Firstly, the OPSi samples were immersed in a 2 % hydro-alcoholic APTES solution (methanol-deionised water (v/v)) for 20 minutes. After rinsing with deionised water (18 MΩ) and drying under with nitrogen flux (N 2 ), the samples were heated at 100 °C in air to remove the solvent traces from the pores. Also and for comparison, some samples without silanization were used in this work.

2-Fabrication of the SERS-active substrate

The spherical GNP deposition was carried out by the immersion of the silanized and unsilanized OPSi surfaces in a colloidal solution of a mono-disperse GNP (5 nm size) for three hours. The small size of GNP could promote their infiltration inside the pores. The samples were then rinsed with deionised water (18 MΩ) and dried with N 2 . Finally, the samples were incubated at 80 °C for 30 min to remove the solvent traces from the pores.

In this study, elaborated samples were incubated in an aqueous solution (for different pH) of a Bovine Serum Albumin (BSA) protein, for 2 days at room temperature. BSA is a protein with a 66.4 kDa molecular weight and it contains 583 amino acids.

In addition, GNP coated OPSi surfaces were immersed in an aqueous Rhodamine 6G (R6G) solution for 20 min.

GNP colloidal solution, APTES, R6G and BSA were purchased from Sigma Aldrich and used without further purification.

3-Instrumentation

UV-Vis-NIR absorption spectra were recorded using a Perkin Elmer spectrophotometer lambda 900. The optical system is equipped with a deuterium lamp and a monochromator (wavelength range 200-2000 nm). The system resolution is of about 1 nm. The AFM scanner (Digital Instruments Nanoscope) was mounted on the sample stage and imaged in tapping mode. The Raman scattering measurements were performed at room temperature on a Micro-Raman system HORIBA Jobin Yvon (LabRAM HR800). The BSA measurements were carried out using an Argon laser excitation 488 nm with a power at the sample surface equal to 47 mW. The excitation was chosen as the closest one to the 5 nm-GNP LSP band. The experiments with R6G were carried out using an He-Ne laser excitation 632.8 nm with a power at the surface less than 1 mW. This wavelength is more suitable for the study of R6G to avoid its fluorescence. The spectrometer resolution is equal to 1 cm -1 . A CCD camera was used to record Raman spectra. Data were collected with an integration time of a few tens of seconds and the intensities were calculated with the unit of counts per second (cps). All Raman and SERS spectra presented in this paper were averaged from five locations on each sample. Some of the Raman spectra were smoothed and we have done the base line.

Results and discussion

In Fig. 1 we report the absorption spectra of the various solutions of colloidal GNP and GNP-BSA conjugate. The GNP's Plasmon band is localized at 523 nm which is related to the collective free electron oscillation on the surface of the GNP (Localized Surface Plasmon Resonance: LSPR). The BSA's absorption band is located at 279 nm. After adding BSA molecules to GNP, we notice a red shift of δλ 1 = 4 nm and a decrease of the integrated intensity of the LSPR band. Such modification could be related to a change of the refractive index of the surrounding medium of the GNP and/or molecular interactions. In addition, the BSA absorption band showed a blue shift of δλ 2 = 3 nm which is an evidence of some GNP-BSA bio-conjugation. It is known that the BSA conformation and its net charge depend on pH value. The pH modification could influence the absorption band position and the interaction types with GNP [33]. In our case, no changes of the LSP frequency were observed for three pH BSA solutions (3.3, 4.9 and 6.9). This result is not consistent with the fact that we should observe a significant shift in the case where the BSA charge is suitable for electrostatic interaction with GNP. In Fig. 2, we present the topographic AFM images of OPSi, after GNP deposition (GNP-OPSi) and after BSA grafting (BSA-GNP-OPSi). The AFM image, in figure 2.b, reveals a small surface roughness variation of the OPSi after GNP deposition. This is probably due to, very low gold nanoparticles concentration (5.5 10 13 mL -1 ) in the colloidal solution. After grafting BSA (figure 2.c), the AFM image did not show an obvious change of the surface roughness. However, we noticed the presence of some agglomeration on the surface (the regions surrounded with dashed circles), which could be attributed to the residual BSA molecules after rinsing. Dashed circles correspond to BSA residue.

Raman and SERS measurements were performed on such silanized OPSi substrate as well as on unsilanized one, for comparison. Silanization step was supposed to give a homogeneous distribution of gold nanoparticles [34]. Since GNP were assumed to be attached on the silanized surface via electrostatic bonds with the amines terminals of APTES molecules.

However, the effect of the functionalization process depends on the metal deposition way and the substrate type (glass, porous materiel, rough substrate…).

Fig. 3 shows the Raman and SERS spectra of BSA at 10 -6 M concentration in the case of silanized substrate. In Fig. 3-a, Raman spectrum of BSA shows the appearance of some vibrational modes, e.g, the amide I which is located at 1660 cm -1 [35]. However, this result was hardly obtained from some regions, which we have attributed to a low adsorbed BSA and/or a strong homogeneity of the adsorbed molecules. The large band between 1400-1500 cm -1 could correspond to the contribution of the ν(CH 2 ) def (1450 cm -1 ) and ν(C-N) (1412 cm -1 ) modes assigned to either BSA and APTES molecules [31]. However, after BSA grafting, the intensity of these bands decreased drastically. We assumed that after rinsing, BSA molecules not linked to APTES molecules were almost removed. Thus, this intensities decrease, can be an indication of the BSA molecules attachment strength.

In Fig. 3 (b-c), we report the SERS spectra before and after BSA grafting. Compared to the spectrum discussed above, the main vibrational modes intensities are clearly amplified. After GNP deposition, this amplification corresponds to the APTES modes and is attributed to a GNP LSPR effect. Also, we observe an increase of the vibrational modes intensities located at 1450 cm -1 and 1350 cm -1 after adding the proteins. However, the SERS effect is not as efficient as for APTES molecules.

Moreover, no clear BSA signature was observed in this case (amide I (1660 cm -1 ), amide II (1550 cm -1 ) and amide III (1260 cm -1 )). To quantify this result, we have calculated the integrated intensity for the 1450 cm -1 mode from SERS (after subtraction of the signal at GNP step) and Raman spectra. Thanks to the contribution of the GNP, the integrated intensity of ν(CH 2 ) mode of the BSA is amplified 2 times. A silanized substrate permits us to obtain a good homogeneity of SERS results but without a good resolution of BSA's specific bands (for instance amide I). BSA modes are located in quite the same frequency region than those of APTES. Thus, it is very difficult to distinguish the vibrational modes of BSA and APTES molecules. To avoid this, one can look for a well adapted functionalized process for instance by choosing a silanization molecule having a different signature than that of BSA protein. This manner could be a rather complicated process. So, we have chosen an unsilanized substrate and GNP were directly deposited on the OPSi substrate. In this case, we can suspect a weak interaction between GNP and hydroxyls terminals of the OPSi substrate surface.

Figures 4-a, 4-b show Raman and SERS spectra of BSA at 10 -6 M concentration at a pH = 6.9

in the range of 1100-1800cm -1 and 2600-3200cm -1 , respectively. SERS spectra show clearly the BSA signature, through the vibrational modes corresponding to amide II (1550 cm -1 ), amide I (1660 cm -1 ), -C=O (1730 cm -1 ) and C-H stretching mode (in the range of 2800-3030 cm -1 ) [36], while Raman spectra show a weak signal of BSA molecules in 2800-3030 cm -1 and no detectable signal of BSA in the range 1100-1800cm -1 . This result proves the efficiency of the SERS substrate based on OPSi to obtain a high sensitivity detection of BSA protein for a pH=6.9. In order to analyze the degree of grafting of the BSA molecules in the pores of the OPSi, we followed the evolution of the C-H mode of BSA, deduced from the Raman spectra recorded in volume. Figure 5 shows SERS spectra of BSA protein at different depths on the cross section of the OPSi. The in-depth Raman measurements shows an enhancement of the signal of BSA modes at the middle (z = 2 -3 µm) of the layer. This indicates that BSA molecules are deeply incorporated into the OPSi layer and it may take places over all the pore walls. It has been reported, that the net charge of BSA and its conformation, strongly depend on its solution pH [37]. The BSA iso-electric point (pI), where the protein is electrically neutral, is between 4.6-4.8 [START_REF] Sigma | Product Information[END_REF]. Thus, the adsorption of BSA molecules on the OPSi surface should be different for various pH (below and above the pI). So, we have studied the BSA pH effect (3, 4.9 and 6.9) on the SERS signal, in the case of unsilanized substrate.

The results are reported in Fig. 6 and show high SERS intensities at a pH close to the pI (4.9).

So, we can deduce that BSA molecules were easily adsorbed on the substrate surface. In addition, we have obtained a better result for a pH=3 than for a pH=6.9. This could be probably related to the net positive charge of the protein at a pH lower than pI, which improves the interactions with GNP via citrate sodium. As a proof of the efficiency of our SERS substrate for a pH= 4.9, we have also reported in Fig. 6 the Raman signal of 10 -6 M BSA molecules. We can see, clearly, the enhancement of all vibrational modes associated to the protein.

In order to estimate the limit of detection (LOD), we have varied the BSA concentration from 10 -6 M down to 10 -8 M and recorded their SERS spectra in the case of unsilanized substrate and for a pH solution equal to 4.9. As it can be seen in Fig. 7, almost all the vibrational bands are visible even at the lowest concentration. However, the band at 1730 cm -1 become indistinguishable, thus we assume that the LOD is up to 10 -8 M. This LOD value is four orders of magnitude lower than those obtained by other authors for BSA molecules [9,27]. This improvement could be related to combined effects of the large surface area of the porous silicon and the small nanoparticles size confined within the pores. Since, the GNP number is low (55 GNP / µm 3 ), it's very difficult to experimentally demonstrate the role of the internal surface of the porous layer. In the other hand and in the goal to verify the efficiency of our SERS substrate, we have used Rhodamine 6G (R6G) molecules. R6G molecules well known to give a high SERS signal, as the other dyes molecules (crystal violet, cyanine...) [25,27]. Fig. 8 shows the SERS spectrum of unsilanized R6G at 10 -8 M concentration. As it can be seen, among the vibrational features, we observed some typical R6G modes at 1128, 1198, 1312, 1356, 1503, 1599 and 1639 cm -1 corresponding to C-H in-plane bend, C-C stretching, aromatic C-C stretching, C-N stretching, aromatic C-C stretching, aromatic C-C stretching and aromatic C-C stretching, respectively [27]. For comparison, R6G normal Raman spectrum correspond to a 10 -6 M

shows just the presence of two vibrational bands located at 1360 and 1508 cm -1 with a very low intensity. In our case, the obtained detection limit of R6G molecules is 10 -8 M. Finally, we can conclude that OPSi is a promising SERS substrate for chemical and biological molecules detection at low concentrations.

Conclusion

SERS based porous silicon substrate represents a viable approach for chemical and biological molecules detection. Spherical gold nanoparticles with small size (5nm) were used to reveal a significant SERS effect for either BSA proteins and R6G molecules at 10 -8 M concentration. We have obtained a better BSA modes resolution in the case of unsilanized OPSi than the silanized one. Although, it has been noted that silanization allows the improvement of target molecule concentration binding on the surface via stronger interaction and it permits also to obtain an homogeneous SERS signal. In our case, the luck of resolution of BSA modes obtained on the silanized substrate could be explained by the interferences between the vibrational modes of APTES and BSA proteins. Furthermore, we have demonstrated that SERS results depend on BSA pH solution. A strong BSA SERS intensity at pH=4.9, corresponding to the iso-electric point of BSA, was obtained.

Figure 1 :

 1 Figure 1: Absorption spectra of colloidal GNP, BSA (10 -5 M concentration) solution and

Figure 2 :

 2 Figure 2: AFM images corresponding to a) OPSi, b) GNP-OPSi and c) BSA-GNP-OPSi.

Figure 3 :

 3 Figure 3: Raman spectrum of BSA (10 -6 M)/APTES/OPSi (a) and SERS spectra of GNP/APTES/OPSi (b) and BSA (10 -6 M)/GNP/APTES/OPSi (c).

1 Figure 4 -

 14 Figure 4-a: Raman spectrum of BSA (10 -6 M)/OPSi (a) and SERS spectrum of

1 Figure 4 -

 14 Figure 4-b: Raman spectrum of BSA (10 -6 M)/OPSi (a) and SERS spectrum of

Figure 5 :

 5 Figure 5: In-depth SERS spectra of BSA/GNP/OPSi in the 2700-3200cm -1 wavelength

Figure 6 :

 6 Figure 6: SERS spectra of BSA (10 -6 M)/GNP/OPSi at pH= 3(a), pH= 4.9 (b) and pH= 6.9 (c) and Raman spectrum of BSA (10 -6 M)/OPSi at pH= 4.9 (d).

1 Figure 7 :

 17 Figure 7: SERS spectra of BSA at pH= 4.9 for 10 -6 M (a) 10 -7 M (b) and 10 -8 M (c).

Figure 8 :

 8 Figure 8: Raman of R6G (10 -6 M)/OPSi (a) and SERS of R6G (10 -8 M)/GNP/OPSi (b).

Acknowledgment

The authors thank the Franco-Tunisian Joint Committee of University Cooperation (CMCU), No. 11G1109 PHC-Utique project for its support and financial of a part of this work.


Anne Bouillard
email: anne.bouillard@ens.fr

Ana Bušić
email: ana.busic@inria.fr

Christelle Rovetta
email: christelle.rovetta@inria.fr

Clones: CLOsed queueing Networks Exact Sampling

Keywords: queueing networks, simulation, coupling from the past

come

INTRODUCTION

Closed queueing networks are largely used in various application domains due to their modeling simplicity and productform stationary distribution in the unlimited capacity case [START_REF] Gordon | Closed queueing systems with exponential servers[END_REF]. This structure is no longer guaranteed when the queues have a finite capacity, due to the effect of blocking. There are different types of blocking, see [START_REF] Onvural | Survey of closed queueing networks with blocking[END_REF] for a detailed discussion. According to the terminology introduced in that survey, our blocking model is RS-RD (repetitive servicerandom destination). Under this blocking policy, the stationary distribution has a product form only if the routing is reversible. When the stationary distribution does not have a product-form, the exact analysis may not be computationally tractable, and we may turn to approximations [START_REF] Balsamo | Queueing networks with blocking: Analysis, solution algorithms and properties[END_REF] or simulation. The main difficulty of the Markov chain Monte Carlo approach is the stopping criterion. Usual tools are the asymptotic variance in the Central Limit Theorem or mixing times [START_REF] Asmussen | Stochastic Simulation: Algorithms and Analysis[END_REF][START_REF] Levin | Markov Chains and Mixing Times[END_REF], unfortunately there is no generic technique for the non-reversible Markov chain case.

In the 1990's, Propp and Wilson introduced a method for sampling a random variable according to the stationary distribution of a finite ergodic Markov chain [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF]: the coupling from the past (CFTP) algorithm. The CFTP algorithm automatically detects and stops when the sample has the correct distribution. The main drawback of the original CFTP is that it considers a coupling from all initial conditions. In the case of closed queueing networks the cardinality of the state space is exponential in the number of queues, which is untractable.

When the network has a product-form distribution, Kijima and Matsui [START_REF] Kijima | Randomized approximation scheme and perfect sampler for closed Jackson networks with multiple servers[END_REF] proposed a perfect sampling algorithm with overall complexity O(K 3 ln(KM )). However, their method strongly relies on the product form representation of the stationary distribution and it cannot be applied to the general case of queues with finite capacity. In a recent paper by Bouillard and al. [START_REF] Bouillard | Perfect sampling for closed queueing networks[END_REF], a new representation of the sets of states has been proposed. This representation is used to derive a bounding chain for the CFTP algorithm for closed queueing networks, that enables exact sampling from the stationary distribution without considering all initial conditions in the CFTP. This method is far more general, as it does not rely on the product-form property.

The complexity of one-step transition of the CFTP algorithm in [START_REF] Bouillard | Perfect sampling for closed queueing networks[END_REF] is O(KM 2 ). In this paper, we derive a new algorithm, based on even more compact representation. The complexity of one-step transition of this new CFTP algorithm is O(KM ). Note that this is the same as the complexity of the computation of the normalizing constant in the product-form case using Buzen's algorithm [START_REF] Buzen | Computational algorithms for closed queueing networks with exponential servers[END_REF]. The overall complexity of our algorithm depends however also on the coupling time of the chain, so it becomes of a practical interest when the network does not have a product form. The main focus of this paper is on algorithmic and implementation issues. We give a detailed description of the matrixbased representations used in our Matlab toolbox Clones.

The paper is organized as follows. In Section 2 we describe the model and give a short overview of results in [START_REF] Bouillard | Perfect sampling for closed queueing networks[END_REF]. We then introduce a new gap-free diagram representation and a new CFTP algorithm with one-step transition complexity of O(KM ). In Section 3 we present our Matlab toolbox Clones and discuss its functionalities and implementation issues. Section 4 is devoted to numerical results. Final remarks and conclusions are contained in Section 5.

BACKGROUND 2.1 Model

Consider a closed network of K ./M/1/C queues with M customers in total. Each queue k ∈ {1, . . . , K} =: Q has a finite capacity C k ≤ M and a service rate ν k . After a customer has been served in queue i, it is directed to queue j with probability pi,j, independently of the current state and past evolution of the network. If queue j is full, then the customer remains in queue i and has to be served again. Its new destination will be chosen independently from j. We assume that the routing matrix P = (pi,j)i,j∈Q is stochastic and irreducible (that is, the network of queues is strongly connected). We set R = {(i, j) | pi,j > 0}.

The state space is denoted by S: it is the set of elements

x = (x1, x2, . . . , xK ) ∈ N K such that K k=1 x k = M and ∀k ∈ Q, 0 ≤ x k ≤ C k . We know that |S| ≤ K+M -1 K-1
, which grows exponentially with K.

For (i, j) ∈ R, ti,j : S → S is the function that describes the routing of a customers from queue i to queue j: ti,j(x) = x + (ejei)1 x i >0 and x j <C j where ei is the vector that has all its coefficients equal to 0 except the i-th, equal to 1 and 1A is the characteristic function of A. This network can be described by an ergodic Markov chain: at each step, a pair (i, j) is randomly chosen, with probability ν i p i,j k∈K ν k . Our objective is to sample the stationary distribution with perfect sampling techniques (as in [START_REF] Propp | Exact sampling with coupled Markov chains and applications to statistical mechanics[END_REF]). Algorithm 1 (PSS) is the direct application of the perfect sampling algorithm to this model. It produces a sample of the stationary distribution in finite expected time.

Algorithm 1: Perfect sampling using sets of state (PSS)

Data: (U-n = (i-n, j-n)) n∈N an i.i.d sequence of r.v Result: x ∈ S 1 begin 2 n ← 1; 3 t ← tU -1 ; 4 while |t(S))| = 1 do 5 n ← 2n; 6 t ← tU -1 • • • • • tU -n ; 7 return x, the unique element of {t(S)}.
This algorithm has a complexity at least linear in |S|, which grows exponentially with K. In [START_REF] Bouillard | Perfect sampling for closed queueing networks[END_REF] a new and compact representation of the state space has been proposed, for which sets of states have a strucured representation, polynomial in K and M . The main idea is to represent states as paths in a directed graph called diagram. In the next subsection, we give a summary of results in [START_REF] Bouillard | Perfect sampling for closed queueing networks[END_REF].

Diagram

Let D = (N, A) be a directed graph where N = {0, . . . , K}× {0, . . . , M } and g : S → P(N 2 ) denote the function which associates a set of arcs to each state:

g(x) = K i=1 (i -1, i-1 k=1 x k ), (i, i k=1 x k ) .
These arcs form a path from node (0, 0) to (K, M ). The

directed graph D = (N, A) is called diagram if there exists S ⊆ S such that A = g(S) := s∈S g(s).
The complete diagram is D = (N, A) where A = g(S), the image of the state space. An example of such a diagram is given in Figure 1

. Note that |A| ≤ |A| ≤ K(M +2)(M +1) 2 . 0, 0 1, 0 2, 0 4, 0 3, 0 1, 1 2, 1 1, 2 3, 1 4, 1 2, 2 3, 2 4, 2 1, 3 2, 3 3, 3 4, 3 5, 4 0 0 0 1 2 Figure 1: Complete diagram with K = 5, M = 3 and C = (2, 1, 3, 1, 2)
. The set of bold (and red) arcs is the image of state x = (0, 0, 2, 0, 1). We can now define the transformation Ti,j corresponding to the service of a customer from queue i to queue j:

Let a = (k -1, m), (k, m ′ ) ∈ A. The value of a is v(a) = m ′ -m
for (i, j) ∈ R, Ti,j(D) = φ • ti,j • ψ(D).
The transformation Ti,j can be performed directly without having to compute ψ, ti,j and φ. It will be shown in Section 3 that Ti,j(D) can be computed in time polynomial in O(KM 2 ) (while computing ti,j for all states, as in PSS, has exponential complexity).

Theorem 1. PSD algortihm terminates in finite expected time and produces an exact sample from the stationary distribution.

Gap-free diagram

In this subsection, we present a sub-class of diagrams that can be represented even more compactly.

A diagram D = (N, A) is called gap-free if ∀k ∈ Q, ∀s ≤ C k , ∀m1 < m2 < m3 ≤ M , ∀m ≤ M , ∀s1 < s2 < s3 ≤ C k , • (k -1, m1), (k, m1 + s) , (k -1, m3), (k, m3 + s) ∈ A ⇒ (k -1, m2), (k, m2 + s) ∈ A; • (k -1, m), (k, m + s1) , (k -1, m), (k, m + s3) ∈ A ⇒ (k -1, m), (k, m + s2) ∈ A; • (k -1, m -s3), (k, m) , (k -1, m -s1), (k, m) ∈ A ⇒ (k -1, m -s2), (k, m) ∈ A.
The second and third implications state that the values of arcs from (to) a node are contiguous. It can be easily seen that the complete diagram is gap-free and that gap-free diagrams have a more compact representation: it suffices to know for each value s and each queue k what are the minimal and maximal arc (from the first implication). Consequently, it can be represented by less than 2K(M + 1) arcs.

In most cases, if D is gap-free, Ti,j(D) is also gap-free. But it is not true in general, as shown on Figure 8. For a diagram D, we denote D gf the smallest gap-free diagram (for arc inclusion) that contains D. This construction is illustrated on Figure 2. From this transformation, we can define a new transition function on gap-free diagrams: for (i, j) ∈ R, The perfect sampling algorithm can be adapted to gap-free diagrams. It suffices to replace every occurrences of T in Algorithm 2 by T gf . This is PSF algortihm.

T gf i,j (D) = [Ti,j(D)] gf . D D gf
Algorithm 2: Perfect sampling using diagrams (PSD)

Data: (U-n = (i-n, j-n)) n∈N an i.i.d sequence of r.v Result: x ∈ S 1 begin 2 n ← 1; 3 T ← TU -1 ; 4 while |ψ(T (D))| = 1 do 5 n ← 2n; 6 T ← TU -1 • • • • • TU -n ; 7 return x, the unique element of {ψ(T (D))}.
Theorem 2. PSF algortihm terminates in finite expected time and produces a sample distributed according to the stationary distribution.

Proof sketch. The algorithm terminates in finite expected time for the same reason as Algorithm 2: the cou-pling sequence exhibited in [START_REF] Bouillard | Perfect sampling for closed queueing networks[END_REF] produces gap-free diagrams only.

Let (Un)

N n=1 ∈ R N . As TU 1 • • • • • TU N (D) ⊆ T gf U 1 • • • • • T gf U N (D)
, and as diagrams representing one state are gapfree, we are ensured that the result of the algorithm is distributed as the stationary distribution.

TOOLBOX IMPLEMENTATION

This section is devoted to the implementation of the transition function of diagrams. As we will see, it is convenient to represent diagrams with matrices. Then a natural choice for implementation is Matlab, as it handles matrix computations (specially when they are sparse) very well. The toolbox Clones is a Matlab toolbox that implements perfect sampling algorithms using different state space representations: state space, diagrams and gap-free diagrams. It is available at http://www.di.ens.fr/∼rovetta/Clones. At the same address there is also a detailed documentation on how to use the Toolbox. In this section, we first introduce the matrix representation of diagrams, then present the transition algorithms and finally describe some additional functionalities of the toolbox.

Matrix representations

Representation of diagrams

Representation of gap-free diagrams

Gap-free diagrams can be more efficiently implemented. Because of the implication ∀k ∈ Q, ∀s ≤ C k , ∀m1 < m2 < m3 ≤ M , (k-1, m1), (k, m1+s) , (k-1, m3), (k, m3+s) ∈ A ⇒ (k -1, m2), (k, m2 + s) ∈ A, each column can be represented by a 2 × (C k + 1) matrix Fk:

• Fk(1, s) = min(m, (k -1, m), (k, m + s) ∈ A) • Fk(2, s) = max(m, (k -1, m), (k, m + s) ∈ A) • with the convention that Fk(1, s) = Fk(2, s) = -1 if {m, (k -1, m), (k, m + s) ∈ A} = ∅.

Diagram transition algorithms

In this section, we present the main algorithm of the Clones toolbox: the transition on diagrams. First, we need to find some classes of paths in the diagram: those that will remain unchanged by the transition and those that will evolve. In the whole section we denote by D the diagram matrix representation of diagram D = (N, A) and describe the transition Ti,j(D).

Finding paths in a diagram

To compute Ti,j(D), we need to know what are the paths of D that will remain unchanged because queue i is empty or queue j is full, and what are the paths that will be modified. These paths can be determine only by the knowledge of their arc value at the i-th and j-th columns. So, we need to define matrices D V i→j that represents the paths of D from columns i to j whose arcs at column i have value in V . It is constructed as follows.

Let V ⊆ {0, 1, . . . , Ci} be a set of values, D V i is the matrix that represents arcs in the column i with values in V :

D V i (m, m ′ ) = Di(m, m ′ )1 m ′ -m∈V .
Paths from column i to column j generated by arcs with value in V in column i are represented by the matrix

D V i→j = [Y1, Y2, . . . , Y |i-j|+1 ]. The matrices Y k are computed as: Case Y1 Y k Y |i-j|+1 i < j D V i diag(1M+1Y k-1 )Dk i > j Dkdiag(Y k+1 1 t M +1 ) D V j where ∀v ∈ {0, 1} N , diag(v)(m, m ′ ) = v(m)1 m=m ′ .
Let 1M be the identity matrix of size M . Note that as diag(v) is a diagonal Boolean matrix, computing diag(1M+1Y k-1 )Dk does not require matrix multiplication. It boils down to copying the lines of Dk corresponding to non-null elements of the diagonal of diag(1M+1Y k-1 ). The most complex operation is to compute 1M+1Y k-1 , which can be done in O(M 2 ). As |j -i| + 1 ≤ K, the complexity of finding those path is in O(KM 2 ).

Transition algorithm

We are now ready to present the transition algorithm. Let D be a diagram and (i, j) ∈ R. We explain how to compute Ti,j(D). First let us explain the transformation on diagram D (see Figure 8). We begin by computing three sets of arcs (not necessarily disjoint), depending on whether they will be modified by the transition:

• Empty is the set of arcs that belong to paths with arc value 0 at column i. Those paths correspond to states where the queue i is empty.

• Full is the set of arcs that belong to paths with arc value Cj at column j. Those paths correspond to states where the queue j is full.

• Transit is the set of arcs that belong to paths with value < Cj at column j and > 0 at column i. Those paths correspond to states that are modified by ti,j.

Second, arcs in Transit are modified in accordance to ti,j:

• In column i, value of arcs are decreased by one (keeping the same source).

• In column k ∈ {min(i, j) + 1, . . . , max(i, j) -1}, arcs are shifted up or down depending on the sign of ij.

• In column j, slopes of arcs are increases by one (keeping the same destination).

If we call Transit ′ this new set, Ti,j(D) = (N, Empty ∪Full∪ Transit ′ ). The transformation using matrix operations is given in Algorithm 3, where the Boolean operators and and or are component-wise and where

+1 -1 D = (N, A) T4,2(D) T ransit T ransit ′ Empty F ull
[←X] m,m ′ = X m,m ′ +1 1 m ′ <M , [→X] m,m ′ = X m,m ′ -1 1 m ′ >0 , [ ↑ X] m,m ′ = X m+1,m ′ 1m<M , [ ↓ X] m,m ′ = X m-1,m ′ 1m>0, [ տ X] = [←[ ↑ X]], [ ց X] = [→[ ↓ X]].
The algorithm is performed on |i-j|+1 matrices. As |i-j|+ 1 ≤ K, and all the matrix operations can be done in O(M 2 ), Algorithm 3 requires O(KM 2 ) elementary operations. It is implemented by fonction Clones T(D, C, [i, j]) in the toolbox Clones.

To compute transitions on gap-free diagrams, the steps of the algorithm are the same. The only difference is that there is no need to compute all the arcs in Empty, Full and Transit, but only the extremal ones in each column, which ;

5 if i < j then 6 D ′ i Transit ← [←D Transit i ]; 7 for k from i + 1 to j -1 do D ′ k Transit ← [ տ D Transit k ] D ′ j Transit ← [ ↑ D Transit j ]; 8 else 9 D ′ i Transit ← [ ↓ D Transit i ]; 10 for k from i + 1 to j -1 do D ′ k Transit ← [ ց D Transit k ] D ′ j Transit ← [→D Transit j ]; 11 D ′ i,j ← D Empty i,j or D Full j,i or D Transit i,j ; 12 return D ′ = [D 1,min(i,j)-1 , D ′ i,j , D max(i,j)+1,K ].
can be computed directly from matrix F. As a consequence, the transformation can be done in O(KM ). This function is implemented in Clones gfT(D, C, M, [i, j]).

Additional functionalities of Clones

The Object-oriented programming can also be used with Clones. Three classes have been defined to handle transitions and perfect sampling using respectively sets of states, diagrams and gap-free diagrams. Let R be the routing matrix. The functions for perfect sampling with sets of states, diagram and gap-free diagrams are respectively Clones PSS(C, M, R, N), Clones PSD(C, M, R, N) and Clones PSF(C, M, R, N). They produce N samples.

PERFECT SAMPLING WITH Clones

In this section we are interested in the running time of the perfect sampling algorithms. Experiments have been performed on a desktop machine.

We chose a small number of queues (K = 5) in order to compare running time of Algorithms PSD and PSF against PSS (that has an exponential complexity). Each queue has a capacity M/2 and the routing matrix is fixed and randomly chosen for the whole experimentation. For each sample, the same sequence of random variables (Un)n is used for the three algorithms. Thus they produce the same sample. For each experiment, we produce 100 samples and are interested in the mean running time. Figure 11 compares the perfect sampling algorithms using respectively diagrams and gap-free diagrams when the number of customers is M = {5, 10, . . . , 200}. We can see that for values from M = 5 to 60, the algorithm with diagrams is slightly faster than the one in gap-free diagrams, whereas for higher values of M , using gap-free diagrams is much more efficient than using diagrams. The latter result was expected, as the algorithmic complexity of one step transition is in O(KM ) for gap-free diagrams instead of O(KM 2 ) for diagrams. The running time when M ≤ 60 is more surprising. We believe that the main reason is that Matlab is optimized for matrix operations. Specially, when it comes to sparse matrices. On the other hand, our implementation of gap-free diagrams cannot benefit from this optimization, and there is room for improvement.

CONCLUSION

In this paper, we presented the toolbox Clones, that performs efficient coupling from the past of closed queueing networks. This toolbox enables the comparison between three different algorithms (CFTP with sets of states, diagrams and gap-free diagrams). Experimentally, PSD and PSF algorithms are very efficient compared to PSS algorithm. For small values of M , the CFTP with diagram implementation terminates faster than with gap-free diagrams. Future work will focus on the improvement of the gap-free one-step transition implementation.

Further investigation also include the study on the coupling Experimentally, the difference between the coupling times of the three implementations is rather small, but obtaining a theoretical bound would prove the efficiency of our approach. Last, we will also explore possible extensions to other classes of closed queueing networks.

  and represents the number of customers in queue k. Graphically, the value of an arc corresponds to its slope. Note that 0 ≤ v(a) ≤ C k . Function g can be used to define the transformation from a set of states to a diagram and conversely: for a diagram D = (N, A) and a set S ⊆ S, we have φ(S) = (N, g(S)) and ψ(D) = g -1 (A).

Figure 2 :

 2 Figure 2: Transformation of a diagram with a gap into a gap-free diagram. Left: D has a gap (no arc ((1, 1), (2, 2)); Right: D gf is constructed by adding this arc to D.

  Let D = (N, A) be a diagram and k ∈ Q. In the toolbox, diagrams are implemented with incidence matrices. More precisely, Dk is a Boolean matrix of size (M + 1) × (M + 1) representing the k-th column of D: Dk(m, m ′ ) = 1 ((k-1,m),(k,m ′ ))∈A . Because a ∈ A implies that 0 ≤ v(a) = mm ′ ≤ C k , matrices Dk are Boolean upper-triangular and Dk(m, m ′ ) = 0 if m ′m > C k . Figure 3 gives an example.

Figure 3 :

 3 Figure 3: Matrix representation D2 of the second column of D.

Figure 4 :

 4 Figure 4: Diagram representation of the complete diagram in Figure 1.

Figure 5 :

 5 Figure 5: Gap-free representation of a column. The gap-free representation of the diagram of Figure 1 is given in Figure 5. Using the notations defined for D, a gap-free diagram is represented by matrix F = F1,K of size 2 × ( K k=1 C k + K). Figure 6 shows F, the gap-free representation of the complete diagram of Figure 1. It is a 2 × 14-matrix.

Figure 6 :

 6 Figure 6: Gap-free representation of the complete diagram in Figure 1.

  Figure 7 depicts D {2} 3→5 on the diagram.

Figure 7 :

 7 Figure 7: Diagram D and D {2} 3→5 in bold (red) paths. It corresponds to paths from column 3 to column 5 generated by arcs of values 2.

Figure 8 :

 8 Figure 8: Steps for computing T4,2(D). Diagram D is obtained from the complete diagram by some transition; D = T4,5 • T4,5 • T2,1(D).

Algorithm 3 :

 3 Transition algorithm Data: D the matrix representation of a diagram D, (i, j) ∈ R Result: D ′ the matrix representation of Ti,j(D).

  conversion of a diagram D into a set of states, ψ(D), is performed by Clones Psi(D). Function Clones CardPsi(D) returns the cardinal of ψ(D). Function Clones StateSpace(C, M) uses Clones Psi(D) to compute the entire state space. Diagrams can be plotted with Clones plotD(D, C, mode); mode allows to write parameters or |ψ(D)| in the graphic (Figure 9). It is also possible to have diagram animations.

Figure 9 :

 9 Figure 9: Plot of a complete diagram.

Figure 10

 10 Figure10compares the running times for the three algorithms with M ∈ {2, 4, . . . , 60}. It can be observed that the running time of Algorithm 1 grows exponentially fast, and the two other algorithms (with diagrams and gap-free diagrams) terminate much faster.

  Figure 11: Mean running time -100 simulations

0, 0 1,0 2, 0 3, 0 4, 0 1, 1 1, 2 1, 3 2, 1 2, 2 2, 3 3, 1 3, 2 3, 3 4, 1 4, 2 4, 3 5, 3 0, 0 1, 0 2, 0 3, 0 4, 0 1, 1 1, 2 1, 3 2, 1 2, 2 2, 3 3, 1 3, 2 3, 3 4, 1 4, 2 4, 3 5, 3

1, 0 2, 0 3, 0 4, 01, 1 1, 2 1, 3 2, 1 2, 2 2, 3 3, 1 3, 2 3, 3 4, 1 4, 2 4, 31, 0 2, 0 3, 0 4, 01, 1 1, 2 1, 3 2, 1 2, 2 2, 3 3, 1 3, 2 3, 3 4, 1 4, 2 4, 3 1, 0 2, 0 3, 0 4, 0 1, 1 1, 2 1, 3 2, 1 2, 2 2, 3 3, 1 3, 2 3, 3 4, 1 4, 2 4, 3 1, 0 2, 0 3, 0 4, 0 1, 1 1, 2 1, 3 2, 1 2, 2 2, 3 3, 1 3, 2 3, 3 4, 1 4, 2 4, 3

ACKNOWLEDGMENTS

This research is supported by the French National Research Agency grant ANR-12-MONU-0019.


V Ilakovac

S Boulfaat

Y Joly

A M Flank

P Lagarde

M Žitnik

M Fialin

H Berger

L Forró

Role of sulfur in BaVS 3 probed by S K-edge absorption spectroscopy

Keywords: Metal-insulator transitions and other electronic transitions, X-ray absorption spectroscopy, Charge-density-wave system

We show that the quasi-1D behavior of BaVS 3 can be understood analyzing the X-ray absorption near edge spectra at the sulfur K edge. Linear dichroïsm experiments, analyzed with the help of ab initio calculations, reveal two strong and polarization dependent pre-edge features, induced by the band character of the 3d vanadium levels. They are related to crystal field split t 2g and e g states. When the temperature is lowered, the t 2g feature shifts progressively to higher energy, and its intensity increases for the polarization along the c-axis, stacking direction of the V-S face sharing octahedra. This behavior points to the depletion of sulfur states and thus the lack of S 3p -V 3d hybridization in the direction of V-S chains.

Introduction

Understanding the interaction of itinerant and localized electrons is one of main challenges in solid state sciences. It determines properties of a large class of correlated systems with very rich phase diagrams, involving a subtle interplay between charge, orbital, and spin degrees of freedom [1]. The complex electronic structure of BaVS 3 , in spite of its apparently simple quasi-1D geometry, ranks it among model compounds for the study of the competition of different ordering parameters.

BaVS 3 structure is constituted of chains of V-S face sharing octahedra running in c-direction and separated by Ba atoms [2]. The distance between neighboring V atoms along the V-S chains is less than two times * Corresponding author: vita.ilakovac-casses@upmc.fr shorter then the one between the chains. It could thus be considered as a quasi-1D compound. Nevertheless, its conductivity anisotropy is extremely low compared to usual quasi-1D compounds [3]. It undergoes three low temperature second order transitions. In stoichiometric BaVS 3 , they occur at 240 K, 69 K and 30 K. The first is related to a structural change, a zigzag deformation of the V-S chains, reducing the crystal symmetry from hexagonal to orthorhombic, without altering its paramagnetic and conducting properties [2]. The metalinsulator (MI) transition at T MI = 69 K is accompanied by the doubling of the unit cell in the c direction [4]. The pre-transitional diffuse lines, observed at the wavevector q MI = 0.5 c * in x-ray diagrams, reveal the formation of the charge density waves (CDW) in V-S chains already at 170 K [5]. It is surprising that only half of active electrons directly participate to the 2k F Peierls instability. Spin degrees of freedom are affected at the MI transition [3], but freeze only below the Néel temperature, T N = 30 K, establishing an incommensurate antiferromagnetic (AF) order in the (a, b) plane, as determined by the inelastic neutron scattering measurements on powder samples [6]. The third component of the AF incommensurability is brought to light by single crystal resonant X-ray diffraction at the V L 3 edge, refining the ordering wave vector q AF = (0.226, 0.226, 0.016) [7,8]. Moreover, these measurements reveal that the spins are polarized along the a direction and ordered alternatively in the b direction of the lowest temperature (monoclinic) unit cell. A conical spin arrangement could explain the incommensurate order along V-S chains [9].

Local density approximation (LDA) calculations reveal an interplay between two different types of t 2g electron states at the Fermi level: two narrow e(t 2g ) bands, and one dispersive band with mainly d z 2 character, related to delocalized states directed along the V-S chains [10,11,12]. They are occupied by a single 3d 1 electron of the formally V 4+ ion. However, according to the LDA, the delocalized d z 2 states are nearly completely filled, which overestimates the itinerant character and fails to explain the gap opening in the insulating phase. LSDA + U calculations, including a static exchange and an antiferromagnetic ordering in the (a,b) plane, find a 0.15 eV gap but do not explain the paramagnetism in the metallic phase [13]. The experimental 1:1 filling ratio of the d z 2 /e(t 2g ) states and the gap opening in the insulating phase are correctly reproduced when the correlation and exchange effects are included via dynamical mean field theory (DMFT) [14].

Already in earliest studies, it was shown that the sulfur content is determinant for BaVS 3 properties. Massenet et al. [16] report lowering of the T S down to 150 K and ferromagnetic behavior below 16-17 K in BaVS 3-δ with δ = 0.12. More systematic study of magnetization, conductivity and NMR measurements are performed with a formal δ = 0.05, 0.10, 0.15 and 0.20 [17]. The deficiency of about (nominal) 3 % turns BaVS 3-δ to paramagnetic semiconductor already at room temperature. A short range 3D charge order (CO) at the wave vector q CO = (2/3, 0, ≈ 1/5) O exists already at 300 K. Quasi-1D Peierls instability is suppressed. Instead, a long range order develops at q CO below 110-220 K [18,19]. The work of Yamasaki et al. [17] confirms that below 15-18 K the sulfur deficient BaVS 3-δ becomes a ferromagnet.

X-ray absorption near edge spectroscopy (XANES) is an experimental tool which is able to give relevant information on both structural and electronic aspects of a studied system. The full power of this spectroscopy is used when the polarization analysis is applied on a single crystal material, as the selection rules permit to probe particular electronic states. The purpose of this paper is to show that such experiments combined with the ab-initio simulations reveal subtle details of the electronic structure of BaVS 3 in different phases of the compound. The originality of the work is in performing the measurements at the S K edge, a priory probing S 3p states, while gaining insight into relative (de)localisation of V 3d states, driven by the V 3d -S 3p overlap.

The details of the realized S K edge polarization and temperature dependent XANES experiment on BaVS 3 and accompanying calculations are presented in the section 2. They are followed by the presentation of results and their interpretations in section 3 and 4, respectively.

Experiment and calculation details

Single crystal BaVS 3 samples are grown by the tellurium flux method [20]. Strong polarization dependence of the S K edge x-ray absorption spectra at room temperature (measured in the fluorescence mode) was c and ⊥ c. The intensity of the fluorescence was measured in a total yield mode by a detector at 90 • scattering angle. The inflection point of the S edge of the pure sulfur was at 2472 eV with the energy resolution of 0.3 eV. The normalization and the self-absorption correction of the spectra are performed by a code using the procedure described in ref. [21,22].

Calculations of the room temperature spectra are performed by the Finite Difference Method Near Edge Spectroscopy (FDMNES) code [23] in the LSDA+U frame. Prior to the evaluation of the cross section, the Schrödinger equation is resolved by the finite difference method in the excited state of a 6 Å cluster embedded in a hexagonal phase BaVS 3 .

Results

Room temperature BaVS 3 S K edge XANES spectra in Fig. 1 show a very strong pre-edge doublet at 2468-2472 eV, a large structured edge feature at 2472-2484 eV and two extended structures at higher energy. The difference between the c and ⊥ c spectra is particularly strong in the pre-edge and the edge region, and becomes less important, but still existent at higher en- ergy. We keep the pre-edge terminology because these features, often present in the S or O K-edges in metallic oxides or sulfides appear at lower energy than the corresponding edge in pure compounds.

FDMNES calculations reproduce all the features in the experimental spectra (see Fig. 2). There are however small differences in relative positions and shapes. In particular, the extended structures are slightly shifted to lower energy and the edge structure at 2480 eV is not as large as the experimental one.

The pre-edge and the edge region (2468-2484 eV) is enlarged in Fig. 3. The polarization dependence is here much more evident compared to the spectra collected at the V L edge [START_REF] Ilakovac | [END_REF]. It is worth to note that V L edge xray absorption spectroscopy probes V 3d states directly via dipolar transition rules, while S K edge XANES gives only an indirect insight, through V 3d -S 3p hybridization. The polarization dependence is particularly important in the region of the second pre-edge peak (2469.2 eV), where σ e is much stronger compared to σ e ⊥ (e stands for experiment). Moreover, the energy difference between the two pre-edge features, is not the same for the two geometries: ∆ e = 1.2 eV, while ∆ e ⊥ = 1.8 eV. The difference of the two spectra (σ eσ e ⊥ ), the linear dichroïsm, is shown in the same figure.

Fig. 4 shows the simulated spectra in the same energy range. The energy difference of the pre-edge peaks is ≈ 50 % larger compared to the experiment, as ∆ c = 2.0 eV and ∆ c ⊥ = 2.4 eV (c stands for calculation). Including electron-electron interaction of 5 eV on Vand 2 eV on S-sites improves relative intensities. It induces a partial spectral weight transfer from 2472 eV to 2476 eV and modifies the pre-peak region for c. However, it does not modify the relative shift of the two pre-edge features: ∆ is not related to the correlation gap. Decreasing the effective core-hole screening to 0.5 elec- trons (transition state method with fractional occupation [START_REF] Stöhr | NEXAFS Spectroscopy[END_REF]) increases the two peaks in the pre-edge region, in particular the first peak for c. It is interesting to note that the calculated c spectra are much more dependent on the cluster size compared to the ⊥ c spectra, indicating that they probe more delocalized states.

The difference of the calculated spectra, σ cσ c ⊥ , is presented by the green line in Fig. 4. It matches well the experimental linear dichroïsm except in two points. First, at the very beginning of the pre-edge, the experiment shows no dichroïsm, while the calculated σ c is more intense than σ c ⊥ . Second, the experiment shows positive dichroïsm at 2482 eV, while the calculations predict slightly negative value. The first discrepancy can be explained by the inability of the LDA based methods to describe correctly the relative occupation of the states at the Fermi level. The second difference can be related to a slight difference of the amount of the sulfur oxide at the surface of the two samples.

The temperature dependence of the two pre-edge structures is shown in Fig. 5. For both orientations, the first pre-edge structure shifts progressively to higher energy as the temperature is lowered. From 300 K to 40 K, the total shift is ≈ 80 meV. Moreover, for c, the first pre-edge structure increases in intensity for ≈ 8% as the temperature is lowered from 300 K to 100 K and remains constant below.

Discussion

Pre-edge structures in ligand K-edge spectra of transition metal compounds are often related to mixing of metal-ligand unoccupied states close to the Fermi level.

In the case of BaVS 3 , V 3d states overlap with S 3p states.

To understand the polarization dependence, we use Cartesian coordinates (x, y, z) with z along the c crystal axis. At the S K edge, quadrupolar transitions are negligible. Thus, when the polarization is along x, y, z, we probe respectively unoccupied 3p x , 3p y , 3p z sulfur states. Moreover, in the hexagonal phase of BaVS 3 , the 3-fold axis makes that in the dipole approximation, the absorption signal is isotrope when the polarization is parallel to the basal (x,y) plane, so we can distinguish only z and ⊥ z orientations. Fig. 6 illustrates the relation between S 3p and V 3d orbitals in V-S octahedra. The point group of vanadium is D 3d . Two of its 3d orbitals corresponding to e g representation have lobes pointing to sulfur atoms are shifted to higher energy by the crystal field splitting (upper panel). Other two, e(t 2g ), point to the faces of the V-S octahedra, similarly to d z 2 states, and are shifted to lower energy (lower panel).

The S 3p and V 3d projected density of states (DOS) calculated by FDMNES is shown in Fig. 7. Similarly to the DOS calculations published previously [10,13,15] S 3p states are strong in the occupied part, while in the unoccupied part V 3d predominates. V 3d unoccupied DOS is increased at E-E F ≈ 0 eV, 1 eV, and 1.8-2.8 eV. The first increase corresponds to d z 2 and e(t 2g ) states at the Fermi level, the second to the d z 2 free-electron parabola and the third to e g states. S 3p DOS is increased at energies where V 3d states have important density showing that their hybridization is strong. In e g region, the S 3p z DOS is important at 1.8 eV, while S 3p x,y has its maximum at 2.7 eV. This can explain the difference in the pre-edge energy splitting, ∆ e (∆ e ⊥ ) of 1.2 eV (1.8 eV), respectively, even if the calculated values are overestimated. The crystal field splitting is differently probed by z and ⊥ z polarization. On the other hand, the low S 3p x,y DOS at E-E F ≈ 1 eV indicates that d z 2 mixes mostly with the S 3p z states. The S 3p z -V 3d z 2 hybridization is thus the main factor of the quasi-1D behavior of BaVS 3 .

Contrary to the polarization dependence which is very strong, the temperature variation of the spectra in Fig. 5 is relatively small. For both geometries, the first pre-edge feature shifts to ≈ 80 meV higher energy, when the temperature is lowered from 300 K down to 40 K. This value is large compared to the 42 -70 meV charge gap opening in V 3d states below T MI [3,26,27,28,29,30]. The shift is progressive, as for 200 K (not shown), its value is in between these measured at 300 K and 100 K. Indeed, at T S = 240 K, BaVS 3 exhibits a hexagonal to orthorhombic second-Figure 6: Schematic presentation of V 3d and S 3p electronic states in the room temperature phase. Five V 3d states, two e g (upper panel), d z 2 and two e(t 2g ) (lower panel), are centered on V atoms. S 3p x,y (yellow lobes) and S 3p z (orange lobes) are centered on one of six S atoms at the edge of V-S octahedra, stacked along the z-direction. Below T S , V atoms are displaced out of the center of the octahedra, alternatively in positive (negative) x-direction. Contrary to the hexagonal phase, were all S atoms are equivalent by symmetry, in the orthorhombic phase, there are two types of sulfur atoms, S 1 , in the (x,z) plane (orange points), and two S 2 (yellow points), per octahedra. order transition. The temperature evolution of cell parameters indicates that the orthorhombic distortion progressively increases when the temperature is lowered [31]. It does not affect remarkably the conducting properties immediately at T S [3], but it induces a degeneracy break of the two e(t 2g ) states. One of them shifts to 100 meV higher energy (at 100 K), while the other hybridizes strongly with d z 2 states [15]. The depletion of V 3d states with strong DOS the Fermi level, reflects in the S K edge XANES as a high energy shift of the first pre-edge peak.

The first pre-edge structure as well increases in intensity for c, when the temperature is lowered. This effect is related to the decrease of the S 3p z -V 3d z 2 hybridization, induced by the zig-zag deformation of vanadium chains. Together with the formation of charge density waves in chains in the pretransitional range (2k F fluctuations below 170 K), it can induce semiconducting behavior already at 130 K, the temperature which is far above the 69 K MI transition [32].

Conclusion

BaVS 3 has attracted strong scientific interest since more than forty years, first for its mixed 1D and 3D behavior, and now as a model compound for the study of the localized/itinerant states' interplay. Its properties still hide mysteries, in particular the question of the semiconducting behavior onset above the metalinsulator transition, and the extreme dependence of its properties on the sulfur content. This work combines S K edge x-ray absorption near edge spectroscopy (XANES) and ab initio calculations in order to provide original data on the relation between sulfur and vanadium electronic states inside V-S chains. S K-edge XANES measurements reveal two strong and polarization dependent pre-edge structures, related to S 3p and V 3d crystal field split, t 2g and e g , states. The first feature shifts progressively to ≈ 80 meV higher energy when the temperature is lowered from 300 K (metallic state) to 40 K (insulating phase). Moreover, when the incoming light polarization is parallel to the c-crystal axis, its intensity increases. This behavior is interpreted by the progressive depletion of S 3p z -V 3d z 2 hybridization due to the zig-zag deformation of V-S octahedra chains.

Figure 1 :

 1 Figure 1: S K edge XANES of BaVS 3 measured for c (σ e , orange cercles) and ⊥ c (σ e ⊥ , black points). Superscript e stands for experiment.

Figure 2 :

 2 Figure 2: Calculated S K edge cross section of BaVS 3 in Mbarn per unit cell for c (σ c , orange line) and ⊥ c (σ c ⊥ , black dotted line). Superscript c stands for calculation.

Figure 3 :

 3 Figure 3: Polarization dependence in the near S K edge region. The difference of the intensities, σ eσ e ⊥ , the linear dichroïsm, is shown by green squares.

Figure 4 :

 4 Figure 4: Calculated near S K edge spectra for c (orange line) and ⊥ c (black dotted line). Linear dichroïsm, σ cσ c ⊥ , is presented by the green dotted line.

Figure 5 :

 5 Figure 5: Temperature dependence of the two pre-edge features for the two orientations, c and ⊥ c, shifted vertically for clarity.

Figure 7 :

 7 Figure 7: S 3p and V 3d projected DOS of BaVS 3 . Sulfur DOS is multiplied by five in order to be more evident in the unoccupied part.

[32] Recent structural refinements, reported by A. Arakcheeva et al. in this Proceedings, point to a new structural transition at 130 K which can be responsible for the semicondutcing behavior below this temperature.

Acknowledgments

Discussions with N. Trcera, Ph. Sainctavit, R. Nemeusat and D. Cabaret are gratefully acknowledged.