ST Christophe Besse

Stephane Descombes

Guillaume Dujardin

Ingrid Lacroix-Violet

Éphane Descombes

Energy preserving methods for nonlinear Schrödinger equations

Keywords: AMS Classification: 35Q41, 81Q05, 65M70 Keywords. Nonlinear Schrödinger equation, Gross-Pitaevskii equation, numerical methods, relaxation methods

teaching and research institutions in France or abroad, or from public or private research centers.

Introduction

The nonlinear Schrödinger equation (NLSE) is a fairly general dispersive partial differential equation arising in many areas of physics and chemistry [START_REF] Abdullaev | Optical Solitons[END_REF][START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF][START_REF] Dauxois | Physics of Solitons[END_REF][START_REF]Lev Pítajevskíj and Sandro Stringari[END_REF][START_REF] Catherine | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF]. One of the most important application of the NLSE is for laser beam propagation in nonlinear and/or quantum optics and there it is also known as parabolic/paraxial approximation of the Helmholtz or time-independent Maxwell equations [START_REF] Abdullaev | Optical Solitons[END_REF][START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF][START_REF] Dauxois | Physics of Solitons[END_REF][START_REF] Catherine | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF]. In the context of the modeling of the Bose-Einstein condensation (BEC), the nonlinear Schrödinger equation is known as the Gross-Pitaevskii equation (GPE), which is a widespread model that describes the averaged dynamics of the condensate [START_REF]Lev Pítajevskíj and Sandro Stringari[END_REF][START_REF] Antoine | Modeling and computation of Bose-Einstein condensates: stationary states, nucleation, dynamics, stochasticity[END_REF]. Depending on the physical situation that one considers, several terms in the right-hand side of the NLSE appear. Our goal is to develop numerical methods for the time integration of a fairly general NLSE including realistic physical situations, that have high-order in time and have good qualitative properties (preservation of mass, energy, etc) over finite times. We develop our analysis in spatial dimension d P t1, 2, 3u because it fits the physical framework, even if most methods and results naturally extend to higher dimensions. The space variable x will sometimes lie in R d , sometimes on the d-dimensional torus T d δ " pR{pδZqq d (for some δ ą 0). In this paper, we consider a NLSE of the form [START_REF] Abdullaev | Optical Solitons[END_REF] iB t ϕpt, xq "

ˆ´1 2 ∆ `V pxq `β|ϕ| 2σ pt, xq `λ `U ˚|ϕpt, ¨q| 2 ˘pxq ´Ω.R ˙ϕpt, xq,
where ϕ is an unknown function from R ˆRd or R ˆTd δ to C, ∆ is the Laplace operator, V is some real-valued potential function, β P R is a parameter that measures the local nonlinearity strength, λ P R is a parameter that measures the nonlocal nonlinearity strength with convolution kernel U , Ω P R d is a vector encoding the direction and the speed of a rotation, and R is a rotation operator that is important in the modeling of rotating BEC (for example, R " x ^p´i∇q when d " 3). The NLSE is supplemented with an initial datum ϕ in . The results presented in this paper extend to more general power law nonlinearities such as

iB t ϕpt, xq " ˜´1 2 ∆ `V pxq `K ÿ k"1
β k |ϕ| 2σ k pt, xq `λ `U ˚|ϕpt, ¨q| 2 ˘pxq ´Ω.R ¸ϕpt, xq, but we restrict ourselves to K " 1 for the sake of simplicity. Equation ( 1) is hamiltonian for the energy functional

(2)

Epϕq " ż R d ˆ1 4 }∇ϕ} 2 `1 2 V |ϕ| 2 `β 2σ `2 |ϕ| 2σ`2 `λ 4 pU ˚|ϕ| 2 q|ϕ| 2 ´Ω 2 ϕRϕ ˙dx,
provided U is a real-valued convolution kernel, symmetric with respect to the origin 1 . In practice the convolution kernel U may for example correspond to a Poisson equation (U pxq " 1{p4π|x|q in dimension d " 3) or it may represent dipole-dipole interactions (see [START_REF] Bao | A simple and efficient numerical method for computing the dynamics of rotating Bose-Einstein condensates via rotating Lagrangian coordinates[END_REF]).

The main goal of this paper is the analysis of numerical methods for the time integration of (1) that preserve the energy [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF] or a discretized analogue of it. In particular, we are interested in the order (in time) of such methods. A well known method is the Crank-Nicolson method introduced in [START_REF] Crank | A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type[END_REF] for parabolic problems (see for example a posteriori error estimates in [START_REF] Akrivis | A posteriori error estimates for the Crank-Nicolson method for parabolic equations[END_REF]) and applied in [START_REF] Delfour | Finite-difference solutions of a nonlinear Schrödinger equation[END_REF] to Schrödinger equations. For all nonlinearities, these methods are fully implicit. However, they have second order in time (see [START_REF] Sanz-Serna | Methods for the numerical solution of the nonlinear Schroedinger equation[END_REF] for the case of the cubic NLS equation and [START_REF] Wang | Crank-Nicolson difference scheme for the coupled nonlinear Schrödinger equations with the Riesz space fractional derivative[END_REF] for the case of a system with possibly fractional derivatives) and preserve discrete analogues of the energy (2) as well as the total mass (squared L 2 -norm) of the solutions. Unfortunately, since they are fully implicit, these methods are costly. To work around this problem, the methods introduced and analyzed in this paper belong to the family of relaxation methods.

Relaxation methods for Schrödinger equations were introduced in [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF][START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF]. They have been applied to different NLS equation for example in the context of plasma physics [START_REF] Oelz | Analysis of a relaxation scheme for a nonlinear Schrödinger equation occurring in plasma physics[END_REF]. For cubic nonlinearities (σ " 1 in (1)), they are linearly implicit hence very popular [START_REF] Antoine | Numerical solution of time-dependent nonlinear schrdinger equations using domain truncation techniques coupled with relaxation scheme[END_REF][START_REF] Dahlby | Plane wave stability of some conservative schemes for the cubic Schrödinger equation[END_REF][START_REF] Antoine | Modeling and computation of Bose-Einstein condensates: stationary states, nucleation, dynamics, stochasticity[END_REF][START_REF] Gazeau | Probability and pathwise order of convergence of a semidiscrete scheme for the stochastic Manakov equation[END_REF][START_REF] Henning | Numerical comparison of mass-conservative schemes for the grosspitaevskii equation[END_REF]. They preserve the L 2 -norm and a discrete analogue of (2). It is well-known that they have numerical order 2 but up to our knowledge there is no proof of order 2 in the literature. This paper presents two new results with respect to relaxation methods for [START_REF] Abdullaev | Optical Solitons[END_REF]. First, we prove rigorously that the classical relaxation method, applied to the classical cubic NLS equation (i.e. (1) with V " 0, β " 1, σ " 1, λ " 0 and Ω " 0) is of order 2. Second, we present a generalized relaxation method that allows to deal with general power law nonlinearities (σ ‰ 1 in (1)) and with the full GPE. The generalized relaxation methods that we introduce in this context are implicit (actually explicit for σ ď 4), have numerical order 2 and we show that they preserve an energy which is also a discretized analogue of (2). This paper is organized as follows. In Section 2 we recall the definition of the Crank-Nicolson method and give a short explanation of its energy preserving property. Section 3 is devoted to relaxation methods applied to (1): In a first part we recall the method introduced in [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF] for the cubic Schrödinger equation and in a second part we give a proof of the optimal order of convergence for an initial datum belonging to H s`4 pR d q, d in t1, 2, 3u and s ą d{2. In section 4 we propose a generalized relaxation method that allows to deal with general nonlinearites and we prove that this method is also an energy preserving method. Section 5 deals with numerical results in different physical models showing the efficiency of the methods.

Preservation of energy and Crank-Nicolson scheme

A usual way to prove the conservation of the energy (2) consists in multiplying the equation (1) by B t ϕpt, xq, where z denotes the conjugate value of a complex z, integrating over R d and taking the real part of the result. This computation relies on the identity which holds for all smooth functions ϕ of time with values into a space of sufficiently integrable functions

(3) Re ż R d |ϕ| 2σ ϕB t ϕ dx " 1 2σ `2 d dt ż R d |ϕ| 2σ`2 dx, @σ ě 0.
A possible way to derive numerical schemes that preserve an energy functional is therefore to mimic this identity at the discrete level.

In 1981, Delfour, Fortin and Payre [START_REF] Delfour | Finite-difference solutions of a nonlinear Schrödinger equation[END_REF], following an idea of Strauss and Vasquez [START_REF] Strauss | Numerical solution of a nonlinear Klein-Gordon equation[END_REF], proposed a way to deal with the nonlinear term |ϕ| 2σ ϕ for the Crank-Nicolson scheme. This method generalizes the second order mid-point scheme for the linear Schrödinger equation

i ϕ n`1 ´ϕn δt " ´1 2 ∆ ϕ n`1 `ϕn 2 ,
where ϕ n pxq denotes an approximation of ϕpt n , xq with the discrete time t n " nδt defined with the time step δt. Their approach can be explained as follows. If one looks for a real-valued function g : C 2 Ñ R such that the scheme takes the form

(4) i ϕ n`1 ´ϕn δt " ˆ´1 2 ∆ `V `βgpϕ n , ϕ n`1 q `λ ˆU ˚ˆ|ϕ n`1 | 2 `|ϕ n | 2 2 ˙˙´Ω.R ˙ϕn`1 `ϕn 2 ,
then, multiplying this relation by iϕ n`1 ´ϕn , integrating overs R d and taking the real part, as we did in the time-continuous setting above, yields to 0 in the left-hand side and several terms in the right-hand side. Amongst these terms, those involving g are equal to

β ż R d gpϕ n , ϕ n`1 q `|ϕ n`1 | 2 ´|ϕ n | 2 ˘dx, since g is real-valued. Let us denote by G the function v Þ Ñ |v| 2σ`2 {p2σ `2q.
A sufficient condition for the method (4) to preserve an energy of the form ( 2) is therefore to have

gpϕ n , ϕ n`1 q `|ϕ n`1 | 2 ´|ϕ n | 2 ˘" Gpϕ n`1 q ´Gpϕ n q.
This is exactly the definition of g chosen in [START_REF] Delfour | Finite-difference solutions of a nonlinear Schrödinger equation[END_REF].

In the following, the Crank-Nicolson method for the GPE (1) is therefore defined using the formula

i ϕ n`1 ´ϕn δt (5) " ˆ´1 2 ∆ `V `β σ `1 |ϕ n`1 | 2σ`2 ´|ϕ n | 2σ`2 |ϕ n`1 | 2 ´|ϕ n | 2 `λ ˆU ˚ˆ|ϕ n`1 | 2 `|ϕ n | 2 2 ˙˙´Ω.R ˙ϕn`1 `ϕn 2 .
We shall use the notation ϕ n`1 " Φ CN δt pϕ n q, for the Crank-Nicolson method (5). In the expression above, the term corresponding to the nonlinearity should be understood as

(6) β σ `1 |ϕ n`1 | 2σ`2 ´|ϕ n | 2σ`2 |ϕ n`1 | 2 ´|ϕ n | 2 " β σ `1 σ ÿ k"0 |ϕ n`1 | 2k |ϕ n | 2pσ´kq ,
so that it is indeed non-singular and it is consistent with the non-linear term β|ϕ| 2σ . The Crank-Nicolson method is fully implicit. It is known to have order two for the cubic NLS equation [START_REF] Sanz-Serna | Methods for the numerical solution of the nonlinear Schroedinger equation[END_REF]. Moreover it preserves exactly the L 2 -norm of the solution as well as the following energy:

(7) E CN pϕq " Epϕq,
with E defined by [START_REF] Ablowitz | Solitons and the Inverse Scattering Transform[END_REF].

In Section 4, we shall use similar ideas to derive energy-preserving relaxation methods for general power laws nonlinearities. Before doing so, we first deal with the classical relaxation method in Section 3.

The classical relaxation method

3.1. An energy preserving method. In [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF], Besse introduced the usual relaxation method [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF] applied to the nonlinear Schrödinger equation (1) with V " 0, λ " 0, Ω " 0 and σ " 1 that is known as the cubic nonlinear Schrödigner equation [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] iB The relaxation method then consists in discretizing both equations respectively at discrete times t n and t n`1{2 and to solve iteratively

(10) $ ' & ' % Υ n`1{2 `Υn´1{2 2 " |ϕ n | 2 , i ϕ n`1 ´ϕn δt " ˆ´1 2 ∆ `βΥ n`1{2 ˙ϕn`1 `ϕn 2 ,
to compute approximations ϕ n of ϕpnδtq. This system is usually initialized with Υ ´1{2 " |ϕp´δt{2q| 2 or by second order approximation of |ϕp´δt{2q| 2 . This method is linearly implicit (recall that σ " 1). Moreover, it is known to preserve exactly the L 2 -norm and the discrete energy [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF]:

(11) E rlx pϕ, Υq " 1 4 ż R d }∇ϕ} 2 dx `β 2 ż R d Υ|ϕ| 2 dx ´β 4 ż R d Υ 2 dx. Indeed Re ˆΥn`1{2 ϕ n`1 `ϕn 2 ϕ n`1 ´ϕn ˙" Υ n`1{2 2 ´|ϕ n`1 | 2 ´|ϕ n | 2 ¯. But Υ n`1{2 ´|ϕ n`1 | 2 ´|ϕ n | 2 ¯" Υ n`1{2 ´|ϕ n`1 | 2 ´|ϕ n | 2 ¯`Υ n´1{2 ´|ϕ n | 2 ´|ϕ n | 2 " ´Υn`1{2 |ϕ n`1 | 2 ´Υn´1{2 |ϕ n | 2 ¯´´Υ n`1{2 ´Υn´1{2 ¯|ϕ n | 2 .
Using the definition of Υ ¨`1{2 in [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF], a simple computation leads to

´Υn`1{2 ´Υn´1{2 ¯|ϕ n | 2 " ´Υn`1{2 ¯2 ´´Υ n´1{2 ¯2 2 .
We therefore conclude that Re ˆΥn`1{2

ϕ n`1 `ϕn 2 ϕ n`1 ´ϕn ˙" ´Υn`1{2 |ϕ n`1 | 2 ´Υn´1{2 |ϕ n | 2 ¯`´Υ n`1{2 ¯2 ´´Υ n´1{2 ¯2 2 ,
which allows to prove the conservation of the discrete energy [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF].

It is interesting to note the consistency of the energy associated to relaxation scheme with the energy (2) for cubic nonlinear Schrödinger equation E rlx pϕ, |ϕ| 2 q " Epϕq.

The relaxation method was proved to converge in [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] but consistency analysis was missing. We present it in the next subsection.

3.2.

Consistency analysis for NLS equation with cubic nonlinearity. The aim of this subsection is to prove that this method has temporal order 2 under fairly general assumptions. This fact is supported by numerical evidences in the literature for years. We provide the first rigorous proof below.

The first equation in [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF] is the discrete equivalent of the continuous constraint Υ " |ϕ| 2 . In particular, this constraint is not an evolution equation. Therefore, we use the ideas introduced in [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] and rewrite the continuous equation (1) (recall that V " 0, λ " 0 and Ω " 0 and σ " 1) as the system

(12) $ ' ' & ' ' % iB t ϕ `1 2 ∆ϕ " βΥϕ, B t Υ " 2Repvϕq, iB t v `1 2
∆v " βpB t Υϕ `ΥB t ϕq.

The discrete system (10) has a discrete augmented equivalent (see [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF][START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF]). Let us denote by v n`1 2 " ϕ n`1 ´ϕn δt the discrete time derivative of ϕ n and define the nonlinearities as

(13) $ ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' % Φ n`1 2 " Υ n`1 2 ˆϕn`1 `ϕn 2 ˙, Ξ n`1 2 " 2Re ˆvn`1 2 ˆϕn`1 `ϕn 2 ˙˙, V n`1 2 " ˆΥn`3 2 `Υn´1 2 2 ˙ˆv n`3 2 `2v n`1 2 `vn´1 2 4 2Re ˆvn`1 2 ˆϕn`1 `ϕn 2 ˙˙ˆϕ n`2 `ϕn`1 `ϕn `ϕn´1 4
˙.

The augmented system writes ( 14)

$ ' ' ' ' ' ' ' & ' ' ' ' ' ' ' % i ϕ n`2 ´ϕn`1 δt `1 2 ∆ ˆϕn`2 `ϕn`1 2 ˙" βΦ n`3 2 , p14.bq Υ n`3 2 ´Υn´1 2 2δt " Ξ n`1 2 , p14.aq i v n`3 2 ´vn´1 2 2δt `1 2 ∆ ˆvn`3 2 `2v n`1 2 `vn´1 2 4 ˙" βV n`1 2 . p14.cq This system allows to compute pϕ n`2 , Υ n`3 2 , v n`3 2 q from X n :" ´Υn´1 2 , Υ n`1 2 , ϕ n´1 , ϕ n , ϕ n`1 , v n´1 2 , v n`1 2 ¯.
We consider the system (14) as the mapping

X n Þ Ñ X n`1 .
The seven variables involved in X n are not independent. If they satisfiy the five relations

$ ' ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' ' % v n´1 2 " ϕ n ´ϕn´1 δt , v n`1 2 " ϕ n`1 ´ϕn δt , Υ n`1 2 `Υn´1 2 " 2|ϕ n | 2 , i ϕ n ´ϕn´1 δt " ˆ´1 2 ∆ `βΥ n´1{2 ˙ϕn `ϕn´1 2 , i ϕ n`1 ´ϕn δt " ˆ´1 2 ∆ `βΥ n`1{2 ˙ϕn`1 `ϕn 2 ,
then the seven variables in X n`1 satisfy the same five relations with n replaced by n `1. This fact is proved in [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF]. We describe now how to build the seven initial data in X 0 from ϕ 0 " ϕ in and Υ ´1 2 so that they satisfy the five relations above:

(15) $ ' ' ' ' ' ' ' ' ' & ' ' ' ' ' ' ' ' ' % Υ 1 2 " 2|ϕ 0 | 2 ´Υ´1 2 , ϕ ´1 " ˆ2 ´iδt ˆ´1 2 ∆ `βΥ ´1{2 ˙˙´1 ˆ2 `iδt ˆ´1 2 ∆ `βΥ ´1{2 ˙˙ϕ 0 , ϕ 1 " ˆ2 `iδt ˆ´1 2 ∆ `βΥ 1{2 ˙˙´1 ˆ2 ´iδt ˆ´1 2 ∆ `βΥ 1{2 ˙˙ϕ 0 , v ´1 2 " pϕ 0 ´ϕ´1 q{δt, v 1 2
" pϕ 1 ´ϕ0 q{δt.

Let us define the operators

A " pi ´δt∆{4q ´1pi `δt∆{4q and B " pi ´δt∆{4q ´1, and the matrix of operators

C " ¨I 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 B 0 0 0 0 0 0 0 B 0 0 0 0 0 0 0 B 0 0 0 0 0 0 0 B 0 0 0 0 0 0 0 B ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' . The mapping X n Þ Ñ X n`1 reads (16) ¨Υn`1 2 Υ n`3 2 ϕ n ϕ n`1 ϕ n`2 v n`1 2 v n`3 2 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' " ¨0 I 0 0 0 0 0 I 0 0 0 0 0 0 0 0 0 I 0 0 0 0 0 0 0 I 0 0 0 0 0 0 A 0 0 0 0 0 0 0 0 I 0 0 0 0 0 A A ´I‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' ¨Υn´1 2 Υ n`1 2 ϕ n´1 ϕ n ϕ n`1 v n´1 2 v n`1 2 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' `δtC ¨0 2Ξ n`1 2 0 0 βΦ n`3 2 0 2βV n`1 2 ‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ ' .
In a more compact form, we define B and M so that the mapping ( 16) reads [START_REF] Delfour | Finite-difference solutions of a nonlinear Schrödinger equation[END_REF] X n`1 " BX n `δtCMpX n , X n`1 q.

We introduce the hypotheses that will allow us to prove our consistency and convergence result for the classical relaxation method [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF] in Theorem 8.

Remark 1. The results below extend to more general cases. In particular, one may treat the case where V is non zero smooth autonomous potential such that the multiplication by V is a bounded linear operator between Sobolev spaces and the case where V is a nonautonomous such operator with sufficient regularity with respect to time.

Hypotheses 2. We fix d P t1, 2, 3u and s ą d{2. We assume ϕ 0 P H s`4 pR d q is given. We denote by T ˚ą 0 the existence time of the maximal solution ϕ of the Cauchy problem (1) (with V " 0, λ " 0, Ω " 0 and σ " 1) in H s`4 pR d q. We assume there exists δt 0 ą 0 such that τ Þ Ñ ϕpτ, ¨q is a smooth map from p´δt 0 , T ˚q to H s`4 pR d q. Moreover, we assume that there exists R 1 , R 2 ą 0 such that for all Υ ´1{2 P H s`4 pR d q with }Υ ´1{2 } H s`4 ď R 1 , the numerical solution X n (with initial datum (15)) is uniquely determined by ( 16) for all δt P p0, δt 0 q and all n P N such that nh ď T , and it satisfies }X n } pH s`2 pR d qq 7 ď R 2 for all such n.

Remark 3. The hypotheses above on the exact solution are fullfilled in several cases. For example, for the exact solution ϕ, it is well known (see [START_REF] Ginibre | On a class of nonlinear schrödinger equations part i, ii[END_REF]) that T ˚" `8 in at least two cases:

' if ϕ in has small H s`4 -norm and β ă 0 ' if β ą 0 and ϕ in P H s`4 pR d q. For the numerical solution, they are fullfilled provided T ˚ă `8 (see [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF]) and also when T ˚" `8 and β ą 0 (see [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF]).

Let t n " nδt denote the discrete times and t Þ Ñ Xptq the vector [START_REF] Gazeau | Probability and pathwise order of convergence of a semidiscrete scheme for the stochastic Manakov equation[END_REF] Xptq "

¨|ϕpt ´δt{2, ¨q| 2 |ϕpt `δt{2, ¨q| 2 ϕpt ´δt, ¨q ϕpt, ¨q ϕpt `δt, ¨q B t ϕpt ´δt{2, ¨q B t ϕpt `δt{2, ¨q‹ ‹ ‹ ‹ ‹ ‹ ‹ ‹ '
.

Using the definition of M (see ( 13) and ( 17)), the fact that H s pR d q, is an algebra since s ą d{2, and the fact that the exact and numerical solutions stay in a bounded set of H s pR d q, it is easy to prove the following lemma.

Lemma 4. Assume Hypotheses 2 is satisfied. There exists C ą 0 such that for all δt P p0, δt 0 q, all n P N such that pn `1qδt ď T

}MpX n , X n`1 q ´MpXpt n q, Xpt n`1 qq} pH s pR d qq 7 ď C `}X n ´Xpt n q} pH s pR d qq 7 `}X n`1 ´Xpt n`1 q} pH s pR d qq 7 ˘.
Note that the constant C depends only on the initial data.

Before starting the proof of our main result of this section (see Theorem 8), we state and prove another lemma.

Lemma 5. There exists a constant b ą 0 such that the operator B defined in ( 16) and ( 17) satisfies for all δt ą 0, and n P N,

(19) ~Bn C~ď b,
where ~¨~is the norm of linear continuous operators from pH s pR d qq 7 to itself.

Proof. The operator B is defined by three diagonal blocks

B 1 " ˆ0 I I 0 ˙, B 2 " ¨0 I 0 0 0 I 0 0 A ', and B 3 " ˆ0 A A A ´I˙.
The first block B 1 is an isometry from pH s pR d qq 2 to itself and so are all its powers. The powers of the second block B 2 read for all n ě 2

B n 2 " ¨0 0 A n´2 0 0 A n´1 0 0 A n '.
Since all the powers of A are of norm less than one, we infer that the norm of B n 2 is less than ? 3. Since the norm of B is less than one, we infer that the norm of

B n 2 ¨B 0 0 0 B 0 0 0 B ', is less than ? 3.
The block B 3 can be diagonalized by blocks as

B 3 " ˆI I ´I A ˙ˆ´I 0 0 A ˙ˆI I ´I A ˙´1 ,
so that for all n P N,

(20) B n 3 " ˆI I ´I A ˙ˆp´Iq n 0 0 A n ˙ˆI I ´I A ˙´1 .
Therefore, the last 2 ˆ2 block of B n C reads

B n 3 ˆB 0 0 B ˙" ˆpp´Iq n A `An qBpI `Aq ´1 pA n ´p´Iq n qBpI `Aq ´1 pp´Iq n`1 A `An`1 qBpI `Aq ´1 pp´Iq n `An`1 qBpI `Aq ´1˙.
Since BpI `Aq ´1 " p1{2qI, we infer

B n 3 ˆB 0 0 B ˙" 1 2
ˆpp´Iq n A `An q pA n ´p´Iq n q pp´Iq n`1 A `An`1 q pp´Iq n `An`1 q

˙.
Since all the powers of A are of norm less than 1, we infer that

@n P N, B n 3 ˆB 0 0 B ˙ ď 4.
This proves the result.

Remark 6. The powers of the operator B are not uniformly bounded. However, the powers of B mutiplied by C are uniformly bounded as shown above. The main reason is that the third matrix in the right hand side of (20) becomes singular at the end of the spectrum of ∆.

Lemma 7. Assume d, s, R 1 , δt 0 and ϕ in are given as in Hypotheses 2. There exists c ą 0 such that for all Υ ´1{2 P H s`4 with }Υ ´1{2 } H s`4 ď R 1 and all δt P p0, δt 0 q, (21)

}X 0 ´Xp0q} pH s`2 pR d qq 7 ď c `}Υ ´1{2 ´|ϕp´δt{2q| 2 } H s`4 pR d q `δt 2 ˘.
Proof. The H s`2 -norm of each of the seven components of the vector X 0 ´Xp0q is estimated separately. The H s`2 -norm of the first component Υ ´1{2 ´|ϕp´δt{2q| 2 is bounded by the H s`4 -norm of the same quantity. For the H s`2 -norm of the second component, we define f ptq " |ϕptq| 2 , which is a smooth function from p´δt 0 , δt 0 q to H s`4 , thanks to Hypotheses 2. We may write using a Taylor formula at 0

f ˆ´δt 2 ˙" |ϕ 0 | 2 ´2Re pϕ 0 B t ϕp0qq δt 2 `ż ´δt{2 0 ˆ´δt 2 ´σ˙f 2 pσqdσ,
and similarly

f ˆδt 2 ˙" |ϕ 0 | 2 `2Re pϕ 0 B t ϕp0qq δt 2 `ż δt{2 0 ˆδt 2 ´σ˙f 2 pσqdσ.
Therefore, one has

Υ 1{2 ´f ˆδt 2 ˙" 2|ϕ 0 | 2 ´Υ´1{2 ´f ˆδt 2 " f ˆ´δt 2 ˙´ż ´δt{2 0 ˆ´δt 2 ´σ˙f 2 pσqdσ ´ż δt{2 0 ˆδt 2 ´σ˙f 2 pσqdσ ´Υ´1{2 .
By triangle inequality, We infer that

(22) }Υ 1{2 ´|ϕpδt{2q| 2 } H s`2 ď }Υ 1{2 ´|ϕpδt{2q| 2 } H s`4 ď c `}Υ ´1{2 ´|ϕp´δt{2q| 2 } H s`4 `δt 2 ˘,
where c " maxp1, sup σPp´δt0q,δt0 }f 2 pσq} H s`2 q only depends on the exact solution of (1). We now estimate the H s`2 -norm of the fifth component of X 0 ´Xp0q. Note that the H s`2 -norm of the third component can be estimated the very same way and that the H s`2 -norm of the fourth component is zero. We start with the identity

ϕ 1 " ϕ 0 ´iδt ˆ´1 2 ∆ `βΥ 1{2 ˙ϕ1 `ϕ0 2 ,
and we denote by rpδtq the consistency error defined by rpδtq " ϕpδtq ´ϕ0 `iδt

ˆ´1 2 ∆ `βΥ 1{2 ˙ϕpδtq `ϕ0 2 .
Using the fact that t Þ Ñ ϕpt, ¨q is a smooth function from p´δt 0 , δt 0 q to H s`4 , we may write another Taylor expansion to obtain rpδtq

" ϕ 0 `δtB t ϕp0q `δt 2 2 B 2 t ϕp0q `ż δt 0 pδt ´σq 2 2 B 3 t ϕpσqdσ ´ϕ0 `iδt ˆ´1 2 ∆ `βΥ 1{2 ˙˜ϕ 0 `δt 2 B t ϕp0q `1 2 ż δt 0 pδt ´σqB 2 t ϕpσqdσ " ´iδt 2 βRe ˆi 2 ϕ 0 ∆ϕ 0 ˙ϕ0 `iδtβ `|ϕ 0 | 2 ´Υ´1{2 ˘ˆϕ 0 ´i δt 2 ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 i δt 2 ˆ´1 2 ∆ `βΥ 1{2 ˙ż δt 0 pδt ´σqB 2 t ϕpσqdσ `ż δt 0 pδt ´σq 2 2 B 3 t ϕpσqdσ " ´iδt 2 βRe ˆi 2 ϕ 0 ∆ϕ 0 ˙ϕ0 `iδtβ `|ϕ 0 | 2 ´|ϕ p´δt{2q | 2 ˘ˆϕ 0 ´i δt 2 ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 iδtβ `|ϕ p´δt{2q | 2 ´Υ´1{2 ˘ˆϕ 0 ´i δt 2 ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 i δt 2 ˆ´1 2 ∆ `βΥ 1{2 ˙ż δt 0 pδt ´σqB 2 t ϕpσqdσ `ż δt 0 pδt ´σq 2 2 B 3 t ϕpσqdσ. Using ˇˇˇϕ ˆ´δt 2 ˙ˇˇˇ2 " |ϕp0q| 2 `δtRe ˆiϕ 0 ˆ´1 2 ∆ϕ 0 `β|ϕ 0 | 2 ϕ 0 ˙˙`ż ´δt{2 0 p´δt{2 ´σqB 2 t p|ϕ| 2 qpσqdσ, we obtain rpδtq " iδtβ `|ϕ p´δt{2q | 2 ´Υ´1{2 ˘ˆϕ 0 ´i δt 2 ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 δt 3 2 βRe ˆi 2 ϕ 0 ∆ϕ 0 ˙ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 ´iδtβ ˆϕ0 ´i δt 2 ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 ˙ż ´δt{2 0 p´δt{2 ´σqB 2 t p|ϕ| 2 qpσqdσ `i δt 2 ˆ´1 2 ∆ `βΥ 1{2 ˙ż δt 0 pδt ´σqB 2 t ϕpσqdσ `ż δt 0 pδt ´σq 2 2 B 3 t ϕpσqdσ.
Note that we have

› › › › › ż δt 0 pδt ´σq 2 2 B 3 t ϕpσqdσ › › › › › H s`2 ď cδt 3 , › › › › › i δt 2 ˆ´1 2 ∆ `βΥ 1{2 ˙ż δt 0 pδt ´σqB 2 t ϕpσqdσ › › › › › H s`2 ď cδt 3 , › › › › › iδtβ ˆϕ0 ´i δt 2 ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 ˙ż ´δt{2 0 p´δt{2 ´σqB 2 t p|ϕ| 2 qpσqdσ › › › › › H s`2 ď cδt 3 , and › › › › δt 3 2 βRe ˆi 2 ϕ 0 ∆ϕ 0 ˙ˆ´1 2 ∆ `β|ϕ 0 | 2 ˙ϕ0 › › › › H s`2 ď cδt 3 ,
where c doesn't depend on δt. Then we infer from the estimates above that

(23) }rpδtq} H s`2 ď cδt `}Υ ´1{2 ´ϕp´δt{2q} H s`2 `δt 2 ˘.
Now, we denote by epδtq the fifth component ϕ 1 ´ϕpδtq of the vector X 0 ´Xp0q and we have

epδtq " ´i δt 2 ˆ´1 2 ∆ `βΥ 1{2 ˙epδtq ´rpδtq.
We want to estimate the H s`2 -norm of epδtq using this relation and the estimate [START_REF]Lev Pítajevskíj and Sandro Stringari[END_REF]. To this aim, we take α P N d with |α| " α 1 `¨¨¨`α d ď s `2 and differentiate the relation above to obtain that

B α x epδtq " iδt 1 4 ∆B α x epδtq ´i δt 2 βB α x pΥ 1{2 epδtqq ´Bα x rpδtq.
Multiplying this relation by B α x epδtq integrating over R d , and taking the real part, we obtain

(24) }B α x epδtq} 2 2 " ´β δt 2 Re ˆi ż R d B α x epδtqB α x pΥ 1{2 epδtqq ˙´Re ˆżR d B α x epδtqB α x rpδtq ˙.
When α " 0 N d , the first term on the right hand side in the equation above vanishes and }epδtq} 2 2 ď }epδtq} 2 }rpδtq} 2 , using Cauchy-Schwartz inequality. We infer, [START_REF] Shen | Vortex solitons under competing nonlocal cubic and local quintic nonlinearities[END_REF] }epδtq} 2 ď }rpδtq} 2 .

Now, when α ‰ 0 N d , the Leibniz' rule applied in [START_REF] Sanz-Serna | Methods for the numerical solution of the nonlinear Schroedinger equation[END_REF] provides us with:

}B α x epδtq} 2 2 " ´β δt 2 ÿ tk:kďαu cpk, αqRe ˆi ż R d B α x epδtqB k x Υ 1{2 B α´k x epδtq ˙´Re ˆżR d B α x epδtqB α x rpδtq ˙,
where cpk, αq are integers. Note that since Υ 1{2 is real valued, the term corresponding to k " 0 N d in the sum vanishes. Then, with Cauchy-Schwarz inequality, we infer

}B α x epδtq} 2 2 ď β δt 2 ÿ tk:kďα,|k|ě1u cpk, αq}B α x epδtq} 2 }B k x Υ 1{2 B α´k x epδtqq} 2 `}B α x epδtq} 2 }B α x rpδtq} 2 . ď β δt 2 ÿ tk:kďα,|k|ě1u cpk, αq}B α x epδtq} 2 }B k x Υ 1{2 } 8 }B α´k x epδtqq} 2 `}B α x epδtq} 2 }B α x rpδtq} 2 .
Since k ď α in the sum above, we have |k| ď |α|. For such k, one has

}B k x Υ 1{2 } 8 ď c}B k x Υ 1{2 } H pd`1q{2 ď c}Υ 1{2 } H pd`1q{2`|k| ď c}Υ 1{2 } H pd`1q{2`s`2 ď c}Υ 1{2 } H s`4
, where c is the Sobolev constant of the injection from H pd`1q{2 pR d q to L 8 pR d q and where we have used the fact that d P t1, 2, 3u. Since the H s`4 -norm of Υ 1{2 is controlled by [START_REF] Oelz | Analysis of a relaxation scheme for a nonlinear Schrödinger equation occurring in plasma physics[END_REF], we have }B α

x epδtq} 2 ď c α }epδtq} H |α|´1 `}B α x rpδtq} 2 , where c α does not depend on δt P p0, δt 0 q. Then, by induction on |α| P t0, . . . , s `2u starting with (25) for α " 0 N d , we have for some positive constant c s`2 :

}B α x epδtq} 2 ď c s`2 }rpδtq} H s`2 ,
for all δt P p0, δt 0 q and all α P N d with |α| ď s `2. Therefore, we have }epδtq} H s`2 is controlled by }rpδtq} H s`2 and the conclusion for the fifth term follows using ( 23):

(26) }ϕ 1 ´ϕpδtq} H s`2 ď cδt `}Υ ´1{2 ´ϕp´δt{2q} H s`4 `δt 2 ˘,
where c does not depend on δt P p0, δt 0 q. It remains to estimate the H s`2 -norm of the sixth and seventh components of the vector X 0 ´Xp0q. We only give the details for the seventh term since the computation is similar and even simpler for the sixth term. Let us denote by ppδtq the consistency error defined as

ppδtq " ϕpδtq ´ϕ0 δt ´Bt ϕ ˆδt 2 ˙.
Using a Taylor expansion, since the exact solution is a smooth function from p´δt 0 , δt 0 q to H s`2 , we have

(27) }ppδtq} H s`2 ď cδt 2 .
Let us denote by qpδtq " v 1{2 ´Bt ϕ ˆδt 2 ˙the seventh component of X 0 ´Xp0q. We have qpδtq ´ppδtq " ϕ 1 ´ϕpδtq δt .

Using estimates ( 26) and ( 27) we have by triangle inequality,

}qpδtq} H s`2 ď }ppδtq} H s`2 `1 δt }ϕ 1 ´ϕpδtq} H s`2 ď cδt 2 `c `}Υ ´1{2 ´ϕp´δt{2q} H s`4 `δt 2 ˘.
This concludes the proof of the lemma.

We prove below that the relaxation method (10) is of order 2.

Theorem 8. Assume d, s, R 1 ą 0, ϕ in P H s`4 pR d q, T ă T ‹ are given and satisfy Hypotheses 2. There exists C ą 0 and δt 0 ą 0 (smaller than the one in the hypotheses) such that for all Υ ´1{2 P H s`4 pR d q with }Υ ´1{2 } H s`4 pR d q ď R 1 , all n P N and all δt P p0, δt 0 q with nδt ď T ,

}X n ´Xpt n q} pH s pR d qq 7 ď C `}Υ ´1{2 ´|ϕp´δt{2q| 2 } H s`4 pR d q `δt 2 ˘.

Proof. First, the initial datum X 0 is computed from ϕ 0 " ϕ in and Υ ´1{2 using [START_REF] Dahlby | Plane wave stability of some conservative schemes for the cubic Schrödinger equation[END_REF]. Therefore, using Lemma 7, there exists a constant c ą 0 such that for all δt P p0, δt 0 q and all }Υ ´1{2 } H s`4 pR d q ď R 1 , we have the estimate ( 21)

}X 0 ´Xp0q} pH s`2 pR d qq 7 ď c `}Υ ´1{2 ´|ϕp´δt{2q| 2 } H s`4 pR d q `δt 2 ˘.
Second, we define the consistency error of the relaxation scheme at time t n " nδt ď T by the formula R n pδtq " BXpt n q `δtCMpXpt n q, Xpt n`1 qq ´Xpt n`1 q.

Substracting this definition from [START_REF] Delfour | Finite-difference solutions of a nonlinear Schrödinger equation[END_REF], we obtain (29) X n`1 ´Xpt n`1 q " BpX n ´Xpt n qq `δtC pMpX n , X n`1 q ´MpXpt n q, Xpt n`1 qqq `Rn pδtq.

From now on, we set for all n P N such that nδt ď T , e n " X n ´Xpt n q. Iterating the relation above, we obtain, as long as nδt ď T ,

e n " B n e 0 `δt n´1 ÿ k"0 B n´k´1 C pMpX k , X k`1 q ´MpXpt k q, Xpt k`1 qqq `n´1 ÿ k"0 B n´k´1 R k pδtq.
This implies

}e n } pH s pR d qq 7 ď }B n CC ´1e 0 } pH s pR d qq 7 `δt n´1 ÿ k"0 }B n´k´1 C pMpX k , X k`1 q ´MpXpt k q, Xpt k`1 qqq } pH s pR d qq 7 `n´1 ÿ k"0 }B n´k´1 CC ´1R k pδtq} pH s pR d qq 7 ď b}C ´1e 0 } pH s pR d qq 7 `δtb n´1 ÿ k"0 } pMpX k , X k`1 q ´MpXpt k q, Xpt k`1 qqq } pH s pR d qq 7 `b n´1 ÿ k"0 }C ´1R k pδtq} pH s pR d qq 7 ď b}e 0 } pH s`2 pR d qq 7 `Cδtb n´1 ÿ k"0 `}e k } pH s pR d qq 7 `}e k`1 } pH s pR d qq 7 b n´1 ÿ k"0 }R k pδtq} pH s`2 pR d qq 7 ,
using Lemmas 4 and 5. Then

p1 ´Cδtbq}e n } pH s pR d qq 7 ď b}e 0 } pH s`2 pR d qq 7 `2Cδtb n´1 ÿ k"0 }e k } pH s pR d qq 7 `b n´1 ÿ k"0 }R k pδtq} pH s`2 pR d qq 7 .
Let δt be small enough to ensure that It is however known to not conserve energy functional. As a generalization of the classical relaxation method [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF], which is designed for the special cubic case (σ " 1), we propose the following method, which allows for a general nonlinearity exponent σ P N ‹ .

We propose to substitute system (31) by ( 32)

# γ σ pt, xq " |ϕpt, xq| 2σ , iB t ϕpt, xq " ´1 2 ∆ϕpt, xq `βγ σ ϕpt, xq.
The modification seems ligth but allows to build an energy preserving scheme. The second equation is approximated to second order at time t n`1{2 by

i ϕ n`1 ´ϕn δt " ˆ´1 2 ∆ `βγ σ n`1{2
˙ϕn`1 `ϕn 2 .

We now want to find an approximation of γ σ " |ϕ| 2σ " γ σ´1 |ϕ| 2 that allow energy conservation following the ideas that were presented in section 2. As for the classical relaxation method, we note that Re ˆγσ

n`1{2 ϕ n`1 `ϕn 2 ϕ n`1 ´ϕn ˙" γ σ n`1{2 2 ´|ϕ n`1 | 2 ´|ϕ n | 2 ¯.
The last term also reads

γ σ n`1{2 ´|ϕ n`1 | 2 ´|ϕ n | 2 ¯" γ σ n`1{2 ´|ϕ n`1 | 2 ´|ϕ n | 2 ¯`γ σ n´1{2 ´|ϕ n | 2 ´|ϕ n | 2 " ´γσ n`1{2 |ϕ n`1 | 2 ´γσ n´1{2 |ϕ n | 2 ¯´´γ σ n`1{2 ´γσ n´1{2 ¯|ϕ n | 2 .
The only choice that allow to preserve energy is to choose

´γσ n`1{2 ´γσ n´1{2 ¯|ϕ n | 2 " σ σ `1 ´γσ`1 n`1{2 ´γσ`1 n´1{2 ¯.
Moreover, we remark that

1 σ `1 γ σ`1 n`1{2 ´γσ`1 n´1{2 δt " 1 σ γ σ n`1{2 ´γσ n´1{2 δt |ϕ n | 2
is a second order approximation of γ σ " γ σ´1 |ϕ| 2 at time t " t n . This method is therefore designed so that it preserves exactly the following energy

E rlx pϕ, γq " 1 4 ż R d }∇ϕ} 2 dx `β 2 ż R d γ σ |ϕ| 2 dx ´β σ 2pσ `1q ż R d γ σ`1 dx.
Since at continuous level γ σ " |ϕ| 2σ , E rlx pϕ, γq reduces to the true energy

1 4 ż R d }∇ϕ} 2 dx `β 2σ `2 ż R d |ϕ| 2σ`2 dx.
Moreover, the generalized relaxation method preserves the L 2 -norm of the solution. We present in Theorem 9 a more general method, which includes possibly non zero other terms in the equation (see Section 4.2).

The first equation of ( 36) is implicit and local. In order to solve it, we propose two different approaches. The first one, for small integer values of σ, consists in using exact formulae, since the equation is polynomial of low degree in these cases. The second approach, is to use a fixed point method, as described in (35). The second equation is explicit, as in the classical relaxation method [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF]. The third equation of ( 36) is implicit and non-local. We propose again two different approaches to solve it. The first approach, presented here for x P T d δ , consists in using another fixed point iteration method: one starts from ϕ 0 n`1 " ϕ n and computes ϕ p`1 n`1 from ϕ p n`1 using the iterative procedure

ˆ1 `i δt 2 ξ 2 2 ˙Fpϕ p`1 n`1 qpξq (37) " ˆ1 ´i δt 2 ξ 2 2 ˙Fpϕ n qpξq ´iδtF ˆ´V `βγ σ n`1{2 `λpU ˚Υn`1{2 q ´Ω.R ¯ˆϕ p n`1 `ϕn 2 ˙˙pξq,
where F stands for the Fourier transform in space as defined in Appendix A, and we set ϕ n`1 " lim pÑ8 ϕ p n`1 . In practice, the iterative procedure ϕ p n`1 Ñ ϕ p`1 n`1 stops whenever the L 2 -norm of the difference between two consecutive steps is below some small tolerance parameter. Note that we decided to implicit the Laplace operator in the iterative procedure (37) in order to ensure that the Sobolev regularity of ϕ p n`1 does no decrease a priori when p increases. Alternatively, the other approach to solve the last equation of (36) consists in following [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF], and using the linearity of the equation. Let us introduce the new unknown ϕ n`1{2 " pϕ n`1 `ϕn q{2, so that the equation reads

ˆI ´i δt 2 ∆ 2 `i δt 2 V `i δt 2 βγ σ n`1{2 `i δt 2 λpU ˚Υn`1{2 q ´i δt 2 Ω ¨R˙ϕ n`1{2 " ϕ n .
This equation can also be written as (38) Lϕ n`1{2 " P ´1ϕ n , where L " ˆI `i δt 2 P ´1 ´V `βγ σ n`1{2 `λpU ˚Υn`1{2 q ´Ω ¨R¯˙, and

P " ˆI ´i δt 2 ∆ 2 ˙.
The precondionning operator P can be easily inverted in Fourier space and the solution of (38) can be obtained by a Krylov method. Note that other choices of preconditionning operator are possible (see [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF]).

In the following, we shall use the notation pϕ n`1 , γ n`1{2 , Υ n`1{2 q " Φ rlx δt pϕ n , γ n´1{2 , Υ n´1{2 q, for the method (36) above. 4.2.1. Energy preservation property for generalized relaxation methods. The generalized relaxation method (36) is designed to preserve exactly the following energy:

E rlx pϕ, γ, Υq " ż R d ˆ1 4 }∇ϕ} 2 `1 2 V |ϕ| 2 `β 2 γ σ |ϕ| 2 ´β σ 2pσ `1q γ σ`1 ˙dx (39) `żR d ˆλ 2 pU ˚Υq|ϕ| 2 ´λ 4 pU ˚ΥqΥ ´Ω 2 ϕRϕ ˙dx,
as we prove in Theorem 9. Note that this energy is consistent with the energy (2) of equation ( 1) in the sense that E rlx pϕ, |ϕ| 2 , |ϕ| 2 q " Epϕq.

Theorem 9. The generalized relaxation method (36) applied to the equation (1) with initial data ϕ 0 " ϕ in , γ ´1{2 and Υ ´1{2 preserves exactly the L 2 norm and the energy functional E rlx defined in (39) in the sense that for all n P N such that a solution of (36) is defined, we have

(40) }ϕ n`1 } 2 " }ϕ 0 } 2 and E rlx pϕ n`1 , γ n`1{2 , Υ n`1{2 q " E rlx pϕ 0 , γ ´1{2 , Υ ´1{2 q.
Proof. Multiplying the last equation of (36) by ϕ n`1 ´ϕn , integrating in space, and taking the real part, one finds that a sum of 5 terms is equal to

(41) ż R d Re `i|ϕ n`1 ´ϕn | 2 ˘dx " 0.
The first term reads

(42) ż R d Re ˆ´1 2
`∆ ϕ n`1 `ϕn 2 ˘`ϕ n`1 ´ϕn ˘˙dx " 1 4 ż R d }∇ϕ n`1 } 2 dx ´1 4 ż R d }∇ϕ n } 2 dx.
The second term reads

(43) ż R d Re ˆ`V ϕ n`1 `ϕn 2 ˘`ϕ n`1 ´ϕn ˘˙dx " 1 2 ż R d V |ϕ n`1 | 2 dx ´1 2 ż R d V |ϕ n | 2 dx.
Using the first equation of (36), the third term reads

ż R d Re ˆβγ σ n`1{2 `ϕn`1 `ϕn 2 ˘`ϕ n`1 ´ϕn ˘˙dx " β ˆ1 2 ż R d γ σ n`1{2 |ϕ n`1 | 2 dx ´σ 2pσ `1q ż R d γ σ`1 n`1{2 dx β ˆ1 2 ż R d γ σ n´1{2 |ϕ n | 2 dx ´σ 2pσ `1q ż R d γ σ`1 n´1{2 dx ˙. (44)
Using the fact that the convolution kernel U is real-valued and symmetric with respect to the origin, and also using the second equation of (36), the fourth term reads

λ ż R d Re ´pU ˚Υn`1{2 q ϕ n `ϕn`1 2 ϕ n`1 ´ϕn ¯dx (45) " λ ´1 2 ż R d pU ˚Υn`1{2 q|ϕ n`1 | 2 dx ´1 4 ż R d pU ˚Υn`1{2 qΥ n`1{2 dx λ´1 2 ż R d pU ˚Υn´1{2 q|ϕ n | 2 dx ´1 4 ż R d pU ˚Υn´1{2 qΥ n´1{2 dx ¯.
Finally, thanks to the fact that the operator R is symmetric, the last term reads

(46) ´Ω ż R d Re ˆ´R ϕ n`1 `ϕn 2 ¯ϕn`1 ´ϕn ˙" ´Ω ż R d pϕ n`1 Rϕ n`1 q `Ω ż R d pϕ n Rϕ n q.
The proof of the preservation of the energy (40) is completed by adding altogether (42), ( 43), ( 44), (44), and (46), and using relation (41). The preservation of the L 2 -norm follows from multiplying the last relation in (36) by ϕ n`1 `ϕn , integrating in space, and taking the imaginary part.

Numerical experiments

In this section, we make some numerical experiments and show that the classical and generalized relaxation methods are efficient methods that preserve mass and energy to machine epsilon. 5.1. One dimensional example: the one-dimensional quintic and septic NLS equation. We present in this subsection some numerical experiments to show the efficiency of the generalized relaxation (34) compared with the Crank-Nicolson scheme [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF] with ϕp0, xq " ϕ in pxq and pt, xq P r0, T s ˆpx , x r q. To deal with the space variable, we discretize the space operators using Fourier spectral approximation and consider periodic boundary conditions. The spatial mesh size is δx ą 0 with δx " px r ´x q{J with J " 2 P , P P N ˚. The time step is δt " T {N for some N P N ‹ . The grid points and the discrete times are

x j :" x `jδx, t n :" nδt, j " 0, 1, ¨¨¨, J, n " 0, 1, ¨¨¨, N.

Let ϕ j,n be the approximation of ϕpt n , x j q satisfying ϕ j,n " 1 J J´1 ÿ k"0 φk,n ω jk J , j " 0, ¨¨¨, J ´1, where φk,n denotes the discrete Fourier transform of sequence pϕ j,n q j given by φn k "

J´1 ÿ q"0 ϕ n q ω ´jk J , k " ´J 2 , ¨¨¨, J 2 ´1,
where ω J " exp p´2iπ{Jq. We also apply the discrete Fourier transform to the approximation γ j,n`1{2 of |ϕpt n`1{2 , x j q| 2 . Let us define the discrete gradient operator

∇ d { p∇ d vq k " iµ k vk , v P C J .
Let us denote by Π d the projection operator

Π d : C 0 prx , x r s, Cq Ñ C J ϕ Þ Ñ pϕpx j qq 0ďjďJ´1 .
We define the discrete r norm on C J as }v} r " ˜δx

J´1 ÿ j"0 |v j | r ¸1{r , v P C J , r ě 1, the mean M pvq " δx J´1 ÿ j"0 v j , v P C J , r ě 1,
and the discrete energies:

E d pvq " 1 4 }∇ d v} 2 2 `β 2pσ `1q › › |v| σ`1 › › 2 2
, and

E rlx,d pv, gq " 1 4 }∇ d v} 2 2 `β 2 M ˆgσ p|v| 2 ´σ σ `1 gq ˙.
As in the continuous case, we have E d pvq " E rlx,d pv, vq for any v P C J . Using these definitions, the energy conservation for the relaxation scheme is built through the following relative error (48) E E,δt " sup nPt0,¨¨¨,N u ˇˇErlx,dpΠdpϕexp0, ¨qq, Π d pϕ ex p´δt{2, ¨qqq ´Erlx,d ppϕ n j q j , pγ n´1{2 j q j q ˇĚ rlx,d pΠ d pϕ ex p0, ¨qq, Π d pϕ ex p´δt{2, ¨qqq .

For the Crank-Nicolson scheme, it is defined by (49) E E,δt " sup nPt0,¨¨¨,N u ˇˇE d pΠ d pϕ ex p0, ¨qqq ´Ed ppϕ n j q j q ˇĚ d pΠ d pϕ ex p0, ¨qqq .

We present on Figure 1 the evolution of E E,δt for various δt when σ " 2 and σ " 3 for both classical and generalized relaxation scheme compared to Crank-Nicolson scheme. The initial datum is chosen to be ϕ in pxq " expp´x 2 q and β " ´1. The time-space domain is r0, 1{2s ˆr´30, 30s. We consider periodic boundary conditions and the interval r´30, 30s is meshed with 2 13 `1 nodes. As it is known, the standard relaxation scheme does not preserve the energy when σ " 2 or σ " 3 and the error curve show a second order convergence. On the other hand, both generalized relaxation and Crank-Nicolson preserve energy to epsilon machine. where U pxq is the nonlocal kernel and α 1 and α 2 are real parameters. When α j ą 0, j " 1, 2, we are considering focusing nonlinearities, whereas when α j ă 0, the nonlinearities are defocusing. This model is generalized in two dimensions in [START_REF] Shen | Vortex solitons under competing nonlocal cubic and local quintic nonlinearities[END_REF] for the study of vortex solitons where B 2

x is replaced by the two-dimensional Laplace operator. Typically, the kernel is regular and is given by (51) U 1 pxq " Concerning the two dimensionsal case, we consider the kernel U 2 , α 1 " 1 and α 2 " ´0.02. The initial datum is chosen as a vortex beam with angular momentum ϕ in px 1 , x 2 q " Ar m e ´r2 {2 e imφ , where r " a x 2 1 `x2 2 , m " 1 is the topological charge, φ is such that x 1 `ix 2 " r exppiφq and A " 5.8 is the amplitude. The computational domain is r´8, 8s 2 and we use 256 Fourier modes. The final time of simulation is T " 10 and the time step is 5 ¨10 ´3. The width of the Gaussian kernel U 2 is determined by µ " 0.4. Once again, the energy is very well preserved, since the relative error energy is bounded by 10 ´11 at the end of the simulation (see Figure 5). We present in Figure 6 the evolution of the solution with respect to time. The 3D representation of the time evolution is presented in figure 6a. A 2D projection of it is displayed in Fig. 6b and the final solution in Fig. 6c. It is interesting to note that the same kind of breathing behaviour observed in the one dimensional setting is also present in 2D experiment.

1 t|x|ďµu p2µq d , x P R d , or (52) U 2 pxq " 1 pµ ? πq d exp `´|x| 2 {µ 2 ˘, x P R d , µ " 0.5 µ " 2.5
Nonlinear parameter α 2 " ´0.5

Nonlinear parameter α " 0.1 where L x3 " ´ipx 1 B x2 ´x2 B x1 q is the x 3 -component of the angular momentum and Ω represents the rotating frequency. The nonlinear parameters β and λ respectively describe the strength of the shortrange two-body interactions in the condensate and the strength of the dipolar interaction modeled with a Coulomb potential. The real-valued external trapping potential V is chosen as V pxq " pγ , where γ x1 ą 0, γ x2 ą 0 are dimensionless constants proportional to the trapping frequencies in the both directions. The normal n " pn 1 , n 2 , n 3 q T represents the dipole axis and we define n K " pn 1 , n 2 q T and B nK " n K ¨∇. The energy is given by (59)

Epϕq "

1 2 ż R 2 ˆ1 2 |∇ϕ| 2 `V pxq|ψ| 2 `β 2 |ϕ| 4 `λ 2 ψ|ϕ| 2 ´ΩϕL z ϕ ˙dx.
The initial datum is computed as a minizer of the energy on the intersection of the unit sphere of L 2 with the energy space ϕ in " arg min

}φ} 2 "1
Epφq.

We use a preconditioned nonlinear conjugate gradient method developed in [START_REF] Antoine | Efficient spectral computation of the stationary states of rotating Bose-Einstein condensates by preconditioned nonlinear conjugate gradient methods[END_REF]. Contrary to the previous subsection where the convolution kernel was regular, the computation of the nonlocal term ψ with the Coulomb potential is known to be a costly task. We have chosen to apply the technique developed in [START_REF] Vico | Fast convolution with free-space Green's functions[END_REF], which allows fast convolution using truncated Green's functions. The parameters of the simulation are Ω " 0.97, γ x1 " γ x2 " 1, λ " 175, β " p250 ´λq a 5{π. The computational domain is p´16, 16q 2 with 2 8 " 256 Fourier modes in each direction, and the computational time is set to T " 15 with the time step δt " 10 ´3. Initially, the dipole axis is n " p1, 0, 0q T and we change it to n " pcospπ{3q, sinpπ{3q, 0q t . We plot on Figure 7 the evolution of the solution from time t " 0 to time t " 20. We add the direction of the dipole axis on each frame. The axis are removed to make the presentation clearer. As it can be seen on Figure 8, the energy is very well preserved as expected.

t " 0 t " 2.5 t " 5.0 t " 7.5 t " 10 t " 12.5

t " 15 t " 17.5 t " 20 In this paper we have given for the first time a proof of the second order of the relaxation method introduce in [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF] for the cubic nonlinear Schrödinger equation. We also have extended the previous method to deal with general power law nonlinearites showing that this method is still an energy preserving method. Note that, since Crank-Nicolson methods are one-step methods, one may use composition techniques to achieve higher orders (4,6,8, etc) while preserving both an energy and the L 2 -norm, and this is not possible so straightforwardly for relaxation methods. Therefore, in future works we want to focus on energy preserving relaxation methods which have higher order.

Appendix A. Notation

In this section, the symbol F will denote different operators depending on whether x belongs to R d or x belongs to T d δ . When x P R d , F denotes the Fourier transform R d defined for ϕ P L 1 pR d q via the formula Fpϕqpξq " 1 p2πq d{2 ż R d ϕpxqe ´ix.ξ dx, ξ P R d , and its inverse F ´1 is defined for φ P L 1 pR d q through the formula

F ´1pφqpxq " 1 p2πq d{2 ż R d φpξqe `ix.ξ dξ, x P R d .
When x P T d δ " pR{pδZqq d , the F denotes the Fourier transform defined for ϕ P L 1 pT d δ q via the formula

Fpϕqpξq " 1 δ d ż T d δ ϕpxqe ´ix.ξ dx, ξ P 2π δ Z d ,
and its inverse F ´1 is defined for φ P 1 p2πδ ´1Z d q through the formula

F ´1pφqpxq " ÿ kP 2π d Z d
φpkqe `ix.k , x P T d δ .

Figure 2 . 5 Figure 3 .

 253 Figure 2. Solutions |ϕ n | 2 of the generalized relaxation method (55) applied to equation (50) at time T " 30 for α 2 " ´0.5, µ " 0.5 and µ " 2.5.

Figure 4 .Figure 5 .

 45 Figure 4. Evolution of the relative energy error (57) for the generalized relaxation method (55) applied to equation (50)

2 x1 x 1 2 `

 2 (a) Evolution of |ϕpt, x1, x2q| 2 (b) Evolution of |ϕpt, x1, x2q| 2 , slice at x2 " 0 (c) Final state

Figure 6 .

 6 Figure 6. Evolution of |ϕpt, x 1 , x 2 q| 2 for the generalized relaxation method (55) applied to equation (50) for the two-dimensional case

Figure 7 . 6 .

 76 Figure 7. Evolution of the solution |ϕpt, xq| 2 of Eq. (58) modeling rotating dipolar Bose-Einstein condensate

Figure 8 .

 8 Figure 8. Evolution of the energy error with respect to time for the relaxation method (36) applied to equation (58)

  }R k pδtq} pH s`2 pR d qq 7 ď Qδt 3 . }R k pδtq} pH s`2 pR d qq 7 ď 2bQT δt 2 . Then we obtain}e n } pH s pR d qq 7 ď 2b}e 0 } pH s`2 pR d qq 7 `2bQT δt 2 `4Cδtb }e k } pH s pR d qq 7 .Using a discrete Gronwall Lemma (see section 5 in[START_REF] Holte | Discrete Gronwall lemma[END_REF]), we get }e n } pH s pR d qq 7 ď p2b}e 0 } pH s`2 pR d qq 7 `2bQT δt 2 q expp4CT bq.The original relaxation method applied to (30) would consists in adding the variable Υ to (30) and discretizing the following continuous system as discrete times t n and t n`1{2

	#	Υpt, xq " |ϕpt, xq| 2σ ,
	1 ´Cδtb ě k"0 n´1 1 2 . This implies ř n´1 ÿ k"0 This estimate and (21) prove the result. }e This implies that 2b 4. The generalized relaxation method for general nonlinearities 4.1. Generalized relaxation method for NLS equation. We start by considering the simplified equation (1) with zero potential (V " 0), no convolution operator (U " 0) and without rotation (Ω " 0). Assuming σ P N ‹ , we are therefore dealing with the classical nonlinear Schrödinger equation (30) iB (31) iB t ϕpt, xq " ´1 2 ∆ϕpt, xq `βΥϕpt, xq.

n } pH s pR d qq 7 ď 2b}e 0 } pH s`2 pR d qq 7 `4Cδtb n´1 ÿ k"0 }e k } pH s pR d qq 7 `2b n´1 ÿ k"0 }R k pδtq} pH s`2 pR d qq 7 .

Taylor expansions and the fact that one can differentiate the last two lines of

[START_REF] Catherine | The Nonlinear Schrödinger Equation: Self-Focusing and Wave Collapse[END_REF]

with respect to time show that there exists a constant Q such that for all n and δt with nδt ď T t ϕpt, xq " ´1 2 ∆ϕpt, xq `β|ϕpt, xq| 2σ ϕpt, xq, with initial datum ϕ in .

  `α1 U ˚|ϕpt, xq| 2 `α2 |ϕpt, xq| 4 ˙ϕpt, xq, ϕp0, xq " ϕ in ,

				10 ´4	Classical relaxation Generalized relaxation
						Crank Nicolson
		σ " 2,	E E,δt	10 ´8	Slope 2
			10 ´12
			10 ´4 10 ´16	10 ´3	10 ´2	10	´1
						δt
				10	´4
		σ " 3,	E E,δt	10	´8
			10 ´12
			10 ´4 10 ´16	10 ´3	10 ´2	10	´1
						δt
	investigate the interactions of dark solitons under competing nonlocal cubic and local quintic nonlin-
	earities. They consider a 1D optical beam with an amplitude ϕpt, xq with competing nonlocal cubic
	and local quintic nonlinearities. Their model is given by the following nonlocal nonlinear Schrödinger
	equation
	(50)	B t ϕpt, xq " i	ˆ1 2	B 2 x

Figure 1. E E,δt for relaxation schemes and Crank-Nicolson scheme when σ " 2 and σ " 3 5.2. A nonlocal Schrödinger equation with cubic-quintic nonlinearities. In [13], Chen et al.

  Two dimensional rotating dipolar Bose-Einstein condensate. In this final subsection, we present results of simulations with the relaxation scheme for a rotating Bose-Einstein condensate subject to long range dipole-dipole interaction (DDI). The model is (see Eq. (8.15)-(8.17) in[START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF])

		$ ' ' ' ' &	iB t ϕpt, xq "	ˆ´1 2	∆ `V pxq `β|ϕ| 2 pt, xq `λψpt, xq ´ΩL x3 ˙ϕpt, xq, x P R 2 , t ą 0,
	(58)	' ' '	ψpt, xq "	´3 2	pB nKnK	´n2 3 ∆q	ˆ1 2π|x|	˚|ϕ| 2 ˙pt, xq, x P R 2 , t ą 0,
		' %	ϕp0, xq " ϕ 0 pxq, x P R 2 ,

5.3.

For real-valued functions, symmetry with respect to the origin is equivalent to real-valued Fourier transform.

Starting with ϕ 0 " ϕ in and γ ´1{2 approximating |ϕp´δt{2q| 2 , the generalized relaxation method is

˙ϕn`1 `ϕn 2 .

Like for the Crank-Nicolson scheme and equation [START_REF] Antoine | Modeling and computation of Bose-Einstein condensates: stationary states, nucleation, dynamics, stochasticity[END_REF], the first equation also reads

so we also have the other version of generalized relaxation method

Note that, when σ " 1, the generalized relaxation method (34) reduces to the classical relaxation method [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF]. In contrast to the classical relaxation method [START_REF] Besse | Analyse numérique des systèmes de Davey-Stewartson[END_REF], when σ ě 2, the generalized relaxation method (34) is fully implicit on its first stage, and linearly implicit in its second stage. For small values of σ (σ " 1, 2, 3, 4) the first stage of (34) is polynomial of degree σ and hence we can use explicit formulas for the computation of the solution γ n`1{2 . For example, for quintic nonlinearity and σ " 2, we have explicit solutions to the quadratic equation

For higher values of σ, one can use the following iterative fixed-point procedure, starting with

which one stops when }γ n`1{2,p`1 ´γn`1{2,p } L 2 is below some small tolerance parameter, and for this index p, one sets γ n`1{2 " γ n`1{2,p`1 .

Numerically, this generalized relaxation method has order 2, as we will see in the numerical experiments section 5. However, we do not address this theoretical question in this paper. In the next subsection, we show how one can design a generalized relaxation method similar to (34) in order to treat the cases with non-zero potential V , non-zero convolution operator U or non-zero rotation Ω.

4.2.

A generalized relaxation method for GPE. We propose the following generalization of the relaxation method introduced in [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF], which uses two additional unknowns γ and Υ: starting from ϕ 0 " ϕ in , we initialize Υ ´1{2 and γ ´1{2 with approximations of |ϕp´δt{2q| 2 and compute for n P N, pϕ n`1 , Υ n`1{2 , γ n`1{2 q from pϕ n , Υ n´1{2 , γ n´1{2 q using the formulae

the parameter µ allowing to control the width of the kernel, or in other words, the nonlocality strength. At the limit µ Ñ 0, the kernel U j , j " 1, 2, tends to a Dirac distribution and we recover a local nonlinear model. The equation ( 50) is associated to the energy

We reproduce the numerical experiments presented in [START_REF] Chen | Interactions of nonlocal dark solitons under competing cubic-quintic nonlinearities[END_REF][START_REF] Shen | Vortex solitons under competing nonlocal cubic and local quintic nonlinearities[END_REF] and show the ability of the generalized relaxation method to preserve the energy (39) in the form (56) below. The generalized relaxation method for (50) consists in approximating the system of equations ( 54)

and the numerical scheme reads (55)

with ϕ 0 pxq " ϕ in pxqq and Υ ´1{2 pxq " γ ´1{2 pxq is some second order approximation of ϕp´δt{2, xq.

In our numerical experiments, this approximation is obtained by applying the Crank-Nicolson scheme starting from ϕ 0 on reverse time step ´δt{2. The energy associated to ( 54) is

and we have the conservation property (40). We are first interested in the one-dimensional case with kernel U 1 and we choose a defocusing nonlocal nonlinearity by considering α 1 " ´1. Like in the previous subsection, the space variable is discretized using Fourier spectral approximation and we consider periodic boundary conditions. In order to avoid any interaction between the nonlocal kernel and the boundaries, we take a very large domain, typically x P r´256π, 256πs discretized with J " 2 14 `1 nodes. The time step is δt " 5 ¨10 ´3 and the final time is T " 30. The initial datum is made of two solitons at a relative distance of 2x 0 , where x 0 " 1. Following [START_REF] Chen | Interactions of nonlocal dark solitons under competing cubic-quintic nonlinearities[END_REF], we choose ϕ in pxq " tanhpDpx ´x0 qqtanhpDpx `x0 qq, where D is the positive root of the equation cothpDµq Dµ

We present results in Figure 2 and 3 for a defocusing cubic nonlinearity α 2 " ´0.5 and a focusing one where α 2 " 0.1. Two values for µ are proposed µ " 0.5 and µ " 2.5. The behavior of the defocusingdefocusing case is surprising for "strong" nonlocality µ " 2.5 since eventually the two solitons breathe.

The evolution of the relative energy error (57) ˇˇE rlx pϕ n , γ n´1{2 , Υ n´1{2 q ´Erlx pϕ 0 , γ ´1{2 , Υ ´1{2 q ˇĚ rlx pϕ 0 , γ ´1{2 , Υ ´1{2 q is presented on Figure 4. The energy is clearly very well preserved.


Amar Tamra

Marie Deburghgraeve

David Dubuc

Marie-Pierre Rols

Katia Grenier

Microwave dielectric spectroscopy for single cell irreversible electroporation monitoring

Keywords: Dielectric spectroscopy, Electroporation, Biosensor, Single cell analysis, microwave

This paper demonstrates the efficient use of microwave dielectric spectroscopy for single cell irreversible electroporation monitoring. The experimental results point out a high correlation (R 2 higher than 0.94) with biological gold standard technique based on flow cytometry, whereas microwave approach features several key advantages: marker-less, contactless and in-liquid cell(s) monitoring. The developed microwave and microfluidic-based biosensor reveals an increase in the capacitance and the conductance of single cells under investigation subjected to irreversible electroporation reflecting an increase in cell damages, which is in great correlation with results routinely obtained in biology. More interestingly, the post electroporation damages of cell have been experimentally quantified over time, which makes the microwave dielectric spectroscopy a technique of interest for cell's electroporation researches.

I. INTRODUCTION

Accessing the intracellular domain of viable cells has shown a wide interest in the area of biomedical applications. This access can be mediated by several means, such as electroporation. It enables the insertion of proteins into cell membrane, cell fusion, eradication of bacteria, the introduction of small and large molecules into the cell, as well as tumor ablation [START_REF] Kotnic | Electroporation-based applications in biotechnology[END_REF]. Electroporation EP (also known as electropermeabilization), is a physical method that consists in applying well controlled electric pulses (in amplitude and duration) [START_REF] Teissie | Mechanisms of cell membrane electropermeabilization: a minireview of our present (lack of?) knowledge[END_REF], in order to transiently permeabilize the plasma membrane of cells and thus allowing the entry and exit of molecules of interest [START_REF] Neumann | Permeability changes induced by electric impulses in vesicular membranes[END_REF]. In some cases, this approach can cause permanent cell damages and results in cell death: irreversible electroporation (IRE) [START_REF] Rubinsky | Irreversible electroporation: a new ablation modality--clinical implications[END_REF]. When electric pulses are applied to the cell, the electric potential Ψ of the cell membrane is modified: when the sum of the native transmembrane potential ΔΨo and the induced transmembrane potential ΔΨi exceeds a threshold value (-200 mV), the permeabilization occurs [START_REF] Teissie | An experimental evaluation of the critical potential difference inducing cell membrane electropermeabilization[END_REF]. Despite progress in this field, enlarging our knowledge about the mechanism of cell electroporation and kinetics associated with electrochemotherapy remains a pressing need. The assessment of cell permeability and viability after exposure to electroporation is basically done by using conventional microscopy with fluorescent markers, such as calcein AM and propidium iodide (PI), and conventional viability tests such as MTT (Tetrazolium Reduction Assay): measuring the mitochondrial activity of cells. These methods are time consuming, require important volumes of materials, and more importantly are invasive to the cells.

Aside from these techniques, new bio-physical methods, such as Atomic Force Microscopy (AFM) and electrical spectroscopy, have consequently been developed to enrich the cell's physical mechanisms description. The measurement of dielectric properties of biological samples is undeniably expected as of interest to monitor both membrane and intracellular content of the biological cells [START_REF] Pavlidis | High Frequency Wideband Permittivity Measurements of Biological Substance Using Coplanar Waveguides and Application to Cell Suspensions[END_REF]. When using sub-GHz Electrical Impedance Spectroscopy (EIS), the membrane state (permeabilization) is indeed revealed [START_REF] Dean | Electrical Impedance Spectroscopy Study of Biological Tissue[END_REF]. Operating at microwaves (beyond few GHz) becomes then particularly interesting, because (1) the cell's membrane becomes transparent [START_REF] Dalmay | Development of High Frequency Microfluidic Biosensors for Intracellular Analysis[END_REF], and consequently (2) the electromagnetic fields penetrate deeply into the cell and interact directly with the intracellular content.

Here we report our work on the exploitation of our microwave and microfluidic-based biosensor, operating at the single cell level, to monitor the intracellular effects of irreversible electroporation on cells under investigation. We first describe the micro-device architecture combining microfluidic and electrical functionalities. A second part deals with our results on the microwave monitoring of cells submitted to irreversible electroporation; their strong correlation with data obtained with flow cytometry clearly points out the efficacy of our technique in depicting the kinetic effects of cell's electroporation.

II. MATERIALS AND METHODS

This work is based on a micro-biosensor developed by Chen et al. [START_REF] Chen | Microwave biosensor dedicated to the dielectric spectroscopy of a single alive biological cell in its culture medium[END_REF], which takes benefits from the convergence between (1) microfluidic functionalities in charge of the processing of biological material (pure cultured medium or medium with cytotoxic agents or cells suspension) and ( 2) microwave sensor for marker-less, contact-less and in-liquid single cell monitoring. The micro device used to analyze biological cells is presented in (Fig. 1): It is composed of a coplanar waveguide with a capacitive gap in the center. On top of it is placed a microfluidic channel, with a mechanical trap in its center which will help capturing and probing one single cell. The cell's analysis is based on the interaction between electromagnetic fields and the biological sample, which reflects the dielectric properties of cells and makes it possible to assess its state in real time.

In this study, we use THP1 cell line (ATCC#TIB-202): a cell line derived from the peripheral blood of an acute monocytic leukemia patient. The traditional cell culture medium is the Roswell Park Memorial Institute one (RPMI) with 10% of Fetal Bovine Serum (FBS). Electric pulses are delivered to the cells using a Betatech pulse generator with two stainless steel electrodes with a distance of 4 mm. The used electrical parameters are 8 pulses of 1.5 kV/cm at pulse duration of 100 µs with a frequency of 1 Hz. The choice of these parameters was made based on the permeability of cells after EP (Fig. 2) assessed by PI intake: propidium iodide is a small molecule that becomes fluorescent when bound to cellular RNA or DNA. Hence, it only marks cells with permanently or transiently compromised membrane integrity [START_REF] Sukhorukov | Surviving high-intensity field pulses: strategies for improving robustness and performance of electrotransfection and electrofusion[END_REF].

We wanted to evaluate the effect of irreversible electroporation on THP1 cells by studying their response through time. The experimental protocol is composed of 2 major steps. First, cells were pulsed with the generator and reincubated at 37°C and 5% CO2 for 3, 5 and 7 hours after EP. At the end of the incubation time, cells were injected in the RF microdevice or in the flow cytometer in order to compare both results. Fig. 2: Evaluation of THP1 permeability post EP by PI intake.

III. MICROWAVE MEASUREMENTS

A. Irreversible electroporation identification

The microwave measurements are performed from 40 MHz to 40 GHz. The culture medium is first measured alone and considered as our reference medium. When the cells are injected into the microfluidic channel, one cell is captured due to the mechanical trap, directly on top of the sensing zone. Two readouts [START_REF] Chen | Microwave biosensor dedicated to the dielectric spectroscopy of a single alive biological cell in its culture medium[END_REF], the effective capacitance C and conductance G of the sensing volume, are extracted with a specific procedure [START_REF] Chen | Microwave Dielectric Spectroscopy of a Single Biological Cell with Improved Sensitivity up to 40 GHz[END_REF], and compared to the reference medium by calculating contrasts as indicated below:

ΔC = C cell -C reference medium , ΔG = G cell -G reference medium (1)
Fig. 3 shows the spectrums of contrast of capacitance and conductance obtained several hours post EP. As shown in Fig. 3, both contrasts decrease (in absolute values) after exposure to electric pulses reflecting cell death due to irreversible electroporation. This may be explained by the fact that permeable cells have greater exchange with the reference medium, thus the contrast between these two samples is less important. This study allowed us to establish a threshold value between alive and affected cells at the middle of the measured dynamic: ΔC thr =-0.38 fF at 5GHz and ΔG thr =-35 µS at 30GHz. The frequencies on both parameters have been selected to obtain the maximum dynamic as suggested in Fig. 3.

B. Monitoring of irreversible electroporation kinetics

After depicting the different electric signatures of a living and an irreversibly electroporated cell, the next step was to demonstrate the ability of such a technique to monitor the kinetic effect of the IRE. As explained in the materials and methods, the cells were electroporated and studied 3, 5 or 7 hours after EP. At least 10 measurements have been performed for each condition. Once the measurements were done, we proceeded to a cell sorting and a classification based on the thresholds established earlier.

IV. DISCUSSION

The potential of the label-free microwave technique to detect cell death after IRE was evaluated by comparing the results obtained with the well-established flow cytometry technique using PI dye to reveal cell death. The bioexperiments described in part II were conducted with both microwave and flow cytometry techniques. Fig. 5 represents the percentage of dead cells post EP obtained with flow cytometry and microwave measurements (based on both contrasts of capacitance and conductance). A statistical analysis was performed to correlate experiments and it points out strong correlation coefficient R 2 between ∆C/flow cytometry (R 2 =0.94) and ∆G/flow cytometry (R 2 =0.97). These correlations clearly reveal the efficacy and robustness of our technique in depicting the kinetic effects of cell's electroporation.

V. CONCLUSION

In this paper, we experimentally demonstrate that the microwave dielectric spectroscopy of single cells significantly reveals the cell's damages after irreversible electroporation and is able to monitor these damages over time.

Experimental investigations driven by biologists indeed point out that microwave readouts feature a high correlation (R 2 of 0.94 on ∆C and 0.97 on ∆G) with biological gold standard technique based on flow cytometry.

These results convince biologists that microwave dielectric spectroscopy may enrich the technique's arsenal for cell's electroporation researches. The technique could contribute to provide more insights on the effect of electroporation and electrochemotherapy on cancer cells.

Fig. 1 :

 1 Fig.1: Schematic top view of the biosensor composed of a coplanar waveguide with a microfluidic channel placed on top; the cell is blocked due to the mechanical trap.

Fig. 3 :

 3 Fig. 3: Capacitive and conductive contrasts of a living cell and an affected cell several hours after electroporation.

Fig. 4 :

 4 Fig.4: Mean values of the contrast of capacitance obtained different hours after EP. Percentages of viable and affected cells are presented.

Fig. 4

 4 Fig.4presents the mean values of the contrasts obtained after classification as well as the different percentages of cells: viable or affected by the IRE. Compared to the non-treated cells, which feature a viability level (number of viable cells divided by the total cells number) around 73%, IRE induce a decrease of the viability down to 50% after 3 and 5 h post IRE and down to 40% beyond 7 hours after the electrical stress. As one can notice, the mean value of contrasts (same results are obtained on the conductance) is constant after IRE. Only the number of cells in each state varies over time.

Fig. 5 :

 5 Fig.5: Correlation between results obtained with microwave measurements (contrasts of capacitance and conductance) and Flow cytometry.

Affected cells

Living cells