Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Christophe Geuzaine
email: cgeuzaine@uliege.be

Qinglin Tang
email: qinglin_tang@scu.edu.cn

Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates

Keywords: nonlinear Schrödinger equation, rotating Gross-Pitaevskii equation, Bose-Einstein condensate, perfectly matched layers (PML), Fourier pseudospectral method, time-splitting scheme

, some numerical simulations illustrating the strength of the proposed approach.

Introduction

Since its first experimental creation in 1995 [START_REF] Anderson | Observation of Bose-Einstein condensation in a dilute atomic vapor[END_REF][START_REF] Bradley | Evidence of Bose-Einstein condensation in an atomic gas with attractive interaction[END_REF][START_REF] Davis | Bose-Einstein condensation in a gas of sodium atoms[END_REF], the Bose-Einstein condensation (BEC) phenomena provides an incredible glimpse into the macroscopic quantum world and has opened a new era in atomic and molecular physics as well as in condensed matter physics. It has been extensively studied both experimentally and theoretically, and is still a very active research topic [START_REF] Abo-Shaeer | Observation of vortex lattices in Bose-Einstein condensates[END_REF][START_REF] Andersen | Theory of the weakly interacting Bose gas[END_REF][START_REF] Baranov | Theoretical progress in many body physics of dipolar gases[END_REF][START_REF] Bloch | Many body physics with ultracold gases[END_REF][START_REF] Dalfovo | Theory of Bose-Einstein condensation in trapped gases[END_REF][START_REF] Fetter | Rotating trapped Bose-Einstein condensates[END_REF][START_REF] Giorgini | Theory of ultracold atomic Fermi gases[END_REF][START_REF] Leggett | Bose-Einstein condensation in the alkali gases: Some fundamental concepts[END_REF][START_REF] Lieb | Derivation of the Gross-Pitaevskii equation for rotating Bose gases[END_REF][START_REF] Morsch | Dynamics of Bose-Einstein condensates in optical lattices[END_REF][START_REF] Pitaevskii | Bose-Einstein Condensation[END_REF][START_REF] Proukakis | Universal themes of Bose-Einstein condensation[END_REF]. At temperatures T much smaller than the critical temperature T c , the properties of a rotating BEC are well described by a macroscopic complex-valued wave function ψ(x, t) whose evolution is governed by the threedimensional (3D) Gross-Pitaevskii equation (GPE). Solving the d-dimensional (d = 2 or 3) dimensionless GPE with a rotation term [START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] leads to the following boundary-value problem: for a given initial state ψ 0 , find the complex-valued wave function ψ(x, t) solution to

i∂ t ψ(x, t) = - 1 2 ∇ 2 + V (x, t) -ωL z + β|ψ| 2 ψ(x, t), ψ(x, t = 0) = ψ 0 (x), x ∈ R d , t ≥ 0, (1.1)
where x := (x, y, z) (:= (x, y) in 2D) and t are the space and time variables, respectively. Denoting by ∇ the gradient operator, ∇ 2 is then the laplacian operator, and V (x, t) is a real-valued function corresponding to the potential. The real-valued constants β and ω respectively represent the nonlinear interaction strength and the rotating frequency. In addition, L z = i(y∂ x -x∂ y ) is the z-component of the angular momentum [START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF]. When ω = 0, the GPE is also often called the nonlinear cubic Schrödinger equation. Let us remark that many more complex models involving GPEs are widely used in the literature [START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF], including in particular nonlocal dipolar interactions or multi-components gases.

When computing the dynamics of the rotating GPE, one can use various stable and efficient schemes that are known to also respect some dynamical properties (mass/energy conservation, time invariance, dispersion relation, time reversibility,...). Here, we do not want to review all the possible methods and we rather focus on two widely used time discretization schemes, i.e. the (Strang) time-splitting [START_REF] Bao | Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation[END_REF] and relaxation [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] methods which are second-order accurate in time. We refer to [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] for more details about the most popular schemes and their properties. When computing fast rotating BEC through the GPE, quantized vorticity takes place and then highly accurate spatial discretization schemes are required to capture the dynamics of the small quantum defects. Among the available numerical methods, the most standard way [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF] is to use FFT-based pseudospectral approximation schemes in cartesian coordinates. Indeed, such approximations are easy to implement, efficient and highly accurate [START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: dynamics and stochastic simulations[END_REF][START_REF] Antoine | Modeling and Computation of Bose-Einstein Condensates: Stationary States, Nucleation, Dynamics, Stochasticity, Nonlinear Optical and Atomic Systems: at the Interface of Physics and Mathematics[END_REF]. In particular, they can be trivially included into time-splitting schemes by a direct integration in the Fourier domain. The implementation work for the relaxation scheme is a bit more involved, but it remains quite simple as it relies on well-known computational bricks (i.e. FFTs, Krylov subspace solvers and simple matrix-free preconditioners). We also consider in this paper the (linear) finite element method for comparison, which is nevertheless only second-order in space (but higher order polynomials could also be used).

In some situations arising in BECs (as well as in many situations related to the nonlinear Schrödinger equation, e.g. in quantum optics or nonlinear wave propagation), the condensate expands. This can be the case for example when a BEC is created and then the confining potential is no longer active. In other situations related to the nonlinear Schrödinger equation, considering that the wave function propagates out of the a priori fixed computational domain is very standard. In these cases, one needs to use a very large computational box, leading then to extensive computational costs and memory storage limitations. An alternative approach consists in considering a fixed box and then setting suitable boundary conditions at the domain interface. According to their mathematical properties, such boundary conditions are called transparent, artificial or absorbing boundary conditions (see [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for some recent reviews in quantum and relativistic mechanics). While being now widely used in practice, their implementation in a finite-difference or finite element scheme is far from being trivial. In addition, considering such boundary conditions in a FFT-based pseudospectral approximation scheme is not possible since periodic boundary conditions are already imposed at the domain boundary. An alternative is to use Perfectly Matched Layers (PMLs) as introduced in computational electromagnetism by Bérenger [19,[START_REF] Bérenger | Three-dimensional perfectly matched layer for the absorption of electromagnetic waves[END_REF][START_REF] Bérenger | Perfectly matched layer for the FDTD solution of wave-structure interaction problems[END_REF][START_REF] Petropoulos | Reflectionless sponge layers as absorbing boundary conditions for the numerical solution of Maxwell equations in rectangular, cylindrical and spherical coordinates[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF]. This method is now widely used in many areas of computational physics and engineering [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with Perfectly Matched Layers[END_REF][START_REF] Bécache | Perfectly mathched layers for the convected Helmholtz equation[END_REF][START_REF] Bermúdez | An exact bounded perfectly matched layer for timeharmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF][START_REF] Collino | The perfectly matched layer in curvilinear coordinates[END_REF][START_REF] Collino | Application of the perfectly matched absorbing layer model to the linear elastodynamic problem in anisotropic heterogeneous media[END_REF][START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF][START_REF] Hu | On absorbing boundary conditions for linearized Euler equations by a perfectly matched layer[END_REF][START_REF] Hu | A stable, perfectly matched layer for linearized Euler equations in unsplit physical variables[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. A review[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF][START_REF] Zeng | The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[END_REF] since it is easy to implement and accurate. Concerning the application of the PML method to linear and nonlinear Schrödinger-type equations, we refer to [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Collino | Perfectly matched absorbing layers for the paraxial equations[END_REF][START_REF]Perfectly matched layers for coupled nonlinear Schrödinger equations with mixed derivatives[END_REF][START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[END_REF][START_REF] Nissen | An optimized perfectly matched layer for the Schrödinger Equation[END_REF][START_REF] Scrinzi | On the non-equivalence of perfect matched layers and exterior complex scaling[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] where various developments are available in the framework of the finite-difference/finite element methods. (Let us remark that the PML method is different from the exterior complex scaling methods [START_REF] Mccurdy | Solving the three-body Coulomb breakup problem using exterior complex scaling[END_REF][START_REF] Muga | Complex absorbing potentials[END_REF][START_REF] Scrinzi | On the non-equivalence of perfect matched layers and exterior complex scaling[END_REF][START_REF] Weinmueller | Perfect absorption in Schrödinger-like problems using non-equidistant complex grids[END_REF].) To the best of the authors' knowledge, there is no contribution until now related to the application of PMLs in the context of FFT-based methods, for both the time-splitting and relaxation schemes applied to Schrödinger equations and GPEs. The aim of this paper is to address this question, to show that this can be done without great effort and that the resulting scheme is indeed numerically efficient, accurate and that the reflection at the boundary is relatively small.

The paper is organized as follows. In Section 2, we introduce the PML formulation in cartesian coordinates for the GPE without rotating term, which means that ω = 0 in system (1.1). We introduce some notations and give six types of PML absorbing functions that will be tested later. In Section 3, we consider first the time discretization of the PML formulation of the GPE based on the Strang splitting scheme as well as relaxation scheme. Then, we detail how to discretize the formulation in space, first by a finite element method but most importantly next by an FFT-based pseudospectral approximation scheme. In Section 4, we develop a thorough numerical study of the approximation schemes and in particular we focus on the choice of the PML absorbing function and its tuning parameters in the case of the highly accurate pseudospectral schemes. The direct extension to the rotating GPE is non trivial as it is explained in Section 5, since a gradient term is involved into the GPE arising in system (1.1). We then propose a relatively simple solution by considering an equivalent formulation in rotating Lagrangian coordinates which allows to extend the previous PML framework with pseudospectral approximation to the rotating case. Some numerical simulations illustrate the efficiency and accuracy of the method. Finally, we conclude in Section 6.

PMLs formulation in cartesian coordinates for the GPE without rotation (ω = 0)

Let us assume that the initial profile ψ 0 is compactly supported into a bounded domain of physical interest D Phy with boundary Γ Phy := ∂D Phy . To solve the GPE (1.1), we consider now that the equation is set in the truncated domain D Phy

i∂ t ψ(x, t) = - 1 2 ∇ 2 + V (x, t) -ωL z + β|ψ| 2 ψ(x, t), x ∈ D Phy . (2.1)
Therefore, to get a complete initial boundary value problem, additional boundary conditions (BCs) are required at the domain boundary Γ Phy := ∂D Phy . Most particularly, when the potential V (x, t) is not trapping, then outgoing waves may emerge and thus suitable BCs are needed to be reflection-less or nonreflecting [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF]. Our analysis starts in this Section with the non-rotating GPE (ω = 0), the extension to the rotating case will be next studied in Section 5.

Let us first introduce the concept of PML [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] by considering the toy model for the 1D linear Schrödinger equation, i.e.

i∂ t ψ(ν, t) = - 1 2 ∂ 2 x ψ(ν, t), ν ∈ R. (2.2)
We assume that D Phy := [-L * ν , L * ν ] and consider the right traveling outgoing waves only. The PML technique can be carried out by stretching the real coordinate into the complex plane [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF]:

ν = ν + e iϑν ν L * ν σ(s)ds, σ(s) ≡ 0 if s < L * ν . (2.3)
Here, ϑ ν ∈ (0, π 2 ) is a real constant and σ(s) is a real-valued function called absorbing function. Plugging (2.3) into (2.2), we obtain the corresponding PML equation

i∂ t ψ(ν, t) + 1 2S(ν) ∂ ν 1 S(ν) ∂ ν ψ(ν, t) = 0, ν ∈ R, (2.4)
where S(ν) = 1 + e iϑν σ(ν). Constant ϑ ν and function σ(ν) are chosen such that ψ is damped and decays fast enough in ν ∈ [L * ν , +∞). Hence, it suffices to truncate the semi-infinite interval [L * ν , +∞) as a finite right PML region [L * ν , L ν ] and to impose a homogeneous Dirichlet or periodic BC on ψ at ν = L ν . A similar treatment can be applied to the left semi-infinite interval (-∞, -L * ν ]. Denote the computational domain as

D := D Phy ∪D PML := [-L ν , L ν ], with D PML := [-L ν , -L * ν ]∪[L * ν , L ν ],
and set Γ := ∂D := {-L ν , L ν }. Then, the Schrödinger equation (2.2) can be approximated by the following PML equation within the domain D:

i∂ t ψ(ν, t) + 1 2S(ν) ∂ ν 1 S(ν) ∂ ν ψ(ν, t) = 0, ν ∈ D, (2.5)
ψ(ν, t) = 0, ν = ±L ν , or ψ(-L ν , t) = ψ(L ν , t), t > 0, (2.6)
with

S(ν) = 1, |ν| < L * ν , 1 + e iϑν σ(|ν| -L ν ), L * ν ≤ |ν| < L ν .
(2.7)

Here, ϑ ν is a given constant which is fixed as ϑ ν = π 4 [START_REF] Antoine | Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] hereafter, and σ(ν) is a given absorbing function. In this paper, we consider the following six types of absorbing functions which are well studied and successfully applied e.g. to scattering problems [START_REF] Bermúdez | An exact bounded perfectly matched layer for timeharmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF] Type 1 :

σ 0 (ν + δ ν ) 2 , Type 3 : -σ 0 /ν, Type 5 : -σ 0 /ν -σ 0 /δ ν , Type 2 : σ 0 (ν + δ ν ) 3 , Type 4 : σ 0 /ν 2 , Type 6 : σ 0 /ν 2 -σ 0 /δ 2 ν , (2.8)
with δ ν = L ν -L * ν being the thickness of the PML layer. Notice that S(ν) is discontinuous at the interface ν = L * ν if the Type 3 or Type 4 absorbing functions are chosen. We are going to compare these six types of absorbing functions and find the (range of) optimal absorbing strength σ opt 0 for each type of profile with different thickness δ ν and discretization schemes for the full GPE.

Let us derive the PML equation of the GPE (2.1) with ω = 0 in Cartesian coordinates. Here, we only show the 2D case, the extension to the 3D case is straightforward and is omitted for brevity. Without loss of generality, we assume that the computational domain of interest is a rectangle, i.e.,

D Phy = [-L * x , L * x ] × [-L * y , L * y ].
Analogous to the derivation of (2.5), placing a PML region D PML =: [-L x , L x ] × [-L y , L y ] \ D Phy surrounding D Phy , stretching the coordinates in both the x-and y-directions, i.e., taking the following substitutions

∂ x -→ ∂ x S x (x) , ∂ y -→ ∂ y S y (y) , (2.9)
we obtain the PML equation for the non-rotating GPE (2.1) in the computational domain

D := [-L x , L x ] × [-L y , L y ] = D Phy ∪ D PML : i∂ t ψ = - 1 2 1 S x (x) ∂ x 1 S x (x) ∂ x + 1 S y (y) ∂ y 1 S y (y) ∂ y ψ + V + β|ψ| 2 ψ, x ∈ D, t ≥ 0. (2.10)
The homogeneous Dirichlet or periodic BCs can then be applied at the boundary Γ := ∂D. To simplify the presentation, hereafter, we let S ν (ν) = S(ν) (ν = x, y), with S(ν) reading as (2.7) and denote the linear differential operator in (2.10) as

L = - 1 2 1 S x (x) ∂ x 1 S x (x) ∂ x + 1 S y (y) ∂ y 1 S y (y) ∂ y .
(2.11)

Discretization (ω = 0)

In this section, we propose different approaches to discretize equation (2.10). To this end, we choose ∆t > 0 as the time step size and denote the time steps by t n = n∆t, for n ≥ 0.

Time discretization

A first well-known approach to deal with the nonlinearity in the GPE (2.10) is to apply the time-splitting technique [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][START_REF] Bao | Numerical solution of the Gross-Pitaevskii equation for Bose-Einstein condensation[END_REF], i.e., one solves

i ∂ t ψ(x, t) = Lψ(x, t), x ∈ D, t n ≤ t ≤ t n+1 , (3.1)
for the time step of length ∆t, followed by solving

i ∂ t ψ(x, t) = V (x, t) + β |ψ(x, t)| 2 ψ(x, t), x ∈ D, t n ≤ t ≤ t n+1 , (3.2)
for the same time step. Equation (3.1) will be further discretized in the next subsection. For t

∈ [t n , t n+1 ], Eq. (3.2) leaves the density |ψ(x, t)| 2 invariant, i.e., |ψ(x, t)| 2 = |ψ(x, t n )| 2 := |ψ n (x)| 2 . Therefore, (3.2)
collapses to be a linear ODE, which can be integrated analytically as

ψ(x, t) = e -i[V(x,tn,t)+β |ψ n (x)| 2 (t-tn)] ψ(x, t n ), x ∈ D, t ∈ [t n , t n+1 ].
Here, for any

a, b ∈ R, function V(x, a, b) reads V(x, a, b) =: b a V (x, τ )dτ.
Another popular approach is to use the relaxation technique introduced by Besse in [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] iδ

+ t ψ n = Lψ n+ 1 2 + (V n+ 1 2 + ϕ n+ 1 2 )ψ n+ 1 2 , in D, (3.3)
with

ϕ n+ 1 2 = 2β|ψ n | 2 -ϕ n-1 2 , for n ≥ 0.
Here, we set:

ϕ -1 2 = β|ψ 0 (x)| 2 , ψ n (x) = ψ(x, t n ), δ + t ψ n = (ψ n+1 -ψ n )/∆t and f n+ 1 2 = (f n+1 + f n )/2 (f = V and ψ).
The relaxation scheme owns the important property that the resulting equation (3.3) is linear and therefore does not need any nonlinear solver unlike the standard Crank-Nicolson scheme, which greatly reduces the computational time and memory cost.

Space discretization

We now introduce two approaches to further discretize the equations (3.1)-(3.3) in space. The first one is based on the Finite Element Method (FEM) [START_REF] Dular | A general environment for the treatment of discrete problems and its application to the finite element method[END_REF] while the other uses a Fourier Pseudospectral (FP) scheme [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF].

Let us start with the FEM. To this end, we simply impose homogeneous Dirichlet BCs for (2.10) at Γ, and hence for Eqs. (3.1) & (3.3). The weak formulation for (3.1) reads as follows: find ψ(x, •) ∈ H 1 0 (D) such that for ∀φ ∈ H 1 0 (D), we have

i d dt D ψ(x, t)φS x S y dx = 1 2 D Q∇ψ(x, t) • ∇φ dx. (3.4)
On the other hand, the weak formulation for (3.3) reads: find

ψ n+1 ∈ H 1 0 (D) such that for ∀φ ∈ H 1 0 (D) i D δ + t ψ n φS x S y dx = 1 2 D Q∇ψ n+ 1 2 • ∇φ dx + V n+ 1 2 + ϕ n+ 1 2 ψ n+ 1 2 φ S x S y dx. (3.5)
In the above two equations, the 2 × 2 diagonal tensor Q is defined as

Q = diag([S y (S x ) -1 , S x (S y ) -1 ]).
To further discretize (3.4) and (3.5), we define T h as a family of partitions (e.g. triangulations or quadrangulations) of the domain D with mesh size h, i.e. D h := ∪ T ∈T h T . The notation T designates an element of the partition T h and Γ h is the partition of the boundary Γ. Moreover, let P k (k = 1 latter in the simulations) be the space of polynomials with degree less or equal to k ≥ 1 and define the finite element spaces

V h = ψ h ∈ C 0 (D) | ψ h | T ∈ P k , ∀ T ∈ T h , V h 0 = ψ h ∈ V h | ψ h | Γ h = 0 . (3.6)
The FEM approximation of (3.4) for the splitting scheme reads:

iA h d dt ψ h (t) = B h ψ h (t), (3.7)
where ψ h (t) is the nodal FEM vector and A h and B h are respectively the generalized mass and stiffness matrices corresponding to the following bilinear forms:

∀(u h , v h ) ∈ V h 0 × V h 0 (A h u h , v h ) = D h S x S y u h v h dx, (B h u h , v h ) = 1 2 D h Q∇u h • ∇v h dx. (3.8)
The resulting system (3.7) is a linear ODE system, which can be further solved with a well-adapted highorder time discrete scheme. In the present paper, we restrict ourselves to the standard Crank-Nicolson scheme. Combined with (3.2), a fully discrete second-order Strang time-splitting FEM (TSFEM) scheme reads:

ψ (1) = ψ n h * e -i V(tn,tn+∆t/2)+β∆t |ψ n h | 2 /2 , (3.9)
i ∆t A h - 1 2 B h ψ (2) = i ∆t A h + 1 2 B h ψ (1) , (3.10)
ψ n+1 h = ψ (2) * e -i V(tn+∆t/2,tn+1)+β∆t |ψ (2) | 2 /2 . (3.11)
Here, V(s, q) is the vector containing the elements V(x, s, q) for all x at the grid nodes and * means the component-by-component product of two complex-valued vectors. Let us remark that the second step of the TSFEM algorithm can be mostly pre-computed once before entering into the time loop by considering a LU factorization since A h and B h are (symmetrical) matrices that do not depend on n. Therefore, each solution of the second time step only requires the computation of the forward/backward substitution which has a linear cost with respect to the finite element space dimension. Finally, the extension to higher-order TS schemes is direct.

As for (3.5), the FEM approximation leads to the following relaxation FEM (ReFEM) scheme

i ∆t A h - 1 2 B h - 1 2 C n h ψ n+1 h = i ∆t A h + 1 2 B h + 1 2 C n h ψ n h , (3.12)
where C n h is the matrix related to the finite element approximation of the bilinear form:

∀(u h , v h ) ∈ V h 0 × V h 0 (C n h u h , v h ) = D h V n+ 1 2 + ϕ n+ 1 2 S x S y u h v h dx. (3.13)
Let us remark that C n h must be recomputed for each n, and the corresponding linear system must be rebuilt and solved, which is finally computationally highly expensive.

The second spatial discretization scheme is based on the Fourier pseudospectral method. To this end, periodic BCs on the boundary Γ must be imposed for equation (2.10), and hence to both (3.1) and (3.3). Let L and M be two even integer numbers. We define h x = 2Lx L and h y = 2Ly M as the mesh sizes in the xand y-directions, respectively. To simplify the notations, we introduce the indices and grid points sets as well as the basis functions as

T LM = ( , m) ∈ N 2 | 0 ≤ ≤ L, 0 ≤ m ≤ M , T LM = (p, q) ∈ Z 2 | -L/2 ≤ p ≤ L/2 -1, -M/2 ≤ q ≤ M/2 -1 , G xy = {(x , y m ) =: (-L x + h x , -L y + m h y ), ( , m) ∈ T LM } , W pq (x) = e iµ x p (
x+Lx) e iµ y q (y+Ly) , µ x p = πp/L x , µ y q = πq/L y , (p, q) ∈ T LM .

Let f n m (f = ψ, ϕ, V , S x , S y ) be the approximation of f (x , y m , t n ) for ( , m) ∈ T LM and n ≥ 0. We denote by f n (f = ψ, ϕ, V, S x , S y ,) as the vector with components {f n m , ( , m) ∈ T LM }. The Fourier pseudospectral discretization of a function ψ is given by

ψ(x, y, t) = L/2-1 p=-L/2 M/2-1 q=-M/2 ψ pq (t)W pq (x, y), (3.14)
where the Fourier coefficients ψ pq (t) read as

ψ pq (t) = 1 LM L-1 j=0 M -1 k=0 ψ m (t)e -iµ x p (x +Lx) e -iµ y q (ym+Ly) . (3.15)
Hence, the Fourier pseudospectral discretizations of ∂ x ψ and ∂ y ψ are respectively given by: for ( , m) ∈ T LM ,

∂ x ψ m (t) = L/2-1 p=-L/2 M/2-1 q=-M/2 iµ x p ψ pq (t)W pq (x , y m ), (3.16
)

∂ y ψ m (t) = L/2-1 p=-L/2 M/2-1 q=-M/2
iµ y q ψ pq (t)W pq (x , y m ).

(3.17)

Furthermore, we define the operators L , S -1

x , S -1 y , V n+ 1 2 and ϕ n+ 1 2 as follows

L = - 1 2 S -1 x ∂ x S -1 x ∂ x + S -1 y ∂ y S -1 y ∂ y , (3.18)
S -1 x ψ m = ψ m /S x m , V n+ 1 2 ψ m = 1 2 V n m + V n+1 m ψ m , (3.19)
S -1 y Ψ m = ψ m /S y m , ϕ n+ 1 2 ψ m = 1 2 ϕ n+1 m + ϕ n m ψ m . (3.20)
The first scheme brings together the time-splitting and spectral methods, e.g. a second-order Time-Splitting Fourier Pseudospectral method (TSFP) as follows: for ( , m) ∈ T LM

ψ (1) = ψ n * e -i V(tn,tn+∆t/2)+β∆t |ψ n LM | 2 /2 , (3.21)
i I ∆t - 1 2 L ψ (2) = i I ∆t + 1 2 L ψ (1) , (3.22)
ψ n+1 = ψ (2) * e -i V(tn+∆t/2,tn+1)+β∆t |ψ (2) | 2 /2 . (3.23)
Here V m (t 1 , t 2 ) = V(x , y m , t 1 , t 2 ), (ψ n LM ) pq and (ψ (2) ) pq are the discrete Fourier series coefficients of the vectors ψ n LM and ψ (2) , respectively. The second scheme brings together the relaxation technique and spectral method, i.e. the relaxation Fourier Pseudo-spectral methods (ReFP) reads:

i I ∆t - 1 2 L - 1 2 V n+ 1 2 - 1 2 ϕ n+ 1 2 ψ n+1 = i I ∆t + 1 2 L + 1 2 V n+ 1 2 + 1 2 ϕ n+ 1 2 ψ n . (3.24)
For solving the linear systems (3.22) and (3.24), a direct solution cannot be used since all the matrix operators are related to variable coefficients PDEs, in particular because of the presence of the PML layer. An alternative approach [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: dynamics and stochastic simulations[END_REF] consists in using a preconditioned Krylov subspace iteration solver. In the present paper, we use the GMRES [START_REF] Saad | Iterative Methods for Sparse Linear Systems[END_REF][START_REF] Saad | GMRES -a Generalized Minimal RESidual algorithm for solving nonsymmetric linear systems[END_REF], accelerated by suitable preconditioners [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: computation of stationary solutions[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations II: dynamics and stochastic simulations[END_REF], e.g.

(i I ∆t -1 2 L 0 ) -1
, where L 0 is the laplacian operator (for σ = 0) in the extended domain D. This operator can be trivially computed by FFT since it is a constant coefficients linear operator and can then serve as a matrix-free preconditioner to (3.24).

Numerical results (ω = 0)

In this section we first carry out, on a linear case, a detailed performance comparison of the TSFEM (3.9)-(3.11), the ReFEM (3.12), the TSFP (3.21)-(3.23) and the ReFP (3.24) methods proposed in the previous section. Then, we further compare the TSFP and ReFP via a nonlinear case. Finally, we apply the TSFP method to simulate the dynamics of a ground states that released from its trapping potential. To simplify the comparison, we only consider two-dimensional examples. The physical domain is fixed to be a square centered at the origin i.e., D Phys = [-L, L] 2 . Moreover, the size of the PML domain and the strength of the absorbing functions in the x-and y-directions are chosen to be identical, i.e., δ x = δ y = δ and σ x 0 = σ y 0 = σ for some provided constants δ and σ. Hence, the computational domain is set to be D = [-L -δ, L + δ] 2 . Furthermore, for the FEM-based methods, the computational domain is partitioned by uniform rectangles with mesh sizes h x and h y in the x-and y-directions, respectively. For all the four methods, we choose h x = h y = h and denote ψ n num as the vector consisting of the numerical solution of ψ(x, t n ) on given grids that lie in the physical domain D Phys . Meanwhile, we denote by ψ n ref the reference solution in D Phys . To quantized the error, we define the relative ∞ errors as:

e rel ∞ (t = t n ) = ψ n num -ψ n ref ∞ ψ n ref ∞ , (4.1)
where • ∞ represents the standard ∞ -norm on a vector space.

Linear case

Here, we consider a non-rotating linear case, i.e., we set β = ω = 0. Moreover, the potential V and initial data ψ 0 are chosen respectively as

V (x, t) = -8|x| 2 , ψ 0 (x) = e -8|x| 2 . ( 4.2)
With these parameters, the linear Schrödinger equation admits the following exact solution which spreads out as the time evolves

ψ(x, t) = i i -4 tanh(4t) exp -8|x| 2 sech 2 (4t) + 34i|x| 2 tanh(4t) 1 + 16 tanh 2 (4t) . (4.3)
Unless stated otherwise, we fix L = 2 and ∆t = 10 -4 . In addition, the relative solution ψ n ref in (4.1) is taken as the exact solution (4.3) at the grid points at time t = t n .

Example 4.1. We first compare the different kinds of absorbing functions and we determine the associated optimal absorbing strength σ for each type, and for the four schemes: ReFEM,TSFEM, ReFP and TSFP. To this end, we take h = 1 16 . Fig. 4.1 shows e rel ∞ (0.24) for the six types of absorbing functions with PML sizes δ = 0.25 and δ = 0.5. From the figures, we can see that: (i) For ReFEM and TSFEM, all the 6 types absorbing functions are almost of the same quality. They have different "optimal" regions of σ, which are simply a "shift" from each other. A similar conclusion can be made for both the ReFP and TSFP schemes. (ii) For ReFEM and TSFEM, as δ increases, the relative error is saturated for all σ, which indicates that the discretization error is reached (second-order accuracy in space of linear elements). For this case, it seems better to use a PML of type 4 or 6 as δ gets larger since the range of σ parameters appears to be larger as well. Meanwhile for ReFP and TSFP, the error still decays to its discretization error as δ increases, which is much smaller than those in FEM-based schemes. (iii) ReFEM and TSFEM behave similarly for each fixed absorbing function and parameter δ. This is similar for ReFP and TSFP with small δ. While for a larger value of δ, ReFP is more accurate than TSFP. (iv) Although these four schemes have similar region of "optimal" σ, Fourier spectral-based schemes outperform the FEM-based schemes in terms of accuracy in the optimal region and are memory cheap since they are matrix-free. Let us remark that ReFEM needs the reconstruction of the matrix of the linear system at each time step, which is quite time consuming and therefore not as efficient as the TSFEM. The preconditioned solution of the linear systems is also very efficient through the preconditioned GMRES for the FP approaches. (v) For all the four schemes, the region of optimal σ are approximately the same for the types 1 and 2, types 3 and 4 and types 5 and 6, respectively. For types 1 and 2, larger values of σ are preferred while smaller values σ are better for types 5 and 6. We recommend the use of intermediate values of σ for types 3 and 4. (vi) In terms of both accuracy and efficiency, we conclude that ReFP and TSFP outperform ReFEM and TSFEM. Example 4.2. In this example, we take the same parameters as in Example 4.1. We further consider the stability region at different times. Fig. 4.2 shows the corresponding results for types 2, 4, 6 for δ = 0.5 for various types of absorbing functions. From this figure and additional experiments not shown here for brevity, we can see that: (i) For fixed h and δ, the stability of the optimal region for σ is not fundamentally changed for different times. (ii) For other values of δ, the conclusion is similar. (iii) The results for types 1, 3, 5 is similar to types 2, 4, 6, respectively. (iv) ReFEM behaves the same as for TSFEM. TSFP is clearly better than ReFEM and TSFEM, while ReFP has the best accuracy.

Example 4.3. We now consider different mesh sizes h with fixed time and PML sizes δ. From now on, since the FP-based approaches are better than FEM-based schemes in terms of accuracy and efficiency, we only consider the ReFP and TSFP methods. We fix the mesh size as h = 1 16 and varies the absorbing strength σ and PML sizes δ. Fig. 4.3 shows e rel ∞ (0.25) for ReFP and TSFP with the absorbing functions of types 2, 4, 6. From this figure and additional figures not shown here, we can see that: (i) As δ increases, the error decreases for a fixed value of σ. (ii) The performance of types 1, 3 and 5 are the same as those of types 2, 4 and 6, respectively. (iii) Generally, the errors of types 2, 4 and 6 are smaller than those of types 1, 3 and 5, respectively. (iv) For both ReFP and TSFP, as δ changes, the optimal region of σ is shifted a lot for the absorbing functions of types 1 and 2, while for the other four types, the stability region of σ does not change significantly, especially when δ is large enough. For this reason and together with comment (ii), it is preferable to use an absorbing function of type 4 or 6. Again, ReFP is more accurate than (but comparable with) TSFP. Here, following the last example, we consider the error with respect to different mesh sizes h. We fix the PML sizes δ = 0.25 and vary the mesh size h as well as the absorbing strength σ. Fig. 4.4 illustrates e rel ∞ (0.25) for ReFP and TSFP with absorbing functions of types 2, 4, 6. From this figure and additional figures not shown here, we can see that: (i) As h decreases, the error for the same parameter σ usually decreases for both ReFP and TSFP. (ii) As h decreases, the stability region of σ enlarged significantly. (iii) However, the error is saturating for a fixed thickness δ.

The nonlinear case

In this section, we apply the proposed methods to study the nonlinear Schrödinger equation, i.e., β = 0. Recall that from the previous section the FFT-based schemes outperform the FEM-based methods and absorbing functions of types 4 & 6 work better than other types for each method. Therefore, we only consider the TSFP and ReFP methods with absorbing functions of types 4 & 6 in the following cases.

Example 4.5. We first compare the performance of TSFP and ReFP methods for a manufactured example. To this end, we let β = -1. The trapping potential and initial data are chosen respectively as

V (x) = - 1 2 sech 2 (x -t) sech 2 (y -t) cosh(2(x -t)) + cosh(2(y -t)) , (4.4)
ψ 0 (x) = sech(x) sech(y)e i(x+y) . (

With these parameters, the nonlinear Schrödinger equation admits an outgoing solitary solution which can be solved analytically as ψ(x, t) = sech(x -t) sech(y -t)e i(x+y) . (4.6)

For both methods, we choose L = 32 and ∆t = 10 -3 . The reference solution ψ n ref in (4.1) is taken as the exact solution (4.6) at the grid points at time t = t n . We fix the mesh size as h = 1 8 and vary the absorbing strength σ as well as the size of the PML domain δ. Fig. 4.5 depicts e rel ∞ [START_REF] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF] for the ReFP methods with different PML sizes δ and absorbing functions. Fig. 4. [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with Perfectly Matched Layers[END_REF] shows the e rel ∞ (t) at different times for both methods with PML size δ = 2 and type 6 absorbing function. From these figures and other numerical experiments not reported here for conciseness, we could see that: (i) The performance of the TSFP method is similar as the one of ReFP. (ii) The performances of TSFP & ReFP are similar as for the linear cases shown in the previous sections, i.e., the error decays as the PML size δ increases. However, the optimal region of absorbing strength σ is different from the linear cases. (iii) Absorbing function of type 6 is more accurate than the one of type 4 in the optimal region. Moreover, the optimal region of type 4 function will shift as the PML size δ changes. Therefore, we suggest to use the absorbing function of type 6 for the practical computations. (iv) For both types absorbing functions, usually ReFP is more accurate than TSFP (cf. fig. 4.6). However, due to the nonlinearity, the preconditioner of ReFP (3.24) needs to be rebuilt at each time step, which thus is less efficient than for the TSFP (3.21)- (3.23) where the nonlinear part can be integrated explicitly. (v) Now that the acuracy of TSFP is not far from ReFP, we suggest that the TSFP method is used with the absorbing function of type 6 in practice. Example 4.6. In this example, we apply the TSFP method to simulate the dynamics of the GPE with a strong nonlinearity. To this end, we choose β = 500, ω = 0 and L = 16. First, we prepare the initial data by computing the corresponding ground state of the GPE with trapping potential V (x) = |x| 2 /2. Then, we release the ground states from the trap, i.e., we set V (x) ≡ 0 and simulate the dynamics using TSFP methods. Under this set-up, the ground state is expected to extend and spreads out which would cause problems wherever it is close to the computational domain if a PML technique is not applied. We fix the size of the PML region, the mesh size and discretization time step respectively as δ = 2, h = 1/16, ∆t = 10 -3 .

The absorbing function is chosen as the Type 6 function with absorbing strength σ = 100. To compare, we also report the results obtained by using TSFP without PML (i.e. set σ = 0). 5. Extension to rotating BECs (ω = 0)

Direct PML formulation in cartesian coordinates

Let us consider the full GPE (2.1) with rotating term. Basically, the main difference now is that the z-component of the angular momentum -ωL z = -iω(y∂ x -x∂ y ) is involved into the equation. Extending the previous PML technique in a stable way is not trivial here. One possible direction could be inspired by the case of the convected Helmholtz equation (set in the frequency domain) which writes according to the pressure field p

(1 -M 2 )∂ 2 x p + ∂ 2 y p + 2ikM ∂ x p + k 2 p = f, (5.1)
in D Phys . Here, M = v 0 /c 0 and k = ω/c 0 , are the Mach number (-1 < M < 1) and the wavenumber, respectively. In addition, c 0 is the sound velocity in the fluid, ω > 0 is the pulsation of the wave, f is a source term, and the mean velocity v 0 is subsonic and uniform. Following [START_REF] Bécache | Perfectly mathched layers for the convected Helmholtz equation[END_REF], a well-suited stable PML for the duct problem (along the x-direction) can be introduced as

(1 -M 2 )(α(x)∂ x + iλ(x)) 2 p + ∂ 2 y p + 2ikM (α(x)∂ x + iλ(x))p + k 2 p = f, (5.2)
in D. The function α := (1 + iσ/ω) -1 must satisfy (α) > 0 and (α) < 0, and λ(x) ∈ R. The choice of λ depends on the constants M and k. Extending a similar approach to the rotating GPE is not clear for many reasons: the domain is not a duct but a rectangle, the operator L z is much more complicated since it involves a linear combination of terms x∂ y and y∂ x with variable dependent coefficients (i.e. x & y). For example, a direct PML approach via coordinate stretching leads to the following PML formulation of the rotating GPE: for x ∈ D, t ≥ 0

i∂ t ψ = - 1 2 1 S x ∂ x 1 S x ∂ x + 1 S y ∂ y 1 S y ∂ y ψ + V + β|ψ| 2 ψ -iω y S y S x ∂ x - x S x S y ∂ y ψ, (5.3)
with

S ν = 1 + e iϑν ν              ν Lv σ ν (s)ds, ν > L ν , 0, -L ν ≤ ν ≤ L ν , ν -Lv σ ν (s)ds, ν < -L ν , ν = x, y. (5.4)
However, this PML approach will result in inaccurate computations (not reported here). Blow-up always arises at the interface between PML region (i.e. D PML ) and the physical domain (i.e. D Phys ) for all the proposed absorbing functions and numerical schemes.

PML formulation in rotating Lagrangian coordinates

To get a well-posed PML formulation for the rotating GPE, we consider here an alternative based on a reformulation of the initial problem into a rotating Lagrangian coordinates framework [START_REF] Bao | A simple and efficient numerical method for computing the dynamics of rotating Bose-Einstein condensates via rotating Lagrangian coordinates[END_REF]. For any time t ≥ 0, let R(t) be the orthogonal rotational matrix

R(t) = cos(ωt) sin(ωt) -sin(ωt) cos(ωt) , if d = 2, R(t) =   cos(ωt) sin(ωt) 0 -sin(ωt) cos(ωt) 0 0 0 1   , if d = 3.
(5.5)

It is easy to check that R -1 (t) = R T (t), for any t ≥ 0 and R(0) = I, where I is the identity matrix. For any t ≥ 0, we introduce the rotating Lagrangian coordinates x as [START_REF] Antonelli | On the Cauchy problem for nonlinear Schrödinger equations with rotation[END_REF][START_REF] Bao | A simple and efficient numerical method for computing the dynamics of rotating Bose-Einstein condensates via rotating Lagrangian coordinates[END_REF][START_REF] García-Ripoll | Construction of exact solution by spatial translations in inhomogeneous nonlinear Schrödinger equations[END_REF] x

= R -1 (t)x = R T (t)x ⇔ x = R(t) x, x ∈ R d , (5.6)
and we denote by ψ := ψ( x, t) the wave function in the new coordinates system

ψ( x, t) = ψ (R(t) x, t) , x ∈ R d , t ≥ 0. ( 5.7)
Therefore, the rotating GPE can be rewritten as a GPE without rotation term but with a time-dependent potential

i∂ t ψ( x, t) = - 1 2 ∇ 2 + V (R(t) x, t) + β| ψ| 2 ψ( x, t), (5.8)
setting ∇ 2 = ∂ 2 x + ∂ 2 y . Now, the PML technique can be directly applied to the new GPE (5.8). All the details remain the same by just changing the symbols x, y to x, y, respectively. The TSFP and ReFP (as well as TSFEM and ReFEM) proposed for the non-rotating GPE extend here. For the TSFP, the use of a numerical integration [START_REF] Antonelli | On the Cauchy problem for nonlinear Schrödinger equations with rotation[END_REF][START_REF] Bao | A simple and efficient numerical method for computing the dynamics of rotating Bose-Einstein condensates via rotating Lagrangian coordinates[END_REF][START_REF] García-Ripoll | Construction of exact solution by spatial translations in inhomogeneous nonlinear Schrödinger equations[END_REF] is needed for evaluating the time-dependent potential. Let us remark that a different transformation was applied in [START_REF]Perfectly matched layers for coupled nonlinear Schrödinger equations with mixed derivatives[END_REF] when dealing with cross-derivatives into the nonlinear Schrödinger equation and using a fourth-order spatial finite-difference scheme and a fourth-order additive Runge-Kutta time scheme.

Numerical examples

Example 5.1. Here, we apply the TSFP method to simulate the dynamics of the rotating GPE with strong nonlinearity. To this end, we fix β = 10000 and study two cases: one with rotating frequency ω = 2.5 while the other considers ω = 4. The initial data are prepared by computing the ground states of the GPE with the same values of β & ω under the quadratic-quartic trapping potential:

V (x, y) = 3(x 2 + y 2 ) 2 -4(x 2 + y 2 ) 40 .

(5.9)

The contour plot of the density of the ground states, and hence for the initial data |ψ(x, t = 0)| 2 , for each case is illustrated in the first subfigure in Fig. 5.8 & 5.9, respectively. We then release the ground states from the trapping potential (i.

Conclusion

In this paper, we proposed and implemented in a relatively simple way the PML formulation of the rotating GPE in cartesian coordinates for the time-splitting and relaxation pseudospectral schemes. A thorough numerical study shows that the PML absorbing functions proposed in [START_REF] Bermúdez | An exact bounded perfectly matched layer for timeharmonic scattering problems[END_REF][START_REF] Bermúdez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF] are best suited in terms of accuracy and numerical stability, thanks to the tuning parameters. In addition, a high accuracy of the pseudospectral approximation schemes is observed, outperforming as expected a low-order finite element approximation. In the case of the rotating GPE, a Lagrangian coordinates transformation is used to rewrite the original equation in the framework of non-rotating GPEs with time-dependent potential. Following this approach, the PML formulation is directly extended to the rotating case, showing that the method remains very efficient and accurate for truncating the infinite spatial domain. Future investigations concern the extension of the PML formulations to polar coordinates and to more general GPEs, including dipolar interactions and multi-components problems.

Figure 4 .

 4 Figure 4.1: e rel ∞ (0.24) in Example 4.1 for the schemes ReFEM,TSFEM, ReFP and TSFP (from left to right) with different absorbing strengths σ and PML sizes δ.

Figure 4 . 2 :

 42 Figure 4.2: e rel ∞ (t) at different times t and various absorbing strengths σ for the PML size δ = 0.5 in Example 4.2 for ReFEM, TSFEM, ReFP and TSFP (from left to right) with absorbing functions of types 2, 4 and 6 (from top to bottom)

Figure 4 .

 4 Figure 4.3: e rel ∞ (0.25) vs. σ for TSFP (top) and ReFP (bottom) with absorbing functions of types 2, 4, 6 for hx = hy = 1 16

Figure 4 . 4 :

 44 Figure 4.4: e rel ∞ (0.25) vs. σ for TSFP (top) and ReFP (bottom) for the absorbing function of types 2, 4, 6, with PML size δ = 0.25, and for different mesh sizes hx = hy = h in Example 4.4.

6 Figure 4 . 5 :

 645 Figure 4.5: e rel ∞ (31) for ReFP with different PML sizes δ and absorbing functions of type 4 & 6 in example 4.5.

6 Figure 4 . 6 :

 646 Figure 4.6: e rel ∞ (t) at different time t for ReFP and TSFP with PML sizes δ = 2 and type 6 absorbing functions in example 4.5

Fig. 4 .

 4 [START_REF] Antoine | Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF] shows the contour plots of the density |ψ(x, t)| 2 at different times for the TSFP method with and without PML techniques. From the figure we can see that the waves are reflected back into the physical domain leading to a nonphysical dynamics if no PML is imposed. In contrast, adding a PML region can absorb the outgoing waves well, and the dynamics of the ground states in the physical domain can be well reproduced.

Figure 4 . 7 :

 47 Figure 4.7: Contour plots of |ψ(x, t)| 2 at different times for example 4.6 solved by TSFP without PML (upper) and with PML (Lower). The red solid lines represents the interface of the PML region and Physical domain.

  e. set V (x) ≡ 0) and simulate the dynamics, which is solved by the TSFP method with absorbing function of type 6 in rotating Lagrangian coordinates. For both cases, the physical domain (in Cartesian coordinates) is set to be D Phys = [-16, 16] 2 . To recover the wave function ψ(x, t) in D Phys (which is transformed back from ψ( x, t)), in the rotating Lagrangian coordinates we set the "physical domain" and computational domain as D Phys = [-24, 24] 2 and D = [-27, 27] 2 , respectively. Moreover, the mesh size and time step size are chosen as h = 1/16 and ∆t = 10 -3 . Fig.5.8 (ω = 2.5) and 5.9 (ω = 4) show the contour plots of |ψ(x, t)| 2 at different times in Cartesian coordinates. From the figure, we can see clearly that the TSFP method works very well. Due to the absence of trap, the wave function in both cases will spread out and hence all the vortices will exit and nothing left finally in the physical domain.

Figure 5 . 8 :

 58 Figure 5.8: Contour plots of |ψ(x, t)| 2 at different times for example 5.1 for ω = 2.5. The green solid lines represents the interface of PML region and "Physical domain" D Phy in rotating Lagrange coordinates. The black dashed lines represents the physical domain D Phy in Cartesian coordinates.

Figure 5 . 9 :

 59 Figure 5.9: Contour plots of |ψ(x, t)| 2 at different times for example 5.1 for ω = 4. The green solid lines represents the interface of PML region and "Physical domain" D Phy in rotating Lagrange coordinates. The black dashed lines represents the physical domain D Phy in Cartesian coordinates.

Acknowledgements

The authors acknowledge the support from the Inria associate team BEC2HPC (Bose-Einstein Condensates: Computation and HPC simulation). X. Antoine thanks the French National Research Agency project NABUCO grant No. ANR-17-CE40-0025. Q. Tang thanks the support of the Fundamental Research Funds for the Central Universities (YJ201807).