Absorbing Boundary Conditions for Relativistic ... - Semantic Scholar

Comment

Report 3 Downloads 207 Views

Noname manuscript No. (will be inserted by the editor)

X. Antoine · E. Lorin · J. Sater · F. Fillion-Gourdeau · A. D. Bandrauk

Absorbing Boundary Conditions for Relativistic Quantum Mechanics Equations

the date of receipt and acceptance should be inserted later

Abstract This paper is devoted to the derivation of absorbing boundary conditions for the KleinGordon and Dirac equations modeling quantum and relativistic particles subject to classical electromagnetic fields. Microlocal analysis is the main ingredient in the derivation of these boundary conditions, which are obtained in the form of pseudodifferential equations. Basic numerical schemes are derived and analyzed to illustrate the accuracy of the derived boundary conditions. Keyword words: Microlocal analysis, pseudo-differential operators, absorbing boundary conditions, wave equation, Dirac equation, Klein-Gordon equation, numerical approximation.

Contents 1 2 3 4 5 6

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Absorbing Boundary Conditions for the Klein-Gordon Equation Absorbing Boundary Conditions for the Dirac Equation in 2-d Other techniques . . . . . . . . . . . . . . . . . . . . . . . . . . Numerical Simulations . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

. . . . . .

1 5 17 28 31 38

1 Introduction 1.1 Introductory remarks This paper deals with the derivation of Transparent and Absorbing Boundary Conditions (TBC, ABC) for the time dependent Dirac (TDDE) and time dependent Klein-Gordon equations (TDKGE) modeling X. Antoine Institut Elie Cartan de Lorraine, Universit´ e de Lorraine, F-54506 Vandoeuvre-l` es-Nancy Cedex, France E-mail: [email protected] E. Lorin School of Mathematics and Statistics, Carleton University, Ottawa, Canada, K1S 5B6 Centre de Recherches Math´ ematiques, Universit´ e de Montr´ eal, Montr´ eal, Canada, H3T 1J4 E-mail: [email protected] J. Sater School of Mathematics and Statistics, Carleton University, Ottawa, Canada, K1S 5B6 E-mail: [email protected] F. Fillion-Gourdeau Centre de Recherches Math´ ematiques, Universit´ e de Montr´ eal, Montr´ eal, Canada, H3T 1J4 E-mail: [email protected] A. Bandrauk Laboratoire de Chimie Th´ eorique, Universit´ e de Sherbrooke, Sherbrooke, Canada, J1K 2R1 E-mail: [email protected]

the interaction of classical electromagnetic fields with quantum particles. Both equations are relativistic versions of the Schr¨ odinger equation for spin-0 particles (TDKGE) and spin±1/2 particles (TDDE). TDDE is an order 1 linear system of 4 equations, while the TDKGE is an order 2 linear equation. When dealing with the discretization of complex wave-like equations, such as the linear or nonlinear wave, Schr¨ odinger, Maxwell, Klein-Gordon, Dirac equations, a crucial question is the reduction of the size of the computational domain which allows to simultaneously gain in the data storage and decrease the overall computational complexity. However, when local or nonlocal internal waves reach the boundary of the domain where too simple boundary conditions (e.g. Neumann, Dirichlet) are imposed, they may be reflected and generate spurious waves interacting with physical waves. Avoiding this issue is an old problem that many papers have addressed, such as the famous pioneering Engquist and Majda’s paper [11], but also [17] and more recently [2, 3, 5]. In this paper, we derive two techniques to rigorously derive TBC & ABC for these two equations, using microlocal analysis and pseudodifferential operators in the spirit of [20, 21, 29, 31] and more specifically [25]. Some techniques used in this paper, were originally proposed by [4–6] in the framework of Maxwell’s equations and linear and nonlinear Schr¨ odinger equations. We then propose to adapt these ideas in the framework of relativistic quantum mechanics for laser-particle TDDE and TDKGE. The complexity of the solution to laser-molecule TDDE in realistic situations is well-known. For instance, when a intense femtosecond laser pulse (∼ 1015 s) of frequency ω0 interacts with a molecule, it involves multiscale phenomena, in time and space (zitterbewegung ∼ 10−21 s, [32]), in frequency (high order harmonics generation, often beyond 100ω0 ), etc. As a consequence, fine discretization, high order methods, adaptive mesh refinement, etc [12, 13, 26] are necessary. Other usual numerical issues with real space methods approximating TDDE, such as numerical dispersion, also necessitate special treatments [19]. The price to pay is naturally a huge computational complexity of the involved methods. In this paper, we derive ABCs for TDDE and TDKGE on a circle, although the technique can be extended to more general domains with smooth and convex boundary. 1.2 Microlocal approach General Principle. The goal of this paragraph is to recall some important notions and ideas about the derivation of TBC/ABCs using microlocal analysis in order to make the paper accessible to nonspecialists. For more details, we refer to [20, 21, 29] for the theory and to [5, 6] for applications to the linear and nonlinear Schr¨odinger equations and Maxwell’s equations [3]. The framework that we consider here is as follows. Let (1) be a 2-d partial differential equation (or system) in R2 × [0, T ] P (x, y, t, ∂t , ∂x , ∂y )u(x, y, t) = 0,

u(x, y, 0) = u0 (x, y),

(1)

where P = P (x, y, t, ∂t , ∂x , ∂y ) belongs to the set of order m pseudo-differential operators [31], denoted by OP S m . Time T > 0 is chosen such that system (1) is well-posed. Operator P has a symbol which is denoted by p ∈ S m (which is equal to det(P) where P = σ(P ) if (1) is a system). We now recall some basic definitions and facts about wavefronts, bicharacteristic strips, and operator factorizations. Details can be found in the classical references [1, 20, 21, 31]. For u ∈ D′ (R2 × R+ ), the wavefront WF(u), is defined as the set of singularities in real x = (x, y, t) and co-variable spaces ζ = (ξ, η, τ ) that is n WF(u) = (x, ζ) ∈ R6 − {0} : ζ ∈ Σx (u) , where Σx (u) is the set of singular frequencies defined as n oc Σx (u) = ζ : |b u(ζ)| 6 Ck0 (1 + |ζ|)−k0 , ∀k0 .

Bicharacteristic strip is a key notion in microlocal analysis and more specifically to derive TBC/ABC. For P belonging to OP S m , its leading symbol pm (x, ζ), has a degree of homogeneity equal to m and the bicharacteristic strips are integral curves of the following Hamilton-Jacobi equations ∂s x = ∇ ζ p m ,

∂s ζ = −∇x pm .

First, it is easy to see that the solution to (1) propagates along bicharacteristic strips, whose union is equal to the wavefront of the solution. In addition if (x0 , ζ 0 ) = (x0 , y0 , t0 , ξ0 , η0 , τ0 ) belongs to WF(u) so does (x, ζ) = x(s), y(s), t(s), ξ(s), η(s), τ (s) (for s > 0). If pm (x0 , ζ 0 ) = 0, such strips are called 2

the null-bicharacteristic strips. The wavefront of the solution is the union of all the bicharacteristics emanating from points where the principal symbol is zero. For the derivation of ABCs, a fundamental question is the identification of points (x0 , ζ 0 ) of the cotangent bundle of the domain boundary, which belong to the wavefront (see below). From there, outgoing/incoming bicharacteristic strips will be identified, and P can be factorized micro-locally into pseudo-differential operators, according to outgoing or incoming strips. We also recall that WF(u) is included in the hyperbolic region of P , that is the region where det(pm ) has real roots. These roots belong to WF(u) and, as recalled above, the corresponding null-bicharacteristic strips will remain in WF(u). A general algorithm to derive TBC/ABC consists of – – – –

Determining the symbol of Operator P . Determining the roots of the leading symbol det(pm ). Constructing the bicharacteristic strips and selecting the outgoing ones. Factorizing P as product of operators to specify the outgoing wave field.

The factorization of P is in principle possible thanks to theorems (which is recalled in a general framework) such as [1]. P Theorem 11 Let us consider P ∈ OP S m a differential operator of the form k+l6m ak,l Dxk Dtl . Then, there exists a factorization of the form P = Dx − λ1 (x, t)Dt · · · Dx − λm (x, t)Dt + order (m − 1) terms,

where λj ξ denotes the distinct real roots (in the hyperbolic region) in τ of the principal symbol pm = P k l k+l6m ak,l ξ τ . Here, we have set: Dx,t := −i∂x,t .

Geometrical Transformation. We now specify a little bit more the mixed problem that we consider. We first introduce a bounded domain Ω ⊆ R2 , which is assumed to be convex, with smooth fictitious boundary Γ , and such that it strictly contains supp(u0 ). We denote ΩT := Ω × [0, T ]. We search for ˜ to the following mixed problem B such that the solution u   P (x, y, t, ∂t , ∂x , ∂y )˜ u(x, y, t) = 0, on ΩT , B˜ u = 0, on Γ × [0, T ], (2) u ˜ (x, y, 0) = u0 (x, y), on Ω,

˜ kH is as small as possible. Roughly speaking, in the appropriate Hilbert space H, is such that ku|Ω − u when the error is zero, the BC are said to be transparent (TBC), and absorbing if not (ABC). For TDDE and TDKGE, we will search for B in the form ∂ + Λ+ u = 0 on Γ . Bu = ∂n

For any (x0 , y0 ) ∈ Γ and n((x0 , y0 )) normal vector to Γ at (x0 , y0 ), we study the outgoing and incoming (reflected) bicharacteristic strips. Although the approach which is used here is applicable to any domain with smooth boundary, it is natural (if possible) to choose a fictitious domain Ω with simple geometry. In this paper, we will consider in 2-d circular domains Ω of radius R and center (0, 0). In that case the curvature is constant which simplifies a lot the algebraic, analytical and numerical computations. We then introduce locally on the boundary a system of polar coordinates which are denoted by (r, θ) with co-variables (ρ, ω), and where (er , eθ ) denote the normal and angular direction at the boundary Γ . We also consider the annular region Γr = {(R + r, θ), r ∈ [0, ε)}. We now introduce the cotangent bundle to ΣT = Γ × [0, T ], denoted by T ∗ (ΣT ). Its hyperbolic region, in which the solution will propagate, is defined by n o H(ΣT ) = (θ0 , t0 , ω0 , τ0 ) ∈ T ∗ (ΣT ) : det pm (0, θ0 , t0 , ρ0 , ω0 , τ0 ) = 0 has real roots in ρ0 . Naturally the elliptic region (where the symbol does not vanish) does not contain any relevant information as (“pb u = 0”), and the glancing region is reduced to (θ0 , t0 , ξ0 , t0 ) due to geometrical assumptions on Ω. The Hamiltonian equations with initial condition (0, θ0 , t0 , ρ± 0 , ω0 , τ0 ) write  ∂pm ∂pm ∂pm   r′ (s) = , θ′ (s) = , t′ (s) = , ∂ρ ∂ω ∂τ   ρ′ (s) = − ∂pm , ω ′ (s) = − ∂pm , τ ′ (s) = − ∂pm . ∂r ∂θ ∂t 3

Formally, we recall that these equations can be rewritten in the form Hpm γ(s) = 0, where γ(s) = r(s), θ(s), t(s), ρ(s), ω(s), τ (s) denotes the bicharacteristic curves and Hpm ∈ T T ∗ (ΣT )), which has a basis (∂r , ∂θ , ∂t , ∂ρ , ∂ω , ∂τ ). We then select the outgoing (resp. incoming) bicharacteristic strips, which correspond to r′ (s)> 0 (resp. r′ (s) < 0), and we associate to these strips an operator Λ+ (resp. Λ− ). Moreover in H ΣT , P can be formally factorized in the form ∂ ∂ P = + Λ+ − Λ− + R, R ∈ OP S −∞ . (3) ∂n ∂n In addition, from Theorem (11), P can also be rewritten as follows. P = Pm + order (m − 1) terms, with

Pm := Dr − λ1 (r, θ, t)Dt · · · Dr − λm (r, θ, t)Dt ,

where λj ρ denotes real distinct roots, such that, the positive ones, would correspond to the outgoing waves. TBC is then of the form of a Dirichlet-to-Neumann (DtN) operator ∂n + Λ+ u = 0 on Γ,

where Λ+ is associated to the part of Pm involving the positive roots only. An instance of this approach will be detailed in the following section. In this goal, asymptotic analysis will be used to analytically determine Λ+ . We need to recall the following (general) definition, see [1]. Definition 11 Consider P (x, D) a scalar order m operator, with smooth symbol p(x, ξ) ∈ C ∞ . There exists a sequence {pm−j (x, ξ)}j∈N such that for ξ 6= 0, and any β ∈ N3 n X Dxβ p − pm−j = O |ξ|m−n−1 , ∀n ∈ N, and when |ξ| → +∞ . j=0

In that case, we note p(x, ξ) ∼

∞ X j=0

pm−j (x, ξ) and P (x, ∂x ) ∼

Another useful definition is as follows. Definition 12 A symbol p ∈ S m is said to be classical if p ∼ neous functions of degree m − j, that is for |ξ| > 1 and λ > 1,

∞ X

Pm−j (x, ∂x ), where Pj = Op(pj ).

j=0

∞ X

pm−j , where the pm−j are homoge-

j=0

pm−j (x, λξ) = λm−j pm−j (x, ξ) .

For TBC/ABCs, we can generally prove that Λ+ ∈ OP S m and that its total symbol λ+ satisfies the following expansion [1] λ+ ∼

∞ X

λ+ m−j ,

(4)

j=0

m−j where λ+ for j ∈ N and Λ+ := Op(pm−j ) ∈ OP S m−j . Then, Λ+ will be constructed m−j ∈ S P∞ + m−j as an approximation of Op( j=0 λm−j ). The sequence will be determined by the identification of the symbols of the same order in the left- and right-hand sides of (3) assuming that λ+ is classical, which means that it has an asymptotic expansion in the form (4).

1.3 Organization of the paper The rest of the paper is organized as follows. In Section 2, we derive ABCs for 1-d and 2-d TDKGE using the technique based on bicharacteristic strips described above. Section 3 is devoted to the construction of ABCs for TDDE using a diagonalization technique derived from several works, mainly [4–6]. In Section 4, we propose other possible approaches to construct ABCs for TDDE and TDKGE as well as some connections between ABCs for TDKGE and TDDE. Section 5 is dedicated to numerical illustrations. Finally, we conclude in Section 6. 4

2 Absorbing Boundary Conditions for the Klein-Gordon Equation 2.1 One-dimensional Case 2.1.1 Empirical absorbing conditions This short section is devoted to the derivation of simple empirical ABCs. We want here to emphasize the fact that simple conditions could be derived to limit spurious wave reflections. However, the weakness of this approach is double. First we do not provide any mathematical information on the reflected wave (accuracy, etc) and then the parameters are empirically chosen, not according to a rigorous procedure. The general principle, which is applicable to any kind of wave equation (Schr¨ odinger, wave equation, Maxwell, Klein-Gordon, etc) consists of absorbing the waves that reach the boundary of the computational domain in a buffer zone using a well designed function. We consider the 1-d KleinGordon (KG) and Dirac (D) boundary-value problems   PKG,D ψ(x, t) = 0, x ∈ [−a, a], BKG,D ψ(x, t) = 0, x ∈ {−a, a},  ψ(x, 0) = ψ0 (x) ∈ C(or C2 ) ,

with where

PKG = c2 ∂x2 − ∂t2 − ieV ∂t − ic2 qeAx∂x − ie∂t V − ic2 e∂x Ax + e2 V 2 − c2 e2 A2x − m2 c4 , PD = i∂t − αx − ic∂x − eAx (x, t) − βmc2 + eI2 Vc (x), αx =

01 10

,β=

1 0 0 −1

, I2 =

10 01

.

In the equation, e denotes the particle charge, m its mass and c the speed of light. The electric potential (V, Ax ) depends on (x, t). Coulomb potential Vc depends on space only. The goal is to provide a boundary operator BKG,D , limiting spurious reflections at the fictitious domain boundary {−a, a}. Here is the principle of the proposed approach – Let us denote by ε and γ two positive real numbers, where ε ≪ a. – We solve   PKG,D ψ(x, t) = 0, x ∈ [−a, a], BKG,D ψ(x, t) = 0, x ∈ {−a, a},  ψ(x, 0) = ψ0 (x), x ∈ [−a, a] .

– We multiply solution ψ by a function defined for instance, as follows  x−a 2  1 − e−γ( x−a+ε ) , x ∈ (a − ε, a], x+a f (x, t) = 1 − e−γ( x+a−ε )2 , x ∈ [−a, −a + ε),  0, elsewhere .

(5)

(6)

The two constants ε and γ are empirically chosen and may depend on Ax , V and Vc . Numerical experience shows that it can be an easy-to-implement alternative to more rigorous (but also more computationally costly) boundary conditions. Note that this approach is easily extendable to higher dimensions, although in that case the empirical choice of parameters is even more problematic (as ε, γ may also depend on the geometry of the fictitious boundary). 2.1.2 Absorbing boundary conditions for Klein-Gordon’s equation in 1-d with a constant electric field Let us first consider that we are solving the TDKGE for a constant electric field, which means that V and Ax are both x and t independent. In this case, the TDKGE is given by PKG ϕ = (c2 ∂x2 − ∂t2 − ieV ∂t − ic2 eAx ∂x + e2 V 2 − c2 e2 A2x − m2 c4 )ϕ(t, x) = 0. 5

Let us denote by τ the time Fourier covariable. In this case, Fourier transforming the above equation leads to c2 ∂x2 − ic2 eAx ∂x + (τ 2 + eτ V + e2 V 2 − c2 e2 A2x − m2 c4 ) ϕ(τ, b x) = 0. Computing the roots of this equation yields r 3 ieAx 1 ± −τ 2 − eV τ − e2 V 2 + c2 e2 A2x + m2 c4 , λ1 (τ ) = ± 2 c 4

(7)

and the solution writes as the superposition of two plane waves +

−

ϕ(τ, b x) = A+ eλ1 (τ )x + A− eλ1 (τ )x .

(8)

The left traveling solution is associated with λ− 1 (easy to see for the field free equation). Therefore, to get the outgoing solution, one must have A− = 0. Deriving according to ∂x provides the following relation (∂x − λ+ b x) = 0, (9) 1 )ϕ(τ,

which shows that the equation selects the outgoing solution to the domain boundary. Going back to the physical domain by inverse Fourier transform, one gets the TBC r 3 iceAx (c∂x − ∂t2 + ieV ∂t − e2 V 2 + c2 e2 A2x + m2 c4 − )ϕ(t, x) = 0, on ΣT . (10) 4 2 The above DtN operator is nonlocal in time. In the field free case, one has (c∂x − ∂t )ϕ(t, x) = 0,

(11)

which is just the local classical TBC for the wave equation part going to the left. 2.1.3 Absorbing boundary conditions for Klein-Gordon’s equation in 1-d with a variable electric field The approach which is used to derive ABCs for the Klein-Gordon equation was originally presented in [3–5] for linear and nonlinear Schr¨odinger equations. We focus here on the computational details. We recall that the field-particle Klein-Gordon equation PKG ϕ = 0 writes h i 2 (i∂t − eV (x, t))2 − c2 i∂x + eAx (x, t) − m2 c4 ϕ = 0, in [−a, a], which we expand into

h i PKG ϕ = c2 ∂x2 − ∂t2 − ieV ∂t − ic2 eAx ∂x + OD ϕ = 0,

in [−a, a],

where OD is the zeroth-order operator given by

OD := −ie∂t V − ic2 e∂x Ax + e2 V 2 − c2 e2 A2x − m2 c4 .

(12)

We search for an ABC at x = ±a under the form of approximate DtN operators. The principle of the technique developed in [5] consists of factorizing PKG in two operators in a neighborhood of any (t0 , τ0 ) in the hyperbolic region H ΣT , as follows PKG (t, x, ∂t , ∂x ) = c∂x − Λ− (t, x, ∂t , ∂x ) c∂x + Λ+ (t, x, ∂t , ∂x ) + R(t, x, ∂t , ∂x ),

where R(t, x, ∂t , ∂x ) is a smoothing operator (in OP S −∞ ) and Λ± (t, x, ∂t , ∂x ) are two pseudodifferential operators to be determined, such that Λ+ (t, x, ∂t , ∂x ) (resp. Λ− (t, x, ∂t , ∂x )) has only outgoing (resp. only incoming) bicharacteristic strips at the domain boundary. First, we notice that (omitting (t, x, ∂t , ∂x ) for the sake of conciseness) c∂x − Λ− c∂x + Λ+ = c2 ∂x2 + c Λ+ − Λ− ∂x + cOp ∂x λ+ − Λ− Λ+ . 6

By identification of the operators of different orders in front of ∂x , we obtain c Λ+ − Λ− = −ic2 eAx and

cOp(∂x λ+ ) − Λ− Λ+ = −∂t2 − ieV ∂t + OD . Denoting λ± = σ(Λ± ), we get at the symbol level c(λ+ − λ− ) = −ic2 eAx , c∂x λ+ − λ− #λ+ = τ 2 − eV τ + σ(OD ) .

(13)

P ± Now according to [31], any classical symbol of S 1 can be asymptotically expanded as λ± ∼ +∞ j=0 λ1−j , ± where λ1−j is the symbol of a pseudodifferential operator of homogeneity degree 1 − j, for j ∈ N. In order to determine the expansion terms, it is sufficient to recursively identify the terms of same homogeneity in system (13). For j = 0, one gets + λ1 − λ− 1 = 0, (14) + 2 −λ− λ 1 1 = τ . For the next homogeneity, we obtain for j = 1 − 2 c(λ+ 0 − λ0 ) = −ic eAx , + − + − + + ∂x λ1 − λ1 λ0 − λ0 λ1 + i∂ξ,τ λ− 1 ∂x,t λ1 = −eV τ . For j = 2, we have the following expression  + − λ−  −1 ) = 0,  c(λ−1 + + − + − + ∂x λ0 − λ− 1 λ−1 − λ0 λ0 − λ−1 λ1 1 2 − 2 +  + − +  +i∂ξ,τ λ− 1 ∂x,t λ0 + i∂ξ,τ λ0 ∂x,t λ1 + ∂ξ,τ λ1 ∂x,t λ1 = σ(OD ), 2

and finally for j = 3  c(λ+ − λ−  −2 ) = 0,   ∂ λ−2 + − + − + − + − +  x −1 − λ1 λ−2 − λ0 λ−1 − λ−1 λ0 − λ−2 λ1 − + − − + +i∂ξ,τ λ1 ∂x,t λ−1 + i∂ξ,τ λ0 ∂x,t λ+ 0 + i∂ξ,τ λ−1 ∂x,t λ1   1 i 1  − 2 − 3 + + + 2 2 3 2  λ− + ∂ξ,τ 1 ∂x,t λ0 + ∂ξ,τ λ0 ∂x,t λ1 − ∂ξ,τ λ1 ∂x,t λ1 = 0 . 2 2 6

(15)

(16)

(17)

By using system (14), we have

− λ+ 1 = λ1 , + λ1 = ±(−τ 2 )1/2 = ±iτ .

At this point, we have to select the forward bicharacteristics, which corresponds to the one that satisfies x′ (s) > 0 at x = +a (resp. x′ (s) < 0 at x = −a), where x′ (s) = ∂ξ p2 , ξ ′ (s) = −∂x p2 ,

t′ (s) = ∂τ p2 , τ ′ (s) = −∂t p2 ,

with principal symbol p2 (x, t, ξ, τ ) = −c2 ξ 2 + τ 2 and for (±a, t0 , ξ0 , τ0 ) such that p2 (x0 , t0 , ξ0 , τ0 ) = 0. We get the following roots ξ = ±τ /c in the hyperbolic region (corresponding to propagating waves). Now as x′ (s) = −2c2 ξ, we conclude that the forward bicharacteristics is then the one associated to −τ /c ± at x = +a (and +τ /c for x = −a), that is λ+ 1 = iτ . Fixing the principal symbols λ1 (as previously) ± defines uniquely the asymptotic expansion of λ and also the classical first-order pseudodifferential 7

operators Λ± . For building an accurate boundary condition, we compute the next three terms of the symbolical asymptotic expansion. This leads to  i  +    λ0 = 2(eV − ceAx ),   i ih i 1 i + (18) λ = σ(O σ(O σ(O ) τ −1 , ) + ) − σ(O ) − G G F D x t −1  2 2 2 4  i h h   i 1 i  1  λ+ σ(OD ) − σ(OGx ) + σ(OGt ) − σ(OF ) τ −2 , −2 = 4 ∂x − ∂t − ieV 2 2 4

where the operators OF , OGx and OGt are defined as follows  OGx := e∂x V − ce∂x Ax ,     OGt := e∂t V − ce∂t Ax ,     OF := e2 V 2 − c2 e2 A2x .

(19)

We can now state the following factorization result, from which will be constructed ABCs of arbitrary order. Details of the proof can easily be deduced from [6].

− Proposition 21 There exists two uniquefirst-order pseudodifferential operators Λ+ k and Λk such that in a neighborhood of any (t0 , τ0 ) ∈ H ΣT − PKG = c∂x + Λ+ c∂ − Λ + R2−k , x k k

Pk−1 + + + + where R2−k ∈ OP S 2−k and Λ+ j=0 λ1−j . Here, we have fixed λ1 := iτ and k = Op σk , with σk = + the λ1−j are defined in (18) for j = 1, 2, 3. P∞ By truncating the formal expansion j=0 λ+ 1−j up to the k first terms, a k-th order ABC is obtained through the boundary condition c∂n + Λ+ (20) k ϕ = 0, on ΣT . We recall that the Riemann-Liouville operator is defined as follows (for ν > 0) Z t 1 D−ν f (t) = (t − s)(ν−1) f (s)ds Γ (ν) 0

(21)

and has a symbol which is equal to (iτ )−ν . When Ax = V = 0, which is the field-free TDKGE, we easily see that σ3+ =

1 X

λ+ 1−j = iτ +

j=0

m2 c4 2iτ

so that Λ+ 2 = ∂t +

m2 c4 −1 ∂ . 2 t

For Ax 6= 0, V 6= 0, we have  + Λ 1 = ∂t ,     + iceAx ieV − , Λ 2 = ∂t + 2 2   iceA 1 i i ieV 1  + x  Λ 3 = ∂t + − − OD − OF − OGx + OGt ∂t−1 2 2 2 4 2 2 where OD , OF , OG are defined in Eq. (12) and (19). 8

2.2 Two-dimensional TDKE 2.2.1 Field-free TDKGE in 2-d We consider a circular domain of radius R, on which we plan to provide an ABC. In that case, it is natural to consider polar coordinates in all Ω and conditions have to be imposed at r = R. Another possibility is to keep cartesian coordinates and to search for local coordinates on the circle. In fine, the treatment at the boundary will be the same. However the second approach is more general as it is applicable, in principle, to any smooth boundary, whatever the system of coordinates used inside the domain. For r ∈ [0, ε) with ε > 0 small enough and denoting by χr = (R + r) the radius of the circular domain Ωr with boundary Γr := ∂Ωr , the TDKGE in local polar coordinates in Σr = Γr × [0, T ] writes, h i 1 1 (r,θ) PKG ϕ = − ∂t2 + c2 ∂r2 − c2 ∂r + c2 2 ∂θ2 − m2 c4 ϕ(t, r, θ) = 0, χr χr

(22)

for (R + r, θ) ∈ Γr , and t ∈ R+ . Variables (r, θ, t) have duals denoted by (ρ, ξ, τ ). The principal symbol (r,θ) of PKG is given by p2 (r, θ, t, ρ, ξ, τ ) = τ 2 − c2 ρ2 −

c2 2 ξ . χ2r

We are first interested in the bicharacteristic strips at the domain boundary (where r = 0). We have to select a point, in the hyperbolic region of T ∗ Σr , that will belong to the wavefront of the solution, p2 (0, θ0 , t0 , τ0 , ρ0 , ξ0 ) = τ02 − c2 ρ20 −

c2 2 ξ = 0. χ2r 0

τ02 ξ02 1/2 The corresponding roots are ρ± = ±( − ) . The hyperbolic region H Σr , of the cotangent 0 2 2 c χr bundle T ∗ Σr is then the set of (θ0 , t0 , ξ0 , τ0 ) ∈ T ∗ Σr such that τ02 /c2 − ξ02 /χ2r > 0 (then pm has real roots). In the elliptic region the waves are evanescent. The Hamilton-Jacobi equations write   r′ (s) = ∂ρ pm , ρ′ (s) = −∂r pm , θ′ (s) = ∂ξ pm , ξ ′ (s) = −∂θ pm ,  t′ (s) = ∂ p , τ ′ (s) = −∂ p . τ m t m The outgoing characteristics are defined by the condition r′ (s) > 0. Since r′ (s) = ∂ρ pm = −2ρ, we deduce that the outgoing bicharacteristics corresponds to ρ− = −(

τ 2 ξ 2 1/2 − ) . c2 χ2r

(r,θ) To build the ABC, we now consider a local factorization of the operator PKG in H Σr as (r,θ)

PKG (r, θ, t, ∂t , ∂r , ∂θ ) = (c∂r − Λ− (r, θ, t, ∂θ , ∂t ))(c∂r + Λ+ (r, θ, t, ∂θ , ∂t )) + R,

where Λ± ∈ OP S 1 has a total symbol λ± and a principal symbol equal to ρ± , respectively, and R ∈ OP S −∞ . We now have to explicitly compute Λ± . Since we have (r,θ) PKG ∼ c∂r − Λ− c∂r + Λ+ = c2 ∂r2 + c Λ+ − Λ− ∂r + cOp(∂r λ+ ) − Λ− Λ+ , (r,θ)

the expression (22) of PKG leads to the following system by identification  c2   =−  c(λ+ − λ− ) χr 2  + − + 2  c∂r λ − λ #λ = τ − c ξ 2 − m2 c4 .  χ2r 9

P∞ − + As an element of S 1 , λ± can be expanded as λ± ∼ j=0 λ± 1−j . We also recall that, denoting λ #λ = σ Λ− Λ+ , P∞ P∞ − P∞ + P∞ − + σ Λ− Λ+ = λ− #λ+ = j=0 λ1−j j=0 λ1−j − i j=0 ∂ρ,ξ,τ λ1−j j=0 ∂r,θ,t λ1−j −

P∞ 2 P 1 P∞ 2 (−i)α α − + ∂ λ− ∂ α λ+ . ∂ λ ∂ λ + j=0 ρ,ξ,τ 1−j j=0 r,θ,t 1−j α>3 2 α! ρ,ξ,τ 1−j r,θ,t 1−j

In the following, we will use the notation

(−i)α α α ∂ λ− ∂r,θ,t λ+ . σα Λ− Λ+ = α! ρ,ξ,τ

As before, we consider symbols up to those in S −2 that is, we approximate λ± by identify terms of the same homogeneity order. This leads in λ− #λ+ to – – – –

P3

j=0

λ± 1−j , and we

+ 2: λ− 1 λ1 . + − + − + 1: λ− 1 λ0 + λ0 λ1 − i∂ρ,ξ,τ λ1 ∂r,θ,t λ1 . − 2 + − + − + − + + − + 1 2 0: λ−1 λ1 + λ0 λ0 + λ1 λ−1 − i∂ρ,ξ,τ λ− 1 ∂r,θ,t λ0 − i∂ρ,ξ,τ λ0 ∂r,θ,t λ1 − 2 ∂ρ,ξ,τ λ1 ∂r,θ,t λ1 − + − + − + − + − + − + -1: λ−2 λ1 + λ−1 λ0 + λ0 λ−1 + λ1 λ−2 − i∂ρ,ξ,τ λ1 ∂r,θ,t λ−1 − i∂ρ,ξ,τ λ0 ∂r,θ,t λ0 1 2 1 2 i 3 + − 2 − 2 − 3 + + + − i∂ρ,ξ,τ λ− −1 ∂r,θ,t λ1 − ∂ρ,ξ,τ λ1 ∂r,θ,t λ0 − ∂ρ,ξ,τ λ0 ∂r,θ,t λ1 + ∂ρ,ξ,τ λ1 ∂r,θ,t λ1 . 2 2 6

Order Order Order Order

− + In fact, in the field-free configuration, as λ± = 1−j only depends on r, τ, ξ, we easily verify that σα Λ Λ 0 for α 6= 0. We deduce  −  c(λ+ 1 − λ1 ) = 0, c2 +  −λ− = τ 2 − 2ξ2 , 1 λ1 χr and we get up to leading order:

 + −  λ1 = λ1

2  λ+ 1 = ±(−τ +

c2 2 1/2 ξ ) . χ2r

As we select the forward bicharacteristics

λ+ 1 =

− τ2 +

c2 2 1/2 ξ , χ2r

then, we obtain   and so



− c(λ+ 0 − λ0 )

=−

− + − + c∂r λ+ 1 − λ1 λ0 − λ0 λ1 = 0,

    λ−

= λ+ 0 +

0

This provides

c2 , χr

c , χr

c  +  = −c∂r λ+  −λ+ 1. 1 2λ0 + χr   +   λ0

c c + ∂ λ , − r 1 χr 2λ+ 1 c2 2 1/2   2  c∂r λ+ . 1 = c∂r − τ + 2 ξ χr =

10

Considering the two next orders, we have  −  c(λ+ −1 − λ−1 ) 

and

= 0,

− + − + − + 2 4 c∂r λ+ 0 − λ1 λ−1 − λ0 λ0 − λ−1 λ1 = −m c ,

 −  c(λ+ −2 − λ−2 ) For λ+ −1 , we get



= 0,

− + − + − + − + c∂r λ+ −1 − λ1 λ−2 − λ0 λ−1 − λ−1 λ0 − λ−2 λ1 = 0.

 + +  λ−1 = λ−1 , 1 + − + 2 4  λ− −1 = + m c + c∂r λ0 − λ0 λ0 . 2λ1

A little algebra finally gives us λ+ −2 as

 + +  λ−2 = λ−2 , 1 c + + + c∂r λ+  λ− −2 = −1 . + − λ−1 2λ0 + χr 2λ1

Thus, summing up the four first terms of the asymptotic expansion leads to c2 2 1/2 c c 1 + + − + 2 λ = (−τ + ξ ) + ∂ λ − + + [m2 c4 + c∂r λ+ r 1 0 − λ0 λ0 ] j=0 1−j + 2 χr χr 2λ1 2λ1 c 1 + + + c∂ λ 2λ + = + − λ+ r −1 . −1 0 χr 2λ1

P3

We denote

P3

j=−1

(23)

+ + + λ+ 1−j by λ(−2,1) and Λ(−2,1) = Op(λ(−2,1) ). Thus, we have proven:

− Proposition 22 There exist 2 operators Λ+ (−2,1) , Λ(−2,1) such that in a neighborhood of any (ω0 , t0 .ξ0 , τ0 ) ∈ H Σr

r,θ − c∂ − Λ PKG = c∂x + Λ+ x (−2,1) + R−3 , (−2,1)

where R−3 ∈ OP S −3 . As a consequence, the TBC at Ω0

is approximated by the following ABC

c∂r + Λ+ ϕ = 0 c∂r + Λ+ (−2,1) ϕ = 0 .

In the next section we will explicitely define the pseudo-differential operators associated with the symbols appearing in equation (23). 11

2.2.2 Field-particle TDKGE in 2-d We now consider the full TDKGE including the electromagnetic field: A(x, y, t) and V (x, y, t) which represent the vectorial and scalar potentials satisfying Maxwell’s equations. They are assumed given and regular enough. As before, we fix ~ = 1. The equation writes h

i∂t − eV (x, y, t)

2

− c2

i 2 1 ∇ − eA(x, y, t) − m2 c4 ϕ = 0 . i

Expanding the equation leads to h i − ∂t2 + c2 ∆ − ieV ∂t − ic2 eA · ∇ − ie∂t V − ic2 e∇ · A + e2 V 2 − c2 e2 A2 − m2 c4 ϕ = 0 .

The domain Ω is again a disk of radius R. We then need to rewrite the TDKGE in polar coordinate in Σr , for r ∈ [0, ε) with ε > 0 small enough, that is (details are skipped) h c2 c2 c2 r,θ PKGL ϕ = − ∂t2 + c2 ∂r2 − ∂r + 2 ∂θ2 − ieV ∂t − ic2 eAr ∂r − i eAθ ∂θ χr χr χr i 2 2 ic e c ∂θ Aθ − c2 e2 A2r − c2 e2 A2θ + e2 V 2 − m2 c4 , −ie∂t V − ic2 e∂r Ar − i eAr − χr χr

where

(24)

Ax = Ar cos(θ) − Aθ sin(θ), Ay = Ar sin(θ) + Aθ cos(θ) .

Let us define the following operators  c2    + c2 µAr ∂r , OB := −   χr   c2 µ A θ ∂θ , OF := −  χr    c2 µ c2 µ    OD := −µ∂t V − c2 µ∂r Ar − Ar − ∂θ Aθ + e2 V 2 − c2 e2 A2 − m2 c4 , χr χr

and the functions  µ := ie,     c2   qOB := − − c2 µAr ,    χr c2 µ qOF := − Aθ ,    χr    c2 µ c2 µ    qOD := −µ∂t V − c2 µ∂r Ar − Ar − ∂θ Aθ + e2 V 2 − c2 e2 A2 − m2 c4 . χr χr

This allows us to rewrite equation (24) in a more compact form

c2 2 ∂ − µV ∂ + O + O + O t B F D ϕ = 0. χ2r θ (r,θ) As before we plan to factorize this operator PKGL in the form c∂r − Λ− c∂r + Λ+ . We recall that c∂r − Λ− c∂r + Λ+ = c2 ∂r2 + c Λ+ − Λ− ∂r + Op(∂r λ+ ) − Λ− Λ+ . c2 ∂r2 − ∂t2 +

By identification of the operators in front of ∂r , we have  1    Λ + − Λ − ∂r = O B , c c2  + −  cOp(∂r λ ) − Λ Λ+ = −∂t2 − µV ∂t + ∂θ2 + OF + OD ,  χ2r 12

and so in terms of symbols     λ+ − λ−

=

qO B , c

c2    c∂r λ+ − λ− #λ+ = τ 2 − iµV τ − 2 ξ 2 + iqOF ξ + qOD . χr

It is now possible to identify the sequence (λ± 1−j )j . First one, gets  +  λ1 − λ− 1 = 0, c2 − +  −λ1 λ1 = τ 2 − 2 ξ 2 , χr that leads to

 + −  λ1 = λ1 ,

2  λ+ 1 = ±(−τ +

We deduce that

c2 2 1/2 ξ ) . χ2r

c2 2 1/2 ξ ) , (25) χ2r corresponding to the outgoing bicharacteristic strips, for the same reason as in the field-free TDKGE. Both equations have indeed the same principal symbol. In addition, by using the explicit formula of σα Λ− Λ+ for α = 0, 1, 2, we get first  qO  + λ0 − λ− = B, 0 c  c∂ λ+ − λ− λ+ − λ− λ+ = −iµV τ + iqOF ξ . r 1 1 0 0 1 2 λ+ 1 = (−τ +

This leads to

1 + + − iµV τ + iqOF ξ − c∂r λ1 + 2λ1 1 − λ0 = − + − iµV τ + iqOF ξ − c∂r λ+ 1 − 2λ1 λ+ 0 = −

Going to the two next order provides  +  λ−1 − λ− −1 and



qO B , 2c qO B . 2c

(26)

= 0,

− + − + − + − + − + c∂r λ+ 0 − λ1 λ−1 − λ0 λ0 − λ−1 λ1 + i∂ξ λ1 ∂θ λ0 + i∂τ λ1 ∂t λ0 = qOD ,

 + λ−2 − λ−  −2   

= 0,

+ − − + − + − + − + c∂r λ+ −1 − λ1 λ−2 − λ0 λ−1 − λ−1 λ0 − λ−2 λ1 + i∂ξ λ1 ∂θ λ−1 + 1 2 + − + − + − 2 + i∂τ λ− 1 ∂t λ−1 + i∂ξ λ0 ∂θ λ0 + i∂τ λ0 ∂t λ0 + ∂ρ,ξ,τ λ1 ∂r,θ,t λ0 = 0 . 2 Some direct computations lead to

   

λ+ −1 = Similarly, we deduce

1 + − + − + + − + − qOD + c∂r λ0 + i∂ξ λ1 ∂θ λ0 + i∂τ λ1 ∂t λ0 − λ0 λ0 . 2λ1

1 + − + − + − λ+ −1 [λ0 + λ0 ] + c∂r λ−1 + i∂ξ λ1 ∂θ λ−1 + 2λ+ 1 1 2 − 2 + + − + − + +i∂τ λ− 1 ∂t λ−1 + i∂ξ λ0 ∂θ λ0 + i∂τ λ0 ∂t λ0 + ∂ρ,ξ,τ λ1 ∂r,θ,t λ0 . 2 As before we can state the following theorem λ+ −2 =

13

(27)

(28)

− Theorem 21 There exist 2 operators Λ+ (1−k,1) , Λ(1−k,1) such that in a neighborhood of any (ω0 , t0 , ξ0 , τ0 ) ∈ H Σr r,θ − c∂ − Λ PKGL = c∂x + Λ+ x (1−k,1) + R−k , (1−k,1)

where R−k ∈ OP S −k and

Λ+ (1−k,1) = Op

k X j=0

λ+ 1−j ,

where (λ+ 1−j )j=0,...,3 , are given by (25), (26), (27) and (28). At this point, we plan to find explicit ABCs which are defined by local (= differential) operators. This means that we must approximate locally the functions defining the symbols by polynomial functions in the dual variables. We recall that √ 1 1 1 + x = 1 + x − x2 + O(|x|3 ), 2 8 for |x| < 1. As a consequence for |τ | large (high time-frequency regime) λ+ 1 is expanded in: s s 2 2 ic2 ξ 2 ic4 ξ 4 c ξ c2 ξ 2 + − 4 3 + O(|τ |−5) . λ1 = −τ 2 (1 − 2 2 ) = iτ 1 − 2 2 = iτ − 2 χr τ χr τ 2χr τ 8χr τ Remark 21 In the following, λ± j will be computed in the regime |τ | ≫ 1. + Next, we approximate λ+ 0 which is a function of ∂r λ1 and

1 = 2λ+ 1

1 . First, we easily show that 2λ+ 1

1 i −1 ic2 2 −3 τ − 2 ξ τ + O(τ −5 ). = − 2 4χr c2 2 −1 2iτ − 2 iξ τ χr

Furthermore, we have c∂r λ+ 1 =

ic2 ξ 2 −1 ic4 ξ 4 −3 τ = τ + O(τ −5 ) . χ3r 2χ5r

As a consequence, λ+ 0 can be approximated by qO B 1 + + − iµV τ + iqOF ξ − c∂r λ1 + 2λ1 2c i −1 ic2 2 −3 ic4 4 −5 qO ic2 ξ 2 −1 ic4 ξ 4 −3 = − iµV τ + iqOF ξ − + B τ + 2ξ τ + 4ξ τ τ − τ 2 4χr 8χr χ3r 2χ5r 2c

λ+ 0 =

and λ+ 0 =

qO qO µ V + B − F ξτ −1 + 2 2c 2

c2 µc2 + V ξ 2 τ −2 + O(τ −3 ) . 2χ3r 4χ2r

From there, we derive the expression up to a O(τ −3 ). This choice is motivated by the complexity of the symbolic expression beyond that order, that would lead to complex operators, hard to approximate qO B + . and that could deteriorate stability and efficiency of the overall scheme. Note that λ− 0 = λ0 − c + Let us now approximate λ−1 that is defined by λ+ −1 =

1 + − + − + − + + − qOD + c∂r λ0 + i∂ξ λ1 ∂θ λ0 + i∂τ λ1 ∂t λ0 − λ0 λ0 . 2λ1 14

+ − In that goal, we expand and approximate c∂r λ+ 0 and λ0 λ0 as

c∂r λ+ 0

µc 1 c c = ∂r V + ∂r qOB − ∂r qOF ξτ −1 + 2 2 2 2 2χr

3c2 µc2 µc2 − 2 − ∂r V ξ 2 τ −2 + O(τ −3 ) V + χr χr 2

and − λ+ 0 λ0

where qOH :=

2 2 qO µ2 2 q O µ −1 B = V − 2 − qOF V ξτ + µV qOH + F ξ 2 τ −2 + O(τ −3 ) , 4 4c 2 4

µc2 c2 + V . In addition, we have 2χ3r 4χ2r + i∂ξ λ− 1 ∂θ λ0 =

µ 2

∂θ V +

∂θ qOB c2 −1 c2 ∂θ qOF 2 −2 ξτ − ξ τ + O(τ −3 ) 2c χ2r 2χ2r

and µ ∂t qOB ∂t qOF −1 c2 µ ∂t qOB 2 −2 + i∂τ λ− ∂ λ = − ∂ V − + ξτ − ∂t V + ξ τ + O(τ −3 ) . t t 1 0 2 2 2c 2 2χr 2 2c This provides the expansion λ+ −1

i −1 ∂t qOB µc 1 µ −3 = − τ + O(τ ) − qOD − ∂t V − + ∂r V + ∂r q O B 2 2 2c 2 2 2 µ c2 c ∂ q ∂t qOF −1 µ2 2 q O θ O B B ∂θ V + ξτ −1 + ξτ − V + 2 − ∂r qOF ξτ −1 + 2 4 4c 2 2 2c χ 2 r µ + qOF V ξτ −1 + O(τ −2 ) . 2

We now set  2 iq iµc i iµ i iµ2 2 iqO   B  qOI := OD − ∂r V − ∂r q O B + V − + ∂t V + ∂t q O B , 2 4 4 8 8c2 4 2 4c µ i ∂θ qOB ic iµqOF ic   V − ∂t q O F − ∂θ V + .  qOJ := ∂r qOF − 4 4 4 2 2c 2χ2r

−3 The symbol λ+ : −1 can finally be approximated by truncating the terms higher than τ −1 λ+ + qOJ ξτ −2 + O(τ −3 ) . −1 = qOI τ

The next and last step consists in expanding λ+ −2 =

1 + − + − + − λ+ −1 [λ0 + λ0 ] + c∂r λ−1 + i∂ξ λ1 ∂θ λ−1 + 2λ+ 1 1 2 + − + − + − 2 + +i∂τ λ− 1 ∂t λ−1 + i∂ξ λ0 ∂θ λ0 + i∂τ λ0 ∂t λ0 + ∂ρ,ξ,τ λ1 ∂r,θ,t λ0 . 2

− + We note that λ+ 0 + λ0 = 2λ0 −

qO B . We set c

qOG = 2λ+ 0 −

qO B c

= µV − qOF ξτ −1 +

c2 µc2 + V ξ 2 τ −2 + O(τ −3 ) . χ3r 2χ2r 15

Moreover, we have  c2   ∂θ λ+ = 2 ∂θ qOI ξτ −2 + O(τ −3 ),  i∂ξ λ− 1 −1   χr       + −1  i∂τ λ− − ∂t qOJ ξτ −2 + O(τ −3 ),  1 ∂t λ−1 = −∂t qOI τ      ∂θ qOB −1 iqOF µ − + i∂ λ ∂ V + τ ∂ λ = − ξ θ θ 0 0  2 2 2c    µ ∂θ q O B c 2 µc2 −2 iqOF ∂θ qOF    + 2i ∂ V + + V ξτ + O(τ −3 ), + θ   4 2 2c 2χ3r 4χ2r          i∂τ λ− ∂t λ+ = iqOF µ∂t V + ∂t qOB ξτ −2 + O(τ −3 ) . 0 0 2 2 2c

+ 2 2 −3 Finally, we observe that ∂ρ,ξ,τ λ− ). We now deduce that 1 ∂r,θ,t λ0 /2 = O(τ

λ+ −2 =

qO µ ∂θ qOB −2 qOI qOG − c∂r qOI + ∂t qOI − F ∂θ V + τ . 2 4 2 2c

i

It is now possible to write out the whole approximation

qO qO ic2 2 −1 µ ξ τ + V + B − F ξτ −1 λ+ 1−j = iτ − 2 2χ 2 2c 2 2 r c µc2 + + V ξ 2 τ −2 + qOI τ −1 + qOJ ξτ −2 2χ3r 4χ2r i qO µ ∂θ qOB −2 qOI qOG − c∂r qOI + ∂t qOI − F ∂θ V + τ + O(τ −3 ) + 2 4 2 2c µ qO ic2 qO = iτ + V + B + − 2 ξ 2 − F ξ + qOI τ −1 2 2c 2χ 2 r 2 c µc2 2 i + + 2 V ξ + qOJ ξ + qOI qOG − c∂r qOI + ∂t qOI 3 2χr 4χr 2 ∂θ qOB −2 qO F µ ∂θ V + τ + O(τ −3 ) . − 4 2 2c P3 + + Let us set: λ+ j=0 λ1−j as the four symbols approximation of λ . Going back to operators, (−2,1) := 1 Λ+ (−2,1) is given by µ c2 1 c2 2 1 µc2 2 −1 Λ+ = ∂ + ∂ − − V + O + − ∂ − O ∂ + iO + V ∂θ t B F θ I t (−2,1) 2 2c 2χ2r θ 2 2χ3r 4χ2r 1 ∂θ OB −2 µ i ∂t . −iOJ ∂θ + OI OG − c∂r OI + ∂t OI − OF ∂θ V + 2 4 2 2c P3

j=0

We have proven the following proposition.

Proposition 23 For |τ | sufficiently large, a family of ABCs for the field-particle TDKGE is given by j = 1, 2, 3, 4, c∂r + Λ+ (1−j) ϕ = 0, where for

– j = 0: Λ+ (1,1) = ∂t – j = 1: 1 µ + Λ+ (0,1) = Λ(1,1) + 2 V + 2cOB 1

O· stand for op(qO· )

16

– j = 2: + Λ+ (−1,1) = Λ(0,1) +

−

c2 2 1 −1 ∂ − O ∂ + iO F θ I ∂t 2χ2r θ 2

– j = 3: Λ+ (−2,1)

=

c2 µc2 2 + V ∂θ − iqOJ ∂θ 2χ3r 4χ2r 1 i ∂θ OB −2 µ + OI OG − c∂r OI + ∂t OI − OF ∂θ V + ∂t . 2 4 2 2c

Λ+ (−1,1)

−

−

In practice, the order of the ABC will be chosen in accordance with the order of accuracy of the interior scheme which is used to solve the TDKGE in Ω.

3 Absorbing Boundary Conditions for the Dirac Equation in 2-d In this section, we are interested in the derivation of ABCs for the 2-d TDDE. The strategy which is used is slightly different from the one used for TDKGE, although some close connections exist and will be discussed later. The idea here is to apply a microlocal diagonalization procedure as usually done for hyperbolic systems (see for instance [16]).

3.1 Dirac equation We denote by PD the 2-dimensional Dirac operator such that PD u = 0 in a bounded domain Ω, and as before we derive a boundary condition on ∂Ω to avoid spurious wave reflections. This operator is given by

where

PD := I4 ∂t + αx c∂x − ieAx + αy c∂y − ieAy + i(Vc + eV )I4 − iβmc2 ,

– Vc : Ω → R is an interaction (Coulomb) potential. – V : Ω × R+ → R is a combination of the self-consistent and external electric potentials. – A : Ω ×R+ → R2 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. – u : R+ × Ω → C4 is the Dirac wavefunction, that is a spinor field initially normalized to 1: ku(0, ·)k 2 4 = 1. L (Ω)

– Hermitian Dirac matrices αx , αy , β are defined by  01 02 σx   , σ = , α =  x x  10 σx 02    02 σy 0 −i αy = , , σy = σ 0 i 0  y 2    I 0    β = 02 −I2 , 2 2

and 02 and I2 are respectively the zero and identity matrices in M2 (C). The following relations hold: α2x = α2y = β 2 = I4 and {αx , αy } = {αx , β} = {αy , β} = 04 . – J = (Jx , Jy ), with Jx = ec u, αx u C4 , Jy = ec u, αy u C4 , denotes the current density and ρ denotes P the particle density which is equal to e 4i=1 |ui |2 .

In the following, we will assume that (A, V ) is given at any time. 17

3.2 Field-particle TDDE in 2-d As before, we consider a circular domain ΩR of radius R. The system is rewritten at the boundary in polar coordinates, where r ∈ [0, ε) with ε > 0: I4 ∂r + L(r, θ, t, ∂θ , ∂t ) u = 0,

(29)

with L = L1 + L0 ∈ OP S 1 and Li is an operator of order i(= 0, 1). We denote R2+ := R∗+ × [0, T ] and + ν(x0 ,y0 ) a neighborhood of (x0 , y0 ). Finally, we introduce ΣTR := ∂ΩR ×[0, T ] and ΣTR+r = ∂ΩR+r ×[0, T ], + + with ΩR+r a disc of radius R + r containing ΩR . We then study the TDDE in the crown ΩR+r − ΩR and design boundary conditions such that incoming waves at r = 0 are set to zero or equivalently to vanish waves that leave the crown surrounding ΩR . Some basic computations lead on ΣTR+r to 1 1 1 ˜ L1 := α ˜ x ∂t + α ˜ x β, ˜xα ˜ y ∂θ , L0 := α c χr c with χr =

1 , R+r 

0 0  0 0 α ˜x =   0 e−iθ eiθ 0

 0 e−iθ eiθ 0  , 0 0  0 0

and 

 0 0 0 −ie−iθ  0  0 ieiθ 0 . α ˜y =  −iθ  0 −ie  0 0 ieiθ 0 0 0 We have β˜ = iβmc2 + i(eV + Vc )I4 − ie(Ar α ˜ x + Aθ α ˜ y ), that is  i Vc + eV + mc2 0 0 −ieAr e−iθ − eAθ e−iθ   0 i Vc + eV + mc2 −ieAr eiθ + eAθ eiθ 0 . β˜ =  −iθ −iθ 2   0 −ieAr e − eAθ e i Vc + eV − mc 0 iθ iθ 2 −ieAr e + eAθ e 0 0 i Vc + eV − mc 

Now since



i 0 α ˜xα ˜y =  0 0

 0 0 0 −i 0 0  0 i 0  0 0 −i

and α ˜ x β˜ is given by  −ieAr + eAθ 0 0 ie−iθ Vc + eV − mc2 iθ 2   0 −ieAr − eAθ 0 ,  ie Vc + eV − mc −iθ 2   0 ie Vc + eV + mc −ieAr + eAθ 0 iθ 2 ie Vc + eV + mc 0 0 −ieAr − eAθ 

18

we get  ω iτ −iθ e 0 0   −χ c r     iτ ω iθ  0 e 0    χr c  L1 = σ(L1 ) =    iτ ω −iθ  0 e − 0    c χr    iτ iθ ω  e 0 0 c χr 

and L˜ := iρI4 + L1 (symbol of I4 ∂r + L1 ). The leading symbol p1 is defined by 2 2 2 ˜ θ, t, ξ, ω) = ρ2 − τ + ω p1 (r, θ, t, ρ, ξ, ω) = det L(r, , 2 2 c χr where  iτ −iθ ω e 0 0   iρ − χ c r     ω iτ iθ   0 iρ + e 0   χ c r ˜ .  L=  iτ −iθ ω   0 e iρ − 0   c χr    iτ iθ ω e 0 0 iρ + c χr 

The symbol p1 has 2 double real roots ρ± = ±

τ 2 ω 2 1/2 − c2 χ2r

in the hyperbolic region of the cotangent bundle, defined by o n ω2 τ 2 < 0 ⊆ T ∗ (ΣTR+r ). − H(ΣTR+r ) = (θ, t, ω, τ ) ∈ T ∗ (ΣTR+r ) : χ2r c2 In H(ΣTR+r ) the outgoing bicharacteristic corresponds to r′ (s) > 0, that is r′ (s) =

∂p1 τ 2 ω 2 = 4ρ ρ2 − 2 − 2 (s). ∂ρ c χr

Therefore, r′ (s) > 0 if and only if one has τ 2 ω 2 τ 2 ω 2 ρ > 0 and ρ2 > 2 − 2 or ρ < 0 and ρ2 < 2 − 2 c χr c χr in a conic neighborhood of H(ΣTR ). We denote by N0 the transition matrix such that L1 = N0 σ(Λ1 )N0−1 where 

iρ−  0 σ(Λ1 ) =  0 0

0 iρ− 0 0 19

0 0 iρ+ 0

 0 0  0  iρ+

and 

 iτ e−iθ iτ e−iθ s s − 0 0   ω2 τ 2  ω2 τ 2 ω ω   c  c − − − + 2 2 2 2   χr c χr χr c χr   iθ iθ   iτ e iτ e (30) N0 =  s 0 − s 0  . ω2 τ 2 2 2  ω ω τ ω    c c − 2+ − 2−   2 2 χr c χr χr c χr     0 1 0 1 1 0 1 0

We denote by V0 the operator Op(N0 ) ∈ OP S 0 , with naturally V0 V0−1 = I + R and R ∈ OP S −∞ . We now apply a technique originally proposed by Taylor, and developed in details by Antoine et al. in [3] in electromagnetism. Starting from

we set w := V0−1 u. The equation becomes

where

∂ + L u = 0, ∂r

∂ + Λ1 w0 = R0 w0 , ∂r R0 := −V0−1 L0 V0 − V0−1

∂V0 . ∂r

From there, Taylor’s method can be implemented. The principle consists of block-diagonalizing the original system up to a certain order. More precisely, we search successive approximations of the solution, in the form of a sequence m−1 X ∂w−m + Λ−j w−m = R−m w−m , ∂r j=−1

where w−m = (I + K−m )w−m+1 (m > 1), Λ−m+1 ∈ OP S −m+1 and R−m ∈ OP S −m . In fine, and according to [3, 31], we have Proposition 31 The solution u to (29), is such that w ∼ V u, where the solution w to ∂ +Λ w ∼0 ∂r

∞ is given by Πp=m+1 (I + K−p )w−m , where σ(Λ− ) 0 σ(Λ) = ∈ M4 (C) 0 σ(Λ+ )

and Λ± are two first-order operators of diagonal principal symbol. In addition, the operator V is defined 1 as Πj=−∞ (I + K−j )V0−1 , for K−j ∈ OP S −j . Sign ∼ is to be understood, as equality up to a smoothing (or regularizing) operator. This is equivalent to say that the operator PD = ∂r +L can in fact be approximately block-diagonalized, or equivalently can be rewritten in the form (∂r + Λ) + R via the ”transition operator” V , and where R, V and Λ are defined explicitly in the proposition above. We now explicitely implement the diagonalizing method, and we first search for an operator K−1 ∈ OP S −1 with K−1 = σ(K−1 ), such that w−1 is defined by w−1 = (I + K−1 )w0 20

and w−1 satisfies the equation ∂w−1 + (Λ0 + Λ1 )w−1 = R−1 w−1 . ∂r According to [3], we have   σ0 (R0 )ij , ∀(i, j) ∈ Jn , (K−1 )ij = λ1,j − λ1,i  0, ∀(i, j) ∈ Jd ,

where

  Jn := (1, 3), (1, 4), (2, 3), (2, 4), (3, 1), (3, 2), (4, 1), (4, 2) ,  Jd := (1, 1), (1, 2), (2, 1), (2, 2), (3, 3), (3, 4), (4, 3), (4, 4) ,

and, denoting R0 = σ0 (R0 ) which represents the order 0 part of R0 (R0 )ij , ∀(i, j) ∈ Jd , (Λ0 )ij = 0, ∀(i, j) ∈ Jn .

By construction, see again [3], we can prove that ( − K−m , Λ1 ] + Op σ−m+1 (R−m+1 ) , ∀(i, j) ∈ Jd , (Λ−m+1 )(i,j) = ij 0, ∀(i, j) ∈ Jn , and

Here,

Pm−1

j=−1

m−1 X ∂w−m + Λ−j w−m = R−m w−m . ∂r j=−1

Λ−j is a block diagonal first-order operator with principal symbol   −1 0 0 0 m−1 X  0 −1 0 0  σ1 Λ−j = λ  , 0 0 1 0 j=−1 0 0 01

where we have set: λ := iλ+ . In fine, ABCs are deduced from − ∂w−m + Λ− w−m = 0, ∂r

where Λ− ∈ M2 (C) −

Λ =

(P

m−1 j=−1

Λ−j 0, otherwise,

kl

, for (k, l) ∈ (1, 1), (1, 2), (2, 1), (2, 2) ,

and ∞ w = Πp=m+1 (I + K−p )w−m ,

where for w = (w1 , w2 , w3 , w4 ) ∈ C4 . At this point we have to evaluate R0 . We again refer to [3] for a detailed description of algebraic symbolic computations. We just apply these formulas skipping the computational details. We have ∂N ∂V0 0 = N0−1 . σ0 V0−1 ∂r ∂r We must estimate σ0 (V0−1 LV0 ) now. We first easily prove that σ0 V0−1 LV0 = N0−1 L0 N0 , 21

where



 1 β 0 α−β 0 β − α   −2iθ   e α  0 0    α−β α−β  =  1  α   0 β−α 0 α−β     β e−2iθ 0 0 β −α β−α

N0−1

and

 iτ e−iθ    s α :=   ω  ω2    − − c  χr χ2r iτ e−iθ    s β :=   ω  ω2    c − +  χ χ2 r

r

τ 2 c2 τ 2 c2

,

.

After some computations and simplifications, we can prove that   2 iθ τ r 0 0 − r − riτ e2 2 2 2 2 2   2c ωχ2 − τc2 2c2 ω + τc2 χωr − ω − τc2 χ2 χ2 r   r r   2 −iθ τ iτ e   r 0 − 0 r r   ω2 2 2 τ2 2 2 ω ω τ τ ω 2 2c − 2 − 2c − − 2   χr c χr χ2 c2 χ2 c2 r r . N0−1 =  2  riτ eiθ τ  0 0 r r   2 2 2 2   2c ωχ22 − τc22 2c2 ω − τc2 χωr − ω − τc2 χ2 χ2   r r r   2 iτ e−iθ τ r r   − r 0 0 2c

ω2 χ2 r

2

− τc2

2c2

ω2 χ2 r

2

+ τc2

ω χr

−

ω2 χ2 r

2

− τc2

Now, by derivation, one gets



αr 0 ∂N0  0 −e2iθ βr = 0 0 ∂r 0 0

 βr 0 0 −e2iθ αr   0 0 0 0

where: αr = ∂r α and βr = ∂r β. This finally gives   αr βr α−β 0 α−β 0     βr αr    0 0   β − α β − α −1 ∂N0 .  N0 =  β α r ∂r r   β−α 0 β −α 0     βr αr  0 0 α−β α−β

This enables us to determine R0 then K−1 . First, let us introduce the following functions  ieAr − eAθ    γ1 := − ,   c   −iθ  ie (Vc + eV − mc2 )   γ2 := , c ieA + eA  r θ   γ := − ,   3 c  −iθ   ie (Vc + eV + mc2 )   γ4 := . c 22

Then, the symbol N0−1 L0 N0 is equal to 

αγ1 +αγ4∗ β+γ2 −βγ3 α−β

  0   γ1 α+γ4 |α|2 +γ3 α−γ2  α−β 0

βγ1 +γ4∗ β 2 +γ2 −βγ3 α−β α∗ γ3 e−2iθ +γ4 |α|2 +γ1 α−γ2∗ e−2iθ α−β

0

0

γ3 β ∗ e−2iθ +γ4 αβ ∗ +γ1 α−γ2∗ e−2iθ α−β

0

0

γ1 β+αβγ4∗ +γ3 α−γ2 α−β

0

0

α∗ βγ4 +e−2iθ α∗ γ3 +γ1 β−γ2∗ e−2iθ β−α

|β|2 γ4 +β ∗ γ3 e−2iθ +γ1 β−γ2∗ e−2iθ β−α



  .  

Let us recall that

K−1

 0 0 σ0 (R0 )13 σ0 (R0 )14 1  0 0 σ0 (R0 ) 23 σ0 (R0 ) 24   . =   σ0 (R0 )31 σ0 (R0 )32 0 0 2λ σ0 (R0 ) 41 σ0 (R0 ) 42 0 0 

We now have to determine an explicit expression for K−1 in order to deduce K−1 , and to derive explicit TBC/ABCs. In this goal, it is necessary to approximate K−1 to provide an expression of K−1 involving simple (differential/pseudodifferentiel) operators. A natural approach consists of considering different kinds of regimes. Typically, we will consider high time-frequency regimes to avoid small oscillations, that is such that |τ | ≫ 1. Naturally, other frequency regimes (on ω ≫ 1, ω/τ ≫ 1, ≪ 1, etc) could be considered from there as long as we remain in the hyperbolic region of the boundary cotangent bundle. The asymptotic expansion for large τ for K−1 ’s components is as follows. A simple way to achieve approximations of K−1 consists of expanding N0−1 and N0 for large τ . First as √ ie−iθ τ τ 2 r2 √ 2 ωc − ω 2 c2 − τ 2 r2 − √ 2 2 ω c − τ 2 r2 ωc − ω 2 c2 − τ 2 r2 , ie−iθ τ ω 2 c2 √ = ωc − √ 2 ω 2 c2 − τ 2 r 2 ωc − ω 2 c2 − τ 2 r2

αr =

then for large τ , αr as well as βr are O(τ −1 ) and α − β is O(1). As a consequence N0−1 ∂r N0 is a O(τ −1 ) and has then no contribution to σ0 (R0 ) which is then equal to σ0 (N0−1 L0 N0 ). We easily see now that s 2ice−iθ ω 2 τ 2 2ie−iθ ωc α−β = − , α + β = , τ χ2r c2 rτ and α2 + β 2 = −

2c2 2ω 2 τ 2 − 2 . τ 2 χ2r c

We remark that for large τ eiθ 1 ω 2 c2 1 =− 1 + 2 2 + O(τ −4 ) α−β 2 2r τ

and α=− c

s τ2 c2

τ e−iθ −

ω2 ω + i χ2r χr

= −e−iθ 1 −

iωc + O(τ −2 ). τ χr

Then, we have 1 α iωc = 1− + O(τ −2 ), α−β 2 τ χr

1 β iωc = 1+ + O(τ −2 ). β−α 2 τ χr 23

Note that

(−1,0) N0



iωc −iθ −e 1 −  τ χr   :=  0 −eiθ   0 1

e−iθ 1 +

0 iωc 1+ τ χr 1 0

0

iωc τ χr

0 1

0

iθ

e



  iωc   1− τ χr   1 0

(31)

and 

(−1,0) N0−1

eiθ − 0  2  −iθ   0 −e  2 :=   eiθ  0  2  −iθ  e 0 2

 1 iωc 1+ 2 τ χr    iωc  0 1−  τ χr . 1 iωc   0 1− 2 τ χr    iωc 0 1+ τ χr 0

1 2 1 2

(32)

In fact, here only the order 0 terms in τ are of interest N0−1 = N0−1 )(0,0) + O(τ −1 ),

(0,0)

N0 = N0

+ O(τ −1 ),

where 

N0−1

(0,0)

eiθ 0 − 2   e−iθ  0 −  2 :=  iθ  e  0  2  −iθ e 0 2

0 1 2 0 1 2

 1 2   −iθ  −e 0   0 (0,0) =  , N0 0 1  1  2 0

 0 e−iθ 0 eiθ 0 −eiθ  . 1 0 1 0 1 0

Then for large τ , we have −1 σ^ 0 (R0 ) = − N0

(0,0)

(0,0)

L0 N0

− N0−1

(0,0) ∂N0 (0,0) ∂r

In addition, we have the following expansion ic iω 2 c3 1 = − − 3 2 + O(τ −4 ). 2λ 2τ 4τ r (−1,0)

So, we approximate K−1 by K−1

(−1,0)

K−1

defined by

 0 0 σ0 (R0 )13 σ0 (R0 )14 ic  0 0 σ0 (R0 ) 23 σ0 (R0 ) 24  . =−    0 0 σ0 (R0 )31 σ0 (R0 )32 2τ 0 0 σ0 (R0 ) 41 σ0 (R0 ) 42 

24

(33)

Since we have    σ0 (N0−1 L0 )11         σ0 (N0−1 L0 )14        σ0 (N0−1 L0 )22        σ0 (N0−1 L0 )23

   σ0 (N0−1 L0 )31        σ0 (N0−1 L0 )34        σ0 (N0−1 L0 )42        σ0 (N −1 L0 )43 0

= = = = = = = =

1 c 1 c 1 c 1 c 1 c 1 c 1 c 1 c

eiθ (iAr e − eAθ ) + eiθ i(Vc + eV + mc2 ) , − i(Vc + eV − mc2 ) − (iAr e + eAθ ) ,

ie−iθ (Vc + eV + mc2 ) + e−iθ (ieAr + eAθ ) , i(Vc + eV − mc2 ) − iAr e + eAθ , − eiθ (iAr e − eAθ ) + eiθ i(Vc + eV + mc2 ) , i(Vc + eV − mc2 ) − iAr e − eAθ , ie−iθ (Vc + eV + mc2 ) − e−iθ (ieAr + eAθ ) , i(Vc + eV − mc2 ) − iAr e + eAθ ,

−1 then we deduce that σ0 (N^ 0 L0 N0 ) is equal to

 −ieAr − i(Vc + eV ) 0 imc2 − eAθ 0 2 2 0 eAθ + imc 0 −i(Vc + eV ) − iAr e  .  −eAθ − imc2 0 −ieAr + i(Vc + eV ) 0 c 0 i(Vc + eV ) − ieAr 0 eAθ − imc2 

We can now define the operators associated to the symbols derived above. They are written as (−1,0) the combination of classical differential and Riemann-Liouville operators. We will denote by K−1 (−1,0) = K−1 + O(τ −2 ). We recall that to any symbol (iτ )−n (for approximate operators such that K−1 n ∈ N∗ ) is associated a Riemann-Liouville operator defined by ∂t−n . (−1,0)

Proposition 32 For large τ and neglecting the terms beyond O(τ −3 ), the operator K−1 imated by

is approx-

˜ (−1,0) = C∂ ˜ t−1 , K −1 where 

 0 0 −imc2 + eAθ 0 0 0 0 i(Vc + eV ) + iAr e   C˜ =  . eAθ + imc2 0 0 0 0 −i(Vc + eV ) + ieAr 0 0 As a consequence, the outflowing part of the wave is selected by annihilating up to a certain order the part of w−m which propagates in the opposite direction of the outward normal vector to the boundary [9]. Theorem 31 The ABCs at order 0 and 1 (corresponding to vanish the outgoing waves from the small crown surrounding ΩR ) are defined as follows. Setting w0 = V0−1 u and w−1 = (I + K−1 )w0 , the ABCs write – At order 0: V0−1 u 1,2 = w0 1,2 = 0 on ΣTR . – At order 1: V−1 u 1,2 = w−1 1,2 = 0 on ΣTR , where V−1 = I + K−1 V0−1 . This simply consists of vanishing the incoming waves (inside ΩR ).

No additional condition (beyond the TDDE itself) on outgoing waves is required. To actually determine an explicit ABC I + K−1 V0−1 u = 0, 1,2

25

it is necessary to evaluate V0−1 and K−1 V0−1 . In this goal and as above, we have to expand N0−1 and K−1 N0−1 for τ large using (32) and an asymptotic expansion. From the above study, we deduce that   eiθ 1 0 0  − 2 2   −iθ   −e 0 e−iθ 0 e−iθ 1    0 − 0 iθ iθ (0,0)   0 −e (0,0) 0 e  2 2  = (34) V0−1 =  iθ , V . 0 1 0 1 1 0  e  0 0  1 0 1 0  2 2  eiθ 1  0 0 2 2 We furthermore have



V0−1

(−1,0)

and

(−1,0) V0

eiθ − 0  2  −iθ   0 −e  2 =  eiθ  0  2   eiθ 0 2

 1 ic −1 0 1 + ∂θ ∂t  2 χr   1 ic −1  0 1 − ∂θ ∂t  2 χr   1 ic −1  0 1 − ∂θ ∂t  2 χr   1 ic −1 0 1 + ∂θ ∂t 2 χr

(35)



ic ic −iθ 1 − ∂θ ∂t−1 0 e−iθ 1 + ∂θ ∂t−1  −e χr χr  ic  −1 iθ = 0 −e 1 + ∂θ ∂t 0 eiθ 1 −  χ r  0 1 0 1 0 1 (−1,0)

Now using, (K−1 )(−1,0) and (V0 )(−1,0) , we get an explicit expression of V−1 (−1,0) V−1

1 = 2

A(−1,0) B (−1,0) C (−1,0) D(−1,0)

0



  ic  ∂θ ∂t−1  (36)  χr  1 0

which is:

(37)

with (−1,0)

A−1

=

1 2

(−1,0) B−1

(−1,0)

C−1

=

1 2

− (imc2 − eAθ )∂t−1 − 1 eiθ 0 , 0 (i(Vc + eV ) + ieAr e)∂t−1 − 1 e−iθ 1 = 2

(−1,0) D−1

0 −(imc2 − eAθ )∂t−1 + 1 −1 (i(Vc + eV ) + ieAr e)∂t + 1 0

,

0 1 − (eAθ + imc2 )∂t−1 eiθ , 0 1 + (i(Vc + eV ) − ieAr )∂t−1 e−iθ

1 = 2

0 1 + (eAθ + imc2 )∂t−1 −1 1 − (i(Vc + eV ) − ieAr )∂t 0

.

The process continues with higher order ABCs but are not computed here. The discretization of this condition will be discussed in Section 5. The ABC which are derived so far depends on the 4 components of the unknown and are independent of the derivative with respect to r. The ABC, V0 u 1,2 = 0, then (V0 + K−1 V0 )u 1,2 = 0 are easily analytically evaluated. 26

3.3 Application to 1-d Dirac equation A direct application of the above conditions in a 1-d framework simply necessitates to vanish angular derivatives and set θ to 0. The domain is denoted by (−a, a) and the equation simply writes ∂r + L(r, θ, t, ∂t ) u = 0, where

1 1 ˜ L1 = α ˜ x ∂t L 0 = α ˜ x β, c c 

0 0 α ˜x =  0 1

and

 001 0 1 0 , 1 0 0 000

β˜ = iβmc2 + i(eV + Vc )I4 − ieAx αx . Some calculations show that α ˜ x β˜ is equal to   −ieAx 0 0 i Vc + eV − mc2 2   0 −ieAx 0 .  i Vc + eV − mc 2   0 i V + eV + mc −ieA 0 c x i Vc + eV + mc2 0 0 −ieAx We deduce that the TDDE in 1-d writes  c∂ u + ∂t u4 − ieAx u1 + i(Vc + eV − mc2 )u4   x 1 c∂x u2 + ∂t u3 − ieAx u2 + i(Vc + eV − mc2 )u3 2   c∂x u3 + ∂t u2 + i(Vc + eV + mc2 )u2 − ieAx u3 c∂x u4 + ∂t u1 + i(Vc + eV + mc )u1 − ieAx u4

= 0, = 0, = 0, = 0.

Since we have

(−1,0)

V0

  −1 0 0 −1 0 1 0 1  0 −1 1  0 −1 0 1  −1 (−1,0) = , V0 =  1 0 0 0 1 0 1 2 0 1 1 1 0 10 

the first right (resp. left) ABC writes

 1 0 , 1 0

u1 − u4 = 0, u2 − u3 = 0, (resp. u1 + u4 = 0, u2 + u3 = 0). (−1,0) equal to Similarly, we derive explicit ABCs from V0−1 + K−1 V0−1 1,2 u = 0, with V−1 

 −imc2 ∂t−1 − 1 0 0 1 − imc2 ∂t−1  1 0 (i(Vc + eV ) + ieAx e)∂t−1 − 1 (i(Vc + eV ) + ieAx e)∂t−1 + 1 0   2 −1 2 −1  .  1 − imc ∂t 0 0 1 + imc ∂t 2 0 1 + (i(Vc + eV ) − ieAx )∂t−1 1 − (i(Vc + eV ) − ieAx )∂t−1 0

The second ABC is written as

and

− imc2 ∂t−1 − 1 u1 + − imc2 ∂t−1 + 1 u4 = 0 (i(Vc + eV ) + ieAx )∂t−1 − 1 u2 + (i(Vc + eV ) + ieAx e)∂t−1 + 1 u3 = 0 . 27

4 Other techniques We now present some alternative techniques to the derivation of ABCs for TDDE or TDKGE.

4.1 Volkov approach We here give some elements of derivation of ABCs based on the Volkov wavefunction in 3-d. The principle is based on the fact that the TDDE has explicit solution in R3 called the Volkov wavefunction, when the interaction potential is set to zero. The approach is very similar to the one proposed in [27], for field-molecule time-dependent Schr¨odinger equations (TDSE). Roughly speaking the idea is as follows. We first rewrite the TDDE in the form   ∂t u = (P1 + P2 )u on R3 × [0, T ], (38)  u(·, 0) = u (·) on R2 , i where in 3-d

and

P1 = −αx c∂x − ieAx − αy c∂y − ieAy − αz c∂z − ieAz − iβmc2 , P2 = −(iVc + eV )I4 , αz =

02 σz σz 02

σz =

0 −1 1 0

with

.

We introduce a fictitious domain Ω such that supp(u0 ) ⊆ Ω and supp(Vc ) ⊆ Ω. As usual, we determine an operator B  ˜ = (P1 + P2 )˜ ∂t u u on Ω × [0, T ],     ˜ (·, 0) = u0 (·) on Ω, u     ˜ = 0 on ∂Ω × [0, T ], B·u

˜ . We remark that in Ω c ⊆ R3 , the Dirac equation writes such that u|Ω ∼ u ∂t u = P1 u.

(39)

It turns out that this equation has an exact solution in R3 , see [28], that we call uv . Denoting and assuming that – the incoming wave is a plane electromagnetic wave, that is V = 0. Its amplitude is denoted by A0 . – A := (Ax , Ay , Az )T is supposed to be in the form A(x, t) = A0 f (ωt − k · r), with x = (x, y, z), f is a given function, ω (resp. k) the imposed incoming pulse frequency (resp. 3-d photon momentum). – A := (0, AT ), k := (ω, k), n := k/kkk, e := A/kAk2 and α = (αx , αy , αz ). Volkov wavefunction uv is searched in the form uv (x, t) = e−i(p0 t−p·x) φ(τ ), where p := (p0 , p) is a constant four-components vector chosen such that p20 − kpk2 = m2 and τ = k·r t− . Now, according to [28], we have the following proposition. ω 28

Proposition 41 For c set to 1 and for c an arbitrary constant bispinor, a solution to (38) is given by h uv (x, t) = 1 +

i c eαs A √ eiS(τ ) , 2(p · k) 2p0

where the phase function is given by S(τ ) = −p0 t + p · x +

Z

τ

h e(p · A) p·k

−

e2 kAk2 i dτ, 2(p · k)

τ =t−

k·r , ω

with αs = αx + αy ey + αz ez +

k α × α · (n × e), 2ω

p · k = p0 ω − p · k.

In fine, we should sum this solution over all the momenta (p0 , p) which naturally leads a huge computational complexity. An important effort should be done on developing efficient techniques to evaluate this Volkov wavefunction for the field-particle TDDE. Note that, it is naturally possible to select certain frequency regimes in order to reduce the computational cost due to the overall sum. Beyond, the computational difficulty to numerically evaluate this function, it is unfortunately not possible to directly impose uv as boundary condition, due to the fact this vectorial function is defined as an integral in R3 . In fact, from the mathematical view point, this approach is not totally rigorous as the Volkov wave function is solution to an IVP in R3 and not from an IBVP in Ω c as used in practice. As a consequence, to be valid, the Coulomb potential has to be zero everywhere. In that case (Vc ≡ 0), imposing B · u = uv on ∂Ω would constitute a TBC. In the general case (Vc 6= 0), it is then necessary to include the Coulomb potential in the expression of uv to obtain a relevant ABC. A natural idea consists of splitting the equation as follows. Assuming that the solution is known at time tn , and for some small ∆tn > 0, the solution to  ˜  ∂t u = (P1 + P2 )˜ u for (tn , tn + ∆tn ], u ˜ (·, tn ) = u ˜ 0 (·),

can be approximated, using the Trotter-Kato formula [10], by the solution of  ˜ ˜ for (tn , tn + ∆tn ], ∂t u = P1 u        ˜ (·, tn ) = u ˜ n (·), u   ˜ ˜ for (tn , tn + ∆tn ], ∂t u = P2 u       ˜ n (·). ˜ (·, tn ) = e∆tn P1 u u

(40)

˜ (·, tn ) − e∆tn P1 e∆tn P2 ) u ˜ (·, tn )kH In the appropriate Hilbert space H, we can prove that ke∆tn (P1 +P2 ) u 2 is a O ∆tn [P1 , P2 ] . This gives an approximate Volkov-Coulomb wavefunction for TDDE denoted by Pn un+1 vc . As a consequence, the new mixed problem consists of solving, for all n such that i=1 ∆ti 6 T ,  ˜ ∂t u = (P1 + P2 )˜ u on Ω × (tn , tn+1 ],     ˜ (·, tn ) = un (·) on Ω, u     ˜ u = unvc on ∂Ω × (tn , tn+1 ], where at time tn , the numerical solution is denoted un+1 and un+1 is solution to (40). Higher order vc ABCs are naturally possible. The main interest of this approach comes from the fact that both ˜ = P1 u ˜ , ∂t u ˜ = P2 u ˜, ∂t u 29

can be solved “analytically” (the equation system is a trivial diagonal system). As a consequence, from the analytical point of view, it is possible to derive very accurate ABCs (simply by increasing the order of the operator splitting). However, we face a similar problem as the one for TDSE. To have a small ˜ = P1 u ˜ in error, ∆tn has to be chosen small. The issue comes from the fact that the solution to ∂t u (40) possesses an oscillatory spatial integral of frequency 1/∆tn . As a consequence, a fine numerical computation of this integral, requiring ∆tn ≪ 1, can be computationally costly. The approach which was proposed in [27] for TDSE is still applicable here. This consists in using Filon’s approximation of highly oscillatory integrals (typically the convergence is a positive power of the time step [22, 23]). This approach for TDDE, will be implemented and studied in a forthcoming paper, devoted to accurate numerical discretizations of ABCs for TDDE and TDKGE. 4.2 Basic transformations on TDKGE and TDDE In the last 2 sections we have derived ABC for TDKGE and TDDE using two distinct but close techniques. This section is devoted to some simple transformations on TDKGE. We show that these transformations allow to rewrite these equations in such form that existing techniques can (almost directly) be applied to derive ABCs. The presentation will be field-free, although most of these ideas are applicable with field-particle equations. Let us consider the system  PKG,D ψ(x, t) = 0,     BKG,D ψ(x, t) = 0,     ψ(x, 0) = ψ0 (x) ∈ C (or C4 ), with x = (x, y) and

  PKG = c2 ∂x2 + ∂y2 + ∂z2 − ∂t2 − m2 c4 , 

PD = iI4 ∂t + icαx ∂x + icαy ∂y + icαz ∂z − βmc2 .

Let us now explain the different transformation that can be used – TDKGE to reaction-diffusion equation. We set χ1 ϕ χ := = . χ2 ∂t ϕ In that case, it is easy to show that χ satisfies the following equation ∂t χ = A(∂x2 + ∂y2 + ∂z2 )χ + Bχ, where A=

0 0 c2 0

,

B=

0 1 −m2 c4 0

.

Reaction-diffusion problems are naturally very studied in the literature (in chemistry, fluid dynamics, etc) in particular from a boundary condition point of view [18]. – TDKGE to TDDE. There exists a fundamental connection between field-free TDKGE and TDDE. Indeed, we easily check that i 1 (∂x2 + ∂y2 + ∂z2 − 2 ∂t2 + m2 c4 )I4 = αx ∂x + αy ∂y + αz ∂z + I4 ∂t + mc2 β c c (41) i 2 x αx ∂x + αy ∂y + αz ∂z + I4 ∂t + mc β . c As a consequence, the derivation of TBC/ABCs from TDKGE (resp. TDDE) can be useful to derive TBC/ABCs to TDDE (resp. TDKGE). We roughly describe a possible approach. Assuming that ABCs have been derived for TDDE. The factorization (41), consists of writing TDKGE as I4 PKG = PD PD∗ . 30

Now, we have seen that PD can be block-diagonalized, which formally writes PD = V −1 (I4 ∂r − Λ)V + R, where the transition operator V can be explicitly constructed, and R ∈ OP S −∞ . Similarly, we have seen that PKG can be factorized as follows I4 PKG = I4 (∂r + Λ+ )(∂r − Λ− ) + R . As a consequence we can formally state that ∗ (∂r + Λ+ )(∂r − Λ− )I4 = V −1 (I4 ∂r − Λ)V V −1 (I4 ∂r − Λ)V + R,

where Λ and Λ± are first order operators and V is an operator of leading order 0. From there, it is possible to find connections between these three operators. We do not go further in this direction, but we think that this could be an interesting question to study.

5 Numerical Simulations This section is devoted to some basic 1-d numerical illustrations of the boundary conditions derived above. More advanced numerical simulations, as well as derivation of accurate discretization of 2-d boundary conditions and analysis of the overall schemes will be presented in a forthcoming paper.

5.1 Discretization for TDKGE The discretizations which are proposed in this paper are relatively naive but still reasonablyP accurate. We respectively denote by ∆x the constant space step and ∆tn the time step at time tn = ni=1 ∆ti . As usual ϕnj denotes an approximation of the exact solution ϕ, in (j∆x, tn ), for j ∈ Z and n ∈ N∗ . The initial data is a wavepacket defined by ϕ0 (x) =

ck0 1 p + ik x , exp − 0 4x2 mc2 + m2 c4 + c2 k02

where the wavenumber k0 = 5. In the following, we impose ~ = c = m = e = 1. Case 1 We first consider the field-free Klein-Gordon equation on a bounded domain [−a, a], approximated by an explicit scheme 2 c n−1 n 2 2 4 n n n n ϕn+1 − 2ϕ + ϕ = ∆t − m c ϕ ϕ − 2ϕ + ϕ j j . j j−1 j j ∆x2 j+1 As seen above, an ABC at +a is given by

m2 c4 c∂x ϕ(x, t) + ∂t ϕ(x, t) + 2

Z

t

ϕ(x, s)ds = 0

(42)

0

which is an improvement of the simple condition: ∂t + c∂x ϕ(x, t) = 0 .

(43)

Above and in the following j ∈ {1, · · · , N }, with x1 = −a and xN = +a. We numerically compare the following conditions: – Dirichlet: ϕ(±a, t) = 0, for all t; + – transport-like condition (43) (whose symbol includes λ+ 1 and λ0 ); + + – improved transport-like condition (42) (whose symbol includes λ+ 1 , λ0 and λ−1 ); 31

Dirichlet − Transport 500

2

450 0

400 350

−2

space

300 −4

250 200

−6

150 100

−8

50 50

100

150

200

250 time iteration

300

350

400

450

−10

Fig. 1 Klein-Gordon with Dirichlet (left) and transport (right) boundary conditions

Dirichlet − Improved transport 500

2

450 0

400 350

−2

space

300 −4

250 200

−6

150 100

−8

50 50

100

150

200

250 time iterations

300

350

400

450

−10

Fig. 2 Klein-Gordon with Dirichlet (left) and improved transport (right) boundary conditions

Dirichlet − Empirical 500

2

450 0

400 350

−2

space

300 −4

250 200

−6

150 100

−8

50 50

100

150

200

250 time iterations

300

350

400

450

−10

Fig. 3 Klein-Gordon with Dirichlet (left) and empirical (right) boundary conditions

We represent, as usual, the quantity t, x, log u(x, t) which enlights the reflections at the domain boundary. In the tests, we choose ∆x = ∆tn = ∆t = 0.04 for all n, T = 18 and a = 10. At the left boundary is implemented Dirichlet’s boundary condition. At the right boundary, we implement the transport-like condition Fig. 1, and improved transport-like condition Fig. 2. As expected, the improved transport-like condition (42) derived in this paper is demonstrated to be better than (43) and Dirichlet. Case 2 In the second example, we consider a field-particle Klein-Gordon equation. This time, the scheme (and 32

the derived ABC) is more complex and writes 2 Anj c n+1 n−1 n 2 n n n ϕj − 2ϕj + ϕj = ∆t ϕjn+1 − ϕn−1 ϕ − 2ϕ + ϕ + j+1 j j−1 j ∆x2 2∆t Bjn n ϕj+1 − ϕnj−1 + Cjn ϕnj , + 2∆x where

 n Aj = −ieV (xj , tn ),     Bjn = −ieAx (xj , tn ),     n Cj = −ie∂t V (xj , tn ) − ic2 ∂x Ax (xj , tn ) + e2 V 2 (xj , tn ) − c2 e2 A2x (xj , tn ) − m2 c4 ,

with boundary conditions imposed at x = ±a, according to (18) and (20). At −a, we impose a Dirichlet condition, and at +a, we impose i) a transport-like condition (operator symbol includes only λ+ 1 ) (see Fig. 4), ii) a transport-like condition including an order 0 operator (operator symbol includes λ+ 1 and λ+ ) (see Fig. 5) iii) a transport-like condition including an order 0 and order −1 operators (operator 0 + + symbol includes λ+ 1 , λ0 and λ−1 ) (see Fig. 6). We assume that A(t) = A0 cos(ωt), V (t) = E0 sin(ωt) with A0 = 1, E0 = 0.1, and ω = 1. Again, we fix: a = 10, ∆x = ∆t = 0.4. Dirichlet−Transport 500

0

450 400

−2

space

350 −4

300 250

−6

200 150

−8

100 50 50

100

150

200 time iteration

250

300

350

400

−10

Fig. 4 Field-Klein-Gordon with Dirichlet (left) and transport-like (right) boundary conditions

Dirichlet−Improved transport − 0 500

0

450 400

−2

space

350 −4

300 250

−6

200 150

−8

100 50 50

100

150

200 time iteration

250

300

350

400

−10

Fig. 5 Field-Klein-Gordon with Dirichlet (left) and improved transport-like at order 0 (right) boundary conditions

As expected the absorption is improved by including additional operators in the ABC. Stability for the interior scheme: field-free-TDKGE. The stability analysis for the interior scheme in the field-free case is standard. Denoting by g the amplificator factor, and using the usual notations [30], we get for all θ ∈ [−π, π) g 2 (θ) − 2 1 − β(θ) g(θ) + 1 = 0, 33

Dirichlet−Improved transport − 1 500

0

450 400

−2

space

350 −4

300 250

−6

200 150

−8

100 50 50

100

150

200 time iteration

250

300

350

400

−10

Fig. 6 Field-Klein-Gordon with Dirichlet (left) and improved transport-like at order −1 (right) boundary conditions

with β(θ) = 2c2

∆t2 ∆t2 m2 c4 sin2 (θ/2) − . 2 ∆x 2

Stability is ensured when the roots g1 , g2 of this characteristic equations are less or equal to 1, for all θ ∈ [−π, π). It is easy to prove that the scheme is stable provided that β(θ) 6 2 for all θ, that is ∆t satisfies ∆t 6 r c

∆x m2 c2 2 ∆x 1+ 4

.

Stability for the interior scheme: field-TDKGE. In that case, we expect the time step to be time-dependent, as the equation involves a time-dependent electric field. Assuming, to simplify the analysis, that ∂x Ax = ∂x V = 0, the characteristic equation in g writes: α ¯n = 0, g 2 (θ) − 2αn 1 − βn (θ) g(θ) + αn

with An = Anj , B n = Bjn and C n = Cjn for all j (as the electric field is assumed to be space independent) ∆tn An −1 αn = 1 + 2

and βn (θ) =

2c2 ∆t2n ∆x2 (−ie∂t V (tn ) + e2 V 2 (tn ) − c2 e2 A2x (tn ) − m2 c4 ) 2 ∆x2 4c ∆xeA (t ) sin(θ) x n + sin2 (θ/2) − . 2c2

Again the stability condition can be written as a function of βn using α ¯ n 2 |g1 (θ)g2 (θ)|2 = = 1 αn

(44)

and g1 (θ) + g2 (θ) = 2αn (1 − β(θ)). According to (44), the product of the root modulus is equal to 1. The roots are given by i p h g1,2 (θ) = 1 − βn (θ) αn ± αn (αn − 1)

and must satisfy |g1,2 (θ)| 6 1 for all θ. Numerically, a stability condition can easily be determined at (int) each time step, ∆tn 6 ∆tn . For instance, assuming V = 0, then αn = 1 and the scheme is stable provided βn 6 2, for all θ, that is: βn (θ) 6 |βn (θ)| 6

∆x2 (c2 e2 (Anx )2 + m2 c4 ) ∆xeAnx 2c2 ∆t2n 1 + 6 2. + ∆x2 4c2 2c2 34

We conclude that in that simplified situation ℓ2 -stability is ensured provided that ∆tn 6 ∆t(int) := n

∆x 1 r e2 (An )2 m2 c2 . n c eAx x + 1 + 2 ∆x + ∆x2 2c 4 4

Discretization of the boundary conditions. The discretization which is proposed is done according to the boundary conditions derived in Proposition 21: ϕn+1 = ϕnN − N

Pn ce c∆tn n n k (ϕ − ϕnN −1 ) − i∆tn (VNn − Anx,N )ϕnN − ∆tn KN k=0 ∆tk ϕN . ∆x N 2

We can rewrite this equation

n n n ϕn+1 = LnN ϕnN + MN ϕN −1 − ∆tn KN N

where we have set

Pn−1 k=0

∆tk ϕkN ,

 1 n  n n n n   2O − O − 2iO + 2iO K := D,N F,N Gx ,N Gt ,N , N   4  c∆tn ce n LnN := 1 − − ∆t2n KN − i∆tn (VNn − Anx,N ),  ∆x 2    c∆tn  n  MN := . ∆x

As for all n > 1, ϕnN does not appear in the interior scheme, the stability analysis can be done independently at the boundary. We need to show that 2 ∆x|ϕn+1 N | 6 ∆x (a)

N X j=1

|ϕnj |2 ,

(a)

which is trivially satisfied if ∆tn 6 ∆tn for ∆tn ensuring simultaneously that n n |KN | 6 1, |LnN | 6 1, |MN | 6 1. (−a)

Similarly, the scheme is stable for ∆tn 6 ∆tn

.

Remark 51 Boundary conditions can also be derived following Section 2.1.2 and approximated in the spirit of [8] using Pad´e’s approximants. Set Z := −∂t2 + e2 V 2 − c2 e2 A2x − m2 c4 . The ABC can formally be written at +a √ ceAx c∂x − i Z − i . 2 Pad´e’s approximants for k = 1, · · · , m for m > 1, write −1 m (2k − 1)π iα/2 2 (2k − 1)π , dk = eiα tan2 , am = e m cos k 4m 4m

α ∈ R,

and as a consequence, we have the classical approximation:

m m X X √ m m −1 Z∼ am − am . k k dk Z + dk k=1

k=1

−1 n n ϕn , that is Z + dm Set vkn := Z + dm k vk = ϕ . Then, at x = xN = a with N ∈ N such that k a = N ∆x c∂x ϕn+1 −i N

m X

k=1

ak ϕn+1 + N

m

m

k=1

k=1

ceAnx n i X m m n+1 i X m m n ϕ = 0, ak dk vk + ak dk vk − i 2 2 2 N 35

where vkn+1

∆t2n n−1 2 n+1 n n − v Mkn − dm ) − ieV 1 − ieV vkn 2 − k N N − ∆tn ϕN k 2 = ∆t2 1 + 2 n (Mkn − dm k )

n and MN = m2 c4 − e2 (VNn )2 + e2 c2 (Anx,N )2 . The quantities Anx,N and VNn denote Ax (a, tn ) and V (a, tn ). n+1 Finally, we approximate ∂x ϕn+1 at a by (ϕn+1 − ϕn+1 N N −1 )/∆x, where ϕN −1 is computed by the interior scheme. In order to illustrate this approach we have solved the field-free TDKGE on (−10, 10) with c = m = 1, see Fig. 7. The following numerical data are chosen: N = 400 grid points, m = 100 (for which convergence is reached) and α = −π/4.

Conclusion: stability of the numerical scheme. We can conclude that, from the interior scheme (int) analysis, for ∆tn 6 ∆tn , ∆x

N −1 X j=2

|ϕn+1 |2 6 ∆x j

N −1 X j=2

|ϕnj |2

(a)

(−a)

and from the exterior scheme analysis that, for ∆tn 6 ∆tn and ∆tn 6 ∆tn 2 ∆x|ϕn+1 N |

6

(int)

, ∆tn

Then for ∆tn 6 min ∆tn

(−a)

, ∆tn

(a)

∆x

∆x|ϕn+1 |2 1

∆x|ϕnN |2 , N X j=1

6

∆x|ϕn1 |2

,

.

the scheme is ℓ2 −stable, that is

|ϕn+1 |2 6 ∆x j

N X j=1

|ϕnj |2 .

Naturally, an implicit scheme would require a more technical analysis. Comparison Pade to Dirichlet 400

0

−1 350 −2 300 −3 250

space

−4

−5

200

−6 150 −7 100 −8 50 −9

50

100

150 time iteration

200

250

300

−10

Fig. 7 Klein-Gordon with Dirichlet (left) and Pad´ e (right) boundary conditions

5.2 Discretization for TDDE In this section, we propose a simple numerical illustration of the derived boundary conditions for the one-dimensional TDDE. The quantity which is represented here is again (t, x, log |ψ1 (x, t)|) which enlights the reflections at the domain boundary. We first rewrite PD in the form PD = i∂t − iA∂x + Bmc2 , where A=

c 0 0 −c

,

B=i 36

0 mc2 −mc2 0

.

A simple and natural condition to impose at ±a is (∂t ±c∂x )ψ1,2 = 0. This boundary condition is called a transport-like boundary condition, and corresponds in fact to the order 0 condition in Proposition (31). We can easily check that for a CFL= 1 upwind scheme, Dirichlet’s boundary condition (solution set to zero at the boundary) is equivalent to impose this transport-like condition. Initially, we set: x2 ψ1 (x, 0) = exp − 2 exp(ik0 x), δ

ψ2 (x, 0) =

x2 ck0 p exp − 2 exp(ik0 x), δ mc2 + m2 c4 + c2 k02

where c = 1, k0 = −50, δ = 0.5, m = 1. The numerical domain is (−a, a) with now a = 5. The final time is T = 10 and ∆t = ∆x = 1/60. One-dimensional ABC can easily be deduced from Proposition 32. The numerical scheme we consider is as follows. We denote by (φnj , ψjn ) and approximation of the exact two-spinors (ψ1 (xj , tn ), ψ2 (xj , tn )) at (xj , tn ). The numerical scheme writes  ∆t n   φn+1 (φ − φnj ) + ∆tmc2 ψjn , = φnj + c j ∆x j+1   ψ n+1 = ψ n − c ∆t (ψ n − ψ n ) − ∆tmc2 φn . j j−1 j j ∆x j

Stability for the interior scheme: field-free TDDE. One denotes g1 (θ), g2 (θ) the amplification factors for φ and ψ, with θ ∈ [−π, π), and ξ∆x, where ξ is the dual variable to x. We have  ∆t iθ   g1n+1 (θ) = 1 + c (e + 1) g1n (θ) + ∆tmc2 g2n (θ), ∆x   g n+1 (θ) = 1 − c ∆t (1 − e−iθ ) g n (θ) − ∆tmc2 g n (θ) . 2 1 2 ∆x

A necessary and sufficient condition for stability is |g1 (θ)|, |g2 (θ)| are less than 1 for any θ ∈ [−π, π). That is g1,2 are roots of the following equation with ν = c∆t/∆x, g 2 (θ) − 2 1 − 2ν sin2 (θ/2) g(θ) − 4ν sin2 (θ/2) + 4ν 2 sin2 (θ/2) + 1 + ∆t2 m2 c4 = 0 .

We deduce g1 (θ)g2 (θ) = 1 + ∆t2 m2 c4 − 4ν sin2 (θ/2)(1 − ν) and g1 (θ) + g2 (θ) = 2(1 − 2ν sin2 (θ/2). To ensure the stability, we naturally need |g1 (θ)g2 (θ)| 6 1 for all θ ∈ [−π, π). Then, a sufficient condition for ℓ2 −stability writes c

∆t 6 1. ∆x

Stability at the boundary: field-free TDDE. The implementation of the above boundary conditions for TDDE is straightforward. We can indeed rearrange the scheme. For j 6= 0 and j 6= N  n n n n n ψ1,j+1 + ψ2,j+1 + ψ1,j−1 − ψ2,j−1 ψ1,j  n+1  ψ1,j = + ∆tmc2 2 n n n n n   ψ n+1 = ψ1,j+1 + ψ2,j+1 − ψ1,j−1 + ψ2,j−1 − ∆tmc2 ψ1,j 2,j 2

n + ψ2,j , 2 n + ψ2,j . 2

n n n n At order 0, we impose ψ1,0 − ψ2,0 = 0 at the left boundary, and ψ1,N + ψ2,N = 0 at the right one. Remark that for the very particular interior scheme which is considered here, Dirichlet’s boundary conditions are equivalent to this order 0 condition (transport-like). At the next order, −1, a RiemannLiouville integral is added to the boundary condition. Stability is again trivially satisfied under CFL condition. Figures 8, 9 show the propagation of |ϕ1,2 | imposing order 0 (left) and −1 (right) boundary conditions.

37

Transport (left&right): ψ

1

600

−3

500

−4 −5

space

400 −6 300 −7 200

−8

100

−9

50

100

150

200

250

300

350

400

450

−10

time iteration Transport (left) − Improved transport (right): ψ

1

600

−3

500

−4 −5

space

400 −6 300 −7 200

−8

100

−9

50

100

150

200

250

300

350

400

450

−10

time iteration

Fig. 8 ABC for log(ψ1 ) a) order 0, b) order −1

Transport (left&right): ψ

2

600

−3

500

−4 −5

space

400 −6 300 −7 200

−8

100

−9

50

100

150

200

250

300

350

400

450

−10

time iteration Transport (left) − Improved transport (right): ψ

2

600

−3

500

−4 −5

space

400 −6 300 −7 200

−8

100

−9

50

100

150

200

250

300

350

400

450

−10

time iteration

Fig. 9 ABC for log(ψ2 ) a) order 0, b) order −1

6 Conclusion This aim of this paper was to rigorously derive absorbing boundary conditions for quantum relativistic equations. Both Klein-Gordon and Dirac equations were considered for quantum particles subject to classical electromagnetic fields. Using usual microlocal analysis tools, sequences of more and more accurate absorbing boundary conditions were derived. Simple numerical discretizations and simulations were proposed and analyzed to illustrate the increasing accuracy of these ABCs. Several analytical questions still remain to be addressed: – The well-posedness of the mixed problem  PD,KG wD,KG = 0 on ΩD,KG × [0, T ],     BD,KG wD,KG = 0 on ∂ΩD,KG × [0, T ],     wD,KG (·, 0) = u0D,KG on ΩD,KG , 38

where the index D stands for Dirac and KG for Klein-Gordon, is naturally the very first important question. – The error analysis to estimate kuD,KG −wD,KG kH(ΩD,KG ) in the appropriate Hilbert space H(ΩD,KG ), where uD,KG satisfies  on R2 × [0, T ],  PD,KG uD,KG = 0 u

D,KG (·, 0)

= u0D,KG ,

is also a fundamental question from a theoretical as well as practical point of view. These questions were for instance treated in [9] for Maxwell’s equations. From a computational point of view, several problems have also to be considered – Although it is easy to derive accurate numerical schemes for the TDKGE and TDDE on unbounded domains, the inclusion of boundary conditions, especially for spectral, finite element or implicit finite volume and difference schemes which necessitate the solving of linear systems, may introduce some computational difficulties. As an example, the numerical approximation of TDSE with absorbing boundary conditions derived in a similar way, necessitates a very subtle discretization to maintain the unconditional ℓ2 -stability of a simple Crank-Nicolson type scheme in 1-d [7]. In addition, the conditioning of matrices involved in the linear systems may increase [24]. – A key problem, which is naturally coupled with the ones described above, is the derivation of accurate and efficient discretization of the pseudo-differential operators involved in the derived ABC. Some of them are nonlocal in time and in space such as the Riemann-Liouville operators (21) which complicates their approximation. – The implementation of these conditions in numerical codes for the computation of multidimensional TDKGE and TDDE is naturally the final objective. For instance, production particle/anti-particle paires using intense and short laser pulses is a very active research area [14, 15]. However, the corresponding simulations, based on laser-particle TDDE, necessitate huge computational efforts [12, 13]. Accurate and efficient ABCs could reduce drastically the computational cost of these simulations. These questions will be treated in forthcoming papers and will justify the use of such absorbing boundary conditions. References 1. S. Alinhac and P. G´ erard. Pseudo-differential operators and the Nash-Moser theorem, volume 82 of Graduate Studies in Mathematics. American Mathematical Society, Providence, RI, 2007. Translated from the 1991 French original by Stephen S. Wilson. 2. X. Antoine, A. Arnold, C. Besse, M. Ehrhardt, and A. Schaedle. A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schr¨ odinger equations. Communications in Computational Physics, 4(4):729–796, 2008. 3. X. Antoine and H. Barucq. Microlocal diagonalization of strictly hyperbolic pseudodifferential systems and application to the design of radiation conditions in electromagnetism. SIAM J. Appl. Math., 61(6):1877–1905 (electronic), 2001. 4. X. Antoine and H. Barucq. On the construction of approximate boundary conditions for solving the interior problem of the acoustic scattering transmission problem. In Domain decomposition methods in science and engineering, volume 40 of Lect. Notes Comput. Sci. Eng., pages 133–140. Springer, Berlin, 2005. 5. X. Antoine, H. Barucq, and A. Bendali. Bayliss-Turkel-like radiation conditions on surfaces of arbitrary shape. J. Math. Anal. Appl., 229(1):184–211, 1999. 6. X. Antoine and C. Besse. Construction, structure and asymptotic approximations of a microdifferential transparent boundary condition for the linear Schr¨ odinger equation. J. Math. Pures Appl. (9), 80(7):701– 738, 2001. 7. X. Antoine and C. Besse. Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schr¨ odinger equation. J. Comput. Phys., 188(1):157–175, 2003. 8. X. Antoine, C. Besse, and P. Klein. Absorbing boundary conditions for the one-dimensional Schr¨ odinger equation with an exterior repulsive potential. J. Comput. Phys., 228(2):312–335, 2009. 9. H. Barucq and M. Fontes. Well-posedness and exponential stability of Maxwell-like systems coupled with strongly absorbing layers. J. Math. Pures Appl. (9), 87(3), 2007. 10. R. Dautray and J.-L. Lions. Analyse math´ ematique et calcul num´ erique pour les sciences et les techniques. Vol. 5. INSTN: Collection Enseignement. [INSTN: Teaching Collection]. Masson, Paris, 1988. Spectre des op´ erateurs. [The operator spectrum], With the collaboration of Michel Artola, Michel Cessenat, Jean Michel Combes and Bruno Scheurer, Reprinted from the 1984 edition. 11. B. Engquist and A. Majda. Absorbing boundary conditions for the numerical simulation of waves. Math. Comp., 31(139):629–651, 1977.

39

12. F. Fillion-Gourdeau, E. L. De La Grandmaison, and A. D. Bandrauk. Relativistic ground state of diatomic molecules from the numerical solution of the Dirac equation on parallel computers. Journal of Physics: Conference Series, 341(1), 2012. 13. F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk. Numerical solution of the time-dependent Dirac equation in coordinate space without fermion-doubling. Comput. Phys. Commun., 183(7):1403–1415, 2012. 14. F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk. Relativistic stark resonances in a simple exactly soluble model for a diatomic molecule. Journal of Physics A: Mathematical and Theoretical, 45(21), 2012. 15. F. Fillion-Gourdeau, E. Lorin, and A. D. Bandrauk. Resonantly enhanced pair production in a simple diatomic model. Physical Review Letters, 110(1), 2013. 16. J.-M. Ghidaglia and F. Pascal. The normal flux method at the boundary for multidimensional finite volume approximations in CFD. Eur. J. Mech. B Fluids, 24(1):1–17, 2005. 17. L. Halpern and J. Rauch. Error analysis for absorbing boundary conditions. Numer. Math., 51(4):459–467, 1987. 18. L. Halpern and J. Rauch. Absorbing boundary conditions for diffusion equations. Numer. Math., 71(2):185– 224, 1995. 19. R. Hammer and W. Ptz. Staggered grid leap-frog scheme for the (2 + 1) d dirac equation. Computer Physics Communications, 2013. cited By (since 1996)0; Article in Press. 20. L. H¨ ormander. Linear partial differential operators. Springer Verlag, Berlin, 1976. 21. L. H¨ ormander. The analysis of linear partial differential operators. III. Classics in Mathematics. Springer, Berlin, 2007. Pseudo-differential operators. 22. A. Iserles. On the numerical quadrature of highly-oscillating integrals. I. Fourier transforms. IMA J. Numer. Anal., 24(3):365–391, 2004. 23. A. Iserles. On the numerical quadrature of highly-oscillating integrals. II. Irregular oscillators. IMA J. Numer. Anal., 25(1):25–44, 2005. 24. R. Kechroud, X. Antoine, and A. Soula¨ımani. Numerical accuracy of a Pad´ e-type non-reflecting boundary condition for the finite element solution of acoustic scattering problems at high-frequency. Internat. J. Numer. Methods Engrg., 64(10):1275–1302, 2005. 25. R. Lascar. Propagation des singularit´ es des solutions d’´ equations pseudo-diff´ erentielles quasi homog` enes. Ann. Inst. Fourier (Grenoble), 27(2):vii–viii, 79–123, 1977. 26. E. Lorin and A. Bandrauk. A simple and accurate mixed P 0 -Q1 solver for the Maxwell-Dirac equations. Nonlinear Anal. Real World Appl., 12(1):190–202, 2011. 27. E. Lorin, S. Chelkowski, and A. D. Bandrauk. Mathematical modeling of boundary conditions for lasermolecule time-dependent Schr¨ odinger equations and some aspects of their numerical computation—onedimensional case. Numer. Methods Partial Differential Equations, 25(1):110–136, 2009. 28. J.T. Mendon ¸ca and A. Serbeto. Volkov solutions for relativistic quantum plasmas. Physical Review E Statistical, Nonlinear, and Soft Matter Physics, 83(2), 2011. 29. L. Nirenberg. Lectures on linear partial differential equations. American Mathematical Society, Providence, R.I., 1973. 30. J. C. Strikwerda. Finite difference schemes and partial differential equations. Society for Industrial and Applied Mathematics (SIAM), Philadelphia, PA, second edition, 2004. 31. M. E. Taylor. Partial differential equations III. Nonlinear equations, volume 117 of Applied Mathematical Sciences. Springer, New York, second edition, 2011. 32. B. Thaller. The Dirac equation. Texts and Monographs in Physics. Springer-Verlag, Berlin, 1992.

40

Recommend Documents

Adaptive absorbing boundary conditions for Schrödinger-type ...

Absorbing boundary conditions for the two-dimensional Schrödinger ...

absorbing boundary conditions for the schr odinger equation - CiteSeerX

Absorbing boundary conditions for nonlinear SchrÃ¶dinger equations