Xavier Antoine
email: xavier.antoine@iecn.u-nancy.fr

Christophe Besse
email: christophe.besse@math.univ-lille1.fr

AND Pauline Klein
email: pauline.klein@iecn.u-nancy.fr

ABSORBING BOUNDARY CONDITIONS FOR GENERAL NONLINEAR SCHR ÖDINGER EQUATIONS *

Keywords: absorbing boundary conditions, pseudodifferential operators, nonlinear Schrödinger equation with potential, stable semi-discrete schemes, fixed point algorithm, relaxation scheme AMS subject classifications. 35Q41, 35Q55, 47G30, 26A33, 65M12, 65M60

This paper addresses the construction of different families of absorbing boundary conditions for the one-and two-dimensional Schrödinger equation with a general variable nonlinear potential. Various semi-discrete time schemes are built for the associated initial boundary value problems. Finally, some numerical simulations give a comparison of the various absorbing boundary conditions and associated schemes to analyze their accuracy and efficiency.

1. Introduction. Schrödinger equations have many applications in physics (see e.g. [START_REF] Sulem | The Nonlinear Schrödinger Equation, Self-Focusing and Wave Collapse[END_REF][START_REF] Tappert | The Parabolic Approximation Method[END_REF][START_REF] Levy | Parabolic equation methods for electromagnetic wave propagation[END_REF][START_REF] Hasimoto | A soliton on a vortex filament[END_REF][START_REF] Peskin | Mechanical equilibrium determines the fractal fiber architecture of aortic heart leaflets[END_REF][START_REF] Bao | Numerical solution of the gross-pitaevskii equation for bose-einstein condensation[END_REF]). They can include repulsive or attractive variable potentials and nonlinear versions are also often met in practice. All these situations represent challenging applications where numerical methods are of upmost importance for predicting the system behavior. These problems are physically set in an unbounded domain which requires the introduction of a fictitious boundary for an efficient numerical simulation. On this boundary must be set an admissible boundary condition such that the restriction of the exact solution to the initial problem coincides with the solution to the bounded problem. It is generally extremely complicate and even often impossible to obtain the exact boundary condition which is given through the Dirichlet-to-Neumann (DtN) operator (see [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF] for a complete review). This DtN map can be derived in the special case [START_REF] Zheng | Exact nonreflecting boundary conditions for one-dimensional cubic nonlinear schrödinger equations[END_REF] of one-dimensional integrable systems following the works by Fokas and his co-authors [START_REF] Boutet De Monvel | Analysis of the global relation for the nonlinear schrödinger equation on the half-line[END_REF]. However, in most situations, it is impossible to obtain the exact DtN map explicitly and the associated boundary condition which is called Transparent Boundary Condition (TBC). Therefore, approximations of this operator are generally considered. They lead to artificial or Absorbing Boundary Conditions (ABCs).

The aim of this article is to provide some constructive processes to obtain robust and accurate approximations of the DtN maps for nonlinear one-and two-dimensional Schrödinger equations subject to potentials

i∂ t ψ + ∆ψ + g(x, t, ψ) ψ = 0, (x, t) ∈ R d ×]0; T ], d = 1, 2, ψ(x, 0) = ψ I (x), x ∈ R d , (1.1)
where ψ I is the initial datum and the function g is a combination of a nonlinearity β(ψ) and a potential V (x, t). The maximal time of computation is denoted by T .

In the sequel of the paper, we assume that the initial datum ψ I is zero outside the computational domain Ω. For obvious reasons linked to the numerical solution of such problems, we are led to truncate R d to a spatial bounded computational domain Ω with a fictitious boundary Σ := ∂Ω. Let us introduce the space-time domains Ω T = Ω × [0; T ] and Σ T = Σ × [0; T ]. Considering the fictitious boundary Σ, we are now led to solve the problem      i∂ t ψ + ∆ψ + g(x, t, ψ) ψ = 0, (x, t) ∈ Ω T , ψ(x, 0) = ψ I (x),

x ∈ Ω, B.C.(x, t, ∂ n ψ, ψ), (x, t) ∈ Σ T , (1.2)
where B.C. denotes a Boundary Condition function on Σ T . Many works have been devoted to the derivation of ABCs for linear Schrödinger equations [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF], that is for β(ψ) = 0. The contributions to the nonlinear case are much more recent. Some developments have been made for the one-dimensional NonLinear Schrödinger (NLS) equation, mainly based on pseudo-and paradifferential approaches. In [START_REF] Antoine | Artificial boundary conditions for one-dimensional cubic nonlinear Schrödinger equations[END_REF], a direct and gauge change strategies are developed for the one-dimensional cubic NLS equation. In [START_REF]Absorbing boundary conditions for nonlinear Schrödinger equations[END_REF], the case of the one-dimensional NLS equation with a cubic term is considered following a potential (direct) approach and pseudodifferential operators. These approaches are compared in [START_REF] Antoine | Towards accurate artificial boundary conditions for nonlinear pdes through examples, Cubo[END_REF]. Other applications concerns NLS equations with a derivative term u∂ x u which can be treated by a pseudo-or paradifferential approach. In [START_REF] Szeftel | Design of absorbing boundary conditions for Schrödinger equations in R d[END_REF], an artificial boundary condition for the linear equation is used for the 2D NLS equation. Another point of view has been developed in [START_REF] Zhang | Unified approach to split absorbing boundary conditions for nonlinear Schrödinger equations[END_REF][START_REF]Unified approach to split absorbing boundary conditions for nonlinear Schrödinger equations: Two-dimensional case[END_REF] where the authors propose a splitting approach for the one-and two-dimensional case. Let us also cite [START_REF] Xu | Adaptive absorbing boundary conditions for Schrödinger-type equations: Application to nonlinear and multi-dimensional problems[END_REF] where the authors analyze an adaptive approach for NLS equations. An extension to coupled NLS equations is given in [START_REF] Zhou | Numerical solution to coupled nonlinear Schrödinger equations on unbounded domains[END_REF]. An alternative solution is related to PMLs. In [START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF], the authors numerically analyze the application of PMLs to the simulation of Bose-Einstein condensates. In [START_REF]A perfectly matched layer approach to the nonlinear Schrödinger wave equations[END_REF], a study of PMLs to one-and two-dimensional linear and nonlinear Schrödinger equations is developed. The results show a quite satisfactory behaviour of the PMLs for nonlinear problems.

The aim of the present paper is to introduce and analyze some new families of ABCs for one-and two-dimensional general nonlinear Schrödinger equations based on pseudodifferential operator approaches. The fictitious boundary Σ is supposed to be smooth and convex for the 2D case. We also propose, for the 2D case, an improved version of the approximation by using a relaxation scheme. Finally, numerical simulations show the accuracy of our methods for the computation of moving soliton solutions. These works extends the developments in [START_REF] Antoine | Artificial boundary conditions for one-dimensional cubic nonlinear Schrödinger equations[END_REF][START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] to very general nonlinear Schrödinger equations subject to potentials. The second section is devoted to the construction of families of ABCs for general one-dimensional nonlinear Schrödinger equations. We present in the third section the different semi-discretizations in view of a solution by fixed point procedure. Section 4 deals with the numerical finite element approximation of the problem. Numerical examples are given in Section 5. The sixth Section concerns the construction of ABCs for the 2D NLS equation. The resulting boundary-value problems are studied in Section 7. Section 8 concerns the discretization of the ABCs by using a fixed point procedure and finite element methods. Stability issues are considered there. Section 9 proposes a relaxation scheme for the problem with ABCs. In Section 10, we consider the application of our ABCs/scheme to the solution of the 2D cubic NLS equation with a soliton initial data. We analyze the behavior of the different schemes in terms of efficiency and accuracy and provides numerical simulations. Finally, the last Section concludes the paper.

2. ABCs for one-dimensional NLS equations. Let us consider the nonlinear Schrödinger equation

i∂ t ψ + ∂ 2 x ψ + α(ψ)∂ x ψ + β(ψ)ψ = 0, (x, t) ∈ Ω T , ψ(x, 0) = ψ I (x), x ∈ Ω. (2.1)
The derivation of a TBC can be extremely difficult and arise only in rare and non general cases. The only known TBCs for (2.1) appear in the paper of Zheng [START_REF] Zheng | Exact nonreflecting boundary conditions for one-dimensional cubic nonlinear schrödinger equations[END_REF] for α = 0 and β = ±2|ψ| 2 . For this choice of parameters, the nonlinear Schrödinger equation is the cubic nonlinear Schrödinger equation which is known to be integrable. The derivation of TBCs relies on inverse scattering technique. Therefore, they are specifically built for these special cases and cannot be extended to more general situations. Another work of Zheng dealing with PML is also of interest [START_REF]A perfectly matched layer approach to the nonlinear Schrödinger wave equations[END_REF].

In the context of ABCs, contributions can be found in the papers [START_REF] Antoine | Artificial boundary conditions for one-dimensional cubic nonlinear Schrödinger equations[END_REF][START_REF]Absorbing boundary conditions for nonlinear Schrödinger equations[END_REF][START_REF] Antoine | Towards accurate artificial boundary conditions for nonlinear pdes through examples, Cubo[END_REF]. The derivation of such ABCs relies again on tools of mathematical analysis such as pseudodifferential calculus [START_REF] Taylor | Pseudodifferential operators[END_REF] or paradifferential calculus [START_REF] Bony | Calcul symbolique et propagation des singularités pour les equations aux dérivées partielles non linéaires[END_REF]. The two techniques have their own pros and cons. To get ABCs for a general nonlinear Schrödinger equation with α = 0, we extend here the ABCs derived in [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] for the potential case V (x, t) to a nonlinear term β(ψ) by using the formal substitution: V (x, t) → β(ψ)(x, t). For the sake of clarity, we only present below the nonlinear version of ABC N j , j = 0, 1 (see [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]) that are denoted by NLABC N j , for α = 0, where N denotes the order of the ABC. These strategies are called gauge change strategy for j = 0 and direct strategy for j = 1.

The nonlinear boundary conditions NLABC N 0 are

∂ n ψ + Λ N 0 ψ = 0, on Σ T , (2.2)
with (NLABC 2 0 )

Λ 2 0 (x, t, ∂ t ) ψ = e -iπ/4 e iB(x,t,ψ) ∂ 1/2 t e -iB(x,t,ψ) ψ (2.3) and (NLABC 4 0 ) Λ 4 0 (x, t, ∂ t ) ψ = Λ 2 0 (x, t, ∂ t ) ψ + i ∂ n β(ψ) 4 e iB(x,t,ψ) I t e -iB(x,t,ψ) ψ , (2.4)
setting: B(x, t, ψ) = t 0 β(ψ)(x, s) ds. The integral operators I α/2 t
are defined by the relation

I α/2 t f (t) = 1 Γ(α/2) t 0 (t -s) α/2-1 f (s) ds, for α ∈ N, (2.5)
and the fractional differential operator ∂

1/2 t is given by ∂ 1/2 t f (t) = 1 √ π ∂ t t 0 f (s) √ t -s ds. (2.6)
Used as boundary conditions for Eq. (2.1) with α(ψ) = 0, the NLABC 2 0 boundary conditions lead to well posed problem. However, we can not hope to get a similar result for the NLABC 4 0 conditions. Indeed, when considered with a single variable potential, proposition 4 (Section 2.4, pp.319) in [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] states that ABC 4 0 gives a well posedness result if the potential satisfies an hypothesis on the sign of its normal derivative which has to be constant in time. We can not assume such a property for a general nonlinearity β(ψ) that involves the solution ψ itself. Following the proof for the potential case, we therefore have the following result.

Proposition 2.1. Let ψ I ∈ L 2 (Ω) be a compactly supported initial datum such that Supp(ψ I ) ⊂ Ω. Let us denote by ψ a solution of the initial boundary value problem

     i∂ t ψ + ∂ 2 x ψ + β(ψ)ψ = 0, in Ω T , ∂ n ψ + Λ 2 0 ψ = 0, on Σ T , ψ(x, 0) = ψ I (x), ∀x ∈ Ω.
(2.7)

Then, ψ fullfils the following energy bound

∀t > 0, ||ψ(t)|| L 2 (Ω) ≤ ||ψ I || L 2 (Ω) , (2.8)
In [START_REF] Antoine | Towards accurate artificial boundary conditions for nonlinear pdes through examples, Cubo[END_REF], the second (direct) strategy was applied with principal symbol √ τ and has led to the nonlinear boundary conditions NLABC N

1 ∂ n ψ + Λ N 1 ψ = 0, on Σ T , (2.9)
with (NLABC 2 1 )

Λ 2 1 (x, t, ∂ t ) ψ = e -iπ/4 ∂ 1/2 t ψ (2.10)
and (NLABC 4 1 )

Λ 4 1 (x, t, ∂ t ) ψ = Λ 2 1 (x, t, ∂ t ) ψ -e -iπ/4 β(ψ) 2 I 1/2 t ψ. (2.11)
Nothing limits the choice of the principal symbols and we therefore are allowed to change it for NLABC N 1 . We present the modified nonlinear boundary conditions when the principal symbol is chosen to be τ -β(ψ) (see [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]). In this case, we get the ABCs NLABC N 1m (for the direct strategy) given by

∂ n ψ + Λ N 1m ψ = 0, on Σ T , (2.12)
with (NLABC 2 1m )

Λ 2 1m (x, t, ∂ t ) ψ = -i i∂ t + β(ψ)ψ (2.13) and (NLABC 4 1m ) Λ 4 1m (x, t, ∂ t ) ψ = Λ 2 1m (x, t, ∂ t ) ψ + ∂ n β(ψ) 4 (i∂ t + β(ψ)) -1 ψ. (2.14)
In terms of pseudodifferential conditions, the new nonlinear ABCs NLABC N 1m are more precise than NLABC N 1 . These nonlinear ABCs can also be extended to a very general nonlinearity including variable potentials. It is sufficient to replace the nonlinear term β(ψ) by the function g which includes the potential and to define it by g(x, t, ψ) = β(ψ) + V (x, t).

(2.15)

3. Semi-discretization schemes. The aim of this section is to present semidiscrete time schemes for the initial value problem

i∂ t ψ + ∂ 2 x ψ + g(x, t, ψ)ψ = 0, (x, t) ∈ Ω T , ∂ n ψ + Λ N 0,1m ψ = 0, on Σ T , ψ(x, 0) = ψ I (x), x ∈ Ω, (3.1)
for a maximal time of computation T . We consider an interior Crank-Nicolson scheme for the time discretization of system (3.1). The interval [0; T ] is uniformly discretized using N T intervals. Let ∆t = T /N T be the time step and let us set t n = n∆t. Furthermore, ψ n stands for an approximation of ψ(t n ) and V n = V (x, t n ). A suitable semi-discrete approximation for the system (3.1) is the scheme developed by Durán and Sanz-Serna for the cubic nonlinear Schrödinger equation [START_REF] Durán | The numerical integration of relative equilibrium solutions. the nonlinear Schrödinger equation[END_REF]. One of the most interesting property is the symmetry of this scheme. Adapted to (3.1), it reads

i ψ n+1 -ψ n ∆t + ∂ 2 x ψ n+1 + ψ n 2 + g(x, t n+1/2 , ψ n+1 + ψ n 2 ) ψ n+1 + ψ n 2 = 0, (3.2)
for n = 0, • • • , N T -1.
Let us remark that, for implementation issues, it is useful to set:

ϕ n+1 = ψ n+1 +ψ n 2 = ψ n+1/2
, with ψ -1 = 0 and ψ 0 = ψ I . Then, the time scheme

(3.2) reads 2i ϕ n+1 ∆t + ∂ 2 x ϕ n+1 + g(x, t n+1/2 , ϕ n+1 )ϕ n+1 = 2i ψ n ∆t . (3.3)
We propose below one approximation for each kind of nonlinear ABCs. The approach for NLABC N 0 is based on semi-discrete convolutions for the fractional operators involved in (2.2). Considering NLABC N 1m , we present a scheme based on the approximation of the operators i∂ t + g(x, t, ψ) and (i∂ t + g(x, t, ψ)) -1 through the solution of auxiliary differential equations which can be solved explicitly. t . The second one deals with the nonlinearity. This last issue is shared with NLABC N 1m conditions and we will focus on this question in the section discussing the finite element approximation. In order to emphasize on this point, we need to distinguish the linear and nonlinear contributions taking place in NLABC N 0,1m conditions. We write the normal derivative of ϕ n+1 on Σ T as:

∂ n ϕ n+1 = B L ϕ n+1 + B N L (ϕ n+1
), where B L and B N L respectively stand for the linear and nonlinear contributions.

The approximation of fractional operators can be performed by various ways. We choose in this paper to use the discrete convolutions of [START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF][START_REF] Antoine | Numerical schemes for the simulation of the twodimensional Schrödinger equation using non-reflecting boundary conditions[END_REF][START_REF]Approximation, stability and fast evaluation of exact artificial boundary condition for the one-dimensional heat equation[END_REF]. If {f n } n∈N is a sequence of complex numbers approximating {f (t n )} n∈N , then the approximations of ∂

1/2 t f (t n ), I 1/2
t f (t n ) and I t f (t n ) with respect to the Crank-Nicolson scheme for a time step ∆t are given by the numerical quadrature formulas

∂ 1/2 t f (t n ) ≈ 2 ∆t n k=0 b n-k f k , (3.4)
I 1/2 t f (t n ) ≈ ∆t 2 n k=0 a n-k f k , (3.5)
I t f (t n ) ≈ ∆t 2 n k=0 c n-k f k , (3.6)
where the sequences (

a n ) n∈N , (b n ) n∈N and (c n ) n∈N are such that    (a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , ...) = (1, 1, 1 2 , 1 2 , 3 8 , 3 8 , ...), b k = (-1) k a k , ∀k ≥ 0, (c 0 , c 1 , c 2 , c 3 , ...) = (1, 2, 2, ...). (3.7)
Thanks to these definitions, the discretization of NLABC 4 0 (for a function g) is given on Σ T by

∂ n ϕ n+1 + e -iπ/4 e iG n+1 2 ∆t n+1 k=0 b n+1-k e -iG k ϕ k + i 4 ∂ n g(x, t n+1/2 , ϕ n+1 )e iG n+1 ∆t 2 n+1 k=0 c n+1-k e -iG k ϕ k = 0, (3.8)
with the notation: G n+1 = I t g(x, t n+1/2 , ϕ n+1 ). Let us introduce the following notations

E p = exp(iG p ) = E p-1 exp ∆t 2 g(x, t p-1/2 , ϕ p ) , (3.9)
setting E 0 = 1 with: E p-1 = exp i∆t p-1 l=1 g(x, t l-1/2 , ϕ l ) . Therefore, the approximation (3.8) can be recast as

∂ n ϕ n+1 = B 4 0,L ϕ n+1 + B 4 0,N L (ϕ n+1 ), (3.10)
with: B 4 0,L = -e -iπ/4 2 ∆t and

B 4 0,N L (ϕ n+1 ) = -e -iπ/4 e i ∆t 2 g n+1 E n 2 ∆t n k=0 b n+1-k E k ϕ k - i 4 ∂ n g n+1 ∆t ϕ n+1 2 + e i ∆t 2 g n+1 E n n k=0 E k ϕ k , setting g n+1 = g(x, t n+1/2 , ϕ n+1 ).
The stability property of the Durán Sanz-Serna scheme linked to the nonlinear absorbing boundary conditions can be proved for the second order ABCs. The result is developed in the following proposition and closely follows the proof of the potential case by using the Z-transform.

        2i ϕ n+1 ∆t + ∂ 2 x ϕ n+1 + g(x, t n+1/2 , ϕ n+1 )ϕ n+1 = 2i ψ n ∆t , in Ω ∂ n ϕ n+1 + Λ 2,n+1 0 ϕ n+1 = 0, on Σ, ψ 0 = ψ I , in Ω, (3.11)
for n = 0, ..., N -1, where the semi-discrete operator Λ 2,n+1 0 is defined by

Λ 2,n+1 0 ϕ n+1 = e -iπ/4 e iG n+1 2 ∆t n+1 k=0 b n+1-k e -iG k ϕ k , (3.12)
is L 2 (Ω) stable and the inequality

∀t > 0, ||ψ n || L 2 (Ω) ≤ ||ψ 0 || L 2 (Ω) , (3.13)
holds.

Let us remark that other discretizations of the fractional operators could have been used like for the fast approximations proposed in [START_REF] Jiang | Fast evaluation of nonreflecting boundary conditions for the Schrödinger equation in one dimension[END_REF]. However, we think that estimates like in Proposition 3.1 are out of reach in this case.

Padé based discretizations for NLABC N

1m . While the previous strategy based on discrete convolution operators is accurate and provides a stability result, it may lead to long computational times. For NLABC N 1m , the discretizations of the resulting pseudodifferential operators involved is not easy to obtain. An alternative approach to discrete convolutions (which cannot be applied here) consists in approximating the square-root operator i∂ t + g(x, t, ψ) by using rational functions and auxiliary functions. More specifically, we consider here the m-th order Padé approximants [START_REF] Milinazzo | Rational square-root approximations for parabolic equation algorithms[END_REF]

√ z ≈ R m (z) = α m 0 + m k=1 α m k z z + δ m k = m k=0 α m k - m k=1 α m k δ m k z + δ m k , (3.14)
where the coefficients (α m k ) 0≤k≤m and (δ m k ) 1≤k≤m are given by

α m 0 = 0 , α m k = 1 m cos 2 (2k+1)π 4m , δ m k = tan 2 (2k + 1)π 4m . (3.15)
Formally, i∂ t + g(x, t, ψ) is approximated by

R m (i∂ t + g) = m k=0 α m k - m k=1 α m k δ m k (i∂ t + g + δ m k ) -1 . (3.16)
Applying this process to the equation (2.12), with N = 4, we have the new relation

∂ n ψ -i m k=0 α m k - m k=1 α m k δ m k (i∂ t + g(x, t, ψ) + δ m k ) -1 ψ + 1 4 (i∂ t + g(x, t, ψ)) -1 ∂ n g(x, t, ψ)ψ = 0. (3.17)
To write a suitable form of the equation in view of an efficient numerical treatment, we classically introduce m auxiliary functions v k , for 1 ≤ k ≤ m, and u (see Lindmann [START_REF] Lindmann | Free-space boundary conditions for the time dependent wave equation[END_REF]) solutions of the following ordinary differential equations

(i∂ t + g(x, t, ψ) + δ m k ) φ k = ψ, (3.18)
with the initial condition φ k (x, 0) = 0 and

(i∂ t + g(x, t, ψ))Ψ = ∂ n g(ψ)ψ. (3.19)
The ABC (3.17) at points x = x l , x r therefore becomes

∂ n ψ -i m k=0 α m k ψ - m k=1 α m k δ m k φ k + 1 4 Ψ = 0. (3.20)
The ODEs (3.18) and (3.19) are discretized with a classical Crank-Nicolson scheme which leads to the discretized ABC

∂ n ϕ n+1 = B 4 1m,L ϕ n+1 + B 4 1m,N L (ϕ n+1 ), (3.21) with B 4 1m,L = i m k=0 α m k and B 4 1m,N L (ϕ n+1 ) = -i m k=1 α m k δ m k φ n+1 k +φ n k 2 -1 4 ∂ n g n+1 Ψ n+1 +Ψ n 2
, where g n+1 = g(x, t n+1/2 , ϕ n+1 ) and φ n+1 k is solution to

i ∆t + 1 2 g n+1 + δ m k 2 φ n+1 k = ϕ n+1 + i ∆t - 1 2 g n+1 - δ m k 2 φ n k and Ψ n+1 satisfies: i ∆t + 1 2 g n+1 Ψ n+1 = ϕ n+1 + i ∆t -1 2 g n+1 Ψ n . 4.
Finite element implementation. Finite element methods are now used for approximating the solution to the semi-discrete equations. This method is particularly relevant since we deal with a Dirichlet-to-Neumann formulation of our ABCs. Indeed, the weak formulation of (2.1) is, for any test-function χ belonging to the Sobolev space

H 1 (Ω), 2i ∆t Ω ϕ n+1 χdx + [∂ x ϕ n+1 χ] xr x l - Ω ∂ x ϕ n+1 ∂ x χdx + Ω g(x, t n+1/2 , ϕ n+1 )ϕ n+1 χdx = 2i ∆t Ω ψ n χdx.
Since the approximations of NLABC N 0,1m are given at points x l , x r by means of a generic relation:

∂ n ϕ n+1 = B N L ϕ n+1 + B N N L (ϕ n+1
), it is easy to replace the space derivative in the term [∂ x ϕ n+1 χ] xr x l . The problem being globally nonlinear, the approximation is performed by means of a fixed point procedure. At a given time t n , we set ζ 0 = ϕ n , and we compute the solution ϕ n+1 as the limit of the iterative scheme

2i ∆t Ω ζ p+1 χdx + [∂ x ζ p+1 χ] xr x l - Ω ∂ x ζ p+1 ∂ x χdx = - Ω g(x, t s+1/2 , ζ p )ζ p χdx + 2i ∆t Ω ψ n χdx, (4.1)
with the boundary condition

∂ n ζ p+1 = B N j,L ζ p+1 + B N j,N L (ζ p ), j = 0, 1m. (4.2)
The discretization of the variational system (4.1)-(4.2) is done with a linear finite element method. The mesh is uniform, the mesh size being h for n h grid points on Ω.

Let us denote by S the stiffness matrix and by M the mass matrix of the linear finite element procedure. The boundary matrix term associated with B N j,L is represented by B N j,L . We denote by b N,p j the discrete version of the right-hand side of (4.1) including the contribution of B N j,N L ζ p . At every iteration of the fixed point algorithm, we have to solve the linear system: 2i ∆t M -S + B N j,L ζ p+1 = b N,p j , j = 0, 1m, with the initialization ζ 0 = ψ n . The bold notation ψ n defines the vector of nodal values in the linear finite element basis. The stopping criterion is given by ζ p+1 -ζ p L 2 (Ω) < ε, for an a priori fixed tolerance ε. Some extra definitions for functions φ n+1 k and Ψ n+1 are necessary for NLABC N 1m but lead to a similar algorithm.

5.

Numerical examples for a mixed potential. We present in this section two numerical examples of a mixed potential: the case of a linear combination between a nonlinearity and a space dependent potential. The first example is presented for the mixed potential g(x, t, ψ) = 5x 2 + |u| 2 . The initial data is taken as the value of the soliton for a quadratic media at t = 0:

ψ I (x) = 2a q sech( √ ax) exp i( cx 2 + θ 0 + a) .
The numerical parameters are q = 1, a = 2, c = 15 and θ 0 = 0. The velocity c is chosen sufficiently large to generate a high frequency solution. The boundary points are x l = -10 and x r = 10 and the computational domain Ω = (-10, 10). The discretization parameters are ∆t = 10 -3 , n h = 4000 (h = 5 • 10 -3 ) and m = 50. The parameter ε involved in the fixed point procedure is fixed to 10 -8 . In order to show the effectiveness of our ABCs, we plot the quantity log 10 (|u(x, t)|). This choice guarantees a fair presentation of the quality between different ABCs. We plot the reference solution on Fig. 5.1 and we present the solution generated by NLABC N 0 boundary conditions. The second-order ABC produces magnitude of reflection of the order of 10 -2.5 as it can be seen on Fig. 5.2. It is improved a little by using the fourthorder ABC (see Fig. 5.3). This is generally the case for other kinds of nonlinearities.

Intensive simulations show that the same accuracy is obtained for the Padé-type ABCs NLABC N j , j = 0, 1m, N = 2, 4, but at the price of a noticeable increase of the number of iterations of the fixed point procedure. Finally, a comparison with PMLs (for optimal tuning parameters) show that our ABCs produce a better accuracy at no computational cost.

The second example allows to validate the approach for a repulsive interaction which models expansion of Bose-Einstein condensate [START_REF] Xu | Absorbing boundary conditions for nonlinear schr à ¶dinger equations[END_REF]. In this specific case, the mixed potential is g(x, t, ψ) = -e -0.5 * (x-15) 2 -2|u| 2 and the initial datum is ψ I (x) = e -0.1 * (x-15) 2 . The solution is computed on the computational domain Ω = (0, 30). The discretization parameters are ∆t = 10 -2 , n h = 2000. The quantity log 10 (|u(x, t)|) is plotted on Fig. 5.4-5.5. There is no reflection back, confirming hence the accuracy of the proposed ABCs.

6. ABCs for the two-dimensional NLS equation. In this part, we consider the generalization of the previous results to the two-dimensional case with a nonlinear potential, defined by: g(x, t, ψ) = V (x, t) + β(ψ). Therefore, we consider the twodimensional NonLinear Schrödinger (2D-NLS) system

i∂ t ψ + ∆ψ + g(x, t, ψ)ψ = 0, x ∈ R 2 , t > 0, ψ(x, 0) = ψ I (x), x ∈ R 2 . (6.1)
The function β is a real-valued positive function in practice and will be always a function of |ψ|. We consider the following geometrical situation: Ω is a smooth convex bounded domain of R 2 with boundary Σ. The initial data ψ I is again supposed to be compactly supported in Ω. To get some nonlinear ABCs, we apply the same approach as in the one-dimensional: the potential V is formally replaced by the nonlinear potential g = V + β(u).

Taylor-based ABCs.

A first class of ABCs concerns Taylor expansion based boundary conditions. Refeering to [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF], these ABCs are labelled ABC M 2,T coming from the "direct approach". As in the one-dimensional case, we cannot symmetrize the nonlinear ABCs according to (∂ n g)I t since for a general nonlinear function g we cannot a priori confirm if sg(∂ n g) is constant or not on Σ. However, it is more realistic to consider that the sign of g = V + β(u) is constant over Σ for example if V and β are either positive or negative. In this case, we get the fourth-order nonlinear

ABC ∂ n ψ + e -iπ/4 ∂ 1/2 t ψ + κ 2 ψ -e iπ/4 κ 2 8 + ∆ Σ 2 I 1/2 t ψ -e -iπ/4 sg(V ) 2 |V | I 1/2 t |V | ψ -e -iπ/4 sg(β(ψ)) 2 |β(ψ)| I 1/2 t |β(ψ)| ψ + i ∂ s (κ∂ s ) 2 + κ 3 + ∆ Σ κ 8 - ∂ n V 4 - ∂ n β(ψ) 4 I t ψ = 0, (6.2)
that we denote by NLABC 4 2,T . This extends the ABCs (2.9)-(2.11) to the 2D case. Lower order ABCs, denoted by NLABC M 2,T , for M = 2 and M = 3, can then be easily deduced. In the above relations, s is the curvilinear abscissa along Σ in the counterclockwise direction, κ = κ(s) is the curvature of the surface and ∆ Σ = ∂ 2 s is the Laplace-Beltrami operator over Σ. Finally, n is the outwardly directed unit normal vector to Ω.

In [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF], a second family labelled ABC M 1,T is developed. The situation is more delicate for these boundary conditions. Essentially, the ABCs are complex to code for non radial linear potentials. In our case, even if β is a radial function, then the nonlinear potential V + β(ψ) is not for example if the initial data ψ I is given with a positive wavenumber. Indeed, in this case, ψ is not symmetrical and so β(ψ) too. For this reason, we only study here the adaptation of the second-order ABC, ABC 2 1,T , which is

∂ n ψ + e -iπ/4 e iB ∂ 1/2 t e -iB ψ + κ 2 ψ = 0. (6.3)
We denote (6.3) by NLABC 2 1,T and define the nonlinear phase function B by

B(x, t, ψ) = t 0 β(ψ)(x, σ)dσ. (6.4)
The extension to a general mixed nonlinear function g is direct.

6.2. Padé-based ABCs. ABCs based on Padé approximants can also be designed. Here, we restrict our study to the ABCs called ABC M 2,P obtained by the "direct" method. In [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF], we developed two families of ABCs, the "direct approach" and "gauge change" ABCs, which produce the same accuracy. However, we also show that the "direct approach" ABCs are much less affected by the numerical dispersion than the "gauge change" ABCs, meaning that we get more accurate results with coarser finite element meshes. Furthermore, the adaptation of the "gauge change" approach is not well-adapted to the case of non radial potentials which is the case here. For all these reasons, we focus on the Padé-based family of ABCs ABC M 2,P , and most particularly on the second-order ABC. The boundary condition NLABC 2 2,P , adapted to a nonlinearity, writes down on Σ as

∂ n ψ -i i∂ t + ∆ Σ + V + β(ψ) ψ + κ 2 ψ - κ 2 (i∂ t + ∆ Σ + V + β(ψ)) -1 ∆ Σ ψ = 0. (6.5)
Next we introduce the Padé approximants of order m to approximate the square-root operator

∂ n ψ -i m k=0 a m k ψ + i m k=1 a m k d m k (i∂ t + ∆ Σ + V + β(ψ) + d m k ) -1 ψ + κ 2 ψ - κ 2 (i∂ t + ∆ Σ + V + β(ψ)) -1 ∆ Σ ψ = 0, (6.6)
and then the auxiliary functions (ϕ k ) 1≤k≤m and Ψ defined respectively by

(i∂ t + ∆ Σ + V + β(ψ) + d m k ) ϕ k = ψ on Σ, (6.7)
for 1 ≤ k ≤ m, and

(i∂ t + ∆ Σ + V + β(ψ)) Ψ = ∆ Σ ψ on Σ, (6.8)
with the initial condition ϕ k (x, 0) = 0 for 1 ≤ k ≤ m and ψ(x, 0) = 0, for x ∈ Σ. Finally, the boundary condition (6.6) writes

∂ n ψ -i m k=0 a m k ψ + i m k=1 a m k d m k ϕ k + κ 2 ψ - κ 2 Ψ = 0. ( 6.9)
System (6.9)-(6.7)-(6.8) is the representation of the condition NLABC 2 2,P by m Padé approximants. 2 1,T and NLABC M 2,T . We consider here the second-order ABCs NLABC 2 1,T and NLABC 2 2,T , as well as the symmetrized third-order boundary condition NLABC 3 2,T . Indeed, among the conditions with gauge change, we only consider the second-order ABC. For the direct method and for the potential case, the boundary condition was symmetrized with respect to the term ∂ n V I t , and then by using that ∂ n V has a constant sign with respect to time on Σ. However, we cannot confirm here anything about the sign of ∂ n β(ψ) on Σ implying that we exclude the fourth-order boundary condition. Then, we get the following result by adapting the proof developed for the potential case [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF] to the nonlinear problem.

Study of the boundary value problem in the case of NLABC

Theorem 7.1. Let ψ I ∈ L 2 (Ω) be an initial data with compact support in Ω. Let V ∈ C ∞ (R 2 × R + , R) be a real-valued potential and β a real positive function. We denote by ψ ∈ C([0; T ], H 1 (Ω)) a solution of the problem with initial data and with the ABC NLABC

2 2,T        i∂ t ψ + ∆ψ + V ψ + β(ψ)ψ = 0, in Ω T , ∂ n ψ + e -iπ/4 ∂ 1/2 t ψ + κ 2 ψ = 0, on Σ T , ψ(x, 0) = ψ I (x),
on Ω, (7.1)

or the ABC NLABC 2 1,T        i∂ t ψ + ∆ψ + V ψ + β(ψ)ψ = 0, in Ω T , ∂ n ψ + e -iπ/4 e iG ∂ 1/2 t e -iG ψ + κ 2 ψ = 0, on Σ T , ψ(x, 0) = ψ I (x), on Ω, (7.2)
where the phase function G is given by (6.4) for g instead of β. Then, ψ satisfies the energy inequality

∀t > 0, ψ(t) L 2 (Ω) ≤ ψ I L 2 (Ω) . (7.3)
Moreover, if we consider the system with the symmetrized ABC NLABC 3 2,T

∂ n ψ+e -iπ/4 ∂ 1/2 t ψ+ κ 2 ψ-e iπ/4 κ 2 8 + ∆ Σ 2 I 1/2 t ψ-e -iπ/4 sg(V ) 2 |V | I 1/2 t |V | ψ -e -iπ/4 sg(β(ψ)) 2 |β(ψ)| I 1/2 t |β(ψ)| ψ = 0, (7.4)
and assume that V is a positive potential, and β is a positive real-valued function, then the inequality (7.3) is still valid.

8. Semi-discretization of the ABCs and fixed point.

(ϕ n+1 -ψ n ) ∆t + ∆ϕ n+1 + W n+1 ϕ n+1 + β(ϕ n+1 )ϕ n+1 = 0, with the usual notation ϕ n+1 = ψ n+1/2 = ψ n+1 +ψ n 2 and W n+1 = V n+1 +V n 2
. Then, we write the variational formulation for a test-

function χ ∈ H 1 (Ω) 2i ∆t Ω ϕ n+1 χ dΩ + Σ ∂ n ϕ n+1 χ dΣ - Ω ∇ϕ n+1 • ∇χ dΩ + Ω W n+1 ϕ n+1 χ dΩ + Ω β(ϕ n+1 )ϕ n+1 χ dΩ = 2i ∆t Ω ψ n χ dΩ. (8.1)
To apply a fixed point, we write

∂ n ϕ n+1 = -B n+1 L ϕ n+1 -B n+1 N L (ϕ n+1 ) (8.2)
where B n+1 L is a coefficient, and B n+1 N L is a function of ϕ n+1 which only contents the nonlinear terms in ϕ n+1 . Then, going to the discretization of the equation leads to the application of the fixed point algorithm as the successive solution to

2i ∆t M -B n+1 L -S + M W n+1 ζ p+1 = 2i ∆t Mψ n -M β(ζ p ) ζ p + b n + B n+1 N L (ζ p ) (8.3)
for p ≥ 0, where ζ is initialized by ζ 0 = ϕ n . When the convergence of the fixed point is reached, we get ϕ n+1 = ζ p+1 . In the above notations, M F designates the matrix representation of the generalized mass matrix associated to a function F and the linear FEM.

8.2. Numerical treatment of the ABCs NLABC M 1,T and NLABC M 2,T . We consider the time discretization of the ABCs related to the Taylor expansion approach and precise some of the numerical implementation issues in a FEM background.

Let us consider the discretization of the ABC NLABC 4 2,T obtained by the direct method. We have

∂ n ϕ n+1 + e -iπ/4 2 ∆t b n+1 + κ 2 ϕ n+1 -e iπ/4 ∆t 2 κ 2 2 a n+1 0 + 1 2 a n+1 2 + 1 2 sg(W n+1 ) |W n+1 | a n+1 V + 1 2 sg(β(ϕ n+1 )) |β(ϕ n+1 )| a n+1 β + i ∆t 2 1 2 ∂ s (κd n+1 1 ) + κ 3 + ∆ Σ κ 8 d n+1 0 - ∂ n W n+1 4 d n+1 0 - ∂ n β(ϕ n+1 ) 4 d n+1 0 = 0, (8.4)
where the coefficients b n+1 , a n+1

a n+1 γ = α n+1 ⋆ ∂ γ s ϕ n+1 , γ ∈ {0, 1, 2}, (8.6)
a n+1 V = α n+1 ⋆ |W n+1 | ϕ n+1 , (8.7)
d n+1 γ = δ n+1 ⋆ ∂ γ s ϕ n+1 , γ ∈ {0, 1}, (8.8)
d n+1 V = δ n+1 ⋆ |∂ n W n+1 | ϕ n+1 , (8.9)
and where

a n+1 β = α n+1 ⋆ |β(ϕ n+1 )| ϕ n+1 . (8.10)
Let us set

µ 1 = e -iπ/4 2 ∆t , µ 2 = 1 2 , µ 3 = -e iπ/4 ∆t 2 , µ 4 = i ∆t 2 . (8.11)
To get the weak variational formulation, we explicit in (8.4) the discrete convolutions and we localize the terms ϕ n+1 by using the notations µ 1 , µ 2 , µ 3 , µ 4 introduced in (8.11)

-∂ n ϕ n+1 = µ 1 + µ 2 κ + µ 3 κ 2 2 + 1 2 ∂ 2 s + W n+1 2 + β(ϕ n+1 ) 2 +µ 4 1 2 ∂ s (κ∂ s ) + κ 3 + ∆ Σ κ 8 - ∂ n W n+1 4 - ∂ n β(ϕ n+1 ) 4 ϕ n+1 + µ 1 n k=0 β n+1-k ϕ k + µ 3 n k=0 α n+1-k κ 2 2 + 1 2 ∂ 2 s ϕ k + µ 3 n k=0 α n+1-k sg(W n+1 ) |W n+1 | |W k | ϕ k + µ 3 n k=0 α n+1-k sg(β(ϕ n+1 )) |β(ϕ n+1 )| |β(ϕ k )| ϕ k + µ 4 n k=0 δ n+1-k 1 2 ∂ s (κ∂ s ) + κ 3 + ∆ Σ κ 8 - ∂ n W n+1 4 - ∂ n β(ϕ n+1 ) 4 ϕ k . (8.12)
The boundary term

Σ ∂ n ϕ n+1 χ dΣ is then associated to a matrix term -B n+1 ϕ n+1 - b n
, where

B n+1 ϕ n+1 = µ 1 M Σ ϕ n+1 + µ 2 M Σ κ ϕ n+1 + µ 3 1 2 M Σ κ 2 ϕ n+1 - 1 2 S Σ ϕ n+1 + 1 2 M Σ W n+1 ϕ n+1 + 1 2 M Σ β(ϕ n+1 ) ϕ n+1 + µ 4 - 1 2 S Σ κ + 1 8 M Σ κ 3 +∆Σκ ϕ n+1 - 1 4 M Σ ∂nW n+1 ϕ n+1 - 1 4 M Σ ∂nβ(ϕ n+1 ) ϕ n+1 ,
and

b n = µ 1 n k=1 β n+1-k M Σ ϕ k + µ 3 n k=1 α n+1-k 1 2 M Σ κ 2 - 1 2 S Σ ϕ k + µ 3 n k=1 α n+1-k 1 2 M Σ sg(W n+1 ) √ |W n+1 | √ |W k | + 1 2 M Σ sg(β(ϕ n+1 )) √ |β(ϕ n+1 )| √ |β(ϕ k )| ϕ k + µ 4 n k=1 δ n+1-k - 1 2 S Σ κ + 1 8 M Σ κ 3 +∆Σ - 1 4 M Σ ∂nW n+1 - 1 4 M Σ ∂nβ(ϕ n+1 ) ϕ k . (8.13)
In the above equations, M Σ and S Σ respectively designate the mass and stiffness matrices associated with Σ for linear elements. The extension to generalized mass and stiffness matrices for a function F is specified by M Σ F and S Σ F . To apply a fixed point method, the term related to ϕ n+1 can be decomposed into linear and nonlinear terms with respect to ϕ n+1 as (8.2), where the matrix B n+1 L is given by

B n+1 L = µ 1 M Σ + µ 2 M Σ κ + µ 3 1 2 M Σ κ 2 - 1 2 S Σ + 1 2 M Σ W n+1 + µ 4 - 1 2 S Σ κ + 1 8 M Σ κ 3 +∆Σκ - 1 4 M Σ ∂nW n+1 . (8.14)
Finally, at each iteration, we solve the linear system

2i ∆t M -B n+1 L -S + M W n+1 ζ p+1 = 2i ∆t Mψ n -M β(ζ p ) ζ p + b n + 1 2 M Σ β(ζ p ) ζ p - 1 4 M Σ ∂nβ(ζ p ) ζ p , (8.15)
where ζ p is initialized by ζ 0 = ϕ n , the matrix B n+1 L is given by (8.14), and b n , which only depends on ϕ 0 , . . . , ϕ n , by (8.13). At convergence, we define then ϕ n+1 = ζ p+1 .

Remark 1. In practice, we do not numerically treat the fourth-order boundary condition. Indeed, this condition takes into account the term ∂ n β(u) which is delicate to numerically compute as already shown for the potential case [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF]. However, this should not be a too restrictive condition since the third-order ABC already gave some satisfactory results for the purely potential case.

We can now study the case of condition NLABC 2 1,T which comes from the gauge change. Its discretization is given by 4 2 ∆t e iW n+1 e iW n+1 n+1 k=1

∂ n ϕ n+1 + e -iπ/
β n+1-k e -iW k e -iW k ϕ k + κ 2 ϕ n+1 = 0, (8.16)
where the phase function W n+1 is associated to V and W n+1 to the nonlinearity β(ψ)

W n+1 = t n+1/2 0 V (x, σ)dσ, (8.17
W n+1 = t n+1/2 0 β(ϕ n+1 )(x, σ)dσ. (8.18)
Let us set for p ≥ 1

E p = exp (iW p ) = exp i ∆t 2 β(ϕ p ) E p-1 ,
where E p-1 does no longer depend on ϕ n+1 : E p-1 = exp i∆t p-1 l=1 β(ϕ l ) for p ≥ 1, and E 0 = 1. Expliciting the terms in ϕ n+1 , we obtain

∂ n ϕ n+1 = µ 1 ϕ n+1 + µ 2 κϕ n+1 + µ 1 e iW n+1 e i ∆t 2 β(ϕ n+1 ) E n n k=1 β n+1-k e -iW k E k ϕ k .
Then, the equation solved by the fixed point procedure is

∂ n ζ p+1 = µ 1 ζ p+1 + µ 2 κζ p+1 + µ 1 e iW n+1 e i ∆t 2 β(ζ p ) E n n k=1 β n+1-k e -iW k E k ϕ k ,
where ζ p is initialized by ζ 0 = ϕ n . The matrix form of the boundary term

Σ ∂ n ζ p+1 χ dΣ is then µ 1 M Σ + µ 2 M Σ κ ζ p+1 + µ 1 n k=1 β n+1-k M Σ e i(W n+1 -W k ) e i ∆t 2 β(ζ p ) E n E k ϕ k . (8.19)
We can see that the right hand side b n is a convolution which includes some generalized mass matrices associated with the function e i(W n+1 -W k ) . If the linear part V of the potential is radial, the expression (8. [START_REF] Lindmann | Free-space boundary conditions for the time dependent wave equation[END_REF]) can be partially simplified as

µ 1 M Σ + µ 2 M Σ κ ζ p+1 + µ 1 e iW n+1 n k=1 β n+1-k M Σ e i ∆t 2 β(ζ p ) E n E k e -iW k ϕ k . (8.20)
Globally, for the condition NLABC 2 1,T , the fixed point equation to solve is

2i ∆t M -B Σ -S + M W n+1 ζ p+1 = 2i ∆t Mψ n -M β(ζ p ) ζ p + b n , (8.21)
with the boundary matrix:

B Σ = µ 1 M Σ + µ 2 M Σ κ , and the right-hand side b n = µ 1 n k=1 β n+1-k M Σ e i(W n+1 -W k ) e i ∆t 2 β(ζ p ) E n E k ϕ k .
If the potential or the domain is non radial, then this condition is computationally very expensive since n generalized mass matrices must be reassembled at each time t n . On a whole resolution with N time steps, assembling the terms b n successively requires the assembling of N (N + 1)/2 generalized mass matrices. Therefore, we observe a large computational time for NLABC 2 1,T , compared to the conditions NLABC M 2,T . Furthermore, this is more and more penalizing as ∆t goes to zero. 8.3. Stability of the semi-discrete problems with NLABC 2 1,T and NLABC M 2,T . For the conditions based on Taylor's expansions, we can prove the following result. We do not develop the proof and refer to [START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF] for more details.

Theorem 8.1. Let (ψ n ) 0≤n≤N be a solution of the system

         i ψ n+1 -ψ n ∆t + ∆ϕ n+1 + W n+1 ϕ n+1 + β(ϕ n+1 ) ϕ n+1 = 0, in Ω, ∂ n ϕ n+1 + Λ M,n+1 2,T
ϕ n+1 = 0, on Σ, for M ∈ {2, 3},

ψ 0 = ψ I , in Ω, (8.22)
where the operators Λ 2,n+1 2,T and Λ 3,n+1 2,T are defined by

Λ 2,n+1 2,T ϕ n+1 = e -iπ/4 2 ∆t b n+1 + κ 2 ϕ n+1 , and
Λ 3,n+1 2,T ϕ n+1 = Λ 2,n+1 2,T ϕ n+1 -e iπ/4 ∆t 2 κ 2 2 a n+1 0 + 1 2 a n+1 2 + 1 2 sg(W n+1 ) |W n+1 | a n+1 V + 1 2 sg(β(ϕ n+1 )) |β(ϕ n+1 )| a n+1 β .
For M = 2, we have the inequality ∀n ∈ {0, . . . , N },

ψ n L 2 (Ω) ≤ ψ 0 L 2 (Ω) . (8.23)
Moreover, if V is a time independent positive potential on Σ, and if β is a positive function, then the inequality (8.23) remains valid for M = 3. Now let us study the condition NLABC 2 1,T with gauge change which is an extension of the pure potential case (see [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Klein | Construction et analyse de conditions aux limites artificielles pour des équations de Schrödinger avec potentiels et non linéarités[END_REF] for a detailed proof). Theorem 8.2. Let (ψ n ) 0≤n≤N be a solution of the system with the ABC NLABC 2

1,T              i ψ n+1 -ψ n ∆t + ∆ψ n+1 + W n+1 ψ n+1 + β(ψ n+1 ) ψ n+1 = 0, in Ω, ∂ n ψ n+1 + e -iπ/4 2 ∆t e iW n+1 e iW n+1 b n+1 + κ 2 ψ n+1 = 0, on Σ, ψ 0 = ψ 0 , in Ω, (8.24)
where the phase functions W n+1 and W n+1 are defined by (8.17) and (8.18). Then we get the inequality ∀n ∈ {0, . . . , N }, for the Padé approximation is expressed through the three equations (6.9), (6.7) and (6.8). The discretization of the system is

ψ n L 2 (Ω) ≤ ψ 0 L 2 (Ω) . ( 8
                         ∂ n ϕ n+1 -i m k=0 a m k ϕ n+1 + κ 2 ϕ n+1 + i m k=1 a m k d m k ϕ n+1/2 k - κ 2 Ψ n+1/2 = 0, 2i ∆t + ∆ Σ + W n+1 + β(ϕ n+1 ) + d m k ϕ n+1/2 k -ϕ n+1 = 2i ∆t ϕ n k , 2i ∆t + ∆ Σ + W n+1 + β(ϕ n+1 ) Ψ n+1/2 -∆ Σ ϕ n+1 = 2i ∆t Ψ n , ϕ 0 k = 0 for 1 ≤ k ≤ m, Ψ 0 = 0, on Σ. ( 8
.26) To solve this nonlinear system, we apply a fixed point algorithm, initialized by

ζ 0 = ϕ n , φ 0 k = ϕ n-1/2 k , for 1 ≤ k ≤ m, χ 0 = Ψ n-1/2
. Next, we determine at each step the solutions ζ p+1 , φ p+1 k and χ p+1 to the system

                   ∂ n ζ p+1 -i m k=0 a m k ζ p+1 + κ 2 ζ p+1 + i m k=1 a m k d m k φ p+1 k - κ 2 χ p+1 = 0, 2i ∆t + ∆ Σ + W n+1 + d m k φ p+1 k -ζ p+1 = 2i ∆t ϕ n k -β(ζ p )φ p k , 2i ∆t + ∆ Σ + W n+1 χ p+1 -∆ Σ ζ p+1 = 2i ∆t ψ n -β(ζ p )χ p . (8.27)
The rule is to replace ϕ n+1 (resp. ϕ n+1 k , Ψ n+1 ) by ζ p+1 (resp. φ p+1 k , χ p+1 ) in the linear expressions in ϕ n+1 , and by ζ p (resp. φ p k , χ p ) in the nonlinear expressions in ϕ n+1 or the expressions like ϕ n+1 ϕ n+1 k . Then, we write the variational formulation and discretize system (8.27) through a FEM by injecting the first equation into the variational formulation (8.1) of the interior equation. We get at the matrix level

2i ∆t M -S + M W n+1 + i m k=0 a m k M Σ - 1 2 M Σ κ ζ p+1 -i m k=1 a m k d m k M Σ φ p+1 k + 1 2 M Σ κ χ p+1 = 2i ∆t M ψ n -M β(ζ p ) ζ p , (8.28)
coupled to [START_REF] Zhang | Unified approach to split absorbing boundary conditions for nonlinear Schrödinger equations[END_REF]) and

-M Σ ζ p+1 + 2i ∆t M Σ -S Σ + M Σ W n+1 + d m k M Σ φ p+1 k = 2i ∆t M Σ ϕ n k -M Σ β(ζ p ) φ p k , (8.
S Σ ζ p+1 + 2i ∆t M Σ -S Σ + M Σ W n+1 χ p+1 = 2i ∆t M Σ ψ n -M Σ β(ζ p ) χ p . (8.30)
Globally, the matrix form of the system with the ABC NLABC 2 2,P leads to solve

           A Ω A 0 0 A A Σ B 1 . . . B m B C 1 D 1 0 0 . . . . . . 0 C m 0 D m 0 C 0 D                       ζ Ω,p+1 ζ Σ,p+1 φ p+1 1 . . . φ p+1 m χ p+1            =            E Ω E Σ F 1 . . . F m F            (8.31)
with unknowns ζ p+1 , φ p+1 k , for 1 ≤ k ≤ m, and χ p+1 . The different elements of the system are given by

A = 2i ∆t M -S + M W n+1 + i m k=1 a m k M Σ - 1 2 M Σ κ , B k = -ia m k d m k M Σ , C k = -M Σ , for 1 ≤ k ≤ m, D k = 2i ∆t M Σ -S Σ + M Σ W n+1 + d m k M Σ , for 1 ≤ k ≤ m, C = S Σ , B = 1 2 M Σ κ , D = 2i ∆t M Σ -S Σ + M Σ W n+1 , (8.32)
and

E = 2i ∆t M ψ n -M β(ζ p ) ζ p , F = 2i ∆t M Σ ψ n -M Σ β(ζ p ) χ p , F k = 2i ∆t M Σ ϕ n k -M Σ β(ζ p ) φ p k , for 1 ≤ k ≤ m. (8.33)
In the above notations, the global matrix A and the global vector E have been reordered according to the mesh by splitting the interior and boundary points. In particular, let us consider that we have n Ω vertices in Ω and n Σ on Σ for a global number of points n P = n Ω + n Σ . Then, the notation A Ω , A , A and A Σ are the corresponding matrices of respectively size n Ω × n Ω , n Ω × n Σ , n Σ × n Ω and n Σ × n Σ . coming from a matrix A and vector E such that

A = A Ω A A A Σ , E = E Ω E Σ . (8.34)
9. Semi-discretization and implementation of the relaxation method. Instead of solving the NLS equation by a fixed point, we can use the relaxation scheme

[9]        i ψ n+1 -ψ n ∆t + ∆ψ n+1/2 + Υ n+1/2 ψ n+1/2 = 0, Υ n+3/2 + Υ n+1/2 2 = β(ψ n+1 ), (9.1)
for 0 ≤ n ≤ N , setting Υ n+1/2 = Υ n+1 +Υ n 2
, and by initializing Υ by Υ -1/2 = β(ϕ 0 ). The weak formulation associated with (9.1) writes down, for a test-

function χ ∈ H 1 (Ω), 2i ∆t Ω ϕ n+1 χ dΩ - Ω ∇ϕ n+1 • ∇χ dΩ + Σ ∂ n ϕ n+1 χ dΣ + Ω Υ n+1/2 ϕ n+1 χ dΩ = 2i ∆t Ω ψ n χ dΩ, (9.
2) coupled to the update: Υ n+3/2 = 2β(ϕ n+1 ) -Υ n+1/2 . In the boundary term Σ, we replace ∂ n ϕ n+1 by using one of the previous ABCs and writing Υ = β(ψ). One of the crucial advantages of this method is that it does not require any iteration like a fixed point and so is much faster. Furthermore, the adaptation of a code for variable potential (linear problems) to nonlinear problems is direct and therefore easy to implement. The gain in computational cost is most particularly important for NLABC 2 1,T (even if it is still more computationally consuming than the other ABCs). The relaxation scheme (9.1) is converging [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] for a nonlinearity β(ψ) = λ|ψ| 2σ . However, since the proof is already extremely complex we do not intend here to study the stability of the resulting scheme. Another important property of the relaxation scheme is that it preserves the two invariants, the density R 2 |ψ| 2 dx and the energy R 2 |∇ψ| 2 + q 2 |ψ| 4 dx. One drawback is that for the nonlinearity β(ψ) = λ|ψ| 2σ with λ ∈ R * presented in [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF], it is not proved that we get Υ n+1/2 ≥ 0 for n ≥ 0 while it is true for Υ. However, no counter-example has been found until now. 9.1. Application of the relaxation scheme to the ABC NLABC 2 2,P . For the ABC NLABC 2 2,P , the relaxation scheme writes down

                                         2i ∆t ϕ n+1 + ∆ϕ n+1 + Υ n+1/2 ϕ n+1 = 2i ∆t ψ n , on Ω, ∂ n ϕ n+1 -i m k=0 a m k ϕ n+1 + κ 2 ϕ n+1 + i m k=1 a m k d m k ϕ n+1/2 k - κ 2 ψ n+1/2 = 0, on Σ, 2i ∆t + ∆ Σ + Υ n+1/2 + d m k ϕ n+1/2 k -ϕ n+1 = 2i ∆t ϕ n k , on Σ, 2i ∆t + ∆ Σ + Υ n+1/2 Ψ n+1/2 -∆ Σ ϕ n+1 = 2i ∆t ψ n , on Σ, ϕ 0 k = 0 for 1 ≤ k ≤ m, Ψ 0 = 0, on Σ, Υ n+3/2 = 2β(ϕ n+1 ) -Υ n+1/2 , on Ω. (9.
3)

The associated matrix formulation of the global system is

           A Ω A 0 0 A A Σ B 1 . . . B m B C 1 D 1 0 0 . . . . . . 0 C m 0 D m 0 C 0 D                        ϕ Ω,n+1 ϕ Σ,n+1 ϕ n+1/2 1 . . . ϕ n+1/2 m Ψ n+1/2             =            E Ω E Σ F 1 . . . F m F            (9.4)
where the different matrices are given by

A = 2i ∆t M -S + M Υ n+1/2 + i m k=1 a m k M Σ - 1 2 M Σ κ , B k = -ia m k d m k M Σ , for 1 ≤ k ≤ m, B = 1 2 M Σ κ , C k = -M Σ , for 1 ≤ k ≤ m, C = S Σ , D k = 2i ∆t M Σ -S Σ + M Σ Υ n+1/2 + d m k M Σ , for 1 ≤ k ≤ m, D = 2i ∆t M Σ -S Σ + M Σ Υ n+1/2 , (9.5)
and the vectors in the right hand side are

E = 2i ∆t M ψ n , F k = 2i ∆t M Σ ϕ n k , for 1 ≤ k ≤ m, F = 2i ∆t M Σ Ψ n . (9.6)
The extension to the case of a mixed potential V (x, y, t) + β(ψ) is direct. 9.2. Application of the relaxation scheme to the ABC NLABC 3 2,T . The third-order boundary condition NLABC 3 2,T writes, for a nonlinearity β(ψ) where β is a positive function,

∂ n ψ = -e -iπ/4 ∂ 1/2 t ψ - κ 2 ψ + e iπ/4 κ 2 8 + ∆ Σ 2 + β(ψ) 2 I 1/2 t ψ, on Σ T .
Let us remark that we did not symetrize the potential term in front of I 1/2 t since we do not know for the relaxation scheme is the discretization of Υ remains positive. The application of the relaxation scheme to NLABC 3 2,T gives

                           2i ∆t ϕ n+1 + ∆ϕ n+1 + Υ n+1/2 ϕ n+1 = 0, on Ω, ∂ n ϕ n+1 = -e -iπ/4 2 ∆t n+1 k=1 β n+1-k ϕ k - 1 2 κϕ n+1 + e iπ/4 2 ∆t 2 n+1 k=1 α n+1-k κ 2 2 + 1 2 ∂ 2 s + Υ n+1/2 2 ϕ k , on Σ, Υ n+3/2 = 2f (ϕ n+1 ) -Υ n+1/2 , on Ω. ( 9.7)
The finite element formulation leads to the linear system

2i ∆t M -B Σ -S + M Υ n+1/2 ϕ n+1 = 2i ∆t Mψ n + b n , (9.8)
where the boundary terms B Σ are

B Σ = e -iπ/4 2 ∆t M Σ + 1 2 M Σ κ + e iπ/4 2 ∆t 2 1 2 M Σ κ - 1 2 S Σ + 1 2 M Σ Υ n+1/2 , (9.9)
and the right hand side b n is

b n = e -iπ/4 2 ∆t n k=1 β n+1-k M Σ ϕ k + e iπ/4 2 ∆t 2 n k=1 α n+1-k 1 2 M Σ κ - 1 2 S Σ + 1 2 M Σ Υ n+1/2 ϕ k . (9.10)
The solution to (9.8) is then updated with: Υ n+3/2 = 2f (u n+1 ) -Υ n+1/2 . 9.3. Application of the relaxation to the ABC NLABC 2 1,T . We consider now the ABC NLABC 2 1,T

∂ n ψ = -e -iπ/4 e iB ∂ 1/2 t e -iB ψ - κ 2 ψ, on T ,
with the phase function given by (6.4). The relaxation scheme is then .11) We remark that when computing ψ n+1 , the term W n+1 needs Υ k+1/2 until Υ n+1/2 , which are already known. This is different from the fixed point algorithm where E n+1 still depends on ψ n+1 , implying the introduction of E p . The FEM leads to

                           2i ∆t ψ n+1 + ∆ψ n+1 + Υ n+1/2 ψ n+1 = 2i ∆t ϕ n , on Ω, W n+1 = exp i ∆t 2 Υ n+1/2 exp i∆t n-1 k=0 Υ k+1/2 , on Σ, ∂ n ψ n+1 = -e -iπ/4 2 ∆t e iW n+1 n+1 k=1 β n+1-k e -iW k ψ k - 1 2 κ ψ n+1 , on Σ, Υ n+3/2 = 2f (ψ n+1 ) -Υ n+1/2 , on Ω. ( 9
2i ∆t M -B Σ -S + M Υ n+1/2 ψ n+1 = 2i ∆t Mϕ n + b n , (9.12)
with: 428 000 215 000 1 264 2 × 10 -4 2 × 10 -3

B Σ = e -iπ/4 2 ∆t M Σ + 1 2 M Σ κ ,
n = e -iπ/4 2 ∆t n k=1 β n+1-k M Σ e i(W n+1 -W k ) ψ k . ( 9
Table 10.1 Mesh parameters for the disk of radius 10. 10. A numerical example: the cubic nonlinearity. We numerically compare now the three ABCs NLABC 2 1,T , NLABC 3 2,T and NLABC 2 2,P . These are also compared to the case without potential ABC 3 0,T and ABC 2 0,P . If one considers a gaussian initial data, then, for the cubic (defocusing) case |ψ| 2 , the ABC with a null potential already gives a good accuracy. Indeed, in this case, the effect of the laplacian is predominant and the nonlinearity has almost no effect. For an initial soliton data for which their is a trade-off between the dispersion and nonlinear term, we need to use the new ABCs. In the 2D case, there is no expression of the soliton which must therefore be built numerically. Here, this is done by solving the stationary NLS equation with a shooting method [START_REF] Menza | Numerical Computation of Solitons for Optical Systems[END_REF]: we find localized stationary radial solutions u(r, t) = e iµt ψ(r), r

= x = x 2 + y 2 , µ ∈ R, to: ∂ 2 r ψ + 1 r ∂ r ψ -µψ + q|ψ| 2 ψ = 0, r ∈ [0; R],
with boundary conditions: ψ ′ (0) = 0 and ψ(0) = β. Then, we get a soliton solution by setting ψ 0 (x, y) = ψ(r)e -ik0x , for a vector wave k 0 = k 0 x. This moving soliton solution is computed with a high accuracy. We can see on Figures 10.1 and 10.2 that the solution is slowly decaying. In the following, we will consider as computational domain Ω the disk with radius 10. We refer to Table 10.1 for the characteristics of the meshes (n T : number of triangles T , n P : number of nodes, n E : number of edges) for the P 1 FEM. Table 10.2 Comparison of the computational times (on Matlab) for the fixed point resolution and the relaxation scheme, the discretization is given by n T = 423 000 and ∆t = 2 × 10 -3 . are identical and independent of the method. The same conclusion is true for others ABCs. This means that computational times are the next step for choosing the method. We report the results on table 10.2, for the two methods and for the three ABCs NLABC 2 1,T , NLABC 3 2,T and NLABC 2 2,P . As it is expected, we clearly see the efficiency of the relaxation method compared to the fixed point approach. Therefore, from now on, we always consider the relaxation scheme. 10.2. Cubic nonlinear potential β(ψ) = q|ψ| 2 with initial soliton data. We consider the cubic nonlinearity β(ψ) = |ψ| 2 and the initial soliton data previously built on the disk of radius R = 10, with a wavenumber k 0 = 5 and T = 2. The three ABCs are compared with ABC 3 0,T and ABC 2 0,P (which are given by NLABC 3 0,T and NLABC 2 0,P for the potential g = 0). We represent on Figures 10.5-10.8 the results obtained for the discretization n T = 1 700 000 and ∆t = 2 × 10 -3 . First, we observe that the soliton moves without dispersion before it hits the boundary. Then, we remark that ABC 3 0,T generates a lot of reflection at the boundary at a level of about 10 -1.5 (we obtain the same results for ABC 2 0,P ). The NLABCs are able to correctly reproduce the solution. The less accurate is NLABC 3 2,T because a reflection of about 10 -1.5 is still visible along the boundary x = (-10, 0). This is no longer true for NLABC 2 1,T and NLABC 2 2,P which clearly increase the accuracy. The ABC NLABC 2 2,P from the Padé approximation is slightly more precise than for the Taylor approach with gauge change. In terms of computational times, the most efficient ABCs are NLABC 3 2,T and NLABC 2 2,P which have a cost about the same as for ABCs without potential. This is a very interesting point. For the ABCs with gauge change, the CPU time increases a lot, compared to NLABC 3 2,T and NLABC 2 2,P . This is mainly due to the reassembling processes required in the convolution computation. This is most particularly true for fine meshes. 11. Conclusion. Accurate and efficient increasing order Absorbing Boundary Conditions (ABCs) have been proposed for one-and two-dimensional Schrödinger equations including general potentials and nonlinearities. Time discretization schemes have been introduced in conjunction with finite element methods. Numerical examples show the importance of considering these ABCs for an accurate and realistic computation. Most particularly, for the two-dimensional case, we have proposed three ABCs: NLABC 2 1,T , NLABC 3 2,T and NLABC 2 2,P . Two methods of discretization in time have been introduced: fixed point algorithm and relaxation scheme. Stability issues have also been analyzed, most particularly for the Taylor based approximations. The relaxation schemes preserves the accuracy and yields a faster algorithm. In terms of accuracy and efficiency and for a soliton initial data, NLABC 2 2,P provides the best results compared to NLABC 2 1,T and NLABC 3 2,T for the cubic nonlinearity.

3. 1 .

 1 Discrete convolutions based discretizations for NLABC N 0 . The numerical difficulties related to the NLABC N 0 conditions are twofold. The first one concerns the approximation of fractional operators ∂

Proposition 3 . 1 .

 31 The semi-discrete Durán Sanz-Serna scheme for the initial boundary value problem(3.1) given by 

Fig. 5 . 1 .

 51 Fig. 5.1. Reference solution for the mixed nonlinearity

Fig. 5 . 2 .

 52 Fig. 5.2. Mixed nonlinearity and NLABC 2 0

Fig. 5 . 3 .

 53 Fig. 5.3. Mixed nonlinearity and NLABC 4 0

  Fig. 5.4.Repulsive interaction and NLABC 2 0

8. 1 .

 1 Variational formulation and fixed point algorithm. We use the Duràn-Sanz-Serna scheme 2i

V , a n+1 γ for γ ∈ {0, 1 ,

 1 2} and d n+1 γ for γ = {0, 1} are defined by (8.5), (8.7), (8.6) and (8.8) with the notations b n+1 = β n+1 ⋆ ϕ n+1 , (8.5)

  )

  and where the right hand side b n is b

Fig. 10 . 1 .

 101 Fig. 10.1. Modulus of the initial soliton data for the cubic NLS equation on the disk of radius 10 (normal scale).

Fig. 10 . 2 .

 102 Fig. 10.2. Modulus of the initial soliton data for the cubic NLS equation on the disk of radius 15 (logarithmic scale).

10. 1 .

 1 Choice of the method of resolution. We briefly compare the fixed point and the relaxation schemes in terms of computational time as well as accuracy. For the fixed point, the stopping criterion is:ζ p+1ζ p 2 L 2 (C n P ) < ε, for ε = 10-6 . The linear systems are solved by the direct solver UMFPACK. We represent the results on Figures 10.3-10.4 obtained with the boundary condition NLABC 3 2,T for the cubic nonlinearity |ψ| 2 and the soliton initial data. We draw the cut plane (x, t) and compute |ψ|. The discretization considers n T = 423 000 triangles and ∆t = 2 × 10 -3 , for k 0 = 5 and T = 2. We observe on Figures 10.3-10.4 that the results for NLABC 3

Fig. 10 . 3 .

 103 Fig. 10.3. Solution |ψ| (in the (x, t) plane) computed with NLABC 3 2,T for the fixed point algorithm.

Fig. 10 . 4 .

 104 Fig. 10.4. Solution |ψ| (in the (x, t) plane) computed with NLABC 3 2,T for the relaxation scheme.

Fig. 10 . 5 .

 105 Fig. 10.5. Numerical solution |ψ| (in the (x, t) plane) for ABC 3 0,T (log scale).

Fig. 10 . 6 .

 106 Fig.10.6. Numerical solution|ψ| (in the (x, t) plane) for NLABC2 1,T (log scale).

Fig. 10 . 7 .

 107 Fig. 10.7. Numerical solution |ψ| (in the (x, t) plane) for NLABC 3 2,T (log scale).

Fig. 10 . 8 .

 108 Fig. 10.8. Numerical solution |ψ| (in the (x, t) plane) for NLABC 2 2,P (log scale).

  × 10 -5 4 × 10 -4

	n T 1 700 000 858 000 2 528 6 n P n E	min |T |	max |T |

.13)

After computing ψ n+1 , Υ n+3/2 is updated. Again, the extension to mixed potentials is direct.

The authors are partially supported by the French ANR fundings under the project MicroWave NT09 460489.