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

Emmanuel Lorin
email: elorin@math.carleton.ca

Q Tang

A Friendly Review of Absorbing Boundary Conditions and Perfectly Matched Layers for Classical and Relativistic Quantum Waves Equations

Keywords: classical quantum waves, relativistic quantum waves, Schrödinger equation, Klein-Gordon equation, Dirac equation, absorbing boundary conditions, perfectly matched layers, numerical schemes, simulation

est destinée au dépôt

A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations

Introduction

The aim of this paper is to discuss how to build and discretize Absorbing Boundary Conditions (ABCs) and Perfectly Matched Layers (PMLs) [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Givoli | High-order local non-reflecting boundary conditions: a review[END_REF][START_REF] Hagstrom | New results on absorbing layers and radiation boundary conditions[END_REF] to truncate an unbounded spatial computational domain, when solving partial differential equations arising in atomic, molecular and laser physics. More specifically, the equations that will be considered are the Schrödinger, Klein-Gordon and Dirac equations in their simplest form. The goal here is not to give all the technical mathematical details, but rather to give a comprehensive introduction to the topic, from both the modeling and computational points of view. For the equations above, many questions remain still open and constitute some advanced developments in mathematical analysis and computational physics.

To motivate the need to develop such theoretical analysis, we first consider an example to understand why such boundary conditions are required for a practical computation. Let us introduce the following general scalar Schrödinger system    i∂ t ψ + ∆ψ + V(x, t)ψ = 0 , (x, t) ∈ R d

x × [0; T ], lim ||x||→∞ ψ(x, t) = 0, t ∈ [0; T ], ψ(x, 0) = ψ 0 (x) ,

x ∈ R d x .

(1)

In the above system, ψ := ψ(x, t) is the complex-valued wave function to be computed, the spatial variable is x ∈ R d (d ≥ 1) and the time variable is t > 0, T being a maximal time of computation and i := √ -1. Various possible potentials V could be considered: R) the space of infinitely differentiable functions, is a real-valued potential (e.g. V (x, t) = V c (x) + x • E(t) for modeling the action of external electric field E in length gauge, with an usual regularized interaction potential V c (x) in order to strictly satisfy the mathematical constraints; see [START_REF] Javanainen | Numerical simulations of multiphoton ionization and above-threshold electron spectra[END_REF] for instance in a simple framework; the regularity assumption could in fact be strongly weakened in the following developments, in particular in the domain Ω), • V(x, t) = f (ψ)(x, t) is a nonlinear potential (e.g. a cubic nonlinearity: V(x, t) = β|u(x, t)| 2 ), • or a combination of both situations : V(x, t) = V (x, t) + f (ψ)(x, t).

• V(x, t) = V (x, t) ∈ C ∞ (R x × [0; T ],
We also assume that the initial wave function ψ 0 is compactly supported (that is ψ 0 is localized in space; some developments for non compactly supported initial data are available in [START_REF] Ehrhardt | Discrete transparent boundary conditions for Schrödinger-type equations for non-compactly supported initial data[END_REF][START_REF] Vaibhav | Artificial boundary conditions for certain evolution PDEs with cubic nonlinearity for non-compactly supported initial data[END_REF]).

For numerical considerations, the initial problem must be solved in a bounded domain Ω. If the physical phenomena is not confined within a finite domain, problems indeed arise. Let us for example consider the following one-dimensional case (d = 1): V(x) := V (x) = x is a linear potential and the initial wave function is a gaussian beam ψ 0 (x) = e -x 2 +10ix (almost compactly supported numerically). A simple way to solve Problem [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF] is to introduce a bounded spatial computational domain e.g. Ω :=]-5; 15[ and to set homogeneous Dirichlet boundary conditions : ψ |ΣT = 0, i.e. we impose 0 on Σ T , where the boundary Σ T is simply defined here by Σ T := ∂Ω × [0; T ], with Σ := ∂Ω = {-5; 15} and T = 2. We report on Figure 1(a) the amplitude of the reference numerical solution |ψ| as a function of time computed on a large domain (so that we do not see the perturbation of the boundary conditions) and restricted here to Ω T := Ω × [0; T ], and on Figure 1(b) the numerical solution to the bounded problem on Ω T with homogeneous Dirichlet boundary conditions at Σ T . Both approaches use a second-order finite difference scheme in space and a Crank-Nicolson scheme in time [START_REF] Strikwerda | Finite difference schemes and partial differential equations[END_REF]. As we can observe, since the wave field ψ strikes the right boundary x r := 15 at T ≈ 0.25, then some numerical reflections occur. This is due to the Dirichlet boundary condition which does not mimic the property that the wave is outgoing to the computational domain, and actually acts like a "wall". At T = 1.4, some new reflections appear at the left point x := -5 and so on, the computational domain acting as a quantum wave guide structure.

From this example, it is clear that a suitable exact or approximate boundary condition must be imposed at the endpoints of the computational domain to avoid artificial reflections. We introduce a computational domain Ω T :=]x , x r [×[0; T ] and a fictitious boundary Σ T := ∂Ω = {x l , x r } × [0; T ]. The potential is usually a non supported function inside Ω (that is not localized in space). The kind of boundary condition that we are looking for writes as a relation between the (Dirichlet) trace ψ |Γ and the normal derivative (Neumann) trace ∂ n ψ |Γ of the solution at the boundary Σ T : F (ψ, ∂ n ψ |Γ ) = 0. In general, this relation is built in the form

∂ n ψ + iΛ + ψ = 0, on Σ T .
The operator iΛ + is called the Dirichlet-to-Neumann (DtN) operator.

ABCs for the one-dimensional wave equation

To introduce step by step the ideas related to ABCs for molecular physics, let us start with the one-dimensional wave equation

∂ 2 x ψ - 1 c 2 ∂ 2 t ψ = -Qδ xQ , (2)
in R x ×]0; T ], where ψ = ψ(x, t) is the wave function, x is the spatial coordinate and t is the time variable. We consider a constant sound speed c for the homogeneous isotropic medium R x . The source is supposed to be a point source with amplitude Q, located at point x Q . The free-space solution (that is for Q = 0) of Eq. ( 2) can be written as the sum of a left and a right traveling waves (ψ and ψ r , respectively here), if we assume that the initial data ψ 0 and ψ 1 are compactly supported in a domain

Ω :=] -x ; x r [:=] -R; R[ (R > 0), with boundary Σ := {x ; x r }. More precisely, if we have      ∂ 2 x ψ - 1 c 2 ∂ 2 t ψ = 0, (x, t) ∈ Ω T := R x ×]0; T [, ψ(x, 0) = ψ 0 (x), x ∈ R x , ∂ t ψ(x, 0) = ψ 1 (x), x ∈ R x , (3)
under the assumption that supp(ψ 0,1 ) ⊂ Ω, then the solution ψ can be written as: ψ(x, t) = ψ (x + ct) + ψ r (x -ct). Some calculations show that we have

ψ(x, t) = 1 2 (ψ 0 (x + ct) + ψ 0 (x -ct)) + 1 2 x+ct x-ct ψ 1 (u)du, (4)
that is

       2ψ (x) = ψ 0 (x) + x -∞ ψ 1 (u)du, 2ψ r (x) = ψ 0 (x) - x -∞ ψ 1 (u)du. (5)
Let us remark that, for |x| sufficiently large, that is for all x such that |x| > L, the supports of the two wave functions ψ ,r are disjoints. This implies e.g. that for x = -L we have: ψ(x, t) = ψ (x + ct), and so

∂ x ψ = c -1 ∂ t ψ.
In a symmetrical way, we have for the right traveling wave: x = L, and ψ(x, t) = ψ r (x -ct), and so ∂ x ψ = -c -1 ∂ t ψ. By introducing the outwardly directed unit normal vector n = ±1 to ] -L; L[, we have the unification of the two boundary conditions

∂ n ψ + 1 c ∂ t ψ = 0, (6)
at Σ T . This means that, if one wants to solve numerically the initial value problem, which is a difficult task since it is set in an unbounded domain R x , we may rather consider to solve the following Initial Boundary Value Problem (IBVP): find the approximate field ψ a such that

           ∂ 2 x ψ a - 1 c 2 ∂ 2 t ψ a = 0, (x, t) ∈ Ω T , ψ a (x, 0) = ψ a 0 (x), x ∈ Ω, ∂ t ψ a (x, 0) = ψ a 1 (x), x ∈ Ω, ∂ n ψ a + 1 c ∂ t ψ a = 0, (x, t) ∈ Σ T , (7)
which is set in the bounded spatial domain Ω T . The main point here is that numerically computing ψ a in Ω requires a finite number of grid points unlike working in R x . Furthermore, in the special case of a one-dimensional problem with constant wave speed, we can prove that we have: ψ a = ψ |ΩT which means that the restriction of the solution ψ of the initial system exactly coincides with the solution ψ a to the bounded domain problem. Therefore, the two waves travel through the boundary Γ without being reflected back into the domain, as physically expected. For this reason, the boundary condition given by Eq. ( 6) is said to be a Transparent Boundary Condition (TBC). If one considers the wave operator, an interesting remark is that we get the following exact factorization of the wave operator

∂ 2 x - 1 c 2 ∂ 2 t = (∂ x + 1 c ∂ t )(∂ x - 1 c ∂ t ). ( 8
)
An alternative way for designing the TBCs above is based on the Laplace transform

L (ψ)(x, ω) = ψ(x, ω) = ∞ 0 ψ(x, t)e -ωt dt,
with the co-variable ω = σ + iτ , σ > 0 (let us remark that an approach using the Fourier transform is also possible, see Section 6). Under these notations, one gets the following properties

L (∂ t ψ)(x, ω) = ω ψ(x, ω) -ψ(x, 0), L (∂ 2 t ψ)(x, ω) = ω 2 ψ(x, ω) -ωψ(x, 0) -∂ t ψ(x, 0).
The initial full space problem can then be rewritten as a transmission IBVP between the bounded Neumann problem

         ( 1 c 2 ∂ 2 t -∂ 2 x )ψ int = 0, x ∈ Ω, t > 0, ∂ x ψ int = ∂ x ψ ext , x ∈ Σ, t > 0, ψ int (x, 0) = ψ 0 (x), x ∈ Ω, ∂ t ψ int (x, 0) = ψ 1 (x), x ∈ Ω, (9)
and the exterior Dirichlet problem

         ( 1 c 2 ∂ 2 t -∂ 2 x )ψ ext = 0, x ∈ R x /Ω, t > 0, ψ ext (x, t) = ψ ext (x, t), x = ±L, t > 0, ψ ext (x, 0) = 0, x ∈ R x /Ω, ∂ t ψ ext (x, 0) = 0, x ∈ R x /Ω. (10)
Let us Laplace transform the first equation of system [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF]. Then, for x > x r , one obtains the following ODE

ω 2 c 2 ψext -∂ 2 x ψext = 0. (11)
A simple computation shows that the solution writes down as the superposition of two waves travelling in opposite directions

ψext (x, ω) = A + (ω)e ω c x + A -(ω)e -ω c x .
Since the exterior solution must have a finite energy (square-integrable function), this leads to the condition A + = 0 and to the solution (with x r = L)

ψext (x, ω) = e -ω c (x-xr) (ψ ext (x r , •))(ω).
Then by differentiating, by continuity and using the second equation in ( 9)

∂ x ψint (x, ω)| x=xr = - ω c ψint (x, ω)| x=xr ,
and finally by applying the inverse Laplace transform, one obtains

∂ x ψ int (x, t)| x=L = - 1 c ∂ t ψ int (x, t)| x=xr , leading to the TBC for ψ int at Σ T ∂ n ψ int + 1 c ∂ t ψ int = 0, on Σ T . (12)
Therefore, the solution ψ int exactly corresponds to ψ a in system [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]. This also translates the property that the factorization (8) indeed splits the original wave operator into an outgoing and an incoming absorbing operators according to the choice of the sign ± in front of the c -1 ∂ t operator.

3. ABCs for one-dimensional Schrödinger equations 3.1. The free-space case

Let us now come to the case of the Schrödinger equation. We first consider the simple case of the one-dimensional Schrödinger equation in the free-space since it is a constant coefficient partial differential equation. The techniques applied to the wave operator to derive boundary conditions can still be used, but the resulting boundary conditions have some very specific differences that require more developments, in particular from the computational point of view.

Let us first introduce the formulation of the problem as coupled transmission problems

   (i∂ t + ∂ 2 x )ψ int = 0, x ∈ Ω, t > 0, ∂ x ψ int = ∂ x ψ ext , x ∈ Σ, t > 0, ψ int (x, 0) = ψ 0 (x), x ∈ Ω, (13)
and

         (i∂ t + ∂ 2 x )ψ ext = 0, x ∈ Ω ,r , t > 0, ψ ext = ψ int , x ∈ Σ, t > 0, lim |x|→∞ ψ ext (x, t) = 0, t > 0, ψ ext (x, 0) = 0, x ∈ Ω ,r . (14)
Now, we Laplace transform the first equation of system [START_REF]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] (for example in the right semi-infinite subdomain Ω r :=]x r ; +∞[, setting Ω :=] -∞; x [). Then one gets the ODE: iω ψext + ∂ 2

x ψext = 0, whose solution is explicitely given by ψext (x, ω) = A + (ω)e

+ √ -iω x + A -(ω)e -+ √ -iω x ,
where the principal determination of the square-root is such that ( + √ •) > 0. Since we need to find a finite mass solution, one requires that ψ ext ∈ L 2 (Ω r ), leading to A + = 0. We can write that

ψext (x, ω) = e -+ √ -iω (x-xr) L( ψext (x r , •))(ω),
and by differentiating and continuity

∂ x ψint (x, ω)| x=xr = -+ √ -iω ψint (x, ω)| x=xr = -e -iπ/4 ω ψint (x, ω)| x=xr √ ω .
Now, the direct application of the inverse Laplace transform leads to

∂ x ψ int (x, t)| x=xr = -e -iπ/4 ∂ t ( 1 √ π t 0 ψ int (x, s)| x=xr √ t -s ds) = -e -iπ/4 ∂ 1/2 t ψ int (x r , t),
where ∂

1/2 t
designates the Caputo fractional derivative [START_REF] Kilbas | Theory and applications of fractional differential equations[END_REF] of order 1/2. Symmetrically, the same derivation occurs at the left boundary. A unified writing of the transparent boundary condition is then

(∂ n + e -iπ/4 ∂ 1/2 t )ψ int = 0, on Σ T . (15)
Compared to the TBC [START_REF] Antoine | Towards accurate artificial boundary conditions for nonlinear PDEs through examples[END_REF] for the wave equation, there is a fundamental difference when considering [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF] for the Schrödinger equation. Indeed, the BC involves a nonlocal time operator which must be discretized carefully (see Section 3.3). In practice, we are led to solving the following IBVP set in a finite spatial domain

   i∂ t ψ int + ∂ 2 x ψ int = 0 , (x, t) ∈ Ω T , ∂ n ψ int + e -iπ/4 ∂ 1/2 t ψ int = 0, (x, t) ∈ Σ T , ψ int (x, 0) = ψ 0 (x), x ∈ Ω. ( 16
)
In addition, if ψ 0 ∈ H 1 (Ω) (finite energy Sobolev space), then it can be proved [START_REF]Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF] that there exists a unique smooth solution ψ in a well-adapted mathematical setting. Moreover, ψ satisfies the following mass inequality:

ψ(t) L 2 (Ω) ≤ ψ 0 L 2 (Ω) , ∀t > 0.
This property can be physically interpreted as the usual mass conservation property for the full space equation: ψ(t) L 2 = ψ 0 L 2 , ∀t > 0. Indeed, for the bounded domain problem [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF], some mass of the solution goes out of the computational domain, which is exactly what is expected from the TBC. Let us also remark that, similarly to the wave operator, one gets the factorization:

i∂ t + ∂ 2 x = (∂ x + √ -i∂ t )(∂ x - √ -i∂ t ).
It can be shown that the TBC in [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF] writes: (∂ n + √ -i∂ t )ψ int = 0, for (x, t) ∈ Σ T .

Adding a potential or/and a nonlinearity

Let us now consider the case where a potential is included into the equation. Then, if the potential V(x) = V (x) is general in Ω ,r in system (1) but is time independent (in Ω ,r , V can be general in the computational domain Ω, even nonsmooth), we are led to solving a variable coefficients ODE corresponding to: ∂ 2

x ψext +(V (x)+iω) ψext = 0, after applying a Laplace transform. Most of the time, getting the TBC is out of reach [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF]. If V also depends on t, this would lead to handle the Laplace transform of V (x, t)ψ ext which is a convolution.

In some special cases, the TBC can still be computed. For example, for a linear potential V (x) = x, the expression of the exact operator can be made explicit [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] through the Airy functions in the Laplace domain. However, the concrete use of such boundary conditions is nontrivial. If one considers a potential V := V (t), then a change of gauge φ(x, t) = ψ ext (x, t)e -iV(t) leads (after coming back to the interior solution ψ int ) to the TBC

∂ n ψ int + e -iπ/4 e iV(t) ∂ 1/2 t (e -iV(t) ψ int ) = 0, on Σ T , (17)
with: V(t) := I t V , where I t is the integral of V from 0 to t.

Considering the general situation is much more complicate and TBCs can generally not be built. Instead, one can produce families of approximate Absorbing Boundary Conditions (ABCs) of increasing order that minimize the reflection at the boundary. Various directions exist [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF][START_REF] Ehrhardt | Exact artificial boundary conditions for problems with periodic structures[END_REF], but some, based on the generalization of the point of view of the Laplace transform by pseudodifferential operators theory and symbolical calculus [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Hörmander | Linear partial differential operators[END_REF], are systematic and can be also used for many systems of PDEs (in particular, for the Klein-Gordon and Dirac equations). For example, one may derive the following fourth-order ABC (see [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF])

∂ n ψ int + e -iπ/4 e iV ∂ 1/2 t e -iV ψ int + i sg(∂ n V ) |∂ n V | 2 e iV I t ( |∂ n V | 2 e -iV ψ int ) = 0, (18)
on Σ T , where sg is the sign function and V(x, t) := I t V . We directly see that, when V is x-independent, this ABC simplifies to [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF] and is then exact as a TBC. Other kinds of ABCs can be found in the literature [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF][START_REF] Ehrhardt | Exact artificial boundary conditions for problems with periodic structures[END_REF] for variable potentials.

Constructing TBCs for a nonlinear equation is generally impossible (see [START_REF] Zheng | Exact nonreflecting boundary conditions for one-dimensional cubic nonlinear Schrödinger equations[END_REF] for a special case). To get a nonlinear ABC for system (1) in 1D, the main idea is to use the formal substitution V = f (ψ)(x, t) and V = I t (f (ψ))(x, t) in the BC ( 18) [START_REF] Antoine | Artificial boundary conditions for onedimensional cubic nonlinear Schrödinger equations[END_REF][START_REF]Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF][START_REF] Antoine | Towards accurate artificial boundary conditions for nonlinear PDEs through examples[END_REF][START_REF] Szeftel | Absorbing boundary conditions for one-dimensional nonlinear Schrödinger equations[END_REF]. Another approach, called unified approach and prospected in [START_REF] Zhang | Unified approach to split absorbing boundary conditions for nonlinear Schrödinger equations[END_REF][START_REF]Unified approach to split absorbing boundary conditions for nonlinear Schrödinger equations: Two-dimensional case[END_REF], also uses a splitting operator idea but rather to build ABCs for nonlinear Schrödinger equations. Volkov-based absorbing boundary conditions following similar ideas, and specifically designed for laser-molecule interaction, were also proposed in [START_REF] Lorin | Mathematical modeling of boundary conditions for laser-molecule time-dependent schrdinger equations and some aspects of their numerical computation -one-dimensional case[END_REF].

How to discretize the initial boundary-value problem

When one wants to discretize one of the above IBVP, we must careful approximate the boundary conditions to get a stable scheme [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF]Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF]. Indeed, if one considers for example [START_REF] Bandrauk | Atoms and molecules in intense laser fields: Gauge invariance of theory and models[END_REF] which includes all the difficulties, we remark that the boundary condition involves

• some time convolution integrals (fractional derivatives/integration operators),

• some variable coefficients and nonlinear terms, • and the normal derivative.

When discretizing the boundary-value problem, even with a simple Dirichlet (ψ int = 0) or Neumann (∂ n ψ int = 0) boundary condition on Σ T , building a semi-discrete time scheme is far from being trivial. Indeed, in such a case, we often require that some physical properties remain fulfilled at the semi-discrete level: mass/energy conservation, time reversibility, gauge invariance, preservation of the dispersion relation. We refer to [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] for a complete description of the most standard and well-adapted schemes (Crank-Nicolson, Time-Splitting, relaxation schemes...). When considering an absorbing boundary condition, then one must take care to preserve these properties. This is not a trivial task even for the simplest case, i.e. the one-dimensional free-space linear Schrödinger equation [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF]Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF]. To illustrate this difficulty, let us mention that when using a Crank-Nicolson scheme [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF][START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF]Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF] for a uniform time step ∆t > 0, the half-order time derivative operator must be discretized thanks to the following quadrature rule

∂ 1/2 t f (t n ) ≈ 2 ∆t n k=0 β n-k f k , (19)
with

β k = (-1) k α k , ∀k ≥ 0, (α 0 , α 1 , α 2 , ...) = (1, 1, 1 2 , 1 2 , 3 8 , 3 8
, ...), to provide an un-conditionally stable scheme. The variable coefficients and nonlinear terms, as well as the normal derivative, are usually discretized consistently with the interior time and space scheme (e.g. finite difference, finite element). Let us remark that well-known FFT-based pseudo-spectral approximation schemes in space [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 con-densation[END_REF] which strongly assume that the boundary conditions are periodic (or that we have a Dirichlet or Neumann boundary condition when using sine-or cosine-transform) cannot be used with an ABC. Nevertheless, we will see in Section 5 that a suitable truncation of the computational domain can be achieved using a Perfectly Matched Layer (PML) approach. Finally, let us remark that fast evaluation schemes can be employed for the computation of convolution quadratures [START_REF] Jiang | Fast evaluation of nonreflecting boundary conditions for the Schrödinger equation in one dimension[END_REF][START_REF]Fast evaluation of exact transparent boundary condition for one-dimensional cubic nonlinear Schrödinger equation[END_REF].

To illustrate the accuracy of the boundary condition [START_REF] Bandrauk | Atoms and molecules in intense laser fields: Gauge invariance of theory and models[END_REF] for the truncated IBVP (1), we consider the nonlinear potential V(x, t) = 0.1x 2 -|ψ| 2 , and the initial data is taken as

ψ 0 (x) = 2sech( √ 2x)e i 15 2
x . We first compute the reference solution ψ ref with a relaxation time scheme [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF][START_REF]Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF][START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] on a very large domain, so that we do not see the effect of the boundary and plot the amplitude |ψ ref (x, t)| in the domain Ω T :=]-10; 10[×[0; 2] (see Figure 2(a)), for ∆t = 10 -3 . Then, we discretize [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] the truncated IBVP with the same semi-discrete time scheme (relaxation), setting successively as a function of time, the discrete version of the ABC (18) at three different locations: x r = 10 (Figure 2(b)), x r = 8 (Figure 2(c)) and x r = 6 (Figure 2(d)). The spatial discretization scheme is based on a linear finite element method [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] with N = 4 × 10 3 points. As we can see, the ABC (18) yields an accurate simulation of the outgoing wave to the computational domain at x r .

ABCs for two-dimensional Schrödinger equations

4.1. The free-space case: straight boundary The extension of ABCs to higher-dimensions requires mainly the suitable integration of the effect of the boundary Σ. In the case of a straight boundary for the two-dimensional case, the exact DtN operator Λ + can be built [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF] in the form

∂ n ψ + iΛ + (∂ y , ∂ t )ψ = 0, on Σ T := Σ × [0; T ], ( 20
)
where n is the outwardly directed unit normal vector to the left half-space Ω := {x := (x, y) ∈ R 2 /x < 0}. The expression of Λ + is based on the inverse time Laplace transform, the partial Fourier transform in y and an extension of the factorization formula (see e.g. [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schrödinger equations[END_REF]). The operator Λ + is nonlocal both in time and space and can be tricky to implement in a standard code (e.g. based on finite-difference or finite-element approximations). To simplify the situation, local approximations are usually used. For example, the following first-and second-order ABCs can be derived

ABC 1 (∂ n + e -iπ/4 ∂ 1/2 t )ψ = 0, on Σ T , ABC 2 (∂ n + e -iπ/4 ∂ 1/2 t -e iπ/4 1 2 ∆ Σ I 1/2 t )ψ = 0, on Σ T , (21)
where ∆ Σ := ∂ 2 y is the second-order derivative operator (called Laplace-Beltrami operator) over Σ. These boundary conditions are local in space, which is expected for deriving efficient algorithms (indeed, the associated discrete matrices are then highly sparse) but nevertheless remain nonlocal in time (yielding then to consider all the past history of the solution). Similarly to the one-dimensional case, the fractional operators ∂ can be evaluated by suitable discrete formulae, through fast convolution algorithms or localized representations based on rational approximations [START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF][START_REF] Jiang | Fast evaluation of nonreflecting boundary conditions for the Schrödinger equation in one dimension[END_REF]. From the implementation point of view, the approach is relatively close to the onedimensional case. The operator ∆ Σ could be approximated through finite-difference or finite-element methods, depending on the spatial approximation scheme in Ω.

The free-space case: curved boundary

Since the boundary conditions ABC 1,2 defined by the expressions [START_REF] Bermudez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF] are derived in the half-space case, they can be directly applied to a rectangular domain. However, some errors may arise during the numerical simulations due to the corner reflection. Usually, additional well-adapted corner conditions must be enforced to accurately represent the wave function ψ, leading to difficult developments in mathematical and numerical analysis [START_REF] Hagstrom | A new auxiliary variable formulation of high-order local radiation boundary conditions: corner compatibility conditions and extensions to firstorder systems[END_REF]. For practical computations, ABCs are rather preferred for smooth boundaries (e.g. a square with rounded corners, or a circle). Therefore, the ABC must take into account the fact that the fictitious boundary is curved, incorporating in particular the effects of the curvature κ(s) (as well as its variations ∂ s κ, s being the curvilinear abscissa along Σ) and the curvilinear variations of the solution onto the surface Σ, i.e. through some terms ∂ s ψ, ∂ 2 s ψ,... We do not give any detail here since the derivation of families of local ABCs is extremely technical and we refer to [START_REF] Antoine | Construction, structure and asymptotic approximations of a microdifferential transparent boundary condition for the linear Schrödinger equation[END_REF] for more details. For example, the following fourth-order ABC can be obtained

∂ n ψ + e -iπ/4 ∂ 1/2 t ψ + κ 2 ψ -e iπ/4 κ 2 8 + 1 2 ∆ Σ I 1/2 t ψ +i κ 3 8 + 1 2 ∂ s (κ∂ s ) + ∆ Σ κ 8 I t ψ = 0, on Σ T , (22)
where ∆ Σ := ∂ 2 s . Let us remark that [START_REF] Besse | A relaxation scheme for the nonlinear Schrödinger equation[END_REF] simplifies to ( 21) when κ = 0, i.e. for a flat interface (and then s = y). Furthermore, when Σ is a circle, many terms vanish. The implementation of such boundary conditions can be achieved again by using suitable discrete quadratures for the time operators [START_REF] Antoine | Numerical schemes for the simulation of the twodimensional Schrödinger equation using non-reflecting boundary conditions[END_REF] and some finite-element or finite-difference schemes in space.

Adding a potential or/and a nonlinearity

Finally, like in the one-dimensional situation, variable potentials and nonlinear terms could be incorporated into the ABC. The derivation needs a lot of technical developments [START_REF]Absorbing boundary conditions for general nonlinear Schrödinger equations[END_REF][START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF]. Here, to illustrate the purpose, we write a fourth-order ABC that must be seen as the extension of ( 18) and ( 22) to the 2D potential case

∂ n ψ + e -iπ/4 e iV ∂ 1/2 t e -iV ψ + κ 2 ψ -e iπ/4 e iV ( ∆ Σ 2 + κ 2 8 + i∂ s V∂ s + 1 2 (i∂ 2 s V -(∂ s V) 2 ))I 1/2 t (e -iV ψ) +ie iV ( ∂ s (κ∂ s ) 2 + κ 3 8 + ∆ Σ κ 8 + i∂ s κ∂ s V 2 )I t (e -iV ψ) -i sg(∂ n V ) |∂ n V |e iV I t ( |∂ n V |e -iV ψ) = 0, on Σ T , (23)
with V = f (ψ)(x, t) and V = I t (f (ψ))(x, t). The discretization can be achieved through all the previous developments. The reader is refer to [START_REF]Absorbing boundary conditions for general nonlinear Schrödinger equations[END_REF][START_REF]Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential[END_REF] for more details and numerical examples showing that these ABCs perform very well.

Even if ABCs can be produced for complicate equations, they are nevertheless not always so easy to implement. In addition, for Schrödinger-type equations, one of the drawback is that they are usually nonlocal in time. Localization can be done but at the price of some advanced algorithmic developments. In the following section, we introduce the Perfectly Matched Layer (PML) approach which is easier to implement and appears as a simple computational alternative to ABCs.

PMLs for Schrödinger equations

The idea of Perfectly Matched Layer, PML, (sometimes also called sponge layers) has been originally introduced for electromagnetism by a French engineer J.-P. Bérenger, in the 94' paper [START_REF] Bérenger | A perfectly matched layer for the absorption of electromagnetic waves[END_REF]. Since then, many developments have been directed towards the understanding and improvement of PMLs for many PDE systems arising in Physics and Engineering [START_REF] Hagstrom | New results on absorbing layers and radiation boundary conditions[END_REF]. This methodology has some closed connections to what is commonly called "absorbers" in the atomic physics literature and which refers to empirically designed absorbing layers [START_REF] Chelkowski | Electron-nuclear dynamics of multiphoton h2 dissociative ionization in intense laser fields[END_REF]. Notice however that PML are rigorously designed and provide better absorption as analyzed in [START_REF]Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF]. In addition, they do not need to be adapted to the potential and nonlinearity, unlike complex Absorbing Potential (AP) approaches [START_REF] Muga | Complex absorbing potentials[END_REF][START_REF] Scrinzi | Infinite-range exterior complex scaling as a perfect absorber in time-dependent problems[END_REF] (see also Section 7.3 for a comparison between the PMLs and AP approaches for the Dirac equation). For the Schrödinger equation, the developments are more recent (see e.g. [START_REF]Numerical solution of time-dependent nonlinear Schrödinger equations using domain truncation techniques coupled with relaxation scheme[END_REF][START_REF]A perfectly matched layer approach to the nonlinear Schrödinger wave equations[END_REF] for some contributions). Basically, PMLs are designed for linear wave problems [START_REF] Hagstrom | New results on absorbing layers and radiation boundary conditions[END_REF] (see also Section 6.3 for the closely related Klein-Gordon equation). Their application to nonlinear PDEs is direct but remains formal, without any theoretical background as for the ABCs. Nevertheless, as seen below, they perform very well for such situations and are easy to implement which makes them attractive. In particular, and unlike ABCs, they can also be easily incorporated into an existing solver based on pseudo-spectral FFT-based methods [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF], opening the road to efficient time-splitting schemes.

To give a simplified but sufficiently general presentation, we consider system (1) for the 2D nonlinear case with a cubic nonlinearity and a general potential V . Let us assume that we are interested in the computation of an approximate solution inside the rectangular domain D Phy =]-L x ; L x [×]-L y ; L y [. The idea of PMLs is to surround the physical domain D Phy by an absorbing unphysical layer:

D PML := (] -L * x ; L * x [×] - L * y ; L * y [) \ D Phy , with L * ν = L ν + δ ν (ν =
x, y) (represented by the green region in Figure 3). The objective is to artificially damp the waves entering into the layer D PML . To this end, one introduces the following functions

S ν (ν) = 1, |ν| < L ν , 1 + e iϑν σ ν (|ν| -L * ν ), L ν ≤ ν < L * ν , ν = x, y, (24)
where ϑ ν is a constant and σ ν is the so-called absorbing function. Various choices of functions ϑ ν and σ ν are reported in the literature (see e.g. [START_REF] Bermudez | An optimal perfectly matched layer with unbounded absorbing function for time-harmonic acoustic scattering problems[END_REF][START_REF]A perfectly matched layer approach to the nonlinear Schrödinger wave equations[END_REF]). Among these, the following choices are popular and work well for nonlinear Schrödinger equations:

ϑ ν = π 4 and σ ν = σ ν 0 (ν + δ ν ) 2
, where δ ν = L * ν -L ν is the thickness in the ν-direction of the PML layer and σ ν 0 is a real-valued positive constant representing the absorbing strength. Since the wave is damped into the layer, we generally are free to choose the boundary condition on the outer boundary. Usually, homogeneous Dirichlet or Neumann boundary conditions are set. Nevertheless, since our goal in the presentation is here to use a pseudo-spectral FFT-based approximation scheme, we impose periodic boundary conditions. With such choices, the 2D cubic nonlinear Schrödinger equation with a potential V is truncated in the computational domain D

:=] -L x ; L x [×] - L y ; L y [= D Phy ∪ D PML and approximated by i∂ t ψ = - 1 2 ∂ x S x (x) ∂ x S x (x) + ∂ y S y (y) ∂ y S y (y) ψ + V + β|ψ| 2 ψ, x ∈ D, (25)
+periodic BC for x ∈ ∂D, t ≥ 0.

Since system ( 25)-( 26) is a nonlinear Schrödinger equation with variable coefficients and periodic boundary conditions, it can be solved efficiently by using FFTs incorporated into a time-splitting scheme [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF]. To this end, we first apply the time splitting approach to deal with the nonlinearity, i.e., for x ∈ D, one solves

i ∂ t ψ(x, t) = - 1 2 ∂ x S x (x) ∂ x S x (x) + ∂ y S y (y) ∂ y S y (y) ψ(x, t), t n ≤ t ≤ t n+1 , ( 27
)
for a time step ∆t > 0, followed by solving

i ∂ t ψ(x, t) = V (x, t) + β |ψ(x, t)| 2 ψ(x, t), t n ≤ t ≤ t n+1 , (28)
for ∆t. If t ∈ [t n , t n+1 ], Eq. ( 28) leaves |ψ| invariant in t, i.e., |ψ(x, t)| = |ψ(x, t n )| := |ψ n (x)|. Higher order splitting can naturally be implemented. Therefore, (28) reduces to a linear ODE, which can be integrated analytically as:

ψ(x, t) = e -i[V(x,t n ,t)+β |ψ n (x)| 2 (t-t n )] ψ(x, t n ), x ∈ D, t n ≤ t ≤ t n+1 , (29)
with V(x, t 1 , t 2 ) = t2 t1 V (x, τ )dτ . To solve Eq. ( 27) with periodic boundary conditions by a Fourier pseudo-spectral method, we adapt the approach developed e.g. in [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]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] which combines FFT-based evaluations of the operator appearing in the right-hand side of Eq. ( 27) after the semi-discretization in time by a Crank-Nicolson scheme and the solution to implicit linear systems through a preconditioned Krylov subspace solver [START_REF] Saad | GMRES -A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems[END_REF]. This approach is simple to implement since it basically uses algorithmic bricks that are well managed (even from the parallel implementation view point). The resulting second-order in time and spectrally accurate scheme in space appears therefore to be very efficient and robust, and well-adapted to simulate problems set in infinite domains. We refer to [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF] for more details and variant around this method.

To illustrate the method, we apply the resulting scheme to study the dynamics of a manufactured 2D soliton. To this end, we choose the potential and initial data respectively as V (x, t) = -1 2 sech 2 (x -t)sech 2 (y -t)(cosh(2(x -t)) + cosh(2(y -t))), ψ 0 (x) = sech(x)sech(y)e i(x+y) and let β = -1. With these parameters, we can solve exactly the problem and get the analytical outgoing solution ψ(x, t) = sech(x-t)sech(y-t)e i(x+y) . The other parameters are chosen as: L ν = 8, δ ν = 1 and σ µ = 80 (ν = x, y). The time step and mesh size are fixed as ∆t = 10 -2 and h ν = 1/8, respectively (tuning the PML parameters can sometimes be tricky). Figure 4 shows the contour plots of |ψ| 2 at different times. We can see that the soliton enters into the PML region, is very well damped and then no visible waves are reflected back into the physical domain. This shows that PMLs can be extremely useful when one wants to truncate an infinite domain, even when the problem is nonlinear. Let us remark that it is easy to extend the method to higher-dimensions, and that circular PMLs are also available (see e.g. [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF] for the Schrödinger equation). The Time-Domain Klein-Gordon Equation (TDKGE) is a wave equation for which absorbing boundary layers have been extensively studied. We assume that the initial data ψ 0 and ψ 1 are compactly supported in a domain Ω :=]-x ; x r [:=]-R; R[ (R > 0), with boundary Σ := {x ; x r }. A strategy similar to the one used for the Schrödinger equation is then applied to a spinless relativistic particle of charge e and mass m and subject to an external constant electric field, written here in Coulomb gauge and denoted by (V, A) (scalar and vector potential). In this case, the TDKGE reads

   P KG ψ = 0, (x, t) ∈ Ω T := R x ×]0; T [, ψ(x, 0) = ψ 0 (x), x ∈ R x ∂ t ψ(x, 0) = ψ 1 (x), x ∈ R x , ( 30
)
where the operator P KG is given by

P KG = ∂ 2 x -1/c 2 ∂ 2 t -ieV /c 2 ∂ t -ieA∂ x + e 2 /c 2 V 2 -e 2 A 2 -m 2 c 2 .
Without giving the mathematical details, one gets the following TBC by factorizing

P KG ∂ n - 1 c ∂ 2 t + ieV ∂ t -e 2 V 2 + 3 4 c 2 e 2 A 2 + m 2 c 4 - iceA 2 ψ = 0, on Σ T . (31)
Now, when the external field is not constant and m = 0, it is no more possible to explicitly construct the exact transparent operator. However, the computation of approximate boundary conditions is possible. This step requires the introduction of mathematical tools which are beyond the scope of this paper. In addition, the introduction of a high-frequency regime assumption or the use of Padé's approximants is often required. For the interested reader, we refer to [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF] where it can be shown that a sequence of ABCs of increasing order can be built. For not constant A x and V , the resulting ABCs of orders k = 1, 2, 3, are given on Σ T by

ABC 1 ∂ n + 1 c ∂ t ψ = 0, ABC 2 ∂ n + 1 c ∂ t + ieV 2c - ieA x 2 ψ = 0, ABC 3 ∂ n + 1 c ∂ t + ieV 2c - ieA x 2 - 1 2c -ie∂ t V -ic 2 e∂ x A x + e 2 V 2 -c 2 e 2 A 2 x -m 2 c 4 - 1 4 (e 2 V 2 -c 2 e 2 A 2 x ) - i 2 (e∂ x V -ce∂ x A x ) + i 2 (e∂ t V -ce∂ t A x ) I t ψ = 0.
(32)

ABCs for the two-dimensional Klein-Gordon equation

We consider the two-dimensional TDKGE ( = 1) including the electromagnetic field, that is A(x, y, t) and V (x, y, t) which respectively represent the vectorial and scalar potentials satisfying the Maxwell's equations. The equation writes

i∂ t -eV (x, y, t) 2 -c 2 -i∇ -eA(x, y, t) 2 -m 2 c 4 ψ = 0 .
To simplify the presentation, we consider that the computational domain Ω is a disk of radius R. We rewrite the TDKGE in polar coordinates in Σ r , for r ∈ [0, ) with > 0 small enough, that is (the details are skipped)

c 2 ∂ 2 r -∂ 2 t + c 2 χ 2 r ∂ 2 θ -µV ∂ t + O 1 + O 2 + O 3 ψ = 0. ( 33
)
We set above

A x = A r cos(θ) -A θ sin(θ), A y = A r sin(θ) + A θ cos(θ) and        O 1 := - c 2 χ r + c 2 µA r ∂ r , O 2 := - c 2 µ χ r A θ ∂ θ , O 3 := -µ∂ t V -c 2 µ∂ r A r - c 2 µ χ r A r - c 2 µ χ r ∂ θ A θ + e 2 V 2 -c 2 e 2 A 2 -m 2 c 4 .
Again skipping all the technical details, ABCs for the field-particle TDKGE can be obtained on Σ R T := C(0, R) × [0; T ], with C(0, R) the circle of radius R and center 0,

as ABC 1 ∂ r + 1 c ∂ t ψ = 0, ABC 2 ∂ r + 1 c ∂ t + ie 2 V + 1 2c O 1 ψ = 0, ABC 3 ∂ r + 1 c ∂ t + ie 2 V + 1 2c O 1 + - c 2 2χ 2 r ∂ 2 θ - 1 2 O 2 ∂ θ + iO 4 I t ψ = 0,
where O 4 is a complicated algebraic scalar operator, which is defined page 282 in [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. Notice that the first-order ABC (ABC 1 ) is exactly ( 6) in the r-direction.

We now consider the discretization of the 1D TDKGE with ABCs. A second-order explicit scheme is proposed to approximate the equation defined by the corresponding operator P KG in the truncated spatial domain Ω, that is for interior nodes

ψ n+1 j -2ψ n j + ψ n-1 j = ∆t 2 c 2 ∆x 2 ψ n j+1 -2ψ n j + ψ n j-1 + A n j 2∆t ψ n+1 j -ψ n-1 j + B n j 2∆x ψ n j+1 -ψ n j-1 + C n j ψ n j ,
where

   A n j = -ieV (x j , t n ), B n j = -ieA x (x j , t n ), C n j = -ie∂ t V (x j , t n ) -ic 2 ∂ x A x (x j , t n ) + e 2 V 2 (x j , t n ) -c 2 e 2 A 2 x (x j , t n ) -m 2 c 4 .
The initial data is a wave packet given by

ψ 0 (x) = ck 0 mc 2 + m 2 c 4 + c 2 k 2 0 e -x 2 4 +ik0x ,
with k 0 = 10. The external field is such that A(t) = cos(t), V (t) = 0.1 sin(t) and = m = c = 1. The computational domain is Ω :=] -10; 10[, T = 8 and ∆x = ∆t = 0.4. The results for ABC 1 , ABC 2 and ABC 3 (see Eq. ( 32)) are presented in Figure 5 in the (x, t)-plane (log-scale). The discretization of the time integral is based on similar ideas [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF] as [START_REF] Bao | Mathematical theory and numerical methods for Bose-Einstein con-densation[END_REF]. We clearly see that increasing the order of the ABC reduces the artificial reflection at the boundary.

PMLs for the Klein-Gordon equation

We start this section by introducing the notion of PML for the linear wave equation. The beginning of this section is extracted from [START_REF] Johnson | Notes on Perfectly Matched layers (PMLs)[END_REF]. As for the Schrödinger equation, a thin layer surrounds the physical domain D Phy designed as an absorbing and reflectionless material. The construction of PMLs is based on an analytic continuation of the wave function in the complex plane, i.e. the wave function can be evaluated for complex values of x in the absorbing layer. More precisely, we set x → x + if (x) for some differentiable real-valued function f such that f (x) = σ x (x)/ω, where σ x (x) > 0 and a given frequency ω, in the PML region and zero otherwise. Then the derivative with respect to x becomes:

∂ x → 1/(1 + iσ x /ω)∂ x .
The prefactor allows for an absorption of the wave function of frequency ω in the PML region. In the simplest 1D case, it then reads

∂ t ψ = c∂ x φ -σ x ψ, ∂ t φ = c∂ x ψ -σ x φ,
where φ is an intermediate function.

In higher dimensions, PMLs for wave equations are established first thanks to the Laplace transform in time and Fourier transform in the direction of propagation. Then a modal solution/ansatz decreasing exponentially in the direction of propagation is postulated. A new set of equations adapted to the chosen ansatz is then derived in the external layer, in Laplace and Fourier variables. Finally going back to real time and space, we get the full set of equations to solve in the PML [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF]. Notice that auxiliary functions may be necessary to localize in time the PML. This concept applied in 2D [START_REF] Johnson | Notes on Perfectly Matched layers (PMLs)[END_REF] to the wave equation:

∂ 2 t ψ -c 2 (∂ 2 x + ∂ 2 y
)ψ = 0 reads as follows. For some function σ x and transformation x → 1/(1 + iσ x /ω)∂ x , the PML version is

∂ t ψ -c∇φ -σ x ψ + ψ = 0, with ∂ t φ x = c∂ x ψ -σ x φ x , ∂ t φ y = c∂ y ψ and ∂ t ψ = cσ x ∂ y φ y .
Now, for more general wave-like equations, including the Klein-Gordon equation, of the form

∂ t ψ -∇ 2 ψ + k≥0 c k u 2k+1 + d k ∂ t ψ 2k+1 = 0, with c k ≥ 0 , d k ≥ 0,
a set of PMLs was established in [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF]. This nice paper is a friendly theoretical as well as numerical presentation of PMLs for nonlinear wave equations.

We now test the PML approach for the equation : ∂ 2 t ψ -∂ 2 x ψ + mc 2 ψ = 0, with mc 2 = 1. In this case, we implement the PML in [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF]

   ∂ 2 t ψ -∂ 2 x ψ + mc 2 ψ = -σ(x)∂ t ψ + σ(x)η (2) + ασ(x) -mc 2 ψ + ∂ x σ(x)η (1) , ∂ t η (1) + σ(x) + α η (1) + ∂ x ψ = 0, ∂ t η (2) + αη (2) + (mc 2 + α 2 )ψ = 0,
where σ is chosen as in [START_REF] Appelö | Application of a perfectly matched layer to the nonlinear wave equation[END_REF], i.e. σ(x) = 0 in the interior zone ]0; L x [ and σ(x x . We approximate this system by using an elementary finite-difference scheme. We take α = 1, L x = 10, δ x = L x /10 and σ x = 50. The initial data is taken as ψ 0 (x) = e -(x-Lx/2) 2 and the final time is T = 6. For comparison, we compute the reference solution ψ ref on a large domain. The space and time steps are chosen respectively as ∆x = 4 × 10 -2 and ∆t = 2 × 10 -2 . We first compare in real space in Figure 6

) = σ x (1 -δ -1 x (x -L * x ) 2

ABCs/PMLs for the Dirac equation

In this section, we propose an overview of ABCs and PMLs for the Dirac equation. In 2D, the Dirac equation which models the interaction of a quantum particle with an external field (A, V ) reads P D ψ = 0, with

P D := I 4 ∂ t + α x c∂ x -ieA x + α y c∂ y -ieA y + i(V c + eV )I 4 -iβmc 2 ,
where [START_REF] Bandrauk | Atoms and molecules in intense laser fields: Gauge invariance of theory and models[END_REF][START_REF] Lorin | A simple and accurate mixed p0-q1 solver for the Maxwell-Dirac equations[END_REF] • V c : R 2 x,y → R is an interaction (Coulomb) potential. • V : R 2

x,y × R + → R is a combination of the self-consistent and external electric potentials.

• A : R 2

x,y × R + → R 2 is a combination of the electromagnetic potential generated by the particle charge and by an external potential.

• m is the mass of the particle, e their charge, c the speed of light.

• ψ : R + × R 2

x,y → C 4 is the Dirac wave function. • the Hermitian Dirac matrices α x , α y , β are defined by

                       α x = 0 2 σ x σ x 0 2 , σ x = 0 1 1 0 , α y = 0 2 σ y σ y 0 2 , σ y = 0 -i i 0 , β = β 2 0 2 0 2 -β 2 , β 2 = 1 0 0 -1 ,
and 0 2 and I 2 are respectively the zero and identity matrices in M 2 (C). The following relations hold: α 2 x = α 2 y = β 2 = I 4 and {α x , α y } = {α x , β} = {α y , β} = 0 4 .

• J = (J x , J y ), with J x = ec ψ, α x ψ C 4 , J y = ec ψ, α y ψ C 4 , denotes the current density and ρ stands for the particle density which is equal to: e 4 i=1 |ψ i | 2 . In the following, we assume that (A, V ) is given at any time t. The system under consideration then reads in 2D P D ψ = 0, (x, y, t) ∈ Ω T = R 2

x,y ×]0; T [, ψ(x, y, 0) = ψ 0 (x, y), (x, y) ∈ R 2

x,y .

(34)

ABCs for the one-dimensional Dirac equation

In the one-dimensional case, the Dirac equation is commonly called in mathematics a first-order hyperbolic system. Such equations are well-known in fluid dynamics, mechanical engineering, as well as electromagnetism. This property is a consequence of the fact that the matrices α x,y,z (and β) are diagonalizable in R. There exists an extended literature (see for instance [START_REF] Ghidaglia | The normal flux method at the boundary for multidimensional finite volume approximations in CFD[END_REF]) on the derivation and approximation of boundary conditions modeling outgoing waves. Basically, the derivation of ABCs for the Dirac equation follows some very similar ideas. However, since the matrix β is not diagonalizable in the same basis as α x,y,z , even in 1D, it is only possible to derive ABCs and not TBCs. The strategy to construct ABCs for the Dirac equation is in fact mathematically equivalent to the one developed for the Klein-Gordon or Schrödinger equations. This is again nontrivial and it requires some material going beyond the scope of the present paper. We refer to [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF] for the interested and mathematically inclined reader. Let us recall that, in 1D, the Dirac equation can be written as a 2 × 2 system

σ x ∂ x ψ -imcβ 2 ψ - 1 c ∂ t ψ = 0, (x, t) ∈ Ω T := R x ×]0; T [, ψ(x, 0) = ψ 0 (x), x ∈ R x . (35)
To explain the principle to build ABCs for the Dirac equation, we first consider the massless particle case m = 0 and no external field. Then, the Dirac equation reads

∂ t ψ + α x ∂ x ψ = 0, with ψ = (ψ 1 , ψ 2 ), i.e. ∂ t ψ 2 + c∂ x ψ 1 = 0, ∂ t ψ 1 + c∂ x ψ 2 = 0.
If we introduce the two fields φ 1 = (ψ 1 + ψ 2 )/ √ 2 and φ 2 = (ψ 1 -ψ 2 )/ √ 2, we get a system of uncoupled transport/wave equations (first-order wave equations)

∂ t φ 1 + c∂ x φ 1 = 0, ∂ t φ 2 -c∂ x φ 2 = 0. ( 36
)
Equivalently, setting

α x = Π x D x Π -1
x , then the system (36) reads ∂ t φ + D x ∂ x φ = 0, where φ = (φ 1 , φ 2 ) and

Π x = 1 √ 2 1 -1 1 1 , D x = 1 0 0 -1 . At x = x r (resp.
x ), we have the following TBC

   ∂ t φ int + cD x ∂ x φ int = 0, x ∈ Ω, t > 0, φ int 2 = φ ext 2 , x = x r , t > 0, φ int (x, 0) = φ 0 (x), x ∈ Ω,
where φ 0 = Π -1 x ψ 0 . The exact solution is then reconstructed by setting ψ int = Π x φ int . Now, if m = 0, and more generally when the particle is subject to an external electric field, the derivation requires more efforts. However, from the above remark when m = 0, a first simple ABC reads at

x = x r    ∂ t φ int + cD x ∂ x φ int -imc 2 Π x β 2 φ int = 0, x ∈ Ω, t > 0, φ int 2 = φ ext 2 , x = x r , t > 0, φ int (x, 0) = φ 0 (x), x ∈ Ω.
or equivalently

   ∂ t ψ int + cα x ∂ x ψ int -imc 2 β 2 ψ int = 0, x ∈ Ω, t > 0, Π -1 x ψ int 2 = Π -1 x ψ ext 2 , x = x r , t > 0, ψ int (x, 0) = ψ 0 (x),
x ∈ Ω.

To improve the absorption at the boundary a finer approach is necessary, where this time Π x will be replaced by an complex operator Π. A short overview is presented in the two-dimensional case in the next subsection.

ABCs for the two-dimensional Dirac equation

The derivation of ABCs for the Dirac equation in 2D is mathematically non-trivial, so that we will just focus here on the general ideas. Their construction is based on a diagonalization [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF] of the Dirac operator up to a negative order operator (following a diagonalization theorem for pseudodifferential hyperbolic systems [START_REF] Antoine | Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism[END_REF] that generalizes the idea of factorization). The action of the remaining operator on the solution at the boundary will be expected negligible at least for high frequencies. To simplify the presentation, we assume that the domain is a disk of radius R, Ω R = C(0, R), and we define Σ R T := Ω R × [0; T ]. The Dirac equation is then formally rewritten in polar coordinates (r, θ)

I 4 ∂ r + L(r, θ, t, ∂ θ , ∂ t ) ψ = 0, (37)
where L = L 1 + L 0 is the operator defined through

L 1 := 1 c α x ∂ t + 1 χ r α x α y ∂ θ , L 0 := 1 c α x β, χ r := 1 R + r ,
with the matrices

α x =     0 0 0 e -iθ 0 0 e iθ 0 0 e -iθ 0 0 e iθ 0 0 0     , α y =     0 0 0 -ie -iθ 0 0 ie iθ 0 0 -ie -iθ 0 0 ie iθ 0 0 0     .
Finally, we set β = iβmc 2 + i(eV + V c )I 4 -ie(A r α x + A θ α y ). We suppose that A θ , A r , V and V c are independent of (t, θ). We denote by F (t,θ) the Fourier transform in (t, θ) with co-variables (τ, ξ). Then Fourier transforming (37) leads to

I 4 ∂ r + F (t,θ) (L)(r, θ, t, ξ, τ ) F t,θ (ψ) = 0. (38)
Again, if α x,y and β were diagonalizable in the same basis, it would be easier to impose TBCs. However, we can show that there exists a transition matrix F (t,θ) (Π) to F (t,θ) (L), such that one gets an almost exact operator diagonalization: (∂ r + Λ)φ ∼ 0, where

F (t,θ) (Λ) = F (t,θ) (Λ -) 0 0 σ(Λ + ) ∈ M 4 (C),
Λ + (resp. Λ -) representing the outgoing (resp. incoming) wave operator at the boundary. The ABCs, which simply consists of vanishing the incoming waves (inside Ω R ), are given in the following theorem.

Theorem 7.1. The zeroth-and first-order ABCs are

• zeroth-order: Π -1 0 ψ 1,2 = 0 on, Σ R T , • first-order: Π -1 1 ψ 1,2 = 0 on, Σ R T ,
where Π -1 0 and Π -1 1 are operators defined as

Π -1 0 = 1 2     -e iθ 0 0 v + 0 -e -iθ v -0 e iθ 0 0 v - 0 e iθ v + 0     , (39)
with v ± = 1 ± ic χ r ∂ θ ∂ -1 t , and Π -1 1 = 1 2 A 1 B 1 C 1 D 1 , with A 1 = - b -1 + e iθ 0 0 a -1 -e -iθ , B 1 = 0 b 1 - a 1 + 0 , C 1 = b -1 -e iθ 0 0 a -1 + e -iθ , D 1 = 0 b 1 + a 1 - 0 , setting a 1 ± = 1 ± (i(V c + eV ) ± ieA r )I t , a -1 ± = 1 ± (i(V c + eV ) ∓ ieA r )I t , b 1 ± = (1 ± (imc 2 ± eA θ )I t and b -1 ± = (1 ± (imc 2 ∓ eA θ )I t .

Numerical discretization and example of ABC for the Dirac equation

We propose a simple numerical illustration of the derived ABCs for the 1D Dirac equation. The quantity which is represented here is again (t, x, log |ψ 1 (x, t)|) which shows the reflections at the domain boundary. We first rewrite P D in the form

P D = i∂ t -iA∂ x + Bmc 2 ,
where

A = c 0 0 -c , B = i 0 mc 2 -mc 2 0 .
As explained in Subsection (7.1), a simple and natural condition to impose at x and x r is (∂ t ± c∂ x )ψ 1,2 = 0. This boundary condition is called a transport-like boundary condition, and corresponds in fact to the zeroth-order ABC in Theorem (7.1). The initial conditions are

ψ 1 (x, 0) = e -x 2 δ 2 +ik0x , ψ 2 (x, 0) = ck 0 mc 2 + m 2 c 4 + c 2 k 2 0 ψ 1 (x, 0),
where c = 1, k 0 = 20, δ = 0.5 and m = 1. The spatial computational domain is ]x ; x r [, with x r = -x = 5. The final time is T = 7 and ∆t = ∆x = 1/40 (we have 400 segments in space). The one-dimensional ABCs can be deduced from Theorem 7.1. We consider the following discretization scheme: denoting by (φ n j , ψ n j ) an approximation of the exact two-spinors (ψ 1 (x j , t n ), ψ 2 (x j , t n )) at (x j , t n ), for interior points, we have

     φ n+1 j = φ n j + c ∆t ∆x (φ n j+1 -φ n j ) + ∆tmc 2 ψ n j , ψ n+1 j = ψ n j -c ∆t ∆x (ψ n j -ψ n j-1 ) -∆tmc 2 φ n j .
At order 0, the transport condition reads at the discrete level φ n 0 -ψ n 0 = 0 at the left boundary, and φ n N + ψ n N = 0 at the right one. Let us remark that for the very particular interior scheme which is considered here, Dirichlet boundary conditions are equivalent to this order 0 condition (transport-like). At the next order, a timeintegral is added to the boundary condition. Stability is trivially satisfied under a CFL (Courant-Friedrichs-Lewy) condition [START_REF] Strikwerda | Finite difference schemes and partial differential equations[END_REF].

PMLs for the Dirac equation

In this section, we briefly describe the application of PMLs to the Dirac equation following [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] where a detailed presentation is available. As previously mentioned, the PML approach consists in constructing an absorbing layer surrounding the physical domain, say in 3D, D Phy :=] -L x ; L x [×] -L y ; L y [×] -L z ; L z [. Denoting by Σ(v) a continuous, positive and non-decreasing function (usually a quadratic function like σ ν 0 v 2 ), which is null for negative v, the following nonlinear 3D coordinates stretching is proposed in [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF]: η → η + iω -1 (µ - η (η) + µ + η (η)), with µ ± η (v) = v ±Lη Σ(-L η ± u)du, for η = x, y or z. By definition, the change of variable is only active in the thin layer D Phy surrounding the physical domain and is equivalent to the use of absorbing functions in the PML approach (the interpretation through coordinates stretching has been introduced in computational electromagnetism in [START_REF] Chew | A 3D Perfectly Matched Medium from Modified Maxwell's Equations with Stretched Coordinates[END_REF], and is shortly discussed in Section 6. In [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF], an accurate finite-difference discretization of these PMLs is proposed coupled with a time-splitting approach [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Hammer | Single-cone real-space finite difference scheme for the time-dependent Dirac equation[END_REF]. Homogeneous Dirichlet boundary conditions are set on the outer boundary. Let us remark that an interesting research direction would be to adapt the spectral approaches with splitting scheme introduced in Section 5 and [START_REF] Antoine | Perfectly matched layers and pseudospectral schemes for the dynamics of Gross-Pitaevskii equations with rotating term[END_REF] for the Schrödinger equation, to the Dirac equation, again imposing periodic boundary conditions on the exterior boundary of the computational domain D.

To conclude this section, we propose the comparison between the PMLs and complex AP approaches in 2D. The parameters for the Dirac equation are again: = c = m = 1. The approximation scheme given in [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] is a leap-frog scheme [START_REF] Hammer | Single-cone real-space finite difference scheme for the time-dependent Dirac equation[END_REF] which partially avoids fermion doubling and uses splitting as in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]. The domain is ] -5; 5[ 2 , T = 10, ∆x = ∆y = 4×10 -2 and ∆t = ∆x/2. The initial data is a gaussian: ψ 0 (x) := e -|x-x 0 | a 2 , with a = 0.6. The function Σ(v) is null for v ≤ 0 and equal to σ 0 v 2 otherwise. The absorption parameter σ 0 takes some values between 10 -1 and 10 1 . In the first plot of Figure 8, the abscissa represents the absorption coefficient σ = σ 0 and the ordinate corresponds to the L 2 (D Phy )-error between the exact and approximate solutions. We consider three different sizes of the layer: δ x (= δ y ) = 2 (diamond), δ x = 1.5 (circle) and δ x = 1 (circle). We clearly observe that the PML approach is extremely accurate while being simple to implement and yields more accurate results than the complex AP method. In addition, we also report in Figure 8, the evolution of |ψ| over the time with PMLs, where the chosen parameters are σ 0 = 1 and the layer size is δ x = 2. Again, we observe that the behavior of the field is well-reproduced.

Conclusion and perspectives

We presented a simplified and friendly overview of recent developments on ABCs and PMLs for PDEs arising in atomic, molecular and laser physics. The equations of interest are the Schrödinger, Klein-Gordon and Dirac equations. In addition to the main theoretical developments to understand these methods, we also provide some simulations to illustrate why they are useful and more accurate than standard empirical truncation techniques designed by physical arguments, e.g. the complex AP method. Even if many developments have now been achieved, some mathematical and numerical questions still need to be investigated. Finally, the probably most important step concerns now the use of these methods for simulations of high-dimensional realistic situations that arise in atomic and molecular physics. We expect that this paper can provide a guideline to achieve such a goal.

Figure 1 .

 1 Figure 1. Example of the one-dimensional Schrödinger equation with a linear potential (V (x) = x): truncating the physical domain with an unsuited boundary condition generates some spurious unphysical reflections at the fictitious boundaries x = -5 and xr = 15.

  (a) Numerical reference solution. (b) Numerical approximate solution (xr = 10). (c) Numerical approximate solution (xr = 8). (d) Numerical approximate solution (xr = 6).

Figure 2 .

 2 Figure 2. ABCs for the 1D Schrödinger equation (with a quadratic potential minus a nonlinear cubic term): V(x, t) = 0.1x 2 -|ψ| 2 , for t ∈ [0; 2]. The location of the ABC changes from xr = 10 to xr = 8 and finally xr = 6 without almost affecting the solution (at least at some visible levels of accuracy since the boundary condition is not transparent).

Figure 3 .

 3 Figure 3. Surrounding the finite computational physical domain of interest D Phy (white color) by an absorbing layer D PML to damp the waves entering the green region. The resulting complete computational domain is then D.

Figure 4 .

 4 Figure 4. PML for the 2D Schrödinger equation: contour plots of |ψ(x, t)| 2 at times t = 0, 4, 6, 7, 8, 9. The problem consists in the dynamics of a manufactured 2D soliton.

6 .

 6 ABCs/PMLs for the Klein-Gordon equation 6.1. ABCs for the field-particle Klein-Gordon equation

Figure 5 .

 5 Figure 5. ABCs for the 1D TDKGE: |ψ| (log-scale) for ABC 1 (top), ABC 2 (middle) and ABC 3 (bottom).

  ) 8 in the layer ]L x ; L * x [, where L * x = L x + δ x . In this test, PMLs are implemented in ]L x ; L * x [ and homogeneous Dirichlet boundary conditions are set on the left boundary x := 0 and at x = L *

  (left) the reference and the PML solutions. Notice that the PML region is ]L x ; L * x [=]10; 11[. In this Figure, a line at x = L x shows the location of the interface between the PML and the internal region. We also represent in Figure 6 (right) (in log-scale and as a function of time) the logarithm of the error between the reference and the PML solutions, i.e. (x, t, log |ψ ref (x, t) -ψ(x, t)|). This clearly illustrates the absorption at the right side (for the PML) of the domain, yielding a correct numerical solution, and the artificial reflection at the left side (for the homogeneous Dirichlet boundary condition), resulting in a wrong numerical solution.

Figure 6 .

 6 Figure 6. PMLs for the 1D TDKGE: comparison between the PML and the reference solutions (left) and logarithm of the absolute error between the reference and PML solutions (right).

Figure 7 .

 7 Figure 7. ABCs for the 1D potential-free Dirac equation: amplitude |ψ 1 | (log-scale) for the zeroth-order transport-like (top) and improved first-order transport-like conditions (bottom) obtained in Theorem (7.1).

  3). Let us define σ η (v) = Σ(-L η -v) + Σ(v -L η ), η = x,y or z. This leads to the modified PML Dirac operatorP (PML) D = -ic α x ∂ (PML) x + α y ∂ (PML) y + α z ∂ (PML) z + mc 2 β + V (x)1 4 ,where∂ (PML) x ψ = ∂ x ψ + i 2π σ x (x) R e it ω -iσ x (x)∂ x ψ(x, ω)dω.

Figure 8 .

 8 Figure 8. PML for the 2D Dirac equation: comparison between the PML (red curves) and AP (blue curves) approaches: the PML approach is much more accurate than the complex AP approach. In addition, increasing the size of the layer improves the accuracy of the solution. Contour plots of |ψ(x, t)| at times t = 0, 2, 4, 6, 8 (for σ 0 = 1 and δx = δy = 2) (By courtesy of Professor Olivier Pinaud).

FundingX.

  Antoine and Q. Tang thanks the support of the French ANR grants Bond (ANR-13-BS01-0009-01) and ANR-12-MONU-0007-02 BECASIM ("Modèles Numériques" call). E. Lorin thanks NSERC for the financial support via the Discovery Grant program.

Acknowledgement

The three authors warmly thanks Professor Olivier Pinaud (Colorado State University) for discussions about PMLs for the Dirac equation and for providing the simulations from Figure 8.