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

Emmanuel Lorin
email: elorin@math.carleton.ca

A simple pseudospectral method for the computation of the time-dependent Dirac equation with Perfectly Matched Layers

Keywords: Dirac equation, time-splitting, pseudospectral approximation, Perfectly Matched Layers, high-order accuracy

published or not. The documents may come

A simple pseudospectral method for the computation of the time-dependent Dirac equation with Perfectly Matched Layers

Introduction

This paper is concerned with the numerical computation of the time-dependent Dirac equation for physical problems involving delocalized wavefunctions. More specifically, we are interested in the numerical solution to the Dirac equation on a truncated domain with absorbing boundary layers. To this end, we propose a simple combination of a pseudospectral method with Perfectly Matched Layers (PMLs), allowing to consider delocalized wavefunctions, as observed when quantum relativistic particles are subject to strong fields. Thanks to a relatively new and simple Fourier-based discretization of spatial differential operators, it is possible to impose PML for solving the Dirac equation on a bounded domain. Although overall, the Fourier-based method still imposes periodic boundary conditions, the outgoing/incoming waves are in fact mainly absorbed. We think that due to the simplicity of the proposed method, most of existing Fourier-based codes could easily be modified to include the proposed methodology. The Dirac equation is a relativistic wave equation which has gain much attention these past 15 years due to the development of 2-d material, such as graphene [START_REF] Fillion-Gourdeau | Numerical computation of dynamical Schwinger-like pair production in graphene[END_REF][START_REF] Fillion-Gourdeau | Time-domain quantum interference in graphene[END_REF][START_REF] Katsnelson | Chiral tunnelling and the Klein paradox in graphene[END_REF], intense-laser-molecule interaction [START_REF] Faisal | A four-component Dirac theory of ionization of a Hydrogen molecular ion in a super-intense laser field[END_REF][START_REF] Gonoskov | Probing nonperturbative QED with optimally focused laser pulses[END_REF][START_REF] Mocken | Quantum dynamics of relativistic electrons[END_REF][START_REF] Salamin | Relativistic high-power laser-matter interactions[END_REF], in particular for pair production [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF][START_REF] Kjellsson | Alternative gauge for the description of the light-matter interaction in a relativistic framework[END_REF][START_REF] Kjellsson | Relativistic dynamics for hydrogenlike systems[END_REF], or from heavy ion collisions modeling and simulation for quarkantiquark production [START_REF] Gelis | Quark-antiquark production from classical fields in heavy-ion collisions: 1 + 1 dimensions[END_REF][START_REF] Reinhardt | Theory of positron production in heavy-ion collisions[END_REF][START_REF] Wells | Asymptotic channels and gauge transformations of the time-dependent Dirac equation for extremely relativistic heavy-ion collisions[END_REF]. In the same time, there has been a tremendous progress in the development of efficient and accurate computational methods for the solution to the Dirac equation [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF][START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Becker | Solution of the time-dependent Dirac equation by the finite difference method and application for Ca 20+ + U 91+[END_REF][START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timeindependent Dirac equation for diatomic molecules: B-splines without spurious states[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Fillion-Gourdeau | Galerkin method for unsplit 3-d Dirac equation using atomically/kinetically balanced B-spline basis[END_REF][START_REF] Grant | B-spline methods for radial Dirac equations[END_REF][START_REF] Kjellsson | Alternative gauge for the description of the light-matter interaction in a relativistic framework[END_REF][START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF][START_REF] Mocken | Quantum dynamics of relativistic electrons[END_REF][START_REF] Selstø | Solution of the Dirac equation for hydrogenlike systems exposed to intense electromagnetic pulses[END_REF], modeling in particular molecules subject to intense external laser fields. In this spirit, we have developed [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF], a Schwinger-like pair production procedure from the interaction of intense laser fields with heavy molecules occuring at specific resonances [START_REF] Fillion-Gourdeau | Enhanced Schwinger pair production in many-centre systems[END_REF][START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF].

In this paper, we are more specificially interested in the relativistic interaction of atoms, molecules or wavepackets with intense and short laser pulses. The key point is that the laser field actually delocalizes the wavefunction, and the latter can actually interact with the domain boundary. In order to avoid artificial reflections, it is necessary to impose absorbing boundary conditions [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF], or to include in the equation absorbing complex potentials or perfectly matched layers [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF]. Theoretically, this approach benefits from spectral convergence (for smooth functions) and from the simplicity of Fourier-based methods, more generally pseudospectral methods, on bounded domains reducing the periodic boundary condition effect thanks to artificial wave absorption at the domain boundary. Perfectly Matched Layers are now widely used in many engineering and physics simulations codes [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equation[END_REF][START_REF] Bermúdez | An exact bounded perfectly matched layer for time-harmonic 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] Colonius | Modeling artificial boundary conditions for compressible flow[END_REF][START_REF] Farrell | The perfectly matched layer in numerical simulations of nonlinear and matter waves[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] Nissenand | An optimized perfectly matched layer for the Schrödinger equations[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][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] to model exterior domains and to avoid unphysical reflections at the domain boundary. PML for the Dirac equation were developed in [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] and then approximated using a finite difference method. The derivation of high-order absorbing boundary conditions (ABCs) for the Dirac equation were proposed in [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. We also refer to [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for an overview of PMLs and ABCs for quantum wave equations including the Dirac equation and to [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF][START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF][START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Fischer | A B-spline Galerkin method for the Dirac equation[END_REF][START_REF] Guo | Spectral and pseudospectral approximations using Hermite functions: Application to the Dirac equation[END_REF][START_REF] Hammer | A dispersion and norm preserving finite difference scheme with transparent boundary conditions for the Dirac equation in (1+1)d[END_REF][START_REF] Kullie | Dirac-Fock finite element method (FEM) calculations for some diatomic molecules[END_REF][START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF][START_REF] Mocken | FFT-split-operator code for solving the Dirac equation in 2+1 dimensions[END_REF] for different approaches for solving the Dirac equation in real or Fourier space. We insist on the fact that the purpose of this paper is not to provide new PMLs for the Dirac equation, but rather to develop an efficient pseudospectral method for solving Dirac Initial Boundary Value Problems (IBVP).

The combination of pseudospectral methods and PMLs is possible thanks to the following simple. For example in the x-direction, let us denote by ξ x the dual Fourier variable and by F x (resp. F -1

x ) the Fourier (resp. inverse Fourier) transform in x. We consider a as a given x-dependent function. Then, for any H 1 -function f , it is possible to formally rewrite a(x)∂ x f (x) as F -1

x a(x)iξ x F x (f )(ξ x ) (x) on an unbounded domain (pseudodifferential operator representation [START_REF] Taylor | Pseudodifferential Operators[END_REF]). Let us remark that the latter is a real space function, although the derivative is approximated using the Fourier transform, and the function (x, ξ x ) → ia(x)ξ x is nothing but the symbol of a(x)∂ x . In practice the function a involves the stretching coordinates function modeling the PML on a bounded domain allowing real space non-reflecting conditions at the domain boundary.

The time-dependent Dirac equation under consideration reads [START_REF] Itzykson | Quantum Field Theory[END_REF] i∂ t ψ(t, x) = Hψ(t, x),

where ψ(t, x) is the time and coordinate dependent four-spinor, and H is the Hamiltonian operator. The latter is given by

H = α • [cp -eA(t, x)] + βmc 2 + I 4 V (t, x), (2)
where the momentum operator is p = -i∇. More specifically, the Dirac equation reads

i∂ t ψ(t, x) = α x -ic∂ x -eA x (t, x) +α y -ic∂ y -eA y (t, x) +α z -ic∂ z -eA z (t, x) +βmc 2 + I 4 V (t, x) ψ(t, x),
where ψ(t, x) ∈ L 2 (R 3 ) ⊗ C 4 is the time and coordinate (x = (x, y, z)) dependent fourspinor. In (3), A(t, x) represents the three space components of the electromagnetic vector potential, V (t, x) = eA 0 (t, x) + V nuc. (x) is the sum of the scalar and interaction potentials, e is the electric charge (with e = -|e| for an electron), I 4 is the 4 × 4 unit matrix and α = (α ν ) ν=x,y,z , β are the Dirac matrices. In this work, the Dirac representation is used, where

α ν = 0 σ ν σ ν 0 , β = I 2 0 0 -I 2 . (3)
The σ ν are the usual 2 × 2 Pauli matrices defined as

σ x = 0 1 1 0 , σ y = 0 -i i 0 and σ z = 1 0 0 -1 , (4)
while I 2 is the 2 × 2 unit matrix. Note that the light velocity c and fermion mass m are kept explicit in Eq. ( 2), allowing to adapt the method easily to natural or atomic units (a.u.). The paper is organized as follows. We recall the basics of PMLs for the Dirac equation in Section 2. Section 3 is dedicated to the derivation of the pseudospectral approximation applied to IBVP for the time-dependent Dirac equation on a bounded domain with PML. Numerical experiments are presented in Section 4. We finally conclude in Section 5.

PML for the Dirac equation

The time-dependent Dirac equation is considered on a bounded (truncated) physical domain denoted by D Phy . We i) add a layer D PML surrounding D Phy , and ii) stretch the coordinates in all the directions. The overall computational domain is next defined by: D = D Phy ∪ D PML . We refer to [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for the construction of PMLs for quantum wave equations and more specifically to [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] for the derivation and analysis of PMLs for the Dirac equation.

The starting point is the stretching of the real coordinates in the complex plane. In [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF], the author uses the following variables

ν = ν + i ω ν L * ν σ(s)ds,
where ω is the dual Fourier variable to t and the so-called absorbing function σ is such that σ(ν) = 0 if s < L * ν , for some imposed L * ν < L ν . Then, the partial derivative ∂ ν , with ν = x, y, z, is shown to be formally transformed into

∂ ν ψ(t, •) → ∂ t ∂ t + σ ν (ν) ∂ ν ψ(t, •) = 1 2π R ωe iωt ω -i σ ν (ν) ∂ ν F t ψ(ω, •)dω, (5)
where F t (resp. F ν ) denotes the Fourier transform with respect to t (resp. ν) and

σ ν (ν) = σ(|ν| -L ν ), L * ν |ν| < L ν , 0, |ν| < L * ν .
Let us remark that the transformation (5) can also be formally rewritten as

∂ ν ψ(t, •) → 1 1 + σ ν (ν)∂ -1 t ∂ ν ψ(t, •).
In the present paper, we rather consider a more simple change of variables [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] involving only the space variable

ν = ν + e iθ ν L * ν σ(s)ds,
where θ ∈ (0, π/2). We then define S ν (ν) := 1 + e iθν σ(ν),

with ν = x, y, z. In the following, the partial derivatives are simply transformed into

∂ ν → 1 S ν (ν) ∂ ν = 1 1 + e iθν σ ν (ν) ∂ ν , (7)
where σ vanishes in D Phy and S ν is equal to 1. From now on, let us consider the transformation [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF], and the associated new Hamiltonian

H PML = α • [cT • p -eA(t, x)] + βmc 2 + I 4 V (t, x), (8)
where

T (x, y, z) = T x (x), T y (y), T z (z) T is such that T ν (ν) = 1 S ν (ν) with ν = x, y, x . (9)
Several types of functions can be selected. An exhaustive study of the absorbing functions σ is proposed in [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF] for (nonlinear) Schrödinger equations Type I: σ 0 (ν + δ ν ) 2 , Type II: σ 0 (ν + δ ν ) 3 , Type III: -σ 0 /ν, Type IV:

σ 0 /ν 2 , Type V: -σ 0 /ν -σ 0 /δ ν , Type VI: σ 0 /ν 2 -σ 0 /δ 2 ν ,
for ν = x, y, z and δ ν = L ν -L * ν > 0. The main difficulty from the pseudospectral point of view is the space-dependence of the coefficients T ν which prevents the direct application of the Fourier transform on the equation. A simple trick will however allow for combining the efficiency of the pseudospectral method and the computation of the non-constant coefficient Dirac equation.

PSeudospectral-PML (PS-PML) method

A splitting of the Dirac Hamiltonian is used below. Although operator splitting is not fundamentally required in the proposed methodology, it allows to simplify the implementation of the method while keeping a good accuracy. The directional splitting provides a very efficient parallelization (per slice) of the algorithm. On the other hand, as the Dirac matrices do not commute, the splitting naturally induces errors (including possible numerical dispersion). As a consequence, the chosen splitting should have an order at least equal to the order of the time integration solvers. Based on (2), we define the operators

A = -icT x (x)α x ∂ x , (10)
B = -icT y (y)α y ∂ y , (11)
C = -icT z (z)α z ∂ z , ( 12
) D = βmc 2 + I 4 V (t, x) -eα • A(t, x). (13)
From time t n to t n+1 , we then successively solve [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF] (the x-dependence in the wavefunction argument is removed for notational convenience)

i∂ t ψ (1) (t) = Aψ (1) (t), ψ (1) (t n ) = ψ n , t ∈ [t n , t n 1 ) ( 14
)
i∂ t ψ (2) (t) = Bψ (2) (t), ψ (2) (t n ) = ψ (1) (t n 1 ), t ∈ [t n , t n 2 ) ( 15
)
i∂ t ψ (3) (t) = Cψ (3) (t), ψ (3) (t n ) = ψ (2) (t n 2 ), t ∈ [t n , t n 3 ) ( 16
)
i∂ t ψ (4) (t) = Dψ (4) (t), ψ (4) (t n ) = ψ (3) (t n 3 ), t ∈ [t n , t n+1 ) (17) and ψ n+1 = ψ (4) (t n+1 ), ( 18
)
where

t n i -t n = ∆t for i ∈ {1, • • • , 3}.
In the next subsections, we detail the space-time approximations of equations ( 14)- [START_REF] Bechouche | On the asymptotic analysis of the Dirac-Maxwell system in the nonrelativistic limit[END_REF]. Cylindrical coordinates could be used as well (see e.g. [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF]).

Pseudospectral approximation in space

We consider for convenience the 3-dimensional system in cartesian coordinates, for ν = x, y, z in the domain

D = [-a x , a x , ] × [-a y , a y ] × [-a z , a z ], i∂ t ψ(t) = -icα ν T ν (ν)∂ ν ψ(t), ψ(t n ) = ψ n , t ∈ [t n , t n+1 )
and we use the same notations as in [START_REF] Antoine | High-order IMEX-spectral schemes for computing the dynamics of systems of nonlinear Schrödinger/Gross -Pitaevskii equations[END_REF]. We first diagonalize

α ν = Π ν ΛΠ † ν , where Λ =     1 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 -1     , and Π † ν = Π T ν .
The matrices Π ν are defined as follows

Π x = 1 √ 2     0 1 1 0 1 0 0 -1 1 0 0 1 0 1 -1 0     , (19)
Π y = 1 √ 2     0 -i -i 0 1 0 0 1 -i 0 0 i 0 1 -1 0     , (20)
Π z = 1 √ 2     1 0 0 -1 0 -1 -1 0 1 0 0 1 0 1 -1 0     . ( 21
)
We set φ := Π † ν ψ, which then satisfies

i∂ t φ(t) = -icΛT ν (ν)∂ ν φ(t), φ(t n ) = Π † ν ψ n , t ∈ [t n , t n+1 ).
We denote the set of grid-points by

D Nx,Ny,Nz = x k 1 ,k 2 ,k 3 = (x k 1 , y k 2 , z k 3 ) (k 1 ,k 2 ,k 3 )∈O NxNy Nz with O NxNyNz = k 1 = 0, • • • , N x -1; k 2 = 0, • • • , N y -1; k 3 = 0, • • • , N z -1 .
Then, let us introduce the following mesh sizes

x k 1 +1 -x k 1 = h x = 2a x /N x , y k 2 +1 -y k 2 = h y = 2a y /N y , z k 3 +1 -z k 3 = h z = 2a z /N z . (22)
The corresponding discrete wavenumbers are defined by ξ p,q,r := (ξ p , ξ q , ξ r ), where

ξ p = pπ/a x with p ∈ {-N x /2, • • • , N x /2 -1}, ξ q = qπ/a y with q ∈ {-N y /2, • • • , N y /2 -1} and ξ r = rπ/a z with r ∈ {-N z /2, • • • , N z /2 -1}.
In the sequel of the paper, we denote by φ ( ) , with ∈ {1, 2, 3, 4}, the th component of the spectral approximation of φ = Π † ν ψ. We also use the notation φ ( )

k 1 (t, y, z) = φ ( ) (t, x k 1 , y, z), φ ( ) k 2 (t, x, z) = φ ( ) (t, x, y k 2 , z) and φ ( ) k 3 (t, x, y) = φ ( ) (t, x, y, z k 3 ). The partial Fourier coefficients are such that                        φ ( ) p (t, y, z) = Nx-1 k 1 =0 φ ( ) k 1 (t, y, z)e -iξp(x k 1 +ax) , φ ( ) q (t, x, z) = Ny-1 k 2 =0 φ ( ) k 2 (t, x, z)e -iξq(y k 2 +ay) , φ ( ) r (t, x, y) = Nz-1 k 3 =0 φ ( ) k 3 (t, x, y)e -iξr(z k 3 +az) .
We can then introduce the inverse partial Fourier pseudospectral approximations, in the x-, y-and z-directions, respectively,

                         φ ( ) k 1 (t, y, z) = 1 N x Nx/2-1 p=-Nx/2 φ ( ) p (t, y, z)e iξp(x k 1 +ax) , φ ( ) k 2 (t, x, z) = 1 N y Ny/2-1 q=-Ny/2 φ ( ) q (t, x, z)e iξq(y k 2 +ay) , φ ( ) k 3 (t, x, y) = 1 N z Nz/2-1 r=-Nz/2 φ ( ) r (t, x, y)e iξr(z k 3 +az) .
We then denote

φ := { φ } 1 4 , (23)
and define the approximate first-order partial derivatives such that for any ν

∂ ν φ ( ) (t n , x k 1 ,k 2 ,k 3 ) ≈ [[∂ ν ]]φ ( ) k 1 ,k 2 ,k 3 and [[∂ x ]]φ ( ) k 1 ,k 2 ,k 3 := 1 N x Nx/2-1 p=-Nx/2 iξ p φ ( ) k 2 ,k 3 p e iξp(x k 1 +ax) , [[∂ y ]]φ ( ) k 1 ,k 2 ,k 3 := 1 N y Ny/2-1 q=-Ny/2 iξ q φ ( ) k 1 ,k 3 q e iξq(x k 2 +ay) , [[∂ z ]]φ ( ) k 1 ,k 2 ,k 3 := 1 N z Nz/2-1 r=-Nz/2 iξ r φ ( ) k 1 ,k 2 r e iξr(x k 3 +az) . (24)
In the following, the index h will be used (e.g. in

φ n h = {φ n k 1 ,k 2 ,k 3 } k 1 ,k 2 ,k 3 )
to denote a spectral approximation to a given wavefunction (e.g. φ n ). This discretization not only allows to select the spatial step as large as wanted, but it also preserves the very high spatial accuracy, the parallel computing structure and the scalability of the split method developed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF].

Time-Splitting PSeudospectral method with PML (TSSP-PML)

As proposed before, we use a splitting of the Time-Dependent Dirac Equation into four time-dependent systems [START_REF] Balay | Efficient management of parallelism in object oriented numerical software libraries[END_REF], [START_REF] Bandrauk | Improved exponential split operator method for solving the time-dependent Schrödinger equation[END_REF], [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF] and [START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF]. We also define x h = (x h , y h , z h ) as the nodes of a real-space grid. At iteration time n, the approximate wave function is denoted by ψ n h . From a stability point of view, an implicit scheme should be implemented. The TSSP (Time-Splitting pseudoSPectral) algorithm then reads 1. First step: integration of the generalized transport equation in Fourier space in the x-direction. One sets φ n h := Π † x ψ n h , and the system

∂ t φ + cT x (x)Λ∂ x φ = 0, φ(t n , •) = Π † x ψ(t n , •),
is approximately solved by

φ n 1 h = φ n h -c∆tT x (x h )Λ[[∂ x ]]φ n 1 h ,
where [[•]] is defined in [START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF]. We then deduce ψ n 1 h = Π x φ n 1 h . This gives an approximation of Eq. ( 14), where the corresponding operator is denoted by P (x) h (∆t). 2. Second step: integration of the generalized transport equation in Fourier space in the y-direction

∂ t φ + cT y (y)Λ∂ y φ = 0.
We set χ n 1 h := Π † y ψ n 1 h and we approximately solve (15) by using the following scheme

φ n 2 h = χ n 1 h -c∆tT ( y y h )Λ[[∂ y ]]χ n 2 h .
We then obtain ψ n 2 h = Π y φ n 2 h . 3. Third step: setting χ n 2 h := Π † z ψ n 2 h , we apply the same procedure in the z-direction, and compute ψ n 3 h . The corresponding operator is denoted by P (z) h (∆t). 4. Fourth step: integration of the source from t n to t n+1

ψ n+1 h = T exp -i t n+1 tn dτ βmc 2 -eα • A h (τ ) × exp -ie t n+1 tn dτ V h (τ ) ψ n 3 h , (25)
where A h is the approximate electric potential. The corresponding operator is denoted by Q n h (∆t). This splitting scheme is implicit, and can be compactly expressed as a composition of the four operators

ψ n+1 h = Q n h (∆t)P (z) h (∆t)P (y) h (∆t)P (x) h (∆t)ψ n h .
In practice, a second-order splitting scheme should be implemented, such as

ψ n+1 h = P (x
Thanks to the operator [[∂ ν ]], it is possible to implement the real-space approximation of the overall wave function on a bounded domain imposing in particular PMLs, resulting in the TSSP-PML scheme.

Mathematical analysis

In this section, we discuss the stability and convergence of the pseudospectral method developed in Subsection 3.2. The proof mainly relies on [START_REF] Bardos | Stability and spectral convergence of Fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method[END_REF][START_REF] Goodman | On the stability of the unsmoothed Fourier method for hyperbolic equations[END_REF]. To simplify the presentation, we assume that the computational domain is D = T(0, 2π) 3 . The convergence analysis contains two main parts. We first consider the one-dimensional case for φ 0 = φ(0, •) ∈ H s (T(0, 2π)) 4 such that Λ = diag(1, 1, -1, -1) ∈ R 4×4 , and

∂ t φ(t, x) + cT (x)Λφ(t, x) = 0, x ∈ T(0, 2π), T x ∈ C 1 [0, 2π] ,
where T x (x) = 1/S x (x) is the inverse of an absorbing function S x as defined in [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF], which is assumed to be differentiable with bounded derivative. Using the same notations as [START_REF] Bardos | Stability and spectral convergence of Fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method[END_REF], we here define

P N [φ] := |k| N 1 2π T φ(x)e -ikx dx e ikx , Ψ N [φ] = φ ,
where φ is given by [START_REF] Brabec | Intense few-cycle laser fields: Frontiers of nonlinear optics[END_REF]. Thus, we have

P N [φ] = Ψ N [φ] + A N [φ]
, where the aliasing term

A N [φ](x) = |k| N |l| 1 1 2π T φ(x)e -i(k+N l)x dx e ikx
which is responsible for a loss of stability/accuracy is such that [55]

A N [φ] (H s ) 4 N s-r φ (H r ) 4 , r > s > 1/2 .
This inequality illustrates that the loss of accuracy is directly related to the smoothness. Denoting φ N (x, t) = |k| N u k (t)e ikx , from (29) the scheme can be rewritten as

∂ t φ N (x, t) = -cT x (x)Λ[[∂ x ]]φ N (x, t) = -cT x (x)Λψ N ∂ x φ = -cT x (x)Λ∂ x ψ N φ .
Let us remark that we can obtain the first step of the scheme in Subsection 3.2 as

∂ t φ N (x, t) + cT x (x)Λ∂ x P N [φ N = 0 .
The stability analysis then reads as follows. First, we have

φ N (•, t) 2 (L 2 ) 4 = -2Λ cT x (x)φ † N (x, t)∂ x P N [φ N dx -2Λ cT x (x)φ N ∂ x A N [φ N ]dx with 2Λ cT x (x)φ † N (x, t)∂ x P N [φ N dx c T x ∞ φ N 2 (L 2 ) 4 .
Finally, according to [START_REF] Bardos | Stability and spectral convergence of Fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method[END_REF][START_REF] Goodman | On the stability of the unsmoothed Fourier method for hyperbolic equations[END_REF], we conclude that there exists a constant K s > 0 such that

φ N (•, t) (L 2 ) 4 N e cKs T x L ∞ t φ N (•, 0) (L 2 ) 4 .
We then have Lemma 3.1. For s > 2, there exists a constant K s > 0 such that:

φ N (•, t) -Ψ N [φ] (L 2 ) 4 e cKs T x L ∞ t N 1-s φ(•, 0) (H s ) 4 + N 2-s max φ(•, τ ) (H s ) 4 . ( 27
)
We now consider the three-dimensional split system. For any transition matrix Π ν defined in [START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF], Π † ν ψ = φ , the above estimate ( 27) is valid for the first to third steps of the scheme (3.2). Finally, the contribution of the fourth step ( 25) is straightforward in the error estimate. We then have the following proposition. Proposition 3.1. Let us assume that ψ 0 = ψ(0, •) ∈ H s (T(0, 2π)) 4 and that for ν =

x, y, z, T ν := 1/S ν ∈ W 1,∞ (0, 2π). Then, the numerical scheme derived in Subsection (3.2) is such that

ψ N (•, t) -Ψ N [ψ] (L 2 ) 4 exp cK s T x L ∞ + A ∞ t × N 1-s ψ(•, 0) (H s ) 4 + N 2-s max ψ(•, τ ) (H s ) 4 , (28)
for some constant K s > 0.

This proposition shows the spectral convergence of the derived pseudospectral method, in the sense that the convergence increases with the order of regularity of the solution. It also precises the dependence on the regularity of the absorption functions S ν , and more specifically the regularity of T ν but also on ii) T ν L ∞ . Assuming that T ν is smooth, we get

d dν 1 S ν (ν) = - S ν (ν) S 2 ν (ν) = - e iθ σ ν (ν) 1 + e iθ σ ν (ν) 2 ,
where σ is differentiable or piecewise differentiable. More specifically,

• for the Type-I PML, σ ν is differentiable

σ ν (ν) = 2sgn(ν)σ 0 (|ν| -L ν + δ ν ), L * ν |ν| < L ν , 0, |ν| < L * ν , where δ ν = L ν -L *
ν . Thus, we obtain that d dν

1 S ν (ν) = 2σ 0 L ν -|ν| 1 + e iθ σ ν (ν) 2 2σ 0 δ ν .
• for the Type-II PML, σ ν is also differentiable and

d dν 1 S ν (ν) 3δ 2 ν .
• for the Type-IV/VI PML, σ ν is only piecewise differentiable

σ ν (ν) =    -2sgn(ν) σ 0 (|ν| -L ν ) 3 , L * ν < |ν| < L ν , 0, |ν| < L * ν . Therefore, for L * ν < |ν| < L ν , one gets d dν 1 S ν (ν) = 2σ 0 L ν -|ν| 3 1 + e iθ σ ν (ν) 2 .
Let us note that for |ν| close to L ν we have

d dν 1 S ν (ν) 2σ 0 (L ν -|ν|) , while for |ν| close to L * ν one gets d dν 1 S ν (ν) 2σ 0 (L ν -L * ν ) .
For the Type-IV-VI (as well as Type-III-V) PML, the piecewise differentiability can play a role in the loss of precision of the overall method (see inequality [START_REF] Chew | Perfectly matched layers for elastodynamics: A new absorbing boundary condition[END_REF]). However, this was not numerically observed in our tests. In [START_REF] Bardos | Stability and spectral convergence of Fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method[END_REF], the authors propose to perform a de-aliasing thanks to a 2/3-smoothing method by using the factor function σ (see 2.12a in [START_REF] Bardos | Stability and spectral convergence of Fourier method for nonlinear problems: on the shortcomings of the 2/3 de-aliasing method[END_REF]). More precisely, for one-dimensional evolution equations

∂ t φ N (x, t) = -cT x (x)Λ[[∂ x ]]Sφ N (x, t) ,
where Sφ N = |k| 2N/3 σ k u k (t)e ikx , we expect that, for some s > 1,

φ N (•, t) -φ (L 2 S ) 4 e Ks T x L ∞ t N -s φ(•, 0) (H s ) 4 + N 2-s max φ(•, τ ) (H s ) 4 .
We have proposed an implicit scheme to solve the first 3 steps to Subsection 3.2, and which is based on the solution to systems of the form:

φ n+1 N (x, t) = φ n N (x, t) -cT x (x)Λ∂ x Ψ N φ n+1 .
In the following, we also propose an alternative scheme which avoids the numerical computation to linear systems due to the implicit time-discretization. For example, for the First step, we can simply instead solve

∂ t φ + cT x (x)Λ∂ x φ = 0, φ(t n , •) = Π † x ψ(t n , •), (29)
by using

φ n 1 h := 1 S x (x h ) ψ n 1 h + 1 - 1 S x (x h ) φ n h , (30)
where ψ n 1 h is solution to

ψ n 1 h = φ n h -c∆tΛ[[∂ x ]]φ n h . (31)
In other words, φ n 1 h is an approximation to

φ(t n 1 , x) = 1 S x (x) F -1 x 1 -ic∆tξ x F x (φ)(t n 1 , ξ x ) + 1 - 1 S x (x) φ(t n , x).
This approximation is justified as follows. From time t n to t n 1 , we integrate (29) and get

φ(t n 1 , x) = φ(t n , x) - c∆t S x (x) F -1 x iξ x F x (φ)(t n , ξ) + O(∆t 2 ) = 1 S x (x) F -1 x 1 -ic∆tξ x F x (φ)(t n , ξ) + 1 - 1 S x (x) φ(t n , x) + O(∆t 2 ) .
As a consequence, we trivially have Proposition 3.2. The approximation (30)-( 31) is consistent with (29).

Numerical simulations

We propose some numerical experiments in 2-d, illustrating the efficiency of the combined PML-TSSP method with a second-order operator splitting. We are mainly interested in the evolution of wavepacket either subject to a static potential or to an external dynamic field. In our simulations, we have used several libraries which are listed below and consider the atomic units. We compare the effects of different absorption functions with various parameters. We can however not conclude regarding the optimal choice, although Type-IV functions are often among the most efficient. Let us note however that from a theoretical point of view (see Subsection 3.3) the accuracy of the pseudospectral method depends on the regularity of the absorption functions S ν , and more specifically the regularity of 1/S ν .

Technical details about the code and parallel computing aspects

The Fourier transforms are performed with the sequential and parallel version 3.3.4 of fftw. Standard linear algebra libraries, gsl (version 1.9) and blas were also used. Finally openMPI 1.6.3 was used for message passing with non-blocking communication. The code was implemented in C++ and the compiler is the version 4.7.0 of gcc. The tests are performed by using a C++-code with MPI-library. The method is implemented on the cluser mammouth-parallel II from the RQCHP. The total processing power of this machine is 333 400 GFlops, and possesses 39648 cores: 3216 processors AMD Opteron 12 cores at 2.1 GHz, and 88 processors AMD Opteron 12 cores 2.2 GHz. The total memory is 57.6 TB and the computer-networking communications is Infiniband QDR (4 GB/sec). Regarding the communication between processes, by default mammouth II will select the closest processors in the same node, then the processors from the nodes in the same topological ring, and finally through different switch-levels. Notice that for a low number of processors, the communication between the processes within the same node is done thanks to a shared memory. The pseudospectral method for solving the Dirac equation solver without PML is detailed in [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF], and in Section 4. Let us recall the main features of the code in 2-d. In the (x, z) coordinates, we proceed by alternating the directions. We first decompose the domain by layer (slice) in the z-direction, then

• we successively perform the evolution (FFT) sequentially in the x-direction, and by layer in the z-direction. Each processor manages one layer in z. A perfect scaling for this step is expected as it does not require any transmission between nodes.

• for all x, we perform the evolution in the z-direction using the parallel FFT (fftw). The performance of this step is then fully dependent on the parallelization of the one-dimensional FFT.

Notice that the presence of the PML does not deteriorate the efficiency of the overall method compared to usual FFT-methods. In particular, the computational complexity is simple to established, and is given by C after N T time iterations

C = O N T f N x N z log(N x N z ) .
The scalability and performance of the TSSP method without PML is fully studied in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], where in particular it is shown a good speed-up of this method.

Numerical experiments

The second-order TSSP-PML method ( 26) is implemented, where the discrete operators [[•]] given by [START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF] are used in the directions ν = x, z. We then choose S ν as follows, for

ν = x, z S ν (ν) = 1, |ν| < L * ν , 1 + e iθ σ(|ν| -L v ), L * ν |ν| < L ν . ( 32
)
Test 1 : wavepacket subject to an external laser field. In this first test, we consider a wavepacket

ψ(x, y, z, 0) = N [1, 0, 0, 0] T exp -(x 2 + y 2 + z 2 )/∆ 2 × exp(ik x x), (33)
subject to an external potential and propagating with a fixed wavenumber k = (5, 0), and ∆ = 128 in [START_REF] Piazza | Extremely highintensity laser interactions with fundamental quantum systems[END_REF]. The time step is given by ∆t = 4.56 × 10 -4 a.u. The physical domain is D = [-8, 8] 2 and is discretized with N x × N z = 256 2 grid points. We impose a linearly polarized electric field such that A x is identically null, and A z (t, x, z) is as in [START_REF] Drake | Application of discrete-basis-set methods to the Dirac equation[END_REF], with A 0 = 100/137, ω = 100, n = 1, n = 19 and T f = 1.824a.u.. As an example, we report in Fig. 1 the electric potential A z at time T = 0.456a.u. and 1.824a.u. and the initial data. In Figs. 2, we compare the evolution of the density d(T, x, z) = 4 i=1 |ψ i (T, x, z)| 2 as a function of time, with i) periodic boundary conditions with PML of ii) type III and iii) type IV, with σ 0 = 10 -2 , θ = π/4 as well as iv) type V and VI with δ ν = L ν -L * ν . The snapshots correspond to times t = 0.456a.u., 1.14a.u., 1.368a.u. and T = 1.824a.u., for respectively 1000, 2500, 3000 and 4000 time iterations. The best absorption is obtained with the Type-VI absorbing function, although the results are relatively close with all the absorbing functions. We report in Fig. 3, the maximum of the density (in logscale) as a function of the time iterations for 6 different absorbing functions, showing that they all provide relatively good properties for truncating the computational domain.

Test 2 : 3-nucleus system. In this test, a Gaussian wavepacket centered at (0, 0) Fig. 4 (Left), is injected in a 3-nucleus interaction potential Fig 4

(Middle) V (x, z) = - Z A (x -x A ) 2 + (z -z A ) 2 + 2 - Z B (x -x B ) 2 + (z -z B ) 2 + 2 - Z C (x -x C ) 2 + (z -z C ) 2 + 2 (34)
with

Z A = Z B = Z C = 10a.u., (x A , z A ) = (-2, 0), (x B , z B ) = (2, 0) and (x C , z C ) = (0, -2).
The overall computational domain is D = [-6, 6] 2 , with N x × N y = 128 2 grid points and ∆t = 3 × 10 -3 a.u. Notice that this time step does not allow for the capture of very small time scale effects, such as the zitterbewegung [START_REF] Thaller | The Dirac equation[END_REF]. As we are here interested in larger scales effects, more specifically the absorption layer effect, we allow a larger time step. For practical application, we however suggest to use, as usual, a time step 1/2mc 2 (≈ 2.67×10 -5 ) to get a better precision. We report on Fig. 4 (Right) the total density at time T = 3×10 -1 a.u, i.e. d(T, x, z), for (x, z) ∈ D, where d(T, x, z) as the Dirac wave function. Without external excitation the wavepacket is then trapped by the potential. We then plug an external laser field which drives the wavepacket: a linearly polarized electric field is imposed such that A x is identically null, and A z (t, x, z) is defined as follows

= 4 i=1 |ψ i (T, x, z)| 2 , setting ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 )
A z (t, x, z) = A 0 cos ω(x/c -t) f ω(x/c -t) , (35)
with A 0 = 1000/137a.u., ω = 1, and where the envelope function is linear, then constant, and finally linear. More specifically f is defined as

f (t) =          ω nπ t, t ∈ [0, nπ/ω], 1, t ∈ [nπ/ω, (n + n )π/ω], ω nπ t - (2n + n )π ω , t ∈ [(n + n )π/ω, (2n + n )π/ω],
where 2n + n denotes the total number of half-cycles. In the test, we took n = n = 1. We report the z-component of the external laser potential in Fig. 5 at final time. We also include a non-null wavenumber in the initial wavepacket, k = (10, 0), out of the ions influence. The wavepacket will eventually be absorbed by the PML unlike the solution without PML. At times t = 3 × 10 -2 a.u., t = 6 × 10 -2 a.u., t = 3 × 10 -1 a.u. and t = 3.6 × 10 -1 a.u., we compare 3 solutions Fig. 6, 7, 8, 9: 1) the solution of reference computed on [-18, 18] × [-18, 18], 2) the solution without PML, and 3) the solution with Type V-PML and σ 0 = 10 -2 , θ = π/4, and δ ν = 1.2. At time t = 6 × 10 -2 a.u., the 3 solutions look identical, as the wavefunction has not reached the boundary of the computational domain [-6, 6] × [-6, 6]. At t = T = 3×10 -1 a.u., we however clearly see that the wavefunction of reference has mainly left the region [-6, 6] × [-6, 6], while the PML solution was almost totally absorbed and the solution without PML, was maintained in the computational domain due to non-absorbing boundary layers. This is also illustrated on Fig. 10, where the 2 -norm as function of time in the D is represented for the 3 solutions. We see that, unlike the solution without PML, the PML-solution makes decrease the 2 -norm, thanks to the absorbing layers. We also report in Fig. 10 (Top) the solution of reference at time T = 0.3a.u. in a larger domain [-18, 18] × [-18 , 18]. By comparison, we see that on this example, the solution of reference has also a decreasing 2 -norm in the zone D ( t, d Ref |D (t, •) 2 ), see Fig. 10 where the 2norm is represented in logscale as a function of time (Bottom). We also report the solution of reference in the domain [-18.18] × [-18, 18] at time T = 0.6a.u.

Test 3 : evolution of wavepacket subject to a weakly nonlinear potential. We consider the two-dimensional Dirac equation with a nonlinear potential

V (x, z, ψ) = - Z A (x -x A ) 2 + (z -z A ) 2 + 1 - Z B (x -x B ) 2 + (z -z B ) 2 + 1 + |ψ| 2 , (36)
with Z A = Z B = 2, (x A , z A ) = (-1, 4) and (x B , z B ) = [START_REF] Ackad | Numerical solution of the Dirac equation by a mapped Fourier grid method[END_REF][START_REF] Antoine | High-order IMEX-spectral schemes for computing the dynamics of systems of nonlinear Schrödinger/Gross -Pitaevskii equations[END_REF]. The initial density with wave vector k = (2, 10) is represented in Fig. 11. More specifically, the initial data is given by where ∆ = 128 in [START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF] and N is a normalization coefficient, such that ψ(•, 0) (L 2 (D)) 4 = 1.

ψ(x, y, z, 0) = N [1, 0, 0, 0] T exp -(x 2 + y 2 + z 2 )/∆ 2 × exp(ik x x), (37)
The overall computational domain is D = [-8, 8] 2 , with N x × N y = 256 2 grid points and ∆t = 4.56 × 10 -4 a.u.. We report on Fig. 12 the total density at time T = 0.912a.u. in logscale, i.e. log(d(T, x, z), for (x, z) ∈ D, where d(T, x, z) = 4 i=1 |ψ i (T, x, z)| 2 , setting ψ = (ψ 1 , ψ 2 , ψ 3 , ψ 4 ) as the Dirac wave function. At that time, the wavepacket passed through the boundary of the computational domain at z = ±8a.u.. The nonlinearity is numerically treated explicitly. The use of the logscale allows to fairly report the accuracy of the PML. In [START_REF] Desclaux | Computational approaches of relativistic models in quantum chemistry[END_REF], we take L ν = 8a.u., L * ν = 0.8L ν a.u., θ = π/4 (similarly to the Schrödinger equation [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF][START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF]), and we pick the Type-IV PML, i.e. σ : ν → σ 0 /ν 2 (singular profile [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods. Application to rotating Bose-Einstein condensates[END_REF]). We represent the solution with different values of σ 0 . These results show that the PML properly absorbs except for σ 0 = 5, with reflection magnitude as weak as 10 -6 for σ 0 = 10 -2 . Let us remark that other PMLs could be considered as well [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF], but recall that we "only" intend here to prove the feasibility of the method for various PMLs. For completeness and by comparisons with results obtained in Fig. 12, we also report in Fig. 13 the solution to the Dirac equation in D with zero Dirichlet boundary condition by using the real-space quantum lattice Boltztmann method proposed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF] and illustrating the total wave reflection at the domain boundary when PMLs or ABCs are not used. This simple example illustrates the fact without using PML and or low-order ABC, a wavefunction is totally reflected. This is a standard issue when solving the Dirac equation, in particular when studying laser-molecule interaction. We report in Fig. 16 (Top), the maximum of the density as a function of time for σ 0 = 10 -3 , 10 -2 , 10 -1 , 1, 5 in logscale with the PML of Type IV for θ = π/4. We clearly see that σ 0 = 10 -2 provides the best absorption properties. We next compare in Fig. 14, for σ 0 = 10 -2 , the reflection for θ = π/16, θ = π/8, θ = π/4 and θ = π/2 showing the importance of properly selecting θ; the best absorption is obtained with θ = π/16. We report in Fig. 16 (Bottom), the maximum of the density as a function of the time iteration for different values of θ = π/16, π/8, π/4, π/2 with Type IV-PML. We again see an optimal performance with θ = π/16 but θ = π/8 or θ = π/4 also provides some good results, meaning that the PML is relatively stable with θ. In the last test, we compare in Fig. 15, the efficiency of the PML for four different absorbing functions (Types I to IV), with respectively σ(ν) = σ 0 (ν + δ ν ) 2 , σ(ν) = σ 0 (ν + δ ν ) 3 , σ(ν) = -σ 0 /ν and σ(ν) = σ 0 /ν 2 . We see that the best results are obtained for the Type-IV PML, a similar quality being also obtained for the Types V and VI PMLs (not reported here).

These tests clearly show the relevance of the combination of the TSSP method with PMLs, as we simultaneously benefit from the efficiency and accuracy of FFTs and the high absorption feature of PMLs to avoid reflections at the boundary.

Test 4 : convergence. In this last test, we propose a simple convergence benchmark. A wavepacket with wavenumber k = (5, 0) is propagating in vacuum. We plot the initial density, and final solution with periodic conditions on D = [-8, 8] 2 in Fig. 17 (Left, Middle). We report the supremum of the solution for N x = N y = 32, 256, 320, 512, 640, 768 with the PML of Type-III, for σ 0 = 10 -2 , θ = π/4, fixing the time step to ∆t = 1.8 × 10 -3 a.u.. We report in Fig. 17 This test illustrates that artificial wave-reflecting occurs with real space methods (here a quantum lattice Boltzmann method [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]) are implemented without PML or ABC.

accurate convergence up to a certain precision. Notice that what is reported is the maximum of the solution over the physical domain, including the PML. In other words the smaller the maximum, the better the absorption. The saturation error comes from the limit of the PML accuracy and error in time (second-order splitting). This shows the high accuracy (spectral-like) of the TSSP with PML.

Conclusion

In this paper, we have proposed a simple time-splitting pseudospectral method which allows for the numerical computation of IBVP for the Dirac equation with efficient PML. Some numerical experiments have shown that the method preserves the accuracy of the PML with static or time-dependent external potentials. Different absorbing functions were considered and tested. The implementation strategy is developed for the TSSP approximation, but can be directly extended to other numerical schemes including implicit ones. We think that Fourier-based codes solving the time-dependent Dirac equation can easily be adapted to include the absorbing layers, hence drastically reducing the negative effect of periodic conditions and simultaneously avoiding artificial wave reflections. We next plan to analyze mathematically the stability and convergence of the method as well as the PML accuracy for different types of absorbing functions.

Figure 1 :

 1 Figure 1: Test 1 : initial density (Left), A z (T, x, y) at times t = 0.456a.u. (Middle) and 1.824a.u. (Right).

Figure 2 :Figure 3 :

 23 Figure 2: Test 1 : density at times t = 0.456a.u. (first column), t = 1.14a.u. (second column), t = 1.368a.u. (third column) t = T = 1.824a.u. (fourth column): with periodic boundary conditions (without PML) (first line) The second (respectively third) line reports the results for the density with TSSP-PML, with a PML of Type-III (respectively Type-IV), setting σ 0 = 10 -2 , θ = π/4. The fourth (respectively fifth) line gives the same kind of results but for the Type-V (respectively Type-VI) PML, with σ 0 = 10 -2 , θ = π/4 and δ ν = L ν -L * ν .

Figure 4 :

 4 Figure 4: Test 2 : (Left) Initial density. (Middle) 3-nucleus potential. (Right) Density at time T = 3 × 10 -1 a.u. without external laser field.

Figure 5 :

 5 Figure 5: Test 2 : z-component of the external laser potential: A z (T, x, z) for T = 3 × 10 -1 a.u.

Figure 6 :

 6 Figure 6: Test 2 : From left to right: No PML, reference and Type V-PML solutions at time t = 3×10 -2 a.u..

Figure 7 :

 7 Figure 7: Test 2 : From left to right: No PML, reference and Type V-PML solutions at time t = 6×10 -2 a.u..

Figure 8 :

 8 Figure 8: Test 2 : From left to right: No PML, reference and Type V-PML solutions at time t = 3×10 -1 a.u..

Figure 9 :

 9 Figure 9: Test 2 : From left to right: No PML, reference and Type V-PML solutions at time t = 3.6 × 10 -1 a.u..

Figure 10 :

 10 Figure 10: Test 2 : (Top) Solution of reference at time T = 6 × 10 -2 a.u. in the domain [-18, 18] × [-18, 18]. (Bottom) 2 -norm in logscale of the 3 solutions (No PML, Reference, Type V-PML) as a function of time in D.

  (Right), the maximum of the density at final time T in logscale log sup (x,z)∈[-8,8] 2 |d(T, x, z)| as a function of the number of gridpoints. We observe an

Figure 11 :

 11 Figure 11: Test 3 : initial density in logscale (Top) and static potential (Bottom).

Figure 12 :

 12 Figure 12: Test 3 : logarithm of the density at time T = 0.912 (from top-left to bottom-right) for the PML with parameters σ 0 = 1, σ 0 = 10 -1 , σ 0 = 10 -2 , σ 0 = 10 -3 , σ 0 = 5, and finally with periodic boundary conditions (S = 1).

Figure 13 :

 13 Figure 13: Test 3 : logarithm of the density at time T = 0.912a.u. with null Dirichlet boundary condition.This test illustrates that artificial wave-reflecting occurs with real space methods (here a quantum lattice Boltzmann method[START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]) are implemented without PML or ABC.

Figure 14 :

 14 Figure 14: Test 3 : logarithm of the density at time T = 0.912a.u. with σ 0 = 10 -2 and from top-left to right-bottom: θ = π/16, θ = π/8, θ = π/4 and θ = π/2.

Figure 15 :

 15 Figure 15: Test 3 : logarithm of the density at time T = 0.912a.u. for σ 0 = 10 -2 θ = π/4, δ ν = L ν -L * ν , and from top-left to right-bottom: absorption functions σ of type I to IV.

Figure 16 :

 16 Figure 16: Test 3 : PML of type IV (Top) maximum of the density as a function of the time iteration with 5 different values of σ and for θ = π/4; (Bottom) maximum of the density as a function of the time iteration for θ = π/16, π/8, π/4, π/2 and σ 0 = 10 -2 .

Figure 17 :

 17 Figure 17: Test 3 : initial data (Left), periodic boundary conditions (Middle), convergence with the Type-III PML (Right), for σ 0 = 10 -2 , θ = π/4 and N x = N y = 32, 256, 320, 512, 640, 768.

in intense two-color laser fields. Phys. Rev. A, 54(4):3254-3260, Oct 1996.

Acknowledgments. X. ANTOINE was partially supported by the French National Research Agency project NABUCO, grant ANR-17-CE40-0025. E. LORIN received support from NSERC through the Discovery Grant program.

References.


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

Emmanuel Lorin
email: elorin@math.carleton.ca

Asymptotic convergence rates of SWR methods for Schrödinger equations with an arbitrary number of subdomains

Keywords: Schrödinger equation, Schwarz Waveform Relaxation, domain decomposition method, asymptotic convergence rate

may come

Asymptotic convergence rates of SWR methods for Schrödinger equations with an arbitrary number of subdomains

Introduction

We are interested in this paper in the analysis of the rate of convergence of some SWR methods by using an arbitrary number of subdomains. This study is an extension of existing results about the convergence of SWR algorithms on 2 subdomains [START_REF] Antoine | Asymptotic estimates of the convergence of classical schwarz waveform relaxation domain decomposition methods for two-dimensional stationary quantum waves[END_REF][START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF][START_REF] Antoine | On the rate of convergence of schwarz waveform relaxation methods for the time-dependent schroedinger equation[END_REF]. We show that the convergence rates established for 2 subdomains are actually still accurate estimates for an arbitrary number of sufficiently large subdomains and bounded potentials. In this paper, we will mainly focus on the computation of contraction factors from Lipschitz continuous mappings involved in the proof of convergence of SWR methods. As a consequence, we will not introduce technical details about the full proof of convergence. Instead, we refer to several papers where the reader could find all the details of these proofs depending on the equation under consideration. For linear advection and diffusion reaction equations, we refer to [START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF]. The analysis and derivation for the Schrödinger equation in the time-dependent case is presented in [START_REF] Antoine | On the rate of convergence of schwarz waveform relaxation methods for the time-dependent schroedinger equation[END_REF][START_REF] Antoine | Domain decomposition method and highorder absorbing boundary conditions for the numerical simulation of the time dependent Schrödinger equation with ionization and recombination by intense electric field[END_REF][START_REF] Besse | Schwarz waveform relaxation method for one dimensional Schrödinger equation with general potential[END_REF][START_REF] Dolean | An introduction to domain decomposition methods[END_REF], while the stationary equation is studied in [START_REF] Antoine | Asymptotic estimates of the convergence of classical schwarz waveform relaxation domain decomposition methods for two-dimensional stationary quantum waves[END_REF][START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF].

The SWR algorithms are well-established methods for the parallel solution of evolution partial differential equations by allowing computations of the underlying PDE on small subdomains with very good speed-up. The convergence of the overall SWR methods is strongly dependent on the choice of the transmission conditions between the subdomain interfaces. Typically, Dirichlet transmission conditions will usually provide very slow convergence, while transparent transmissions boundary conditions provide a very fast convergence. In this paper, we do not discuss the space discretization of the algorithm but only focus on the convergence of the continuous in-space algorithms. We refer to [1, [START_REF] Dolean | An introduction to domain decomposition methods[END_REF][START_REF] Gander | Optimized Schwarz waveform relaxation methods for advection reaction diffusion problems[END_REF][START_REF] Gander | Dirichlet-Neumann and Neumann-Neumann waveform relaxation algorithms for parabolic problems[END_REF][START_REF] Gander | Overlapping Schwarz for linear and nonlinear parabolic problems[END_REF][START_REF] Gander | Optimal Schwarz waveform relaxation methods for the one-dimensional wave equation[END_REF][START_REF] Gander | Optimized Schwarz methods[END_REF][START_REF] Gander | Optimal convergence for overlapping and nonoverlapping Schwarz waveform relaxation[END_REF][START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF][START_REF]A time-dependent Dirichlet-Neumann method for the heat equation[END_REF] for details about the SWR methods.

This paper is organized as follows. In section 2, we analyze the rate of convergence of the CSWR (Classical SWR) and OSWR (Optimized SWR) methods for the stationary Schrödinger equation solved by the imaginary-time method. The analysis is extended in Section 3, and some numerical experiments are proposed in Section 4.

Stationary states problems

We study the convergence of the SWR methods by using an arbitrary number m 2 of subdomains, for computing the point spectrum of the Schrödinger Hamiltonian using the imaginary-time method [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF] Antoine | GPELab, a Matlab toolbox to solve Gross-Pitaevskii equations I: Computation of stationary solutions[END_REF][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 | Modeling and computation of Bose-Einstein condensates: stationary states, nucleation, dynamics, stochasticity[END_REF][START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein condensation[END_REF][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 | 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]. More specifically, we intend to determine the ground state to the following Schrödinger Hamitonian: -+ V (x), where the potential V is supposed to be smooth and bounded with bounded derivatives. The imaginary-time method reads, also called Normalized Gradient Flow method (NGF), for

t 0 = 0 < t 1 < • • • < t n < t n+1 < • • •                        ∂ t φ(t, x) = φ(t, x) -V (x)φ(t, x)φ(t, x), x ∈ Ω, t n < t < t n+1 , φ(t, x) = 0, x ∈ ∂Ω, t n < t < t n+1 , φ(t n+1 , x) := φ(, t + n+1 , x) = φ(t - n+1 , x) ||φ(•, t - n+1 )|| L 2 (R) , φ(0, x) = ϕ 0 (x), x ∈ R, with ||ϕ 0 || 2 L 2 (R) = 1, (1)
where ϕ 0 is a given initial guess, usually built from an ansatz. The procedure is repeated until the convergence is reached, that is when the following stopping criterion is satisfied

||φ(t n+1 , •) -φ(t n , •)|| L ∞ (R) δ, (2)
for δ > 0 small enough. We propose the following decomposition (with possible overlap):

Ω = ∪ m i=1 Ω i , such that Ω i = (ξ - i , ξ + i ), for 2 i m -1, Ω 1 = (-a, ξ + 1 )
and Ω m = (ξ - m , +a) (see Fig. 1). Moreover, the overlapping size is:

ξ + i -ξ - i+1 = ε > 0, and |Ω i | = L + ε,
where L is assumed to be much larger than ε. We also have ξ ± i+1 -ξ ± i = L. We study the convergence rate of the CSWR and OSWR algorithms. The convergence rate is defined as the slope of the logarithm residual history with respect to the Schwarz iteration number, i.e. {(k, log(E (k) )) : k ∈ N}, where (resp. T (k cvg ) ) denotes the NGF converged solution (resp. time) in Ω i at Schwarz iteration k, and δ Sc is a small parameter. More specifically, we intend to prove that the convergence rate is independent of the number of subdomains as already numerically mentioned in [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF] and Section 4. Finally, φ

E (k) := m i=1 φ cvg,(k) i (ξ - i+1 ,ξ + i ) -φ cvg,(k) i+1 (ξ - i+1 ,ξ + i ) ∞,Γε L 2 (0,T (k cvg ) ) δ Sc , (3)
(k) i is the local solution in Ω i for any i ∈ {1, • • • , m}, at Schwarz iteration k ∈ N.

Potential-free equation

We first consider the potential-free Schrödinger equation in imaginary-time, with

P (∂ t , ∂ x ) = ∂ t -∂ 2
x . The NGF algorithm consists in solving for any n ∈ N, from t n to t n+1 -:

   ∂ t -∂ 2 x φ(t, x) = 0, (t, x) ∈ (t n , t - n+1 ) × Ω, φ(t, -a) = 0, t ∈ (t n , t - n+1 ), φ(t, a) = 0, t ∈ (t n , t - n+1 ).
Then, we normalize the solution

φ(t n+1 , •) = φ(t n+1 -, •)/ φ(t n+1 -) L 2 (Ω)
. The procedure is repeated until convergence following [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF]. Let us start by studying the rate of convergence of the CSWR, then considering the OSWR algorithm.

CSWR algorithm

The CSWR algorithm reads as follows: for k

1 and i ∈ {2, • • • , m -1} from time t n to t n+1 -, solve          P (∂ t , ∂ x )φ (k) i = 0, (t, x) ∈ (t n , t - n+1 ) × Ω i , φ (k) i (0, •) = ϕ 0 , x ∈ Ω i , φ (k) i (t, ξ + i ) = φ (k-1) i+1 (t, ξ + i ), t ∈ (t n , t - n+1 ), φ (k) i (t, ξ - i ) = φ (k-1) i-1 (t, ξ - i ), t ∈ (t n , t - n+1 ).
In Ω 1 , we get

         P (∂ t , ∂ x )φ (k) 1 = 0, (t, x) ∈ (t n , t - n+1 ) × Ω 1 , φ (k) 1 (0, •) = ϕ 0 , x ∈ Ω 1 , φ (k) 1 (t, -a) = 0, t ∈ (t n , t - n+1 ), φ (k) 1 (t, ξ + 1 ) = φ (k-1) 2 (t, ξ + 2 ), t ∈ (t n , t - n+1 ),
and in Ω m

         P (∂ t , ∂ x )φ (k) m = 0, (t, x) ∈ (t n , t - n+1 ) × Ω m , φ (k) m (0, •) = ϕ 0 , x ∈ Ω m , φ (k) i (t, ξ - m ) = φ (k-1) m-1 (t, ξ - m ), t ∈ (t n , t - n+1 ), φ (k) m (t, +a) = 0, t ∈ (t n , t - n+1 ).
To analyze the convergence of this DDM, we set e i := φ exact|Ω i -φ i and we introduce h ± i ∈ H 3/4 0 (0, T ) = φ ∈ H 3/4 (0, T ) : φ(0) = 0 , for i ∈ {1, • • • , m}. We then consider the following system, where we denote by τ the dual variable to t in Fourier space, and by • the Fourier transform with respect to t,

iτ -∂ 2 x e i (τ, x) = 0, (τ, x) ∈ R × Ω i , e i (τ, ξ ± ) = h ± i (τ ), τ ∈ R .
We set α(τ ) := e iπ/4 √ τ , and we get, for any τ ∈ R,

e i (τ, x) = A i (τ )e α(τ )x + B i (τ )e -α(τ )x and for i ∈ {2, • • • , m -1} e i (τ, ξ ± i ) = A i (τ )e α(τ )ξ ± i + B i (τ )e α(τ )ξ ± i = h ± i (τ ) .
By using the boundary conditions, we find,

for i ∈ {2, • • • , m -1},        A i (τ ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) h + i (τ )e -α(τ )ξ - i -h - i (τ )e -α(τ )ξ + i , B i (τ ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) h - i (τ )e α(τ )ξ - i -h + i (τ )e α(τ )ξ + i .
Similarly, we have

       A 1 (τ ) = h + 1 (τ ) e α(τ )(L+ε/2-a) -e -α(τ )(L+ε/2+a) , B 1 (τ ) = h + 1 (τ ) e α(τ )(L+ε/2+a) -e -α(τ )(L+ε/2-a) ,
and

       A m (τ ) = h - m (τ ) e α(τ )(a-L-ε/2) -e α(τ )(L+ε/2+a) , B m (τ ) = - h - m (τ ) e -α(τ )(L+ε/2+a) -e α(τ )(L+ε/2-a) .
We introduce a mapping G (C) from H 3/4 (R)

2m to itself, defined as follows

G (C) : h + i (t), h - i+1 (t) 1 i m = e i (t, ξ - i+1 ), e i+1 (t, ξ + i ) 1 i m .
Thus, we deduce that

G (C) : h + i (τ ), h - i+1 (τ ) 1 i m = e i (τ, ξ - i+1 ), e i+1 (τ, ξ + i ) 1 i m , where, for i ∈ {2, • • • , m -1}, e i (τ, ξ - i+1 ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) h - i (τ )(e εα(τ ) -e -εα(τ ) ) + h + i (τ )(e Lα(τ ) -e -Lα(τ ) )
and

e i+1 (τ, ξ + i ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) h + i+1 (τ )(e εα(τ ) -e -εα(τ ) ) + h - i+1 (τ )(e Lα(τ ) -e -Lα(τ ) ) .
Following the same strategy as for the two subdomains case, see [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF][START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF], we compute

G (C),2 { h + i , h - i+1 } i . Let us set h +,(2) i (τ ) := e i (τ, ξ - i+1 ) and h -, (2)
i+1 (τ ) := e i+1 (τ,

ξ + i ), we get for i ∈ {2, • • • , m -1} e (2) i (τ, ξ - i+1 ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) e i (τ, ξ - i+1 ) e α(τ )L -e -α(τ )L + e i (τ, ξ + i-1 ) e α(τ )ε -e -α(τ )ε and e (2) i+1 (τ, ξ + i ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) e i+1 (τ, ξ + i ) e α(τ )L -e -α(τ )L + e i+1 (τ, ξ - i+2 ) e α(τ )ε -e -α(τ )ε .
We then obtain e

(2)

i+1 (τ, ξ + i ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) 2 h - i+1 (τ ) e α(τ )L -e -α(τ )L 2 + h - i+1 (τ ) e α(τ )ε -e -α(τ )ε 2 + h - i+1 (τ ) e α(τ )L -e -α(τ )L 2 + 2 h + i+1 (τ ) e α(τ )ε -e -α(τ )ε e α(τ )L -e -α(τ )L = e -2εα(τ ) 1 -e -2α(τ )(L+ε) 2 h - i+1 (τ ) 1 -e -2α(τ )L 2 + h - i+1 (τ ) e α(τ )(ε-2L) -e -α(τ )(ε-2L) 2 +2 h + i+1 (τ ) e α(τ )ε -e -α(τ )ε e -α(τ )L -e -3α(τ )L and e (2) i (τ, ξ - i+1 ) = 1 e α(τ )(L+ε) -e -α(τ )(L+ε) 2 h + i (τ ) e α(τ )L -e -α(τ )L 2 + h + i (τ ) e α(τ )ε -e -α(τ )ε 2 +2 h - i (τ ) e α(τ )ε -e -α(τ )ε e α(τ )L -e -α(τ )L = e -2α(τ )ε 1 -e -2α(τ )(L+ε) 2 h + i (τ ) 1 -e -2α(τ )L 2 + h + i (τ ) e α(τ )(ε-2L) -e -α(τ )(ε-2L) 2 +2 h - i (τ ) e α(τ )ε -e -α(τ )ε e -α(τ )L -e -3α(τ )L .
For m 2 subdomains, we have

G (C),2 { h + i (τ ), h - i+1 (τ )} i = e -2α(τ )ε h + i (τ ), h - i+1 (τ ) i + O e -α(τ )L .
In particular, for τ < 0, we obtain

G (C),2 { h + i (τ ), h - i+1 (τ )} i e -ε √ 2|τ | h + i (τ ), h - i+1 (τ ) i + Ce -L √ 2|τ | .
We then get

G (C),2 { h + i , h - i+1 } i H 3/4 (tn,t - n+1 ) 2m C { h + i , h - i+1 } i ,
with C a positive constant lower that 1. We also have

d 2 e i dx 2 (τ, x) = α 2 (τ ) e i (τ, x) . Then, one gets e i L 2 (0,T ;H 2 (Ω i )) R Ω i |α(τ )| 4 |α(τ )||e α(τ )(L+ε) -e -α(τ )(L+ε) | | h + i (τ )|(|e -α(τ )(ξ - i x) | + |e -α(τ )(x-ξ + i ) |)+ | h - i (τ )|(|e -α(τ )(ξ + i -x) | + |e -α(τ )(x-ξ - i ) |) dx Je ne sais pas trop ce que tu veux ecrire Π m i=1 H 2,1 (Ω i ×(0,T ))
This justifies the fact that, as numerically observed in [START_REF] Halpern | Optimized and quasi-optimal Schwarz waveform relaxation for the one-dimensional Schrödinger equation[END_REF] and in Section 4, the convergence rate is independent of the number of subdomains of length L. However, the overall convergence is linearly dependent of the number of subdomains through the summation over the m subdomains in (3). As a consequence, we expect that the logscale slope of the residual history, i.e. the convergence rate, to be independent of m. Finally, the overall error {(k, E (k) )} is still shifted in logscale, by a positive constant linearly dependent on log(m). We conclude by the following of the CSWR method with subdomains of length L is of the form

C (C) ε = sup τ ∈R e -ε √ 2|τ | + O e -L √ 2|τ | .

OSWR algorithm

We now study the OSWR algorithm for m 2 subdomains. To this end, we define the transparent boundary operator

∂ x ± ∂ 1/2
t . The OSWR algorithm reads as follows: for k 1 and for

i ∈ {2, • • • , m -1}, solve          P (∂ t , ∂ x )φ (k) i = 0, (t, x) ∈ (t n , t - n+1 ) × Ω i , φ (k) i (0, •) = ϕ 0 , x ∈ Ω i , (∂ x + ∂ 1/2 t )φ (k) i (t, ξ + i ) = (∂ x + ∂ 1/2 t )φ (k-1) i+1 (t, ξ + i ), t ∈ (t n , t - n+1 ), (∂ x -∂ 1/2 t )φ (k) i (t, ξ - i ) = (∂ x -∂ 1/2 t )φ (k-1) i-1 (t, ξ - i ), t ∈ (t n , t - n+1 ). In Ω 1 , we get          P (∂ t , ∂ x )φ (k) 1 = 0, (t, x) ∈ (t n , t - n+1 ) × Ω 1 , φ (k) 1 (0, •) = ϕ 0 , x ∈ Ω 1 , φ (k) 1 (t, -a) = 0, t ∈ (t n , t - n+1 ), (∂ x + ∂ 1/2 t )φ (k) 1 (t, ξ + 1 ) = (∂ x + ∂ 1/2 t )φ (k-1) 2 (t, ξ + 1 ), t ∈ (t n , t - n+1 ), and in Ω m          P (∂ t , ∂ x )φ (k) m = 0, (t, x) ∈ (t n , t - n+1 ) × Ω m , φ (k) m (0, •) = ϕ 0 , x ∈ Ω m , (∂ x -∂ 1/2 t )φ (k) i (t, ξ - m ) = (∂ x -∂ 1/2 t )φ (k-1) m-1 (t, ξ - m ), t ∈ (t n , t - n+1 ), φ (k) m (t, +a) = 0, t ∈ (t n , t - n+1 ). We then consider iτ -∂ 2 x e i (τ, x) = 0, (τ, x) ∈ R × Ω i , (∂ x ± α(τ )) e i (τ, ξ ± ) = h ± i (τ ), τ ∈ R .
By defining α(τ ) := e iπ/4 √ τ , we obtain

e i (τ, x) = A i (τ )e α(τ )x + B i (τ )e -α(τ )x .
By again using the boundary conditions, we find for

i ∈ {2, • • • , m -1} A i (τ ) = h + i (τ )e -α(τ )ξ + i 2α(τ ) , B i (τ ) = - h - i (τ )e α(τ )ξ - i 2α(τ ) .
We introduce a mapping G (O) defined as follows

G (O) : h + i (t), h - i+1 (t) 1 i m = (∂ x + ∂ 1/2 t )e i+1 (t, ξ + i ), (∂ x -∂ 1/2 t )e i (t, ξ - i+1 ) 1 i m .
Thus, we can write that

G (O) : h + i (τ ), h i+1 (τ ) 1 i m = (∂ x + α(τ )) e i+1 (τ, ξ + i ), (∂ x -α(τ )) e i (τ, ξ - i+1 ) 1 i m = 2α(τ ) A i+1 (τ )e α(τ )ξ + i , -B i (τ )e -α(τ )ξ - i+1 1 i m . 7 where for i ∈ {2, • • • , m -1}. Let us introduce h +,(2) i (τ ) = 2α(τ )A i+1 (τ )e α(τ )ξ + i , h -,(2) i+1 (τ ) = -2α(τ )B i (τ )e -α(τ )ξ - i+1 ,
so that one gets: e

(2)

i (τ, x) = A (2) i (τ )e α(τ )x + B (2)
i (τ )e -α(τ )x . From the boundary conditions, we find,

for i ∈ {2, • • • , m -1},          A (2) i (τ ) = h +,(2) i (τ )e -α(τ )ξ + i 2α(τ ) = A i+1 (τ ) , B (2)
i (τ ) = - h -,(2) i (τ )e α(τ )ξ - i 2α(τ ) = B i-1 (τ ) .
In addition, we have

G (O),2 : h + i (τ ), h - i+1 (τ ) 1 i m = (∂ x + α(τ )) e (2)
i+1 (τ, ξ + i ), (∂ x -α(τ )) e (2) i (τ, ξ - i+1 ) 1 i m = 2α(τ ) A (2) i+1 (τ )e α(τ )ξ + i , -B (2)
i (τ )e -α(τ )ξ - i+1 1 i m = h + i+2 (τ )e -α(τ )(ξ + i+2 -ξ + i ) , h - i-1 (τ )e -α(τ )(ξ - i+1 -ξ - i-1 ) 1 i m = h + i+2 (τ )e -2α(τ )(L+ε) , h - i-1 (τ )e -2α(τ )(L+ε) 1 i m
We can finally state the following result.

Proposition 2.2. The convergence rate C

(O) ε of the OSWR method with subdomains of length L is O e -2α(τ )(L+ε) . Remark 2.1. Using a simple unitary transformation, it is trivial to extend the above results to the case of a linear equation with time-dependent potential V (t) in L 1 loc (R). Indeed, from

i∂ t + ∂ 2 x + V (t) φ(t, x) = 0,
it is sufficient to define the new unknown φ(t, x) = e -t 0 V (s) φ(t, x) by gauge change, φ satisfying then the potential-free Schrödinger equation.

The space variable potential case

We now assume that the potential is space-dependent. Then, the argument used in Remark 2.1 is longer valid. As it was studied in [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], we can no more get a simple expression of the exact solution on each subdomain, and we then have to use approximations. Let us set :

P (t, x, D) = ∂ t -∂ 2 x -V (x)
, where V is smooth, bounded with bounded derivative.

CSWR algorithm

Using the same notations as above, we directly consider the error equation for the CSWR algorithm:

iτ -∂ 2 x + V (x) e i (τ, x) = 0, (τ, x) ∈ R × Ω i e i (τ, ξ ± ) = h ± i (τ ), τ ∈ R (4)
We also assume that V is positive, and that V ∞ exists. From Nirenberg's factorization theorem [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF][START_REF] Nirenberg | Lectures on linear partial differential equations[END_REF]

P (t, x, ∂ t , ∂ x ) = (∂ x + iΛ -)(∂ x + iΛ + ) + R, (5)
where R ∈ OPS -∞ is a smoothing pseudodifferential operator. The operators Λ ± are pseudodifferential operators of order 1/2 (in time) and order zero in x. Furthermore, their total symbols λ ± := σ(Λ ± ) can be expanded in S

1/2 S as λ ± ∼ +∞ j=0 λ ± 1/2-j/2 , (6)
where λ ± 1/2-j/2 are symbols corresponding to operators of order 1/2 -j/2, see [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF]. We denote by e i an approximate solution to [START_REF] Antoine | Robust and efficient preconditioned Krylov spectral solvers for computing the ground states of fast rotating and strongly interacting Bose-Einstein condensates[END_REF] in Ω i of the following form (neglecting the scattering effects)

e i (τ, x) = A i (τ ) exp -i x ξ + i λ + (y, τ )dy + B i (τ ) exp -i x ξ - i λ -(y, τ )dy . By construction, λ -= -λ + if we select λ + 1/2 = -λ - 1/2
(see [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF]). Now, by using the boundary conditions, we find that,

for i ∈ {2, • • • , m -1},                  A i (τ ) = h + i (τ ) -h - i (τ ) exp -i ξ + i ξ - i λ -(y, τ )dy 1 -exp -2i ξ + i ξ - i λ -(y, τ )dy , B i (τ ) = h - i (τ ) -h + i (τ ) exp -i ξ + i ξ - i λ -(y, τ )dy 1 -exp -2i ξ + i ξ - i λ -(y, τ )dy .
Let us introduce

G (C) : h + i (t), h i+1 (t) 1 i m = e i (t, ξ - i+1 ), e i+1 (t, ξ + i ) 1 i m ,
thus we have

G (C) : h + i (τ ), h i+1 (τ ) 1 i m = e i (τ, ξ - i+1 ), e i+1 (τ, ξ + i ) 1 i m .
Next, some computations yield

e i (τ, x) = h + i (τ ) -h - i (τ ) exp -i ξ + i ξ - i λ -(y, τ )dy exp i x ξ + i λ -(y, τ )dy 1 -exp -2i ξ + i ξ - i λ -(y, τ )dy + h - i (τ ) -h + i (τ ) exp -i ξ + i ξ - i λ -(y, τ )dy exp -i x ξ - i λ -(y, τ )dy 1 -exp -2i ξ + i ξ - i λ -(y, τ )dy and e i+1 (τ, ξ + i ) = h + i+1 (τ ) -h - i+1 (τ ) exp -i ξ + i+1 ξ - i+1 λ -(y, τ )dy exp -i ξ + i+1 ξ + i λ -(y, τ )dy 1 -exp -2i ξ + i+1 ξ - i+1 λ -(y, τ )dy + h - i+1 (τ ) -h + i+1 (τ ) exp -i ξ + i+1 ξ - i+1 λ -(y, τ )dy exp -i ξ + i ξ - i+1 λ -(y, τ )dy 1 -exp -2i ξ + i+1 ξ - i+1 λ -(y, τ )dy = h + i+1 (τ ) -h - i+1 (τ ) exp -i L+ε 0 λ -(y + ξ - i+1 , τ )dy exp -i L 0 λ -(y + ξ + i , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i+1 , τ )dy + h - i+1 (τ ) -h + i+1 (τ ) exp -i L+ε 0 λ -(y + ξ - i+1 , τ )dy exp -i ε 0 λ -(y + ξ - i+1 , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i+1 , τ )dy .
Then, one gets

e i (τ, ξ - i+1 ) = h + i (τ ) -h - i (τ ) exp -i ξ + i ξ - i λ -(y, τ )dy exp i ξ - i+1 ξ + i λ -(y, τ )dy 1 -exp -2i ξ + i ξ - i λ -(y, τ )dy + h - i (τ ) -h + i (τ ) exp -i ξ + i ξ - i λ -(y, τ )dy exp -i ξ - i+1 ξ - i λ -(y, τ )dy 1 -exp -2i ξ + i ξ - i λ -(y, τ )dy = h + i (τ ) -h - i (τ ) exp -i L+ε 0 λ -(y + ξ - i , τ )dy exp -i ε 0 λ -(y + ξ - i+1 , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i , τ )dy + h - i (τ ) -h + i (τ ) exp -i L+ε 0 λ -(y + ξ - i , τ )dy exp -i L 0 λ -(y + ξ - i , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i , τ )dy Let us set h +,(2) i (τ ) := e i (τ, ξ - i+1 ) and h -, (2)
i+1 (τ ) := e i+1 (τ, ξ + i ). Then we are led to e

(2)

i+1 (τ, ξ + i ) = e i+1 (τ, ξ - i+2 ) -e i+1 (τ, ξ + i ) exp -i L+ε 0 λ -(y + ξ - i+1 , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i+1 , τ )dy × exp -i L 0 λ -(y + ξ + i , τ )dy + e i+1 (τ, ξ + i ) -e i+1 (τ, ξ - i+2 ) exp -i L+ε 0 λ -(y + ξ - i+1 , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i+1 , τ )dy × exp -i ε 0 λ -(y + ξ - i+1 , τ )dy = h - i+1 (τ ) exp -2iελ -(τ, ξ - i+1 ) + R 1 τ, ε, L, max 1 i m h i and e (2) i (τ, ξ - i+1 ) = e i (τ, ξ - i+1 ) -e i (τ, ξ + i-1 ) exp -i L+ε 0 λ -(y + ξ - i , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i , τ )dy × exp -i ε 0 λ -(y + ξ - i+1 , τ )dy + e i (τ, ξ + i-1 ) -e i (τ, ξ - i+1 ) exp -i L+ε 0 λ -(y + ξ - i , τ )dy 1 -exp -2i L+ε 0 λ -(y + ξ - i , τ )dy × exp -i L 0 λ -(y + ξ - i , τ )dy = h + i (τ ) exp -2iελ -(τ, ξ - i+1 ) + R 2 τ, ε, L, max 1 i m h i ,
where R i = O max 1 i m h i e -α(τ )L . We then have

G (C),2 i { h + i (τ ), h - i+1 (τ )} i = C i (ε, τ, V ) h + i (τ ), h - i+1 (τ ) i + O max 1 i m h i e -α(τ )L .
We yet refer [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], where for large τ , we can rigorously expand the symbol λ ± . We do not proceed to this laborious expansions in this paper. In first approximation [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], the local convergence rate between Ω i , Ω i+1 is given by

C i (ε, τ ) ≈ exp -ε 2|τ | + V (ξ - i+1 )) .
In general, the above procedure does not allow to directly construct a contraction factor, that is, subdomain-independent coefficients C i . This is a consequence to the fact that the convergence rate is dependent on the values of the potential at the subdomain interfaces. However we argue that the contraction factor (related to the convergence rate of the CSWR method)

C (C) ε is such that sup τ exp -ε( 2|τ | + max i∈{1,••• ,m} V (ξ - i+1 )) C (C) ε sup τ exp -ε( 2|τ | + min i∈{1,••• ,m} V (ξ - i+1 )) sup τ exp -ε 2|τ | .
The following lemma actually justifies that the convergence rate obtained in the case of m subdomains is close to the one for two subdomains and numerically observed in [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF] and in Section 4.

Lemma 2.1. Let us assume that f :

x ∈ R N → f (x) ∈ R N , is such that f (x) = γx + ε, with γ < 1 and ε = o(γ). We define the sequence x n+1 = f (x n ) + ε, with x 0 ∈ R N given. Then, we have x n+1 γ n+1 x 0 + o(γ).
For ε > 0 and V ∞ small enough, we directly get

C (C) ε ≈ sup τ exp -2iεα(τ ) .
For the sake of simplicity, we do not detail more the computations. We conclude that, as in the potential-free case the contraction factor for m sufficiently large subdomains (L 1) is close to the one for the two-subdomains case.

Proposition 2.3. The convergence rate C (C) ε

of the CSWR method with m subdomains of length L is the same as for two-subdomains (see [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF]) up to a term of the order of O e -L √ 2|τ | .

Although, the convergence rate C (C) ε is determined, the overall convergence is also proportional to the number of subdomains due to the

H 3/4 (t n , t n+1 ) 2m -norm of {h + i , h - i+1 } 1,••• ,m .
In particular, the transmission from one subdomain to the next one naturally slows down the overall convergence of the SWR method, but without changing the slope of the residual history in logscale (see [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF]).

OSWR algorithm

We now consider the OSWR algorithm. More specifically, from time t n to t n+1 -, the OSWR algorithm with potential reads as follows:

         P (∂ t , ∂ x )φ (k) i = 0, (t, x) ∈ (t n , t - n+1 ) × Ω i , φ (k) i (0, •) = ϕ 0 , x ∈ Ω i , ∂ x + iΛ + (τ, x) φ (k) i (t, ξ + i ) = ∂ x + iΛ + (τ, x) φ (k-1) i+1 (t, ξ + i ), t ∈ (t n , t - n+1 ), ∂ x + iΛ -(τ, x) φ (k) i (t, ξ - i ) = ∂ x + iΛ -(τ, x) φ (k-1) i-1 (t, ξ - i ), t ∈ (t n , t - n+1 ). (7)
The operators Λ ± are coming from ( 5) and the corresponding symbols are denoted by λ ± . The latter can be constructed as an asymptotic series of the form [START_REF] Antoine | Asymptotic estimates of the convergence of classical schwarz waveform relaxation domain decomposition methods for two-dimensional stationary quantum waves[END_REF]. In practice, the series +∞ j=0 λ ± 1/2-j/2 is truncated, and for p ∈ N, we can define λ ± p := p j=0 λ ± 1/2-j/2 and the corresponding operators Λ ± p . In [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF], the SWR convergence rate is established for the transmission operator ∂ t ± iΛ ± p . In this paper, we will only consider Λ ± . The error equation in Fourier (resp. real) space in time (resp. space) is

iτ -V (x) -∂ 2 x e i (τ, x) = 0, (τ, x) ∈ R × Ω i , ∂ x + iΛ ± (τ, x) e i (τ, ξ ± ) = h ± i (τ ), τ ∈ R.
The analysis of convergence is identical to one presented in Subsubsection 2.2. Basically, we define

G (O) : h + i (t), h - i+1 (t) 1 i m = ∂ x + iΛ + (t, x) e i+1 (t, ξ + i ), ∂ x + iΛ -(t, x) e i (t, ξ - i+1 ) 1 i m
Thus, we obtain

G (O) : h + i (τ ), h i+1 (τ ) 1 i m = ∂ x + iλ + (τ, x) e i+1 (τ, ξ + i ), ∂ x + iλ -(τ, x) e i (τ, ξ - i+1 ) 1 i m .
We define approximate solutions to [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF] on each Ω i , neglecting again the scattering effets

e i (τ, x) = A i (τ ) exp -i x 0 λ + (y, τ )dy + B i (τ ) exp i x 0 λ -(y, τ )dy .
Then by construction λ + = -λ -, we get

A i (τ ) = h + i (τ )e ξ + i 0 λ + (τ,y)dy 2iλ + (τ, x) , B i (τ ) = - h - i (τ )e ξ - i 0 λ -(τ,y)dy 2iλ + (τ, x) ,
leading to

G (O),2 : h + i (τ ), h - i+1 (τ ) 1 i m = (∂ x + iλ + (τ, x)) e (2)
i+1 (τ, ξ + i ), (∂ x -iλ + (τ, x)) e (2) i (τ, ξ - i+1 ) 1 i m = h + i+2 (τ )e -i ξ + i+2 ξ + i λ + (τ,y) , h - i-1 (τ )e -i ξ - i+1 ξ - i-1 λ + (τ,y)dy 1 i m = h + i+2 (τ )e -2i L+ε 0 λ + (τ,y+ξ + i ) , h - i-1 (τ )e -2i L+ε 0 λ + (τ,y+ξ - i-1 )dy 1 i m .
We conclude that, as for the CSWR method, the OSWR method has a convergence rate independent, up to a multiplicative constant, of the number of subdomains of length L. Details of the convergence over 2 subdomains for the OSWR and quasi-OSWR can be found in Section 2.3 of [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF].

Scalability

We notice that using a large number of subdomains does not allow for an acceleration of the convergence rate, which is the same as for the 2 subdomain case up to a multiplication coefficient. We however benefit from i) an embarrassingly parallel algorithm, ii) small scale local computations on each subdomain.

Extension to time-dependent problems

The principle for analyzing the rate of convergence in the time-dependent equation is roughly speaking identical to the stationary case (see also [START_REF] Antoine | On the rate of convergence of schwarz waveform relaxation methods for the time-dependent schroedinger equation[END_REF]). Basically, we have to replace t (resp. τ ) by it (resp. iτ ) in the equations. We denote by P (∂ t , ∂ x ) = i∂ t + ∂ 2

x the Schrödinger operator, and consider the IBVP on Ω = (-a, a)

       i∂ t + ∂ 2 x φ(t, x) = 0, (t, x) ∈ (0, T ) × Ω, φ(0, •) = ϕ 0 , x ∈ Ω, φ(t, -a) = 0, t ∈ (0, T ), φ(t, a) = 0, t ∈ (0, T ),
where T > 0 and ϕ 0 ∈ L 2 (R) are given. The convergence rate is defined as the slope of the logarithm residual history according to the Schwarz iteration number, that is {(k, log(E (k) )) : k ∈ N}, where this time

E (k) := m i=1 φ (k) i (ξ - i+1 ,ξ + i ) -φ (k) i+1 (ξ - i+1 ,ξ + i ) ∞,Γε L 2 (0,T ) δ Sc , (8)
δ Sc being a small parameter.

First, the CSWR algorithm reads as follows, for k 1 and for i ∈ {2, • • • , m -1}:

         P (∂ t , ∂ x )φ (k) i = 0, (t, x) ∈ (0, T ) × Ω i , φ (k) i (0, •) = ϕ 0 , x ∈ Ω i , φ (k) i (t, ξ + i ) = φ (k-1) i+1 (t, ξ + i ), t ∈ (0, T ), φ (k) i (t, ξ - i ) = φ (k-1) i-1 (t, ξ - i ), t ∈ (0, T ).
In Ω 1 , we get

         P (∂ t , ∂ x )φ (k) 1 = 0, (t, x) ∈ (0, T ) × Ω 1 , φ (k) 1 (0, •) = ϕ 0 , x ∈ Ω 1 , φ (k) 1 (t, -a) = 0, t ∈ (0, T ), φ (k) 1 (t, ξ + 1 ) = φ (k-1) 2 (t, ξ + 2 ), t ∈ (0, T ),
and in Ω m

         P (∂ t , ∂ x )φ (k) m = 0, (t, x) ∈ (0, T ) × Ω m , φ (k) m (0, •) = ϕ 0 , x ∈ Ω m , φ (k) i (t, ξ - m ) = φ (k-1) m-1 (t, ξ - m ), t ∈ (0, T ), φ (k) m (t, +a) = 0, t ∈ (0, T ).
In order to analyze the convergence of the SWR methods, it is sufficient to repace τ by -iτ in all the equations from Section 2. In particular, for m subdomains, we can write that

G (C),2 { h + i (τ ), h - i+1 (τ )} i = e -2β(τ )ε h + i (τ ), h - i+1 (τ ) i + O e -α(τ )L ,
where β(τ ) := √ τ .

A similar study is possible with OSWR algorithm for m subdomains. First, define the transparent operator ∂ x ± e -iπ/4 ∂ 1/2 t at the interfaces. The OSWR algorithm reads as follows: for k 1 and for

i ∈ {2, • • • , m -1}          P (∂ t , ∂ x )φ (k) i = 0, (t, x) ∈ (t n , t - n+1 ) × Ω i , φ (k) i (0, •) = ϕ 0 , x ∈ Ω i , (∂ x + e -iπ/4 ∂ 1/2 t )φ (k) i (t, ξ + i ) = (∂ x + e -iπ/4 ∂ 1/2 t )φ (k-1) i+1 (t, ξ + i ), t ∈ (t n , t - n+1 ), (∂ x -e -iπ/4 ∂ 1/2 t )φ (k) i (t, ξ - i ) = (∂ x -e -iπ/4 ∂ 1/2 t )φ (k-1) i-1 (t, ξ - i ), t ∈ (t n , t - n+1 ).
In Ω 1 , we get

         P (∂ t , ∂ x )φ (k) 1 = 0, (t, x) ∈ (t n , t - n+1 ) × Ω 1 , φ (k) 1 (0, •) = ϕ 0 , x ∈ Ω 1 , φ (k) 1 (t, -a) = 0, t ∈ (t n , t - n+1 ), (∂ x + e -iπ/4 ∂ 1/2 t )φ (k) 1 (t, ξ + 1 ) = (∂ x + e -iπ/4 ∂ 1/2 t )φ (k-1) 2 (t, ξ + 1 ), t ∈ (t n , t - n+1 ), and in Ω m          P (∂ t , ∂ x )φ (k) m = 0, (t, x) ∈ (t n , t - n+1 ) × Ω m , φ (k) m (0, •) = ϕ 0 , x ∈ Ω m , (∂ x -e -iπ/4 ∂ 1/2 t )φ (k) i (t, ξ - m ) = (∂ x -e -iπ/4 ∂ 1/2 t )φ (k-1) m-1 (t, ξ - m ), t ∈ (t n , t - n+1 ), φ (k) m (t, +a) = 0, t ∈ (t n , t - n+1 ). We then consider iτ -∂ 2 x e i (τ, x) = 0, (τ, x) ∈ R × Ω i , (∂ x ± β(τ )) e i (τ, ξ ± ) = h ± i (τ ), τ ∈ R,
where we have set β(τ ) := e iπ/4 √ τ , and we get: e i (τ, x) = A i (τ )e β(τ )x + B i (τ )e -β(τ )x . From the boundary conditions, we find, for i

∈ {2, • • • , m -1}, A i (τ ) = h + i (τ )e -β(τ )ξ + i 2β(τ ) , B i (τ ) = - h - i (τ )e β(τ )ξ - i 2β(τ ) .
We introduce a mapping G (O) defined as follows

G (O) : h + i (t), h - i+1 (t) 1 i m = (∂ x + e -iπ/4 ∂ 1/2 t )e i+1 (t, ξ + i ), (∂ x -e -iπ/4 ∂ 1/2 t )e i (t, ξ - i+1 ) 1 i m .
We again find the same convergence rate as in [START_REF] Antoine | On the rate of convergence of schwarz waveform relaxation methods for the time-dependent schroedinger equation[END_REF] for 2 subdomain, but the overall error is still shifted in logscale, by a positive constant linearly dependent on log(m).

Numerical experiments

Eigenvalue problem

We present a simple test case illustrating the theoretical results presented in Section 2. We take V (x) = 5x 2 /2, and Ω = (-2, 2) and we solve the Schrödinger equation in imaginary-time using a three-points finite difference scheme, with meshsize ∆x = 1/64, and for a time step ∆t = 0.2. Details about the solver can be found in [START_REF] Antoine | An analysis of Schwarz waveform relaxation domain decomposition methods for the imaginary-time linear Schrödinger and Gross-Pitaevskii equations[END_REF]. We apply the CSWR method on m = 2, 4, 8 subdomains and compare the convergence rate as a function of the Schwarz iterations. Initially, we take φ 0 (x) = φ 0 (x)/ φ 0 2 , with φ 0 (x) = e -4x 2 . The m -1 overlapping zones have a length equal to ε = ∆x, corresponding to 2 nodes. We observe that the asymptotic residual history is numerically independent of the number of subdomains, see Fig. 2.

We also report on Fig. 3 the converged solution (x, φ cvg,(k (cvg) ) (x) = φ g (x)), x ∈ (-2, 2) , the initial guess (x, φ 0 (x)), x ∈ (-2, 2) , as well as the NGF converged solution after 20 Schwarz iterations with 8 subdomains, {(x, φ cvg,(20) (x)), x ∈ (-2, 2)}. We however numerically observe that for 8 subdomains the asymptotic rates of convergence seems to be relatively different than 2 and 4. This can be explained by the fact that for a large number of subdomains, the size of each of these subdomains is relatively small, which induces an inaccuracy in the theoretical slope of convergence. In addition, the overall convergence rate depends on the value of the potential at the subdomain interfaces. In the next example, the overall domain is larger, leading to larger subdomain, then better accuracy of the convergence rate.

Time-dependent equation

In this example the potential is chosen as V (x) = -5 exp(-2x 2 ), with Ω = (-8, 8) and final time T = 10. We solve the Schrödinger equation in real-time with the three-points finite difference scheme, for ∆x = 1/32 and ∆t = 10 -2 . We refer to [START_REF] Antoine | On the rate of convergence of schwarz waveform relaxation methods for the time-dependent schroedinger equation[END_REF] for the details concerning the solver. The CSWR method is applied for m = 2, 4, 8 and 16 subdomains. We compare the convergence rate as a function of the Schwarz iterations. Initially, we take φ 0 (x) = φ 0 (x)/ φ 0 2 , with φ 0 (x) = π -1/4 e -5(x+2) 2 +ik 0 x , for k 0 = 5. The m -1 overlapping zones have a length equal to ε = ∆x, corresponding to 2 nodes. We observe on Fig. 4 that, as proven above, the asymptotic residual history is independent of the number of subdomains. We also report in Fig. 5 the amplitude of the converged solution (x, |φ N,(k (cvg) ) (x)|, x ∈ (-8, 8) , the initial guess (x, |φ 0 (x)|), x ∈ (-8, 8) , as well as the TDSE converged solution after 10 Schwarz iterations with 16 subdomains, (x, |φ N, (7) (x)|), x ∈ (-8, 8) , where N is such that T N = T = 10.

Conclusion

In this paper, we presented an asymptotic analysis of the convergence rate of the multidomains CSWR and OSWR methods for the 1D linear Schrödinger equation in imaginaryand real-time. Asymptotic estimates show that the already existing estimates stated for the two-domains configuration extend here to m domains. This is illustrated through some numerical examples.

[1] M. Al-Khaleel, A.E. Ruehli, and M.J. Gander. Optimized waveform relaxation methods for longitudinal partitioning of transmission lines. IEEE Transactions on Circuits and Systems, 56:1732-1743, 2009.

Figure 1 :

 1 Figure 1: Domain decomposition with possible overlapping.

Proposition 2 . 1 .

 21 The convergence rate C (C) ε

Figure 2 :

 2 Figure 2: Comparison of the asymptotic convergence rates for 2, 4 and 8 subdomains, including a zoom in the asymptotic regime.

Figure 3 :

 3 Figure 3: Initial guess, potential and converged ground state.

Figure 4 :

 4 Figure 4: Comparison of the asymptotic convergence rates for 2,, 8, 16 and 32 subdomains, including a zoom in the asymptotic regime.

Figure 5 :

 5 Figure 5: Amplitudes of the initial guess, converged solution and solution after 7 Schwarz iterations on 16 subdomains.