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

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

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

Matthias Ehrhardt
email: ehrhardt@math.uni-wuppertal.de

Absorbing Boundary Conditions for Solving N-Dimensional Stationary Schrödinger Equations with Unbounded Potentials and Nonlinearities

Keywords: Absorbing boundary conditions, stationary Schrödinger equations, unbounded domain, spatially dependent potential, ground states computation PACS: 03.65.-w, 42.82.Et, 43.30.Ma 2000 MSC: 35J10, 65M60, 65N30

come   L'archive ouverte pluridisciplinaire

Introduction

The solution of the Schrödinger equation occurs in many applications in physics, chemistry and engineering (e.g. quantum transport, condensed matter physics, quantum chemistry, optics, underwater acoustics, . . . ). The considered problem can appear in different forms: time-dependent or stationary equation, linear or nonlinear equation, inclusion of a variable potential among others. One of the main difficulty when solving the Schrödinger equation, and most particularly from a numerical point of view, is to impose suitable and physically admissible boundary conditions to solve numerically a bounded domain equation modeling an equation originally posed on an unbounded domain. Concerning the timedomain problem, many efforts have been achieved these last years. We refer the interested reader e.g. to the recent review paper [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF] and the references therein for further details.

In this paper, we begin to focus on the solution to the one-dimensional stationary Schrödinger equation. For a given potential V possibly extending to infinity, eventually nonlinear (V := V (x, ϕ), we want to solve the following equation

-α d 2 dx 2 + V ϕ = Eϕ, x ∈ R, (1)
or rewritten as

d 2 dx 2 + 1 α E -V ϕ = 0, x ∈ R, (2)
with some parameter α that allows for some flexibility. More precisely, we study the extension of the recently derived time-domain boundary conditions [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] to the computation of stationary states: we determine here the pair (ϕ, E), for a given linear or nonlinear potential V . This eigenvalue problem is also known as the computation of ground states. The energy of the system is then the eigenvalue E and the associated stationary state is the eigenfunction ϕ. In particular, we seek the fundamental stationary state which is linked to the smallest eigenvalue. In practice, higher order states are also of interest. After a careful numerical validation of the 1D case, we develop some ABCs for the generalized N -dimensional case

-α∆u + V (x)u + f (u)u = Eu
in R N for an unbounded potential V and general nonlinearities. Another related problem which is not treated here is linear and nonlinear scattering. We refer to our extended version of the present paper [START_REF] Antoine | Modeling boundary conditions for solving stationary Schrödinger equations[END_REF] where a thorough study of scattering problems with the ABCs is developed.

Let us note here that numerical approaches based on solving the stationary Schrödinger with ABCs is not the only practical possibility. Indeed, other solutions may be used. The most widely used approach with applications for example in Bose-Einstein condensation is to bound the infinite space with a sufficiently large computational domain and then to impose a Dirichlet boundary condition on this fictitious boundary. The fact that we may fix the field at zero at the boundary is related to the property that ground states decay exponentially fast in the far-field. Then, spectrally accurate solutions can be efficiently obtained following the works by Bao et al. [START_REF] Bao | Efficient and spectrally accurate numerical methods for computing ground and first excited states in Bose-Einstein condensates[END_REF][START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF]. Other solutions include whole space solutions based on spectral methods (see [START_REF] Bao | A Generalized-Laguerre-Fourier-Hermite pseudospectral method for computing the dynamics of rotating Bose-Einstein condensates[END_REF][START_REF] Bao | A generalized-Laguerre-Hermite pseudospectral method for computing symmetric and central vortex states in Bose-Einstein condensates[END_REF]). Let us also note the following other strategy [START_REF] Ledoux | MATSLISE: A MATLAB package for the numerical solution of Sturm-Liouville and Schrodinger equations[END_REF][START_REF] Ledoux | Solution of the Schrodinger equation over an infinite integration interval by perturbation methods, revisited[END_REF][START_REF] Ledoux | Study of Study of Special Algorithms for solving Sturm-Liouville and Schrodinger Equations[END_REF] for solving the linear Schrödinger equation with potential. If you approximate the potential one each subinterval by a constant, linear or quadratic function, then you get one each subinterval a problem that can be solved explicitly (harmonic oscillator, Airy equation) and these sub-problems are solved sequentially by some transfer algorithm. For constant potentials this approach was worked out in [START_REF] Ledoux | MATSLISE: A MATLAB package for the numerical solution of Sturm-Liouville and Schrodinger equations[END_REF][START_REF] Ledoux | Solution of the Schrodinger equation over an infinite integration interval by perturbation methods, revisited[END_REF][START_REF] Ledoux | Study of Study of Special Algorithms for solving Sturm-Liouville and Schrodinger Equations[END_REF] (see the MATLISE solver available at http://users.ugent.be/∼vledoux/). Obviously, this strategy can be generalized to linear/quadratic potentials.

For the stationary Schrödinger equation [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF], boundary conditions for solving linear scattering problems with a constant potential outside a finite domain have been proposed e.g. by Ben Abdallah, Degond and Markowich [START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF], by Arnold [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF] for a fully discrete Schrödinger equation and in a two-dimensional quantum waveguide by Lent and Kirkner [START_REF] Lent | The quantum transmitting boundary method[END_REF][START_REF] Kirkner | The numerical simulation of electron transmission through a two-dimensional quantum device by the finite element method[END_REF]. Other recent contributions include [START_REF] Heinen | Radiation boundary conditions for the numerical solution of the three-dimensional time-dependent Schrödinger equation with a localized interaction[END_REF][START_REF] Heinen | Numerical calculation of strong-field laser-atom interaction: An approach with perfect reflection-free radiation boundary conditions[END_REF] where the case of finite range potentials is treated. The case of bound states can be found for specific one-dimensional linear Schrödinger equations in [START_REF] Ehrhardt | Solutions to the discrete Airy equation: Application to parabolic equation calculations[END_REF][START_REF] Ehrhardt | Fast calculation of energy and mass preserving solutions of Schrödinger-Poisson systems on unbounded domains[END_REF][START_REF] Moyer | Numerov extension of transparent boundary conditions for the Schrödinger equation in one dimension[END_REF][START_REF] Moyer | Numerical solution of the stationary state Schrödinger equation using transparent boundary conditions[END_REF][START_REF] Ehrhardt | Exact artificial boundary conditions for problems with periodic structures[END_REF]. These boundary conditions are needed e.g. to improve existing simulation tools for semiconductors that allows to investigate certain stationary (and also transient) behavior of the devices, like conductance, capacity, current-voltage curves. Often the physical relevant effects take place only in a small subregion of the device, and the novel absorbing boundary conditions offer the possibility to confine the computations to this small domain. We refer the reader to [START_REF] Burgnies | Mécanismes de conduction en régime ballistique dans les dispositifs électroniques quantiques[END_REF][START_REF] Frensley | Boundary conditions for open quantum driven far from equilibrium[END_REF][START_REF] Markowich | Semiconductor Equations[END_REF][START_REF] Schmeiser | Computational methods for semiclassical and quantum transport in semiconductor devices[END_REF] for more application details.

The goal of this work is to propose and numerically validate some new boundary conditions for modeling linear and nonlinear variable unbounded potentials stationary onedimensional Schrödinger equations with application to ground-state computation. Finally, we extend these absorbing boundary conditions to higher dimensional problems. The paper is organized as follows. In Section 2, we explain how to obtain the stationary boundary conditions from the time-dependent case. Sections 3 and 4 are respectively devoted to their applications to linear and nonlinear eigenstate computation. Section 5 gives an extension of the boundary conditions to the N -dimensional space for nonlinear stationary Schrödinger equations with unbounded potentials. Finally, Section 6 draws a conclusion and give an outlook for possible future research directions.

Absorbing boundary conditions: from the time-domain to the stationary case: the one-dimensional case

In order to derive some absorbing boundary conditions (ABCs) for the stationary Schrödinger equation ( 2), let us first start with the time-domain situation. In case of the timedependent Schrödinger equation with a linear or nonlinear potential

V i∂ t u + ∂ 2 x u + V u = 0, ∀(x, t) ∈ R × R + , u(x, 0) = u 0 (x), x ∈ R, (3)
the following second-and fourth-order ABCs on the boundary

Σ × R + ABC 2 2 : ∂ n u = i i∂ t + V u, (4)
ABC 4 2 :

∂ n u = i i∂ t + V u - ∂ n V 4 (i∂ t + V ) -1 u, (5)
were derived recently in [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]. Here, for an operator A, √ A denotes the square-root operator of A [START_REF] Taylor | Pseudodifferential Operators[END_REF] with respect to its spectral decomposition. The fictitious boundary Σ is given by the two interval endpoints x ℓ and x r . The outwardly directed unit normal vector to the bounded computational domain Ω =]x ℓ ; x r [ is denoted by n.

To obtain some ABCs for the stationary equations (1) or (2), we consider these equations supplied with a new potential: V := -V /α. Moreover, we are seeking some time-harmonic solutions u(x, t) := ϕ(x)e -i E α t . Since

i∂ t u = E α ϕ(x) e -i E α t ,
we have

i∂ t + V u = 1 √ α √ E -V (ϕe -i E α t )
These considerations yield some stationary ABCs that we designate by SABC M ('S' stands for stationary and M denotes the order of the boundary condition) :

SABC 2 : ∂ n ϕ = i 1 √ α √ E -V ϕ, on Σ, (6)
SABC 4 :

∂ n ϕ = i 1 √ α √ E -V ϕ + 1 4 ∂ n V E -V ϕ, on Σ. ( 7
)
The second-and fourth-order ABCs for the time-dependent Schrödinger equation ( 4), [START_REF] Bao | Efficient numerical methods for computing ground states and dynamics of dipolar Bose-Einstein condensates[END_REF] were developed under a high frequency assumption [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]. This relation can be translated to the stationary case in terms of links between E and V . The new relations will be given for the eigenvalues problems in the next dedicated sections.

Remark 1. For the time-dependent case [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF], we constructed two families of ABCs, denoted by ABC M 1 and ABC M 2 . These ABCs all coincide if the potential is time-independent. In the stationary case, all the potentials fall into this category and thus the ABCs are equivalent. Hence, we get the unique class of stationary ABCs, SABC M (without subscript index). For convenience, the form of the boundary conditions ( 6)-( 7) is based on ABC M 2 (we refer to [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] for more technical details).

Application to the computation of stationary states: the linear case

Let us consider the Hamiltonian

H H = -α d 2 dx 2 + V (x), x ∈ R, (8)
defined through α and V . The task here is to determine the pair (φ E , E) solution to the eigenvalue problem:

Hφ E = Eφ E , x ∈ R. (9)
This problem can also be formulated as follows: find the eigenvalues (E n ) n∈N (energies) and the associated real-valued eigenfunctions (φ n ) n∈N (eigenstates or ground states) as solutions of: Hφ n = E n φ n , x ∈ R. To fix the eigenfunction, it is necessary to impose a normalization condition: φ E L 2 (R) = 1. Let us begin with the case where the potential does not depend on the eigenfunction (called linear case here). The nonlinear case will be treated later in Section 4.

Square-root ABCs

Before discussing the difficulties related to the ABCs, let us consider the numerical solution of our problem with a homogeneous Dirichlet boundary condition. The variational formulation of (9) reads

-α[∂ n φ E ψ] xr x ℓ + α Ω ∂ x φ E ∂ x ψ dx + Ω V φ E ψ dx = E Ω φ E ψ dx, (10)
for some test-functions ψ ∈ H 1 0 (Ω) [START_REF] Chazarain | Introduction to the Theory of Linear Partial Differential Equations[END_REF]. Let S 0 , M 0 and M 0 V be respectively the stiffness matrix, mass and generalized mass matrices associated with the potential V for P 1 finite element and a homogeneous Dirichlet boundary condition (these matrices are some elements of M n h -1 (R)). The discrete problem can be classically formulated as the following generalized eigenvalue problem: find the pair (E, φ E ) as solution to

αS 0 + M 0 V φ E = EM 0 φ E , M 0 φ E 2 = 1, (11)
which is a generalized eigenvalue problem with an equality constraint. Here, φ E is a vector in R n h -1 which is normalized by: M 0 φ E 2 = 1 ( • 2 being the usual Euclidian norm in R n h -1 ). The global algorithm complexity is essentially the sum of the complexities for building the sparse finite element matrices and for computing the eigenvalue problem. In this paper, we use Matlab's eigs function which provides the p smallest positive eigenvalues corresponding to the generalized eigenvalue problem. This function automatically normalizes the eigenvectors in the Euclidian norm hence fulfilling the normalization constraint in [START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF]. eigs is associated with the software ARPACK. In the case where the potential is not always positive, we use the property that the smallest eigenvalue E 0 is larger than the minimum of the potential V min and solve [START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF] by a translation of -V min . Finally, the solution to [START_REF] Ben Abdallah | On a one-dimensional Schrödinger-Poisson scattering model[END_REF] generates the sequence of the p first eigenvalues (E 0 n ) 0≤n≤p-1 , eigenvectors (φ 0 n ) 0≤n≤p-1 and finite element eigenfunctions (φ 0 n ) 0≤n≤p-1 associated with the Dirichlet boundary condition. Since this eigenvalue problem is linear with respect to E, we can solve it without using e.g. a fixed point algorithm, unlike the case of including a square-root ABC as it is explained below. For this reason, the solution is called "direct" in the sequel of the paper.

Let us consider now the SABC 2 boundary condition

∂ n φ E = i √ α √ E -V φ E , on Σ. ( 12
)
The main difficulty with this boundary condition is its nonlinear dependence on E. As a consequence, we cannot isolate the terms (E, φ E ) in the right-hand side of (10) in a linear way, that is under the form Eφ E . More precisely, the nonlinear eigenvalue problem to solve is αS

+ M V + B M (E M ) φ M = E M Mφ M , Mφ M 2 = 1, (13)
using the matrix notations of the scattering problem. We precise that both the eigenvalues and eigenfunctions depend on the chosen boundary condition SABC M by the notation: (E M , φ M ). The first p eigencomponents are indexed as follows:

(E M n , φ M n ), with 0 ≤ n ≤ p-1.
The nonlinear dependence on the boundary term is given by the presence of B M (E M ). To solve the eigenvalue problem with SABC M , we have to apply an iterative scheme like a fixed point method (with a prescribed tolerance ǫ) and update E M at each iteration step j.

This procedure implies that we have to a priori choose an eigenvalue of index n (denoted by E M n ) that we wish to calculate. This is an important drawback since we have to a priori compute successively all the eigenvalues and associated eigenvectors. In fact, it appears that eigs is also able to provide an approximation of the first p eigenvalues (E M,j n ) 0≤n≤p-1 of (E M n ) 0≤n≤p-1 and the corresponding eigenvectors φ M n . As a consequence, we also have to recompute the boundary terms arising in B M (E M,j n ). Hence, the fixed point algorithm reads

αS + M V + B M (E M,j n ) φ M,j+1 = E M,j+1 Mφ M,j+1 , Mφ M,j+1 2 = 1, (14)
each linear problem being solved by using the Matlab routine eigs. More generally, for a boundary condition with a nonlinear dependence on the energy E, we use an associated fixed point algorithm. Even if we iterate through a fixed point algorithm, it appears that the algorithm also simultaneously gives some approximations of the other eigenvalues and eigenvectors (see the numerical section). This approach is therefore designated by "direct" if we only iterate on one a priori fixed eigenvalue. This algorithm can be applied successively by iteration using the fixed point algorithm and keeping only the computed eigenvalue and eigenvector related to the current iteration. Of course, the resulting algorithm is more expensive but at the same time more accurate. This approach is designated by "loop" in the sequel. Let us remark that there is no difference between both approaches for the Dirichlet problem.

Linearized ABCs

Unlike the case of the Dirichlet problem, we previously saw that the algorithm related to the square-root ABCs is iterative because of the nonlinearity. To avoid this problem, we can linearize SABC 2 and SABC 4 . The principle is based on a Taylor's expansion in the regime E ≪ V . This asymptotic regime is justified in particular for an harmonic potential V (x) = 1 2 x 2 since V grows quickly as soon as we do not place the boundary too close to the origin and we restrict our computations to relatively not too high energies. For the boundary condition SABC 2 [START_REF] Arnold | Mathematical concepts of open quantum boundary conditions[END_REF], this leads to the approximation of SABC 2 by SABC 2 lin given by

∂ n φE = - V ℓ,r √ α φE + 1 2 E √ α V ℓ,r φE . (15)
Next we can isolate the linear part according to E as

∂ n φE = β 2 ℓ,r φE + Eγ 2 ℓ,r φE , with β 2 ℓ,r
and γ 2 ℓ,r defined by [START_REF] Heinen | Radiation boundary conditions for the numerical solution of the three-dimensional time-dependent Schrödinger equation with a localized interaction[END_REF]. Including these ABCs in the weak formulation (10) leads, after discretization by the P 1 finite element method, to the following linear eigenvalue problem (M = 2)

(αS + M V + C M ) φM = ẼM (M + D M ) φM , M φM 2 = 1. (16)
We have defined the two matrices (M = 2)

C M =       αβ M ℓ 0 0 0 0 . . . . . . 0 0 . . . . . . 0 0 . . . . . . 0 0 0 0 αβ M r       , B M =       -αγ M ℓ 0 0 0 0 . . . . . . 0 0 . . . . . . 0 0 . . . . . . 0 0 0 0 -αγ M r       . ( 17
)
Problem ( 16) is directly solved without iteration by using eigs. The computational cost is therefore the same as for a Dirichlet boundary condition. Furthermore, since we do not have to iterate, the algorithm provides simultaneously the first p eigenvalues ( ẼM n ) 0≤n≤p-1 and associated eigenvectors ( φM n ) 0≤n≤p-1 . The resulting algorithm is called direct. In the case of SABC 4 , a similar strategy of linearization of (7) leads to the approximation

∂ n φE = β 4 ℓ,r φE + Eγ 4 ℓ,r φE , (18)
with β 4 ℓ,r and γ 4 ℓ,r respectively given by

β 4 ℓ,r = β 2 ℓ,r - 1 4 ∂ n V |x=x ℓ,r V ℓ,r , γ 4 ℓ,r = γ 2 ℓ,r - 1 4 ∂ n V |x=x ℓ,r V 2 ℓ,r (19)
by using the approximation 1 4

∂ n V E -V ≈ - 1 4
∂ n V V - 1 4
∂ n V V 2 E. ( 20
)
Adapting the functions, our problem can be written as [START_REF] Heinen | Numerical calculation of strong-field laser-atom interaction: An approach with perfect reflection-free radiation boundary conditions[END_REF].

Numerical examples

Example 1 (Harmonic potential). We first consider the well-known (positive) harmonic potential V (x) = 1 2 x 2 , i.e. the equation to solve is

- 1 2 φ ′′ E + 1 2 x 2 φ E = Eφ E , x ∈ R, (21)
with α = 1 2 . The square-integrable normalized solutions of ( 21) are the Hermite functions

φ ex n (x) = π -1/4 √ 2 n n! e x 2 /2 d n dx n e -x 2 , n ≥ 0 ( 22
)
and the corresponding eigenvalues (energies) are:

E ex n = n + 1 2 .
The eigenfunctions φ ex n (x) vanish for |x| → ∞, but this decay is slower and slower as n grows.

Let us recall that, for the case of the square-root ABCs, we have the direct and loop strategies. In the sequel, when we present an error calculation with respect to x r or h, this is clearly obtained by the direct approach since n is fixed. When we compute a range of eigenvalues (curves with n as abscissa), we report the results for both strategies to compare the respective accuracies.

A first numerical test consists in presenting the error on both the energy and eigenfunctions depending on the computational domain size. For the harmonic potential, we always consider a symmetric domain Ω =] -x r ; x r [. For a fixed n, the value of an eigenfunction is closer to zero as x r becomes larger. This means that we should observe the impact of the ABCs compared to the homogeneous Dirichlet boundary condition depending on the location of x r . Figure 1 reports, for the fundamental state n = 0 and in logarithmic scale, the absolute error on the eigenvalue

|∆E| = |E num n -E ex n | and the error in the L 2 -norm of the eigenfunction ||∆φ|| L 2 (Ω) = ||φ num n -φ ex n || L 2 (Ω)
when the right endpoint x r varies between 1 and 7.

Figure 2 presents similar results for n = 4 and x r varying between 3 and 10. The calculations are obtained for the numerical eigenvalues E num n equal to E M n (for SABC M ) or ẼM n (for SABC M lin ), depending on the order M of the ABC and its type (square-root or linearized). In the nonlinear case, corresponding to SABC M , the number of iterations is 50 to reach convergence with ε = 10 -12 . The spatial step size is h = 1 • 10 -3 .

For n fixed, we observe an accuracy improvement for both boundary conditions when x r grows. When x r is close to the origin (for example x r = 1 for n = 0, x r < 3 for n = 4), all the conditions lead to inaccurate results. However, even for these small values of x r , the linearized ABCs already give an approximation of the eigenvalue while this is not the case for the Dirichlet boundary condition as well as for SABC 2,4 . Indeed for x r = 3 and n = 4, the ABCs SABC 2,4 lin give Ẽn with an error equal to 10 -2 when the error for the homogeneous Dirichlet boundary condition is about 10 -1 and 1 for the square-root ABCs. The same remark holds for n = 0.

It seems from these tests that the linearized ABCs are the most robust boundary conditions concerning the size of the computational domain. From a general point of view, the ABCs always provide a better precision, at least about the same as with the Dirichlet boundary condition but often far better. The ABCs of different orders generally give a similar accuracy with however a better accuracy behaviour of the square-root ABCs but at a higher computational cost. They improve the accuracy of the Dirichlet boundary condition from a factor between 10 and 10 3 according to the configuration, before attaining the saturation zone. After a certain value of x r , all the boundary conditions lead to the same accuracy which only depends on the spatial mesh size. For the computation of the eigenfunctions, this value can be estimated to x r = 6 for n = 0 and to x r = 6.5 for n = 4.

We also remark that we must increase x r as n grows to get the same accuracy. To confirm this, we compute the variation of the error when Ω =] -x r ; x r [ is fixed and n varies. We set x r = 4 and for n ∈ [0, 10] we report the error on the eigenvalues E n (Figure 3) for both the "direct" and "loop" approaches.

For all the boundary conditions, we can clearly see that the accuracy decays as n increases. Indeed, the ABCs have been built in the high frequency regime. In our context, this means that we require that:

E n -V r ≪ 0 (23)
holds for a given point x r and for a fixed potential V . As a consequence, this limits the calculation of energies under the condition E n ≪ x 2 r /2 for example in the harmonic case. In the proposed simulation, setting x r = 4 leads to E n ≪ 8, which is coherent with the observations in Figure 3. Another way to interpret this property is that increasing the accuracy and the range of eigenvalues must be a priori guided by relation [START_REF] Frensley | Boundary conditions for open quantum driven far from equilibrium[END_REF]. To visualize this, we show in Figure 4 the potential V as well as the first energies E n . We can read from this figure the abscissa x where E -V becomes negative and we can have a first idea of the choice of the minimal abscissa x r to choose to get a sufficiently large gap between E and V (x r ) according to [START_REF] Frensley | Boundary conditions for open quantum driven far from equilibrium[END_REF].

For example, for the fundamental state n = 0, the energy associated with E 0 is the lowest level red curve. From the intersection with the curve of V (x), we can see that E 0 -V (x) is negative for x ≥ 1 and we can estimate that the difference between E 0 and V (x) will be enough starting from about x ≥ 2. Coming back to Figure 1 confirms these values since choosing x r = 1 provides a possible computation but does not necessarily converge towards E 0 while setting x r = 2 gives a correct approximation of E 0 . We can do the same analysis for n = 4 (fifth red curve from the bottom). We see that E 4 -V (x) is negative from x ≥ 3 and "very negative" after x ≥ 4. These values must be connected with the curves of Figure 2.

On this example (n = 4), we also remark that the linearized ABCs are more accurate than the original square-root ABCs, with a gain of a factor 10 in precision. This remark could also have been made on Figure 2 corresponding to n = 4, most particularly for the computation of the eigenvalue. The precision obtained for x r ≥ 4.5 with the linearized boundary conditions is the same as for the square-root boundary conditions but the linearization yields an accuracy improvement on smaller computational domains while the iterative algorithm for the squareroot conditions does not converge (2 ≤ x r ≤ 3). Moreover, let us note that the spectrum is simultaneously obtained in the linear case without iterating which is a crucial gain compared to the "loop" approach, showing hence the need of linearizing. As a consequence, the ABCs SABC M lin are, for a similar computational cost, to privilege to the Dirichlet boundary condition for accuracy purpose and/or for reducing the computational domain. Let us also finally remark that the gain in terms of accuracy of the "loop" approach is interesting as we can see it in Figure 2(a) but for a relatively higher computational complexity.

We now wish to compare the performances of the linearized and square-root ABCs. The previous curves illustrated the question of accuracy. Generally speaking, the square-root ABCs provide a better accuracy but at a higher computational cost even for the "direct" approach since a fixed point is required. We show in Figure 5 the number of iterations when using SABC 2 and SABC 4 , with respect to x r , for two situations: n = 0 and n = 4. Figures 5(a) and 5(b) must be connected to Figures 1 and2 which are their equivalent in terms of accuracy.

For the first value of x r , we often observe the divergence of the algorithm (the maximal number of iterations of the fixed point algorithm is 20). Again, this is one of the interesting property of the linear ABCs since, if we go back to Figures 1 and2, they also give a rough estimate of the eigenvalue.

For a slightly larger value of x r , the number of iterations stagnates to 5. Finally, when the maximal accuracy is reached, the algorithm needs 2 or 3 iterations. Globally, the computational costs for the square-root ABCs are roughly 5 times the costs for the linear ones and the Dirichlet boundary condition. At the same time, a higher accuracy is obtained by the square-root ABCs.

Finally, we present in Figures 6 (n = 0) and 7 (n = 4) the influence of the discretization on the accuracy for a given computational domain ] -x r ; x r [. We fix x r and report the errors |∆E| and ∆φ L 2 (Ω) depending to the mesh size h, for h between h = 5 • 10 -2 and h = 1 • 10 -4 . The value of x r is chosen such that the saturation of the error has not been reached yet so that we can see an effect of the ABCs compared to the Dirichlet boundary condition.

One remarkable property is that for n = 0 (Figure 6) the accuracy remains increasing with the ABCs by refining the mesh while this is not the case for the Dirichlet boundary condition. Indeed, we cannot gain more accuracy after h = 10 -2 if we do not increase the size of the computational domain. Concerning the ABCs (which are already more accurate than the Dirichlet boundary condition for h = 10 -2 ), we can still improve the solution by refining, most particularly with SABC 4 . Note that this remark holds for both the eigenvalues and eigenvectors.

Example 2 (Pöschl-Teller potential). The potential that we analyze now has the property to lead to negative eigenvalues. A necessary condition to justify the application of the previous approach is that

V (x r ) -E ≥ 0. ( 24
)
Hence, according to n and the rank of the eigenvalue that we are looking for, we have to choose x r sufficiently large so that condition (24) is fulfilled. Since V is negative and even if we have E ≪ V , then linearizing SABC 2,4 by using a Taylor's expansion with respect to E/V is no longer relevant since V can be equal to zero. Let us set V min = min x∈R V (x) and using the property that the Schrödinger equation is linear, we define a new positive potential W = V -V min and F n = E n -V min . Problem ( 8)-( 9) is then equivalent to

-αφ ′′ E + W φ E = F n φ E . ( 25
)
The boundary conditions SABC 2,4 lin are so the linearized versions of SABC 2,4 according to 1/(V -V min ) (and not 1/V ) by using the equivalent assumption:

E -V min ≪ V -V min .
The Pöschl-Teller potential [START_REF] Shao | Arbitrarily precise numerical solutions of the one-dimensional Schrödinger equation[END_REF] is given by

V (x) = - λ(λ + 1) cosh 2 (x) , (26)
and α = 1 in (8). This potential is always negative (see Figure 8(a)). For λ = 9, it leads to nine eigenvalues:

E n = -(9 -n) 2 , 0 ≤ n ≤ 8.
In Figure 8(b), we plot the different energy levels, compared to the potential. To take into account the translation, we rather present V (x) -V min and E n -V min .

For a given eigenstate, we can a priori estimate the size of the computational domain to consider that the high frequency hypothesis is satisfied and that the convergence of the iterative algorithm occurs. We set h = 5 • 10 -4 and analyze, for n fixed, 0 ≤ n ≤ 8, the error on E n for the different ABCs depending on the position of x r . We depict the results in Figure 9 for n = 0, n = 4 and n = 8. The ABCs always improve the accuracy compared to the Dirichlet boundary condition. This is most particularly clear for large n. For n = 8 and ] -5; 5[, the accuracy obtained with the Dirichlet boundary condition is less than 10 -2 and about 5 • 10 -5 for SABC 2,4 . To get a similar precision with the Dirichlet boundary condition, we would have to choose x r = 10 leading therefore to a significant larger computational domain. The effect of the linearized ABCs is variable. For n = 0 the ABCs SABC 2,4 lin are almost as precise as the ABCs SABC 2,4 , but when n increases, the accuracy is similar to the one obtained by using the Dirichlet boundary condition.

For the same potential, we observe in Figure 10 the error on E n for x r fixed. For x r = 2, we notice a factor 10 to 100 between the Dirichlet boundary condition and the ABCs for the states 2 ≤ n ≤ 6. The second-and fourth-order ABCs have a similar accuracy. At x r = 4, all the boundary conditions are equivalent for the first eigenstates but when n grows the ABCs remain accurate while the Dirichlet boundary condition is less precise (n = 6, n = 7). Indeed, they yield an accuracy of the eigenvalue about 10 -3 while the Dirichlet boundary condition gives only 10 -1 (n = 8).

Remark 2. In [START_REF] Antoine | Modeling boundary conditions for solving stationary Schrödinger equations[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 also study the case of the Morse potential as well as the Woods-Saxon potential. This extends our conclusions to these two situations.

Application to the computation of ground states: the nonlinear case

Problem and numerical scheme

We are interested in computing ground states for nonlinear Schrödinger equations. Most particularly, we consider a nonlinear potential which is the sum of a cubic nonlinearity and a harmonic potential. This kind of nonlinearity arises e.g. in Bose-Einstein condensates [START_REF] Gross | Structure of a quantized vortex in boson systems[END_REF][START_REF] Pitaevskii | Vortex lines in an imperfect Bose gas[END_REF][START_REF] Bao | The nonlinear Schrödinger equation and applications in Bose-Einstein condensation and plasma physics[END_REF]. The dimensionless one-dimensional Gross-Pitaevskii equation [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF][START_REF] Bao | Ground-state solution of Bose-Einstein condensate by directly minimizing the energy functional[END_REF][START_REF] Dehghan | A compact split-step finite difference method for solving the nonlinear Schrödinger equations with constant and variable coefficients[END_REF] reads

i ∂ψ ∂t = - 1 2 ∂ 2 x ψ + V ψ + β|ψ| 2 ψ, x ∈ R, (27)
setting V (x) = 1 2 x 2 and where the nonlinearity coefficient β can be negative or positive. We restrict ourselves to this special nonlinearity but all results can be directly extended to other cases. In view of computing the stationary solutions we write: ψ(x, t) = e -iEt φ E (x), where E is the chemical potential of the condensate and φ E is a real-valued function independent of time. Let us note that the stability of exactly this kind of problems was studied analytically in [START_REF] Thelwell | Instabilities of one-dimensional stationary solutions of the cubic nonlinear Schrödinger equation[END_REF][START_REF] Carr | Stationary solutions of the one-dimensional nonlinear Schrödinger equation. I. Case of repulsive nonlinearity[END_REF][START_REF] Carr | Stationary solutions of the one-dimensional nonlinear Schrödinger equation. II. Case of attractive nonlinearity[END_REF] and hence can be checked numerically using our proposed ABCs.

Function φ E is then solution to

-α∂ 2 x φ E + V φ E + β |φ E | 2 φ E = Eφ E , x ∈ R, (28)
where α = 1 2 , under the normalization constraint

φ E L 2 (R) = 1. ( 29
)
Finally, the function φ E of the problem ( 28)-( 29) satisfies the boundary conditions φ ′ E (0) = 0 and φ E (±x) → 0 for x → +∞. The resulting system is a nonlinear eigenvalue problem under constraint. The eigenfunction φ E being known, we can determine the associated eigenvalue E by

E = R α |∂ x φ E | 2 + V φ 2 E + βφ 4 E dx. ( 30
)
The problem ( 28)-( 29) is solved on a symmetric computational domain Ω =]-R; R[, with R > 0 and Σ = {-R; R}. We keep on denoting this domain by Ω =] -x ℓ ; x r [. We introduce (E 0 , φ 0 ) as a solution to the boundary value problem with Dirichlet boundary condition

     -α∂ 2 x φ E + V φ E + β|φ E | 2 φ E = Eφ E , in Ω, φ E = 0, on Σ, φ E L 2 (Ω) = 1. (31)
Analogously, we designate by (E M , φ M ) the solution computed with a M -th order nonlinear ABC obtained from the linear stationary ABCs ( 6)- [START_REF] Bao | A generalized-Laguerre-Hermite pseudospectral method for computing symmetric and central vortex states in Bose-Einstein condensates[END_REF]. To this end, we replace formally the potential V by the new nonlinear potential V + β|φ| 2 to get the second-order ABC

∂ n φ E = i √ α E -V -β|φ E | 2 φ E , on Σ, (32)
and fourth-order ABC

∂ n φ E = i √ α E -V -β|φ E | 2 φ E + 1 4 ∂ n (V + β|φ E | 2 ) E -V -β|φ E | 2 φ E , on Σ. ( 33
)
For the sake of clarity, we keep on designating by SABC M the above M -th order ABC. The interior equation is discretized by the semi-implicit scheme

-α∂ 2 x φ M,j+1 + V φ M,j+1 + β|φ M,j | 2 φ M,j+1 = E M,j+1 φ M,j+1 , (34)
for j ≥ 0 and M = 0, 2, 4. Now and independently of the boundary condition, the algorithm must be iterative since the interior scheme is nonlinear. As a consequence, we systematically use the fixed point method on the n-th eigenvalue E M n and eigenfunction φ M n for solving the eigenvalue problem. The variational formulation reads

-α[∂ n φ M,j+1 n ψ] xr x ℓ + α Ω ∂ x φ M,j+1 n ∂ x ψdx + Ω V φ M,j+1 n ψdx + β Ω |φ M,j n | 2 φ M,j+1 n ψdx = E M,j+1 n Ω φ M,j+1 n ψdx, (35)
for any test-function ψ. In the Dirichlet case, by choosing ψ ∈ H 1 0 (Ω), which makes the first term of the equation vanish, the discrete problem is, for

M = 0,    αS 0 + M 0 V + βM 0 |φ M,j n | 2 φ M,j+1 n = E M,j+1 n M 0 φ M,j+1 n , M 0 φ M,j+1 n 2 = 1. ( 36
)
For the ABCs, we use for ∂ n φ M,j+1 n the fixed point version

∂ n φ M,j+1 n = i √ α E M,j n -V -β|φ M,j n | 2 φ M,j+1 n ( 37
)
for the second-order ABC [START_REF] Bao | The nonlinear Schrödinger equation and applications in Bose-Einstein condensation and plasma physics[END_REF] and

∂ n φ M,j+1 n = i √ α E M,j n -V -β|φ M,j n | 2 + 1 4 ∂ n (V + β|φ M,j n | 2 ) E M,j n -V -β|φ M,j n | 2 φ M,j+1 n ( 38
)
for the fourth-order condition [START_REF] Bao | Computing the ground state solution of Bose-Einstein condensates by a normalized gradient flow[END_REF]. Hence, the term -α[∂ n φ M,j+1 n ϕ] xr x ℓ leads, from a discrete point of view, to a matrix contribution B j M φ M,j+1 for the M -th order ABC, where the matrix coefficients B j M only depend on the values of φ M,j n and E M,j n . By applying the fixed point algorithm on the n-th eigenvalue E M n and eigenvector leads to the iterative scheme

φ M n    αS -αB M,j + M V + βM |φ M,j n | 2 φ M,j+1 n = E M,j+1 n Mφ M,j+1 n , Mφ M,j+1 n 2 = 1. (39)
The matrix coefficients B M,j are given by

(B M,j ) 1,1 = i √ α E M,j n -V ℓ -β|φ M,j n,ℓ | 2 + 1 4 ∂ n (V + β|φ M,j n | 2 ) |x=x ℓ E M,j n -V ℓ -β|φ M,j n,ℓ | 2 (40)
and

(B M,j ) n h +1,n h +1 = i √ α E M,j n -V r -β|φ M,j n,r | 2 + 1 4 ∂ n (V + β|φ M,j n | 2 ) |x=xr E M,j n -V r -β|φ M,j n,r | 2 (41)
for SABC 4 (M = 4). For SABC 2 (M = 2), it is sufficient to retain only the first term of each of the above expressions. We have set here: φ M,j n,ℓ = φ M,j n,|x=x ℓ and φ M,j n,r = φ M,j n,|x=xr . As in the linear case, we can formulate the linearized versions of the second-and fourthorder ABCs. These ABCs are then designated by SABC 2,4 lin . Doing so, we have the secondorder ABC

∂ n φ M,j+1 n = - √ V √ α φ M,j+1 n - β 2 1 √ α √ V |φ M,j n | 2 φ M,j+1 n + 1 2 1 √ α √ V E M,j+1 n φ M,j+1 n (42)
and the fourth-order ABC

∂ n φ M,j+1 n = - √ V √ α φ M,j+1 n - β 2 1 √ α √ V |φ M,j n | 2 φ M,j+1 n + 1 2 1 √ α √ V E M,j+1 n φ M,j+1 n + - 1 4
∂ n (V + β|φ M,j n | 2 ) V + β 4 |φ M,j n | 2 ∂ n (V + β|φ M,j n | 2 ) V 2 φ M,j+1 n -E M,j+1 n ∂ n (V + β|φ M,j n | 2 ) V 2 φ M,j+1 n . (43)
The iterative scheme then reads

   αS -αB M,j + M V + βM |φ M,j n | 2 φ M,j+1 n = E M,j+1 n M + αB j E,M φ M,j+1 n Mφ M,j+1 n 2 = 1. (44)
The matrix coefficients B M,j et B M,j E are given by

(B M,j ) 1,1 = - √ V ℓ √ α - β 2 1 √ α √ V ℓ |φ M,j n,ℓ | 2 - 1 4
∂ n (V ℓ + β|φ M,j n,ℓ | 2 ) V ℓ + β 4 |φ M,j n,ℓ | 2 ∂ n (V ℓ + β|φ M,j n,ℓ | 2 ) V 2 ℓ (45) and (B M,j E ) 1,1 = 1 2 1 √ α √ V ℓ - ∂ n (V ℓ + β|φ M,j n,ℓ | 2 ) V 2 ℓ (46)
for the fourth-order ABC. The expression of the coefficients of index (n h + 1, n h + 1) is the same but taking its value at x = x r . We can easily extract the coefficients associated with the second-order ABC by keeping only the first term of each expression. Unlike the linear situation, there is no gain in terms of computational time since the problem is fully nonlinear.

Numerical results

We consider [START_REF] Shao | Arbitrarily precise numerical solutions of the one-dimensional Schrödinger equation[END_REF] for different values of the parameter β. For each value, we uniquely determine the fundamental state n = 0. To get some reference eigenvalues, we numerically compute them on the domain ] -30; 30[, with a step size h = 10 -4 and SABC 2 (33). This method provides some values reported in Table 1 which are conform with the ones given in [START_REF] Bao | Ground-state solution of Bose-Einstein condensate by directly minimizing the energy functional[END_REF]. Let us note here that we do not give some results for larger values of β because the fixed point algorithm then diverges. It would be necessary at this point to use another numerical algorithm (a Newton method or a continuation method) for solving the problem with an ABC. Finally, we present in the sequel the absolute errors: For the simulations, the initialization of the fixed point algorithm uses the exact harmonic potential solution (β = 0): φ M,0 0 (x) = 1 π 1/4 e -x 2 /2 . The fixed point algorithm tolerance is ε = 10 -12 and the mesh size of the linear finite element method is h = 10 -3 . Figures 11,12 and 13 report the error on both the eigenvalue and eigenfunction at the origin depending on the right endpoint x r , for the values β = -6.2742, β = -2.5097 and β = 3.1371, respectively.

∆E = |E num -E ref | and ∆φ(0) = |φ num (0) -φ ref (0)|,
Generally speaking, for a given case, all the algorithms converge with about the same number of iterations, independently of the boundary condition. We also note that, for negative values of β, the linearized ABCs lead to the same accuracy as the nonlinear ABCs (not reported here) for a similar computational time. In Figures 11 and12, we only present the results for the Dirichlet boundary condition and SABC 2,4 . For β > 0 (Figure 13), the linearized ABCs possess an accuracy at least equal to the one with SABC 2,4 . We do not have any explanation about this fact. For β = -6.2742 (Figure 11), all the algorithms converge in 23 iterations. The ABCs improve the accuracy from a factor 10 compared with the Dirichlet boundary condition for x r = 1.5, and almost 100 when x r = 2, then for x r ≥ 2.5, all the boundary conditions have the same accuracy: 10 -5 . The precision of the secondorder ABC is slightly better than the fourth-order ABC. There is no explanation here about this behavior but we recall that ABCs for the nonlinear case where obtained by the formal argument: "the potential is replaced by the nonlinearity". For β = -2.5097 (Figure 12), the convergence takes 14 iterations. The ABCs again provides a gain of precision compared with the Dirichlet boundary condition for x r between 1.5 and 3.5, with a better accuracy for the second-order ABC (see the points x r = 3 and x r = 2.5 for example). For β = 3.13712, the situation is quite similar but requires 77 iterations to converge. Unlike, the two previous cases, the linearized ABCs give a slightly better accuracy than for the original ABCs.

Extension to N -dimensional stationary problems

In the N -dimensional linear time-dependent case, the ABCs take the following form

ABC 1 2 : ∂ n u -iOp( -τ + ∆ Σ + V ) u = 0 (47)
for the first-order condition and

ABC 2 2 : ∂ n u -iOp( -τ + ∆ Σ + V ) u +HOp((-τ + ∆ Σ + V ) -1 ) i∂ t + V u = 0 (48)
for the second-order one on Σ × R + . The Schrödinger equation under consideration is the following

i∂ t u + ∆u + V (x)u = 0, ∀(x, t) ∈ R N × R + , u(x, 0) = u 0 (x), x ∈ R N , (49)
Here, the computational domain Ω is a bounded set of the R N dimensional space. Its boundary Σ is supposed to be a (N -1)-dimensional convex and compact manifold. Its mean curvature H is defined by: H = trace(R)/(N -1), where R is the curvature tensor of the surface. For example, for a two-dimensional surface, we have: H = κ/2, where κ is the local curvature at a point of the surface. For a 2-sphere of radius R, we get H = 1/R. The operator ∆ Σ is the Laplace-Beltrami operator over the surface. For the sphere it is defined by

∆ Σ f := 1 r 2 sin ϕ ∂ ∂ϕ sin ϕ ∂f ∂ϕ + 1 r 2 sin 2 ϕ ∂ 2 f ∂θ 2 .
for a function f expressed in spherical coordinates (r, ϕ, θ). Other expressions exist for the hypersphere and other manifolds. The asymptotics behind the construction of the ABCs is somewhat different which results in a different definition of the order of an ABC. We do not develop the whole theory 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 in the 2-dimensional case.

To get the ABCs for the N -dimensional stationary case

-α∆u + V (x)u = Eu,
in Ω, we make the substitutions: -τ → E/α and V → -V /α. This leads to the first-order Stationary ABC:

SABC 1 : ∂ n u + 1 √ α V -(α∆ Σ + E) u = 0 (50)
and

SABC 2 : ∂ n u + 1 √ α V -(α∆ Σ + E) u + H(V -E)(V -(α∆ Σ + E)) -1 u = 0 (51)
As in the one-dimensional case, both conditions are nonlinear with respect to the energy E. Furthermore, the square-root involves now the surface Laplace-Beltrami operator ∆ Σ . Here, we propose the formal asymptotics: V ≫ α∆ Σ + E which can be justified by theoretical arguments of operator theory [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]. Using a second-order Taylor expansion, we obtain the approximate linearized SABC SABC 1 lin :

∂ n φE + √ V √ α φE - √ α 2 √ V ∆ Σ φE + √ α 2 √ V E φE = 0 (52)
on Σ. For the second-order SABC, we do not really have to linearize the corrective term which can be considered linearly through the introduction of an auxiliary function Ψ E . More precisely, we have

SABC 2 lin : ∂ n φE + √ V √ α φE - √ α 2 √ V ∆ Σ φE + √ α 2 √ V E φE + H(V -E)Ψ E = 0 (53)
coupled to the surface equation: -α∆ Σ Ψ E + (V -E)Ψ E -φE = 0. With such a trick, the coupled system with unknowns ( φE , Ψ E ) remains linear and is well-adapted to a symmetrical weak formulation for instance. Moreover, each term can be very easily implemented in usual numerical codes based for example on finite difference, finite element or spectral methods. The adaptation to the nonlinear stationary Schrödinger equation (for a smooth nonlinearity like f (u) = β|u| 2 ): -α∆u + V (x)u + f (u)u = Eu, in Ω, can be made by simply replacing V (x) by V (x) + f (u) in the above ABCs.

Conclusion

We have proposed some accurate and physically admissible absorbing boundary conditions for modeling linear and nonlinear stationary Schrödinger equations with variable potentials. Based on numerical schemes, these boundary conditions have been validated for many configurations including linear and nonlinear ground-state computations. Furthermore, the extension to N -dimensional problems that can be used in scattering problems like [39, Section 12.1], is given.

Future extensions would include variable mass Schrödinger equations [START_REF] Betcke | Analysis and efficient solution of stationary Schrödinger equation governing electronic states of quantum dots and rings in magnetic feld[END_REF] among others. It might also be valuable to extent the presented work to systems of Schrödinger equations that arise as so-called multiband effective mass approximations (MEMAs) to model electronic states in modern semiconductor nanostructures, cf. [START_REF] Zisowsky | Discrete transparent boundary conditions for transient k • p-Schrödinger equations with application to quantumheterostructures[END_REF][START_REF] Klindworth | Discrete transparent boundary conditions for multiband effective mass approximations[END_REF][START_REF] Odermatt | Bandstructure calculation using the k • p method for arbitrary potentials with open boundary conditions[END_REF]. Let us finally remark that applications to generalized Schrödinger equations could also be developed by adapting the methods developed in [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF][START_REF] Antoine | Numerical schemes for the simulation of the two-dimensional Schrödinger equation using non-reflecting boundary conditions[END_REF].

Figure 1 :

 1 Figure 1: Example 1: Error (n = 0).

Figure 2 :Figure 3 :

 23 Figure 2: Example 1: Error (n = 4).

Figure 4 :

 4 Figure 4: Example 1: Harmonic potential and the energies E n for 0 ≤ n ≤ 9.

= 4 Figure 5 :

 45 Figure 5: Example 1: Number of iterations needed for the algorithms associated with the different ABCs, with respect to x r for n = 0 and n = 4.

Figure 6 :

 6 Figure 6: Example 1: Error (x r = 3.5 and n = 0) with respect to h.

Figure 7 :

 7 Figure 7: Example 1: Error (x r = 4 and n = 4) with respect to h.

Figure 8 :

 8 Figure 8: Example 2: Pöschl-Teller potential and its first nine energy levels.

n = 8 Figure 9 :Figure 10 : 4 (Figure 11 :

 8910411 Figure 9: Example 2: Error of the eigenvalues for the different ABCs and the Pöschl-Teller potential.

Figure 12 :

 12 Figure 12: Errors ∆E and ∆φ(0) for β = -2.5097.

Figure 13 :

 13 Figure 13: Errors ∆E and ∆φ(0) for β = 3.1371.

Table 1 :

 1 where 'ref' refers to the values in Table 1 and 'num' to the ones computed with the proposed method. Numerical values E ref and φ ref E (0) computed on a larger domain for different β.

	β	φ E (0)	E
	-6.2742 1.265512713848083 -4.956873352670034 -2.5097 0.913230941756339 -0.806257128073956 3.1371 0.645961493829006 1.526594842533555

(Matthias Ehrhardt) 1 Supported partially by the French ANR fundings under the project MicroWave NT09 460489 (http://microwave.math.cnrs.fr/).


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

C Besse

P Klein

Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme

Keywords: Complex Absorbing Potentials, Absorbing Boundary Conditions, Perfectly Matched Layers, nonlinear Schrödinger equation, unbounded domain, relaxation scheme

The aim of this paper is to compare different ways for truncating unbounded domains for solving general nonlinear one-and two-dimensional Schrödinger equations. We propose to analyze Complex Absorbing Potentials, Perfectly Matched Layers and Absorbing Boundary Conditions. The time discretization is made by using a semi-implicit relaxation scheme which avoids any fixed point procedure. The spatial discretization involves finite element methods. We propose some numerical experiments to compare the approaches.

Introduction

We propose and compare some numerical schemes to solve the general Schrödinger equations in unbounded media

i∂ t Ψ ∞ (x, t) + ∆Ψ ∞ (x, t) + V (x, t) Ψ ∞ (x, t) + f (|Ψ ∞ |)Ψ ∞ (x, t) = 0, (x, t) ∈ R N ×]0, T [, Ψ ∞ (x, 0) = Ψ 0 (x), x ∈ R N .
(1) The wave function Ψ ∞ is defined in the unbounded domain R N (N ≥ 1). In view of a numerical computation, different solutions may be used. For example, one usual scheme consists in splitting the Laplacian and potential-nonlinear parts of the equation and next solving the first linear equation e.g. by FFT methods and exactly integrating the second nonlinear one (see e.g. [1]). This kind of scheme is efficient and accurate if the solution remains confined within the computational domain (for instance for solving the Gross-Pitaevskii equations). Indeed, then periodic boundary conditions may be applied for the Fourier solution since the wave vanishes on the boundary of a large enough computational domain. In the same situation, other schemes may be used as for example Crank-Nicolson schemes, Runge-Kutta methods in time and finite difference or finite element in space. Spectral techniques may be applied too (see [START_REF] Bao | The nonlinear Schrödinger equation and applications in Bose-Einstein condensation and plasma physics[END_REF] for more details about some of these techniques). Independently of the numerical discretization, one common problem arises when the solution does not remain inside the computational domain. This is for example the case for the defocusing nonlinear cubic Schrödinger equation, for linear Schrödinger equation with laser ionization of a one-dimensional helium atom [3], for strong field laser atom interaction [4]... Then, it is well-known that Dirichlet or periodic boundary conditions on the boundary of the computational domain are not adapted. Our goal here is not to focus on all the physical situations which can arise but rather to propose some different ways of truncating accurately the computational domain for (1) and compare them numerically.

Concerning truncation methods, the case of the free Schrödinger equation is now mastered and many possible solutions can be developed. We refer to [5] for an overview of the techniques. Considering now the linear Schrödinger equation with potential requires more developments. For example, time dependent but space independent potentials can be considered easily by gauge change and be treated as the free-space case. The situation of a space variable potential is much more complex. In some exceptional cases, explicit exact boundary conditions may be written at the fictitious boundary. However, in most situations, approximate boundary conditions must be derived. These boundary conditions are usually called Absorbing Boundary Conditions (ABCs) since they try to absorb waves striking the nonphysical boundary. We refer to [5] for such examples in the one-dimensional case. In the two-dimensional case, only a few solutions can be found for the free-and potential cases [5]. In the nonlinear case which is much more complicate, the ABCs at hand are often formally built from the linear case with potential. To the best of our knowledge, only a few papers propose some ABCs [6,7,8,9,10,11]. For the two-dimensional case, the only strategies for simulating ABCs have been proposed in [11,12]. We propose and numerically test in this paper some new ABCs for (1) for the one-and two-dimensional cases which are related to the ABCs for potentials developed in [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][START_REF] Antoine | [END_REF]. A very common method in physics is related to Complex Absorbing Potentials (CAPs) [START_REF] Muga | [END_REF]. The idea which is physically natural consists in adding to the linear Schrödinger equation a complex absorbing potential to damp the incoming wave in a surrounding layer. We try to extend this approach here to (1). As we will see, this approach fails to work and generates large reflections. In particular, the choice of the absorbing potential and its parameters is non trivial and extension to a nonlinear problem does not appear as natural. A related technique not analyzed here is the method of Exterior Complex Scaling (see examples in [16,17]). A closely powerful method introduced by Bérenger [18] for electromagnetic waves is the Perfectly Matched Layer (PML) approach. The method introduces dissipation but inside the Laplacian term and not the potential term. We apply this technique here [19] to (1) to show its accuracy. It appears that the accuracy that can be expected from the ABCs and PMLs approaches is about the same, generally showing an error of reflection of the order of 0.1% or less, and is therefore useful for practical computations. For an easier implementation of all the truncation techniques, we use a semi-implicit relaxation scheme [20] which leads to a flexible implementation of the method. In particular, it does not require any iteration like in a fixed point or Newton procedure for the nonlinearity. Since the PML and ABCs approaches are accurate for the one-dimensional case, we next introduce their extension to two-dimensional problems. We detail the discretization issues and analyze their accuracy in the case of the propagation of a soliton in a cubic media.

The plan of the paper is the following. In Section 2, we introduce the ABCs, CAPs and PMLs techniques for the one-dimensional equation (1). We propose some schemes based on relaxation techniques coupled to Finite Element Methods (FEM). Finally, we numerically test and compare the different approaches. In the third Section, we extend our methods and discretizations to the two-dimensional nonlinear Schrödinger equation. Some numerical simulations confirm the accuracy of the methods. Finally, Section 4 gives a conclusion. Let us note here that the codes corresponding to Section 2 and 3 can be downloaded freely at http://microwave.math.cnrs.fr/code/index.html if the reader wants to know more about the implementation issues of all these techniques.

One-dimensional nonlinear problems

Absorbing Boundary Conditions (ABCs)

The first approach that we investigate concerns absorption at the boundary. We consider the time-dependent one-dimensional nonlinear Schrödinger equation with a variable potential and a nonlinear term

i∂ t Ψ(x, t) + ∂ 2 x Ψ(x, t) + V (x, t) Ψ(x, t) + f (|Ψ|)Ψ(x, t) = 0, (x, t) ∈ Ω T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω. ( 2
)
We assume here that Ω :=]x , x r [ represents a bounded computational domain of boundary Γ := ∂Ω := {x , x r } and set Ω T := Ω×]0; T [, Σ t := Σ×]0; T [ (see Figure 1). Furthermore, Ψ 0 is supposed to be a compactly supported initial data inside Ω. If the potential V and nonlinear interaction f are constant outside Ω, we respectively call them localized potential and interaction. Then exact absorption at the boundary Σ can be obtained. To write this boundary condition (also called transparent), let us assume that the potentials are equal to zero outside Ω. Then, it is now standard that the boundary condition is given by

∂ n Ψ(x, t) + e -iπ/4 ∂ 1/2 t Ψ(x, t) = 0, on Σ T , ( 3
)
where n is the outwardly directed unit normal vector to Σ. The operator ∂ 1/2 t is the half-order derivative operator

∂ 1/2 t Ψ(x, t) := ∂ t t 0 1 √ π Ψ(x, ) √ t - d (4)
It can be proved that system (2)-( 3) is well-posed in a mathematical setting and that we have the energy bound

Ψ(•, t) 0,Ω ≤ Ψ 0 (•) 0,Ω , (5)
for any time t > 0, where u(•, t) 0,Ω designates the 2-norm of Ψ over Ω

Ψ(•, t) 0,Ω = ( Ω |Ψ(x, t)| 2 dx) 1/2
which can also be interpreted as the probability of finding Ψ in Ω and translates the absorbing property of the boundary condition. The boundary condition (3) is exact in the sense that there is no reflection back into the computational domain. Mathematically, this implies that the solution to (2)-( 3) is strictly equal to the restriction Ψ ∞ |Ω T , solution to (1). In the case of unbounded potential and interaction, then the situation is much more complex. Essentially, the possibility of writing the exact boundary condition (3) is related to the fact that for localized interactions, the Laplace transform can be used in the left and right exterior domains to write down the boundary condition through the Green function. In the case of nonlocal interactions, this no longer possible. Except in some special situations of potentials (e.g. linear potential) and nonlinearity (essentially integrable systems like the cubic case), it is impossible to get the exact expression of the absorption conditions. Here, we present without any mathematical details which are too cumbersome the boundary conditions that can be set at the boundary. We 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 the mathematical details. Essentially, the derivation is based on high-frequency asymptotic expansions in the Laplace domain using the extended theory of pseudodifferential operators. The resulting boundary conditions are no longer exactly non-reflecting. They are then called Absorbing (and not transparent) Boundary Conditions (ABCs), and we need to precise the order related to their asymptotics with the aim of measuring the a priori accuracy of the boundary condition. The ABC of order two is given by

∂ n Ψ -i i∂ t + V + f (|Ψ|) Ψ = 0, (6)
on Σ T . The square-root operator of i∂ t + V + f (|Ψ|) means that we consider the spectral square-root decomposition of this operator. The resulting operator is nonlocal but will be localized later for the numerical purpose through Padé approximants. Higher-order ABCs can be derived [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] but will not be tested here. We cannot expect that an ABC works well for any potential. In practical computations and for f = 0, one physically admissible assumption is that the potential is repulsive which means that V : R × R + is smooth and that we have

x∂ x V (x, t) > 0, ∀(x, t) ∈ Ω ,r × R + , where Ω :=] -∞; x ] and Ω r := [x r ; ∞[. An example is V (x, t) = β 2 |x| a , with 0 < a ≤ 2.
The assumption on the nonlinearity is not clear most particularly when a potential is added. The only "intuitive" assumption is that the solution Ψ is outgoing to the bounded domain and that no nonlinear or potential effect makes it reflecting back into Ω, which is a priori difficult to check because mathematically hard to write. Finally, we can prove [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] that (5) still holds for the second-order ABC (6) and for a time independent potential V (x) which translates the absorbing property of the boundary condition.

Complex Absorbing Potential (CAP), Exterior Complex Scaling and Perfectly Matched

Layers (PMLs) Another useful and widely studied approach for computing solutions to time-dependent Schrödinger equations with a potential term by using an absorbing domain is first the technique of Complex Absorbing Potential (CAP). Essentially, the idea consists in introducing a complex potential in the exterior domain to absorb the travelling wave. Mathematically, this consists in adding a spatial potential -iW in some exterior layers Ω =]x p , x [ and Ω r =]x r , x rp [ (see Figure 2). Of course, to coincide with the solution to (2), W is required to be zero in Ω and with a positive real part in the layers to damp the incoming wavefield. Another way to analyze this approach is called the Exterior Complex Scaling approach which consists in interpreting the introduction of the complex potential as the complex scaling: x → xe iθ , where θ is a rotation angle which must be correctly chosen. Extension includes the Smooth Exterior Scaling approach [16,17]. From the numerical point of view, the CAP approach is direct to code since we have to solve

i∂ t Ψ(x, t) + ∂ 2 x Ψ(x, t) + V (x, t) Ψ(x, t) -iW (x) Ψ(x, t) + f (|Ψ|)Ψ(x, t) = 0, (x, t) ∈ Ω ext T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω ext , (7) in the extended domain Ω ext := Ω ∪ Ω ∪ Ω r =]x p , x rp [ with boundary Σ ext = {x p , x rp }.
Here, according to [16,21], we consider the quadratic profile

W (x) =      W 0 δ -2 (x -x ) 2 , x p < x < x , 0, x < x < x r , W 0 δ -2 (x -x r ) 2 , x r < x < x rp , (8)
for a real positive value of W 0 . The thickness δ of the layer is δ = |x rp -x r | = |x p -x |. Other choices include exponential type absorbing functions [16]. At the boundary points x p and x rp , a boundary condition must be imposed. Here, we consider the classical homogeneous Dirichlet boundary conditions Ψ(x, t) = 0 at x p and x rp . However, a suitable extension does not seem direct for nonlinear problems as we will see later.

We concentrate now on another closely related approach called the Perfectly Matched Layers (PMLs) method which was introduced by Bérenger [18] for Maxwell equations. The idea consists in introducing a complexification of the derivative operator through damping in the extended domains Ω and Ω r . In the case of the nonlinear Schrödinger equation, this can be written down as

i∂ t Ψ + 1 S(x) ∂ x 1 S(x) ∂ x Ψ + V (x, t)Ψ + f (|Ψ|)Ψ = 0, (9)
in Ω ext (Figure 2). Function S is given by S(x) := 1 + Rσ(x). The layer parameters R and σ must be chosen carefully. Optimization techniques and adaptive discretizations can be developed [5]. Here, we will use the parameter values derived in [19], i.e. R = e iπ/4 and σ is the quadratic function

σ(x) =      σ 0 δ -2 (x -x ) 2 , x p < x < x , 0, x < x < x r , σ 0 δ -2 (x -x r ) 2 , x r < x < x rp . ( 10
)
The distance δ := |x rp -x r | = |x -x p | (that we take equal on both sides here for simplicity) corresponds to the thickness of the left and right absorption regions of the computational domain. Again, we fix the homogeneous Dirichlet boundary conditions: Ψ(x, t) = 0 at x p and x rp . It is interesting to note the close form of CAP and PML approaches even if they lead to different equations to solve.

Unlike the ABCs, both CAP, ECS and PMLs must be adapted and optimized according to each situation. They have the advantage of being easy to code but at the price of an extended domain of computation Ω ext where the potential V must be known. This is not always the case if V is given only numerically for instance instead of analytically. One of the advantages of both CAP and PML methods is that they are local in time while it is not a priori the case for the ABCs which include a nonlocal square-root operator. However, this drawback can be removed by using local approximations as seen later. Finally, all the techniques are derived in the linear situation and their application to nonlinear problems is formal. Therefore, their accuracy must be prospected.

Discretization schemes

To compute the solution of the previous problems, we have to introduce some numerical discretizations. First, we have to deal with the time discretization in the interior domain. One widely used scheme is the Crank-Nicolson scheme which reads

2i Φ n+1 ∆t + ∂ 2 x Φ n+1 + W n+1 Φ n+1 + f (|Ψ n+1 |) + f (|Ψ n |) 2 Φ n+1 = 2i Ψ n ∆t ( 11
)
for 0 ≤ n ≤ N -1, and a time step ∆t = T /N . Here, we set

W n+1 = V n+1 + V n 2 , Φ n+1 = Ψ n+1 + Ψ n 2 . (12)
A more adapted scheme to the computation of soliton solutions is the Duran-Sanz-Serna scheme [22]. Unfortunately, in both cases, the scheme remains nonlinear and a fixed point or a Newton algorithm is required. A few iterations are then necessary increasing the global computational cost of the procedure. Instead of using these numerical methods, we can consider the Besse relaxation scheme in [20]. Applied to the nonlinear equation i∂ t Ψ + ∂ 2

x Ψ + V Ψ + f (|Ψ|)Ψ = 0, on Ω T , the method leads to the solution to

i∂ t Ψ + ∂ 2 x Ψ + V Ψ + ΥΨ = 0, on Ω T , Υ = f (|Ψ|), on Ω T . (13)
Then, system ( 13) is discretized as

       2i Φ n+1 ∆t + ∂ 2 x Φ n+1 + W n+1 Φ n+1 + Υ n+1/2 Φ n+1 = 2i Ψ n ∆t , Υ n+3/2 + Υ n+1/2 2 = f (|Ψ n+1 |), (14)
for 0 ≤ n ≤ N -1. The initialization of Υ is chosen as Υ -1/2 = f (|Ψ 0 |). We can see that this time discretization avoid any additional iterative procedure for the nonlinear term and so is well adapted to nonlinear problems. Furthemore, this scheme is known to preserve many invariants for the nonlinear Schrödinger equation [20].

In the case of the CAP, ECS and PML techniques, the relaxation scheme directly applies on Ω ext . In the ABCs case, we have to discretize correctly the square-root operator. To this aim, we can approximate this operator by using Padé approximants R m and Υ

∂ n Ψ -iR m i∂ t + V + Υ Ψ = 0. The function R m is defined by R m (z) = a m 0 + m k=1 a m k z z + d m k = m k=0 a m k - m k=1 a m k d m k z + d m k , (15)
where the coefficients a m k and d m k are given by a m 0 = 1, and for 1

≤ k ≤ m a m k = 1 m cos 2 (2k-1)π 4m , d m k = tan 2 (2k -1)π 4m . (16)
We can explicit this relation (formally with z = i∂ t + V + Υ) by

∂ n Ψ -i m k=0 a m k Ψ - m k=1 a m k d m k i∂ t + V + Υ + d m k -1 Ψ = 0. ( 17
)
Then we introduce some auxiliary functions (ϕ k ) 1≤k≤m defined by

ϕ k = (i∂ t + V + Υ + d m k ) -1 Ψ
for the square-root approximation. The ABC then becomes

∂ n Ψ -i m k=0 a m k Ψ + i m k=1 a m k d m k ϕ k = 0,
and the time discretization is

∂ n Φ n+1 -i m k=0 a m k Φ n+1 + i m k=1 a m k d m k ϕ n+1/2 k = 0. ( 18
)
The definitions of the auxiliary functions lead to the coupled differential equations

i∂ t ϕ k + (V + V + Υ)ϕ k + d m k ϕ k = Ψ, x = x ,r .
These are discretized by ϕ

n+1/2 k = ϕ n+1 k +ϕ n k 2 with 2i ∆t ϕ n+1/2 k + (W n+1 + Υ n+1/2 )ϕ n+1/2 k + d m k ϕ n+1/2 k = Φ n+1 + 2i ∆t ϕ n k , x = x ,r .
This gives an explicit expression for the auxiliary functions at x = x ,r

ϕ n+1/2 k = 1 2i ∆t + W n+1 + Υ n+1/2 + d m k Φ n+1 + 2i ∆t 2i ∆t + W n+1 + Υ n+1/2 + d m k ϕ n k . (19)
Injecting ( 19) into ( 18), we get an explicit expression of ∂ n Ψ n+1 in terms of Ψ n+1 and other updated functions. To discretize the problem in space, we use the weak formulation associated to ( 14) which writes down for any test function ϕ,

2i ∆t Ω Φ n+1 ϕdΩ - Ω ∂ x Φ n+1 ∂ x ϕ dΩ + Σ ∂ n Φ n+1 ϕ dΣ + Ω (W n+1/2 + Υ n+1/2 )Φ n+1 ϕ dΩ = 2i ∆t Ω Ψ n ϕ dΩ, (20)
with the update of Υ on Ω

Υ n+3/2 = 2f (|Φ n+1 |) -Υ n+1/2 . (21)
We can see here that the nonlinearity is just included as a potential and therefore no fixed point or Newton method is necessary since the nonlinearity is explicit through the update (21). This is also definitively the way the boundary condition (18) on Σ is treated. The ABCs are then just a Fourier-Robin boundary condition which can be easily implemented in the finite element code. The formulation (20) can then be solved through (high-order) finite element methods, resulting in the solution of a linear system at each time step. In the case of a CAP, ECS or PML, the adaptation is direct by integrating [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] over Ω ext and by using the homogeneous Dirichlet boundary condition. In the sequel, we use linear finite element methods. The domain is decomposed into n h spaced elementary segments of size h.

Other spatial discretization methods could be used such as the finite difference or spectral techniques.

Numerical examples 2.4.1. Soliton propagation in a cubic media

The first test case concerns the propagation of a soliton in a cubic nonlinear media (V = 0 and f

(|Ψ|) = |Ψ| 2 ) Ψ exact (x, t) = 2 sech √ 2(x -k 0 t) exp i k 0 2 (x -k 0 t) + i(2 + k 2 0 /4)t (22)
for the wavenumber k 0 = 15. The computational domain is [-10; 10] and the maximal time of computation is T = 2. We take the time step ∆t = 10 -3 and n h = 2000 points for the uniform spatial discretization. For the CAP, the size of the layer is δ = 4 and W 0 = 10 for n h = 2800. The same domain is considered for the PML with σ 0 = 10. For the ABC, we take m = 50 Padé functions. We represent the amplitude of Ψ in normal but also logarithmic scale to observe the small reflections which may appear during the numerical solution and to understand the accuracy improvement. We clearly see on Figures 3(a)-3(b) that the CAP approach gives good results until a wave comes back into the computational domain. On figures 3(c)-3(f), the PML and ABC methods yields some good results with about the same accuracy. Low amplitude waves are reflected back into the computational domain (about 0.1% of the maximal amplitude of the soliton). For the PML, the size of the computational domain is 40% more than for the ABCs. Finally, a very good accuracy is obtained for the relaxation scheme which makes it very attractive.

Numerical simulation of a saturable nonlinear Schrödinger equation

The second test that we present concerns the numerical solution of a saturable Schrödinger equation (SNLSE) used in optics. This equation models pulse propagation in optical fibers made from doped silica [23]. The saturating nonlinear term is exponential and the equation that we solve is

   i∂ t Ψ(x, t) + ∂ 2 x Ψ(x, t) + |Ψ| 2 1 + γ|Ψ| 2 Ψ(x, t) = 0, (x, t) ∈ Ω T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω. ( 23
)
The saturation term γ is related to the saturation intensity of the fiber. We take γ = 0.5 to be close to the values in [23]. The initial datum is the soliton ( 22) (for t = 0). All the simulations parameters are the same as in Section 2.4.1. We can again see that the CAP method leads to wrong results while the PML and ABC approaches provide a good accuracy for a small reflection. The relaxation scheme is efficient and accurate.

Two-dimensional nonlinear problems

Let us consider the two-dimensional nonlinear Schrödinger equation

i∂ t Ψ(x, t) + ∆Ψ(x, t) + V (x, t) Ψ(x, t) + f (|Ψ|)Ψ(x, t) = 0, (x, t) ∈ Ω T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω. ( 24
)
Here, the computational domain Ω is a bounded set of R 2 with a regular convex boundary Σ. We compute the solution on a time interval [0; T ]. The resulting time-space domains are therefore: Ω T := Ω × [0; T ], and Σ T := Σ × [0; T ]. The initial data Ψ 0 is compactly supported in Ω. The operator ∆ is the Laplace operator in two-dimensions:

∆ := ∂ 2 x 1 + ∂ 2 x 2 , with x = (x 1 , x 2 ). The potential V is a function from R 2 × [0; T ] onto R which is smooth outside Ω.
The nonlinear function f is supposed to be smooth outside Ω.

Absorbing Boundary Conditions

The development of ABCs for two-dimensional nonlinear Schrödinger equation is a very recent area of research. Very few results are only available. We refer to [11,12] for examples. In the present paper, we consider the new ABCs developed in [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][START_REF] Antoine | [END_REF]. Again, we do not explain the technical theory for constructing such boundary conditions 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][START_REF] Antoine | [END_REF] for more details. The ABC for the two-dimensional case is

∂ n Ψ -i i∂ t + ∆ Σ + V + f (|Ψ|)Ψ + κ 2 Ψ - κ 2 (i∂ t + ∆ Σ + V + f (|Ψ|)) -1 ∆ Σ Ψ = 0, on Σ T . ( 25
)
This boundary condition can be considered as an extension of ( 6). Vector n is the outwardly directed unit normal vector to Ω. The function κ is the curvature of Σ at a point of the surface. Finally, ∆ Σ is the Laplace-Beltrami operator over the surface. For example, for a disk of radius R, the normal derivative is

∂ n = ∂ r , the curvature is R -1 and ∆ Σ = R -2 ∂ 2 θ
, where the polar coordinate system is (r, θ). We will show during the numerical approximation by the relaxation scheme how to suitably localize this operator.

Complex Absorbing Potential and Perfectly Matched Layers

The extension of the CAP method is direct in the form

i∂ t Ψ(x, t) + ∆Ψ(x, t) + V (x, t) Ψ(x, t) -iW (x) Ψ(x, t) + f (|Ψ|)Ψ(x, t) = 0, (x, t) ∈ Ω ext T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω ext . ( 26
)
The choice of the function W is however not clear most particularly for a nonlinear problem and a general domain Ω. Here, the domain Ω is the rectangle ] -a 1 , a 1 [×] -a 2 , a 2 [, and the extended computational domain

Ω ext is ] -(a 1 + δ), a 1 + δ[×] -(a 2 + δ), a 2 + δ[. The thickness of the layer is parameterized by δ. We use W (x) = W 1 (x)W 2 (x)
, where W j is given by

W j (x) = W 0 δ -2 (|x j | -a j ) 2 1 |x j |≥a j (x). ( 27
)
Function 1 O is the characteristic function of a set O. It is not clear if this choice is optimal. We do not develop the ECS here.

For the PMLs approach, the choice of domain Ω ext is restricted in practice since essentially the PMLs are written according to a special coordinates system related to Ω. For example, for the cartesian coordinates associated with the rectangles Ω and Ω ext , we consider the modified system

       i∂ t Ψ + 1 S 1 S 2 (∂ x 1 ( S 1 S 2 ∂ x 1 Ψ) + ∂ x 2 ( S 2 S 1 ∂ x 2 Ψ)) + V Ψ + f (|Ψ|)Ψ = 0, (x, t) ∈ Ω ext T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω, Ψ(x, t) = 0, x ∈ Σ ext T , (28)
with, for j = 1, 2, S j (x) = 1 + Rσ j (x)

σ j (x) = σ 0 ( |x j | -a j δ ) 2 1 |x j |≥a j (x).
Let us now consider an annulus. The physical computational domain is the disk of radius Ω = D R and the PML medium is the annulus D δ of thickness δ = R * -R. Then, going to the polar coordinates system (r, θ), the PML writes down

         i∂ t Ψ + 1 r ŜS ∂ r Ŝr S ∂ r Ψ + S Ŝr ∂ 2 θ Ψ + V Ψ + f (|Ψ|)Ψ = 0, (x, t) ∈ Ω ext T , Ψ(x, 0) = Ψ 0 (x), x ∈ Ω, Ψ(x, t) = 0, x ∈ Σ ext T , (29)
with

Ω ext = D R * , Σ ext = C R * and            S(r) = 1 + Rσ(r), Ŝ(r) = 1 + R r r R σ(s)ds, σ(r) = σ 0 r -R δ 2 1 r≥R (x).
(30)

Discretization schemes

For the Equation ( 24), the interior relaxation scheme which is based on

i∂ t Ψ + ∆Ψ + V Ψ + ΥΨ = 0, on Ω T , Υ = f (|Ψ|), on Ω T , (31)
is given by

       2i Φ n+1 ∆t + ∆Φ n+1 + W n+1 Φ n+1 + Υ n+1/2 Φ n+1 = 2i Ψ n ∆t , Υ n+3/2 + Υ n+1/2 2 = f (|Ψ n+1 |), (32)
for 0 ≤ n ≤ N -1, where the notations are the same as (12).

In the case of the ABC (25), we discretize the relation by using Υ = f (|Ψ|) and some additional auxiliary functions (ϕ k ) 1≤k≤m and ψ. More precisely, we get the scheme

                                 ∂ 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 + ∆ Σ + W n+1 + Υ n+1/2 + d m k ϕ n+1/2 k -Φ n+1 = 2i ∆t ϕ n k , on Σ, 2i ∆t + ∆ Σ + W n+1 + Υ n+1/2 ψ n+1/2 -∆ Σ Φ n+1 = 2i ∆t ψ n , on Σ, ϕ 0 k = 0 on 1 ≤ k ≤ m, ψ 0 = 0, on Σ, Υ n+3/2 = 2f (|Φ n+1 |) -Υ n+1/2 , on Ω. ( 33
)
In view of a finite element approach, we need to write the associated weak formulation for the first equation of system (32), for a test-function ϕ,

2i ∆t Ω Φ n+1 ϕdΩ - Ω ∇Φ n+1 • ∇ϕ dΩ + Σ ∂ n Φ n+1 ϕ dΣ + Ω W n+1 Φ n+1 ϕ dΩ + Ω Υ n+1/2 Φ n+1 ϕ dΩ = 2i ∆t Ω Ψ n ϕ dΩ, (34)
and updating by

Υ n+3/2 = 2f (|Ψ n+1 |) -Υ n+1/2 . ( 35
)
Then, we inject the expression of the normal derivative ∂ n Φ n+1 by using (33). This leads to a linear system of equations with unknowns Φ n+1 in Ω, coupled to the surface equations in terms of additional unknowns (ϕ k ) 1≤k≤m and ψ. For the CAP and PML, both the time and space discretizations are directly realized in the extended computational domain. Finite difference are used for the spatial discretization.

3.4. The example of the cubic media for the propagation of a soliton 3.4.1. Soliton construction by a shooting method In the two-dimensional case, there is no explicit analytical expression of the soliton. Its construction must be realized numerically. To this aim, we compute a solution to the 2D stationary Schrödinger equation by using a shooting method [24]. Let us consider the nonlinear focusing cubic Schrödinger equation

i∂ t ψ + ∆ψ + |ψ| 2 ψ = 0, on R 2 × R + , (36)
and let us compute a stationary solution under the form

ψ(r, t) = e iµt φ(r), (37)
where r = x = x 2 + y 2 , µ ∈ R and φ is supposed to be spatially localized. Then we have to solve the nonlinear elliptic equation

-µφ + ∆φ + |φ| 2 φ = 0, x ∈ R 2 . ( 38
)
We make the assumption that φ has a radial symmetry. Then going to the radial coordinates, we have to solve a second-order ordinary differential equation on the interval [0; R]

∂ 2 r φ + 1 r ∂ r φ -µφ + |φ| 2 φ = 0, r ∈ [0; R]. (39)
We can next work for µ = 1 by the change of variable

φ(r) = 1 √ µ φ r √ µ .
Furthermore, to get a two differentiable smooth solution we impose ψ (0) = 0 to avoid any singularity at 0. Finally, we have to solve the differential nonlinear system

   ∂ 2 r φ + 1 r ∂ r φ -φ + |φ| 2 φ = 0, 0 < r < R, φ (0) = 0, φ(0) = β, (40)
where we try to find a solution which tends towards zero to infinity. A Taylor expansion at zero shows that we must have φ (0) = β -q|β| 2 β. A shooting method is then applied by taking a radial step of discretization ∆r = 10 -3 . The initial data β is ajusted in such a way that suitable decay of the solution φ is obtained. We therefore have a radial solution φ(r), for r ≤ R, which is then extended to the whole disk of radius R by radial symmetry to get the soliton solution. Then, the soliton is multiplied by a gaussian with wavenumber k 0 = k 0 x as ψ 0 (x, y) = φ(r)e -ik 0 x 1 (41)

to make it moves. At R = 10, we have the approximation |ψ 0 (R)| ≈ 5 × 10 -5 . If we extend the domain to R = 15, we get |ψ 0 (R)| ≈ 3 × 10 -7 . Indeed, the soliton has a slow decay which implies that relatively large computational domains must be chosen to be sure that the initial data is numerically compactly supported. At the computer code level, this implies that we have to work in quadruple precision for R = 15. For stability reasons, the number of significative digits is crucial for β. Finally, we will consider next the disk of radius R = 10 as computational domain Ω.

Remark 1. This construction of the soliton directly extends to other nonlinearities such as f (|ψ|) = |ψ| 2σ , with σ > 0. This allows for example to consider the quintic nonlinearity f (|ψ|) = |ψ| 4 . Let us remark that other possible applications of interest of the shooting techniques (for stationary solutions) as well as relaxation schemes (for nonlinear dynamics) could be for instance the numerical solution of nonlinear Schrödinger equations modelling DC and AC Josephson effects for superfluid Fermi in BCS-BEC crossover at zero temperature [25,26]. Then, specific potentials must be considered for realistic simulations. In the same spirit, applications could be considered in the background of Feshback resonance where specific nonlinear potentials with trapping potential and external driving field are included [27].

Accuracy of the truncation techniques

We now consider as initial data the soliton solution computed by the previous shooting method and then modulated by e ik 0 x 1 with k 0 = 5. For the CAP and PMLs, we take ∆t = 10 -3 for a rectangular physical domain Ω := [-10, 10] × [-10, 10] embedded in the extended domain Ω ext := [-12, 12]× [-12, 12]. This last domain is discretized with a uniform grid composed of 501 × 501 points for the finite difference approximation. For the ABC approach, the domain is the circle of radius 10 for n P = 220000 degrees of freedom of the linear finite element method. Figure 5(a) represents the 3D propagation of the soliton solution |ψ| in the space (x 1 , x 2 , t) with ABCs in log-scale. We can observe some small reflections at the boundary when the soliton hits the left plane. For an easier visualization of the results, we choose to report a slice of the wave field in the plane (x 1 , t) for x 2 = 0 again in log-scale. This gives Figures 5(b)-5(d) for respectively the CAPs, PMLs and ABCs approaches. As in the one-dimensional case, the proposed CAP method gives incorrect results while both PMLs and ABCs solution are physically correct with small reflection back. Furthermore, the relaxation schemes yield again a suitable accuracy for a low computational cost.

Conclusion

We proposed in this paper a numerical comparison of CAPs, PMLs and ABCs techniques for one-and two-dimensional nonlinear Schrödinger equations (cubic and saturation media). The initial data is a soliton which is analytically known for the one-dimensional case and numerically built by a shooting method in the two-dimensional case. The numerical schemes are based on a relaxation scheme in time and finite element or finite difference in space. From the numerical simulations, it appears that the PMLs and ABCs provide a suitable and similar accuracy compared to the CAPs method which leads to wrong results. Finally, all the computer codes are available freely at http://microwave.math.cnrs.fr/code/index.html.

Figure 1 :

 1 Figure 1: Computational domain Ω T and fictitious boundary Σ T = Σ×]0; T [.

Figure 2 :

 2 Figure 2: Computational domain for the CAP and PML approaches.

Figure 3 :

 3 Figure 3: Numerical solutions |ψ| for the cubic case (left: normal scale, right: log scale).

Figure 4 :

 4 Figure 4: Numerical solutions |ψ| for the saturation case (log scale).

Figure 5 :

 5 Figure 5: Numerical solutions |ψ| for the 2D cubic case: Figure 5(a): 3D propagation of the soliton solution on the disk of radius R = 10 with ABCs. Figures 5(b)-5(d): CAPs, PMLs and ABCs solution. The figures are built as slices in the (x 1 , t) plane for x 2 = {0}.

partially by the French ANR fundings under the project MicroWave NT09 460489 (http://microwave.math.cnrs.fr/).