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

E Lorin
email: elorin@math.carleton.ca

Computational performance of simple and efficient sequential and parallel Dirac equation solvers

Keywords: Dirac equation, pseudo-spectral method, high performance computing, wave equation

 for solving the Time Dependent Dirac Equation (TDDE) modeling the interaction of classical electromagnetic fields with quantum relativistic particles.

Introduction

In this work, we are interested in the simulation of the relativistic dynamics of an electron of mass m coupled to an external classical electromagnetic field. The time-dependent Dirac equation under consideration reads [START_REF] Itzykson | Quantum Field Theory[END_REF] i∂ t ψ(t, x) = Hψ(t, x)

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

H = α • [cp -eA(t, x)] + βmc 2 + I 4 V (t, x), (2)
where the momentum operator is p = -i∇. More specifically, the Dirac equation under consideration reads [START_REF] Itzykson | Quantum Field Theory[END_REF] i∂ t ψ(t, x) = α x -ic∂ x -eA x (t, x) + α y -ic∂ y -eA y (t, x)

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

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

σ x = 0 1 1 0 , σ y = 0 -i i 0 and σ z = 1 0 0 -1 , (5)
while I 2 is the 2 × 2 unit matrix. Note that the light velocity c and fermion mass m are kept explicit in Eq. ( 2), allowing to adapt the method easily to natural or atomic units (a.u.). Throughout this work, we consider the single particle Dirac equation which is relevant for calculations describing Quantum Electrodynamics (QED) processes coupled to strong classical fields (e.g. particle-antiparticle pair creation from very high intensity electromagnetic classical fields [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF]). The purpose of this paper is not to derive original numerical solvers for the Dirac equation, or to provide simulations of specific quantum relativistic problems, but to compare the efficiency of two existing simple but accurate methods for solving the Dirac equation. Although pseudo-spectral methods are used for decades for solving the Dirac equation, for instance for Graphene or pair production simulations [START_REF] Ackad | Numerical solution of the Dirac equation by a mapped Fourier grid method[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF][START_REF] Piazza | Extremely high-intensity laser interactions with fundamental quantum systems[END_REF][START_REF] Faisal | A four-component Dirac theory of ionization of a Hydrogen molecular ion in a super-intense laser field[END_REF][START_REF] Gonoskov | Probing nonperturbative QED with optimally focused laser pulses[END_REF][START_REF] Salamin | Relativistic high-power laser-matter interactions[END_REF][START_REF] Selstø | Solution of the Dirac equation for hydrogenlike systems exposed to intense electromagnetic pulses[END_REF][START_REF] Wells | Asymptotic channels and gauge transformations of the time-dependent Dirac equation for extremely relativistic heavy-ion collisions[END_REF], recently a real space method [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF][START_REF] Salomonson | Relativistic all-order pair functions from a discretized single-particle Dirac Hamiltonian[END_REF] was established to address these physical questions. This real space method, which has closed connection with the Quantum Lattice Boltzmann method [START_REF] Fillion-Gourdeau | Formal analogy between the Dirac equation in its Majorana form and the discrete-velocity version of the Boltzmann kinetic equation[END_REF][START_REF] Succi | Lattice Boltzmann equation for quantum mechanics[END_REF] is simply a finite-difference method at CFL = 1 [START_REF] Strikwerda | Finite Difference Schemes and Partial Differential Equations[END_REF], for solving a first order linear hyperbolic system with eigenvalues of equal magnitude. This method which is based on the Method of Characteristics (MC) will be referred in this paper to as a MC-based method. Thanks to its simplicity, the MC-based method is shown to be highly scalable and much more efficient sequentially and in parallel than pseudo-spectral methods for a fixed number of degrees of freedom (dof). Its main weakness is however the strict condition which is imposed on the space step. In comparison, the pseudo-spectral methods are much more flexible regarding the choice of the spatial discretization step, that is the number of dof. In summary, for physical problems involving very small space scales, the MC-based will be shown to be much more efficient than pseudo-spectral methods. This will be the case for instance for problems involving heavy ions, or for intense and short laser-molecule interactions. In the opposite, when a fine spatial resolution is not required pseudo-spectral methods will be shown to be more efficient thanks to the possibility to select much coarser meshes, while keeping a very good accuracy (spectral convergence). Naturally, other types of methods exist for solving the Dirac equation, such as variational methods [START_REF] Fillion-Gourdeau | Galerkin method for unsplit 3-d Dirac equation using atomically/kinetically balanced B-spline basis[END_REF][START_REF] Fischer | A B-spline Galerkin method for the Dirac equation[END_REF][START_REF] Grant | Variational methods for Dirac wave equations[END_REF][START_REF] Grant | B-spline methods for radial Dirac equations[END_REF][START_REF] Kullie | Dirac-Fock finite element method (FEM) calculations for some diatomic molecules[END_REF][START_REF] Yang | Fully numerical relativistic calculations for diatomic molecules using the finite-element method[END_REF] or Krylov-type methods [START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF], but these are not discussed in this paper, which focuses on simple and efficient techniques.

The paper is organized as follows. In Section 2 (resp. Section 3), we recall the basics of the MC-based (resp. basic pseudo-spectral) method for the Dirac equation modeling the interaction of an electron with and external electromagnetic field. We next present in Section 4, an overview of the parallel computing aspects. In Section 5 (resp. Section 6), a series of sequential and parallel experiments is presented to illustrate the strengths and weaknesses of each method in the one-dimensional (resp. multi-dimensional) case. We conclude in Section 7.

Method of characteristic-based TDDE Solver in cartesian coordinates

Operator Splitting (first order)

We here recall the principle of operator splitting for the computation of an approximate solution to Eqs. ( 1) and ( 2) at time t n+1 , and denoted by ψ n+1 (x). The initial condition at time t n is assumed given by ψ(t n , x) = ψ n (x). [START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF] As previously discussed in [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF], this can be done [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] with an operator splitting scheme. We first define the operators

A = -icα x ∂ x (7) B = -icα y ∂ y (8) C = -icα z ∂ z (9) D = βmc 2 + I 4 V (t, x) -eα • A(t, x). (10)
The following splitting in Cartesian coordinates is considered [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF] (the x-dependence in the wavefunction argument for notational convenience):

i∂ t ψ (1) (t) = Aψ (1) (t), ψ (1) (t n ) = ψ n , t ∈ [t n , t n+1 ) ( 11
)
i∂ t ψ (2) (t) = Bψ (2) (t), ψ (2) (t n ) = ψ (1) (t n+1 ), t ∈ [t n , t n+1 ) ( 12
)
i∂ t ψ (3) (t) = Cψ (3) (t), ψ (3) (t n ) = ψ (2) (t n+1 ), t ∈ [t n , t n+1 ) ( 13
)
i∂ t ψ (4) (t) = Dψ (4) (t), ψ (4) (t n ) = ψ (3) (t n+1 ), t ∈ [t n , t n+1 ) (14) and ψ n+1 = ψ (4) (t n+1 ) ( 15
)
where the upper subscript in parenthesis on the wavefunction denotes the splitting step number. Note that this simple splitting scheme leads to an error that scales in O(∆t 2 ), corresponding to a first-order numerical scheme (for more details on the analysis of the method, see [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF]). The method consists of solving each equation independently with an initial condition given by the solution of the previous step. Note also that for every step, the time increment is identical, i.e. ∆t ≡ t n+1 -t n . This splitting approach will be used in this paper with both the MC-based and pseudo-spectral-methods.

Method of characteristics (MC) based method

In this section, we summarize the MC-based method which was proposed in [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]. Eqs. [START_REF] Faisal | A four-component Dirac theory of ionization of a Hydrogen molecular ion in a super-intense laser field[END_REF] to [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF] can actually be solved independently using the method of characteristics (MC). We first diagonalize the Dirac matrix, thus decoupling the spinor components. The resulting equation has a form similar to set of advection equations (linear first-order in time and space derivative). The method of characteristics is then used to find an explicit analytical solution. Finally, the solution is transformed back to the original Dirac matrix representation. The explicit solution to these equations and computational details can be found in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], as well as [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] in the cylindrical coordinates. The final result for the solutions of Eqs. [START_REF] Faisal | A four-component Dirac theory of ionization of a Hydrogen molecular ion in a super-intense laser field[END_REF] to ( 13) is

ψ (1) (t n+1 , x) = 1 2 [I 4 + α x ]ψ n (x -c∆t, y, z) +[I 4 -α x ]ψ n (x + c∆t, y, z) (16)
ψ (2) (t n+1 , x) = 1 2 [I 4 + α y ]ψ (1) (t n+1 , x, y -c∆t, z) +[I 4 -α y ]ψ (1) (t n+1 , x, y + c∆t, z) (17)
ψ (3) (t n+1 , x) = 1 2 [I 4 + α z ]ψ (2) (t n+1 , x, y, z -c∆t) +[I 4 -α z ]ψ (2) (t n+1 , x, y, z + c∆t) (18)
The space domain is then discretized in cubic elements with edges of length a = ∆x = ∆y = ∆z, inside which the wavefunction is constant (Q 0 -type elements). The discretized wavefunction and electromagnetic field are then written as

ψ h (t, i) = N m=1 1 m (i)ψ(t, xm ) (19)
A h (t, i) = N m=1 1 m (i)A(t, xm ) (20)
where N = N x N y N z is the total number of elements, ψ h (t, i) and A h (t, i) are the discretized wavefunction and electromagnetic field (i ≡ (i, j, k) ∈ Z 3 are indexing the volumes), the function 1 m (i) is equal to 1 in volume m indexed by i and is zero outside, while xm is the vector pointing to the center of volume m with the following components:

xm = x min + (i + 1 2 )a, y min + (j + 1 2 )a, z min + (k + 1 2 )a (21)
where x min , y min , z min are the lower domain boundary coordinates. The constraint c∆t = a (with space step ∆t and space step a) has to be satisfied for stability and accuracy reasons [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF]. It then allows to write ψ(x ± c∆t) → ψ h (i ± K), with K ∈ N * . This discretized scheme is then exact in each dimension (up to errors coming from the projection of the grid) for any K ∈ N * such that c∆t = Ka. Finally Eqs. ( 16) to [START_REF] Gonoskov | Probing nonperturbative QED with optimally focused laser pulses[END_REF] become for K = 1

ψ n 1 h (i) = 1 2 [I 4 + α x ]ψ n h (i -1, j, k) + [I 4 -α x ]ψ n h (i + 1, j, k) (22)
ψ n 2 h (i) = 1 2 [I 4 + α y ]ψ n 1 h (i, j -1, k) + [I 4 -α y ]ψ n 1 h (i, j + 1, k) (23)
ψ n 3 h (i) = 1 2 [I 4 + α z ]ψ n 2 h (i, k, k -1) + [I 4 -α z ]ψ n 2 h (i, j, k + 1) (24)
ψ n+1 h (i) = exp -iβmc 2 ∆t -i Ṽ n h (i) + iα • Ãh (i) ψ n 3 h (i) (25)
where

Ãn (x) = e t n+1 tn dτ A(τ, x) ≈ A n (x)∆t Ṽ n (x) = e t n+1 tn dτ V (τ, x) ≈ V n (x)∆t. ( 26
)
The constraint a = c∆t/K is very restrictive and makes sense when A, V (including V nuc. ) or ψ(•, 0) possess wavenumbers of the order of 1/mc. In Eq. ( 14), the solution is simply given by

ψ (4) (t n+1 , x) = T exp -i t n+1 tn dτ βmc 2 -eα • A(τ, x) × exp -ie t n+1 tn dτ V (τ, x) ψ (3) (t n+1 , x), ( 27
)
where T is the time-ordering operator. In practice it is shown in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], that the solution in Eq. ( 25) can be rewritten as:

ψ n+1 h (i) = U (i) exp -i Ṽ n h (i) ψ n 3 h (i) (28)
where U (i) is a matrix given explicitly by

U (i) ≡       c(A) -i mc 2 ∆t A s(A) 0 i Ãh,z (i) A s(A) [i Ãh,x (i)+ Ãh,y (i)] A s(A) 0 c(A) -i mc 2 ∆t A s(A) [i Ãh,x (i)-Ãh,y (i)] A s(A) -i Ãh,z (i) A s(A) i Ãh,z (i) A s(A) [i Ãh,x (i)+ Ãh,y (i)] A s(A) c(A) + i mc 2 ∆t A s(A) 0 [i Ãh,x (i)-Ãh,y (i)] A s(A) -i Ãh,z (i) A s(A) 0 c(A) + i mc 2 ∆t A s(A)       (29)
and where c(A) ≡ cos(A) , s(A) ≡ sin(A) and A ≡ (mc

2 ∆t) 2 + Ãh (i) • Ãh (i). ( 30
)
This method was implemented and tested in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF]. Although very robust, linear and very easy to efficiently implement in parallel (see Sections 5 and 6), the restriction on the spatial step makes it irrelevant in several physical configurations. Notice finally, that the MC-method can also be implemented in cylindrical coordinates with for instance azimuthal symmetry, as proposed in [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF]. In that case, the equation is again split in r and z, and the MC-based method can still be implemented in the z-direction, while in the r-direction a Poisson equation solver is used.

Regarding the mathematical analysis, we refer to [START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF], where the authors study the convergence analysis of the MC-based method coupled with a Maxwell equation Q 1 -finite element solver. For the MC-based method only, it is easy to show that it is diffusion-less in any dimension, and dispersion-less only in one dimension. The dispersion issue in higher dimension was considered in [START_REF] Hammer | A dispersion and norm preserving finite difference scheme with transparent boundary conditions for the Dirac equation in (1+1)d[END_REF], where the authors used staggered grids to reduce the numerical dispersion. Let us remark that the boundary condition problem is not discussed in this paper, but we refer to [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum wave equations[END_REF][START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF], for the interested reader.

Notice that in the MC-based method, the main numerical error comes from the operator splitting in four operators. We recall that splitting

ψ t = (A + B + C + D)ψ (with ψ(•, t n ) = ψ n (•))
, where A, B, C, D are four spatial differential or algebraic operators, involves the combination of four equations ψ t = Aψ, ψ t = Bψ, ψ t = Cψ, ψ t = Dψ. For a second order splitting, we get the following approximation: e ∆t(A+B+C+D) -e ∆tA/4 e ∆tB/2 e ∆tA/4 e ∆tC/2 e ∆tD e ∆tC/2 e ∆tA/4 e ∆tB/2 e ∆tA/4 ψ n

(L 2 ) 4 = O(∆t 3 ) (31)
In other words at each time iteration, an error in ∆t 3 is produced by the above time-splitting.

The splitting errors will also be produced in the pseudo-spectral methods (except for the so-called Unsplit-Pseudo-Spectral method) presented below.

Pseudo-spectral (PS) methods

The main issue with the finite volume/difference approach proposed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF] for discretizing the Dirac equation on the cartesian grid, or in the z-direction in cylindrical coordinates [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF], is due to a combination of a stability/ 2 -conservation constraint (CFL = 1), and a physical constraint (∆t 1/mc 2 ). The consequence is that, in order to avoid numerical diffusion, we need to choose a very small space step (∆z 1/mc). As recalled above, this restriction is only acceptable when the particle is subject to a classical electromagnetic field with very high wavenumbers, or when the interaction potential or the initial data possess small spatial scales. Although, the corresponding approach is very attractive in particular from a parallel computing viewpoint, it then suffers from a need of a very high number of gridpoints or finite volumes, or more generally of degrees of freedom (dof). In this section, we present simple alternatives based on pseudo-spectral methods, which simultaneously allow for i) preserving the high accuracy of the generalized transport equation solver, ii) the high scalability of the overall solver, and finally iii) to release the constraint on the size of the space step that is to significantly reduce of the overall number of degrees of freedom. The methodology which is now detailed is also valid in cylindrical coordinates in the z-direction and in cartesian coordinates in the (r, θ) directions. It basically consists of solving the corresponding hyperbolic systems by using one-dimensional FFTs, generalizing the method applied now to the generic one-directional scalar transport equation in 3d, u t (t, x, y, z) + cu x (t, x, y, z) = 0. The solution to this equation can indeed be simply computed by solving the following ordinary differential equation:

u t (t, ξ, y, z) + cF -1 x iξF x u(t, ξ, y, z) = 0
where F x is the Fourier transform with respect to x, and ξ denotes the co-variable of x. Naturally, this will be performed numerically using a discrete Fourier transform. More generally, we consider for convenience the 3-dimensional system in cartesian coordinates, for γ = x, y, z in the domain Ω = [-a x , a x , ] × [-a y , a y ] × [-a z , a z ]:

i∂ t ψ(t) = -icα γ ∂ γ ψ(t), ψ(t n ) = ψ n , t ∈ [t n , t n+1 )
and we use the same notation as [START_REF] Antoine | High-order IMEX-spectral schemes for computing the dynamics of systems of nonlinear Schrödinger/GrossPitaevskii equations[END_REF]. We first diagonalize α γ = Π γ D γ Π T γ where

D γ =     1 0 0 0 0 1 0 0 0 0 -1 0 0 0 0 -1     .
The matrices Π γ are defined as follows.

Π x = 1 √ 2     0 1 1 0 1 0 0 -1 1 0 0 1 0 1 -1 0     , Π y = 1 √ 2     0 -i -i 0 1 0 0 1 -i 0 0 i 0 1 -1 0     , Π z = 1 √ 2     1 0 0 -1 0 -1 -1 0 1 0 0 1 0 1 -1 0     .
We set φ := Π T γ ψ, which then satisfies

i∂ t φ(t) = -icD γ ∂ γ φ(t), φ(t n ) = Π T γ ψ n , t ∈ [t n , t n+1 )
We denote the grid-point set by

D Nx,Ny,Nz = x k 1 ,k 2 ,k 3 = (x k 1 , y k 2 , z k 3 ) (k 1 ,k 2 ,k 3 )∈O NxNy Nz with O NxNyNz = (k 1 , k 2 , k 3 ) ∈ N 3 , : h = 0, • • • , N x -1; k = 0, • • • , N y -1; = 0, • • • , N z -1 .
Then, we define

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

ξ p = pπ/a x with p ∈ {-N x /2, • • • , N x /2 -1}, ξ q = qπ/a y with q ∈ {-N y /2, • • • , N y /2 -1} and ξ r = rπ/a z with r ∈ {-N z /2, • • • , N z /2 -1}.
Then, we can define the partial Fourier pseudospectral approximations φ ( ) with ∈ {1, 2, 3, 4}, in the x, y-and z-directions

               φ ( ) (t, x, y, z) = 1 N x Nx/2-1 p=-Nx/2 φ ( ) p (t, y, z)e iξp(x+ax) φ ( ) (t, x, y, z) = 1 N y Ny/2-1 q=-Ny/2 φ ( ) q (t, x, z)e iξq(y+ay) φ ( ) (t, x, y, z) = 1 N z Nz/2-1 r=-Nz/2 φ ( ) r (t, x, y)e iξr(z+az)
The Fourier coefficients are defined as follows, where we denote φ ( )

k 1 (t, y, z) = φ ( ) (t, x k 1 , y, z), φ ( ) k 2 (t, x, z) = φ ( ) (t, x, y k 2 , z) and φ ( ) k 3 (t, x, y) = φ ( ) (t, x, y, z k 3 ).        φ ( ) p (t, y, z)e iξp(x+ax) = Nx-1 k 1 =0 φ ( ) k 1 (t, y, z)e -iξp(x k 1 +ax) φ ( ) q (t, x, z)e iξq(y+ay) = Ny-1 k 2 =0 φ ( ) k 2 (t, x, z)e -iξq(y k 2 +ay) φ ( ) r (t, x, y)e iξr(z+az) = Nz-1 k 3 =0 φ ( ) k 3 (t, x, y)e -iξr(z k 3 +az)
We finally define the approximate first-order partial derivatives

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

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

In the following subsections, we propose different pseudo-spectral methods by using the formalism that presented above.

Classical FFT-method

It has to be noticed that the pseudo-spectral method which will be presented below, differs from the classical Fourier-type method, which typically requires a splitting of the equation described in this subsection. The equation is first split in [-a x , a x , ] × [-a y , a y ] × [-a z , a z ] using ( 7), ( 8), ( 9), [START_REF] Piazza | Extremely high-intensity laser interactions with fundamental quantum systems[END_REF]. The first order splitting from t n , t n+1 , reads as follows.

1. First step: integration of the source from t n to t n+1 .

ψ n 1 h = T exp -i t n+1 tn dτ [βmc 2 -eα • A h (τ )] × exp -ie t n+1 tn dτ V h (τ ) ψ n h .
This is performed using the time-ordering operator. 2. Second step: integration of the generalized transport equation in Fourier space. One sets φ n 1 h := Π T x ψ n 1 h , and solve

∂ t φ + Λ x ∂ x φ = 0 basically by discretizing in space: φ(t n+1 , •) = F -1 x e -iξ∆tΛx F x (φ) , where F x is the Fourier transform in the x-direction that is for = 1, • • • , 4 φ ( ),n 2 h = 1 N x Nx/2-1 p=-Nx/2 e -iξp∆tλ (x) Nx-1 k 1 =0 φ (j),n 1 h,k 1 e -iξp(x k 1 +ax) e iξp(x+ax j )
where λ (x) is the th eigenvalues of Λ x .

Similarly one computes

φ ( ),n 3 h = 1 N y Ny/2-1 q=-Ny/2 e -iξq∆tλ (y) Ny-1 k 2 =0 φ ( ),n * 2 h,k 2 e -iξq(y k 2 +ay) e iξq(y+ay j )
where φ

n * 2 h = Π T y Π x φ n 2 h and finally φ ( ),n+1 h = 1 N z Nz/2-1 r=-Nz/2 e -iξr∆tλ (z) Nz-1 k 3 =0 φ ( ),n * 3 h,k 3 e -iξr(z k 3 +az) e iξr(z+az j )
where φ

n * 3 h = Π T z Π y φ n 3 h .
Although this approach is classical and allows for a very accurate and fast approximation of each split equation, it also systematically necessitates a splitting of the TDDE, which to a certain extent limits its interest. In the following this method will be referred as the pseudo-spectral FFT-method. Alternatively, it is possible to directly apply 3d FFTs.

SPS-scheme

We here describe the numerical pseudo-spectral method for solving the split multidimensional TDDE in cartesian coordinates (3). It will be referred as the (split-pseudospectral) SPS-method. As proposed before, this will be performed by splitting the TDDE into four time dependent systems ( 7), ( 8), ( 9) and [START_REF] Piazza | Extremely high-intensity laser interactions with fundamental quantum systems[END_REF].

1. First step: integration of the generalized transport equation in Fourier space in the x-direction. One sets φ n h := Π T x ψ n h , and the system

∂ t φ + cΛ x ∂ x φ = 0, φ(t n , •) = Π T x ψ(t n , •) is approximately solved by φ n 1 h = φ n h -c∆tΛ x [[∂ x ]] φ n h . We then deduce ψ n 1 h = Π x φ n 1 h .
Notice that it is possible to use a higher order discretization in time, such as

   φ n * h = φ n h -c∆tΛ x [[∂ x ]] φ n h φ n 1 h = φ n h -c ∆t 2 Λ x [[∂ x ]] φ n h + [[∂ x ]] φ n * h .
The corresponding operator is denoted by P x (∆t). 2. Second step: integration of the generalized transport equation in Fourier space in the y-direction

∂ t φ + cΛ y ∂ y φ = 0. We set φ n * 1 h := Π T y ψ n 1
h , and we solve

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

ψ n+1 h = T exp -i t n+1 tn dτ βmc 2 -eα • A h (τ ) × exp -ie t n+1 tn dτ V h (τ ) ψ n 3 h . (32)
The corresponding operator is denoted by Q(∆t).

This scheme is fully explicit, and can be compactly expressed as [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], the method of characteristics is replaced by a pseudo-spectral method. In space, the convergence is spectral (FFT). In cartesian coordinates the above methodology is applied in each split direction. In cylindrical coordinates, the Poisson solution is required in the r-direction. Another important interest of this method is that the order of accuracy of the splitting can easily be increased, unlike the MC-method which does not allow for splitting of order higher than 2.

ψ n+1 h = P x (∆t)P y (∆t)P z (∆t)Q(∆t)ψ n h . In practice,

UPS-scheme

The approach developed in this subsection is close to the SPS-method, but without splitting the system. The corresponding method will be referred as the (Unsplit Pseudospectral) UPS-method. By using an unsplit method, we avoid errors coming from the noncommutation of the Dirac matrices. Notice however that the mass operator contribution βmc 2 , cannot be solved analytically anymore, and then still requires a very small time step or higher order discretization in time to accurately approximate its contribution. We first present the methodology in one-dimension. One-dimensional UPS-method.

i∂ t ψ(t, x) = -icσ x ∂ x + mc 2 σ z + V (t, x)I 2 -eA(t, x)σ x ψ(t, x) (33)
We introduce the notations

Λ x = 1 0 0 -1 , Π x = Π T x = 1 √ 2 1 1 1 -1
where

σ x = Π x Λ x Π T x , φ = Π T x ψ and Γ x (t, x) = Π T x mc 2 σ z + V (t, x)I 2 Π T x = V (t, x) mc 2 mc 2 V (t, x) .
The TDDE is then rewritten

∂ t φ + cΛ x ∂ x φ -ieA(t, x)Λ x φ + iΓ x (t, x)φ = 0. (34)
In the UPS-formalism, we consider the semi-discrete in space system:

∂ t φ + cΛ x [[∂ x ]]φ -ieA(t, x)Λ x φ + iΓ x (t, x)φ = 0. ( 35
)
Using the same notation as above, and denoting by ψ h the discrete solution to [START_REF] Selstø | Solution of the Dirac equation for hydrogenlike systems exposed to intense electromagnetic pulses[END_REF] and φ h the one to [START_REF] Succi | Lattice Boltzmann equation for quantum mechanics[END_REF], the numerical approximation reads

φ n+1 h = φ n h -c∆tΛ x [[∂ x ]] φ n h + i∆t eA n h Λ x -Γ n
x,h φ n h . Any higher order approximation of the equation is actually possible, such as the following order-2 scheme.

φ n * h = φ n h -∆tΛ x [[∂ x ]] φ n h + i∆t eA n h Λ x -Γ n x,h φ n h , φ n+1 h = φ n h - ∆t 2 Λ x [[∂ x ]] φ n * h + Λ x [[∂ x ]] φ n h -i ∆t 2 eA n h Λ x -Γ n x,h φ n h + eA n * h Λ x -Γ n * x,h φ n * h .
We now extend this approach in the multi-dimensional case.

Multi-dimensional UPS-method. We now solve from t n to t n+1

i∂ t ψ(t, x) = α x -ic∂ x -eA x (t, x) + α y ic∂ y -eA y (t, x) +α z -ic∂ z -eA z (t, x) + βmc 2 + I 4 V (t, x) ψ(t, x).
Assuming that ψ h is known at time t n , the semi-discrete in space scheme reads in [-a x , a

x , ]× [-a y , a y ] × [-a z , a z ] i∂ t ψ h = α x -ic[[∂ x ]] ψ h -eA x,h (t)ψ h + α y -ic[[∂ y ]] ψ h -eA y,h (t)ψ h +α z -ic[[∂ z ]] ψ h -eA z,h (t)ψ h + βmc 2 + I 4 V h (t)ψ h (36)
which is nothing but a system of differential equations, which can be solved using an explicit or implicit high order scheme. Assuming that ψ n h is known, a first order explicit scheme in time reads

iψ n+1 h = iψ n h + ∆t α x -ic[[∂ x ]] ψ n h -eA n x,h ψ n h + α y -ic[[∂ y ]] ψ n h -eA n y,h ψ n h +α z -ic[[∂ z ]] ψ n h -eA n z,h ψ n h + βmc 2 + I 4 V n h ψ n h
where A n x,h stands for A x,h (t n ), etc. However such a scheme would not ensure stability and accuracy in principle. In practice higher order scheme will be used. An alternative from ( 36) is then to formally integrate analytically this expression by expliciting the [[∂ γ ]] ψ terms (γ = x, y, z). Say from time t n+1 to t n :

ψ h (t n+1 ) = exp ie α x t n+1 tn A x,h (s)ds + α y t n+1 tn A y,h (s)ds + α z t n+1 tn A z,h (s)ds -imc 2 (t n+1 -t n )β -iI 4 t n+1 tn V h (s)ds × c(t n+1 -t n ) α x [[∂ x ]] ψ h (t n ) + α y [[∂ y ]] ψ h (t n ) + α z [[∂ z ]] ψ h (t n ) × t n+1 tn exp -ie α x s tn A x,h (τ )dτ + α y s tn A y,h (τ )dτ + α z s tn A z,h (τ )dτ -imc 2 (s -t n )β -iI 4 s tn V h (τ )dτ
which then requires high order approximations of time integrals.

Parallel computing

The parallel computation of the generalized transport systems which is detailed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF], is highly scalable thanks to a spatial domain decomposition (in any direction x, y, z). At each subdomain interface, the Dirac equation is diagonalized in the direction orthogonal to the subdomain interface. We then impose Dirichlet boundary conditions at these interfaces and on the components corresponding to the outgoing characteristics, and these conditions are sent by message passing. The spatial domain can be decomposed in layers in one specific directions, or alternatively in small cubic subdomains.

It is well known that the parallel implementation of the FFT suffers from a relative lack of scalability. However the pseudo-spectral methods presented in this paper allow for using the high accuracy and relative efficiency of the parallel one-dimensional FFTs in the zdirection, and also keeping the nice scalability properties of [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] in the x-and y-directions where spequential FFTs are performed. More specifically in 3 dimensions, we proceed by alternating the directions. We first decompose the domain by layers in the z-directions (see Fig. 1). Then

• We successively perform the evolution (FFT) sequentially in the x-and y-directions, by layer in the z-direction. Each processor manages one layer in z. We expect a perfect scaling for this step, as it does not require any transmission from one node to another.

• For all x and y, we perform the evolution in the z-direction using the parallel FFT (fftw). The performance of this step is fully dependent on the parallelization of the one-dimensional FFT. The main difference with the MC-based methods proposed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] is that we do not have the space-step restriction due to the CFL-number, but keeping simultaneously a high accuracy and high scalability. Stability and convergence properties are also preserved. As a consequence, the number of degrees of freedom in space can largely be reduced. Notice that the parallelization approach is valid for the classical FFT-method, as well as for the SP-methods presented in this paper. The computational complexity (CC) is relatively simple to established. Assume that the real-space grid is composed by N x × N y × N z volumes. We denote by N T f the number of time steps to reach a final physical time T f . The computation of the MC-solver for the split generalized transport equation [START_REF] Guo | Spectral and pseudospectral approximations using Hermite functions: Application to the Dirac equation[END_REF] (resp. [START_REF] Hammer | A dispersion and norm preserving finite difference scheme with transparent boundary conditions for the Dirac equation in (1+1)d[END_REF], resp. ( 24)) requires

c x N T f N x N y N z operations (resp. c y N T f N x N y N z , resp. c z N T f N x N y N z )
, where the integers c x ≈ c y ≈ c z are typically integers smaller than 10. The solution to the differential system (27) also requires c T f N T f N x N y N z operations, where c T f is dependent on the order of the numerical integrator and of the time-ordering operator. In fine, the total complexity C MC , is then given by

C MC = O N T f N x N y N z
with a small prefactor. In addition this method is pleasingly parallel as discussed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], and efficiency greater than 1 can even be observed due to cache effects. By comparison, the CC of the UPS-method, C UPS , and for the SPS-method, C SPS , are

C UPS,SPS = O N T f N x N y N z log(N x N y N z ) .
Compared to the MC-method, the prefactor is however much higher. The SPS-method also allows for a very efficient parallel computing of the "many" one-dimensional FFTs. Indeed, the UPS/SPS require the computation of one-dimensional FFTs of three-dimensional wavefunctions. Two of three real variables are then fixed, while the FFT is applied in the third one. Notice however that multi-dimensional FFTs can be very efficiently computed in parallel, while one-dimensional parallel FFT versions may suffer from a loss of efficiency. We refer to [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF][START_REF] Guo | Spectral and pseudospectral approximations using Hermite functions: Application to the Dirac equation[END_REF] for the presentation and mathematical analysis of pseudo-spectral methods in different regimes. In the relativistic regime, we also refer to [START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Mocken | Quantum dynamics of relativistic electrons[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF] for efficient pseudo-spectral methods.

In the following, the first three one-dimensional tests are performed by using matlab. The fourth one-dimensional test, as well as all the multi-dimensional tests are performed by using a C++-code with MPI-library. The parallel tests are realized on the computer mammouth-parallel II from Compute Canada, with 30984 cores, 24 cores/node, 32 G/node & FAT nodes, opteron 2.1 GHz.

One-dimensional experiments

We now compare the performance of the pseudo-spectral methods with the MC-based method. It is important to recall that pseudo-spectral methods are expected to be more attractive than the MC-based method, only if the space variations of either the electromagnetic field (EM), the interaction potential, or the initial data are much larger than 1/mc. Typically, denoting by k ∞ the largest wavenumber of the EM and of V nuc. , the space steps for respectively the PS-methods and MC-method, are selected as follows:

∆x PS k ∞ /M, ∆x MC min k ∞ /M, 1/mc
where M is a non-zero integer. Notice that the interaction of the intense EM with quantum particle generates very high frequency photons, and may potentially necessitate even finer spatial resolution.

More specifically, we compare the efficiency of the FFT-, SPS-, UPS-and MC-schemes which are presented in Section 3 and Appendix A. Recall that the strength of the UPS-scheme compared to the SPS-and FFT-scheme is that it does not require the split of the Dirac Hamiltonian. On the other side, the FFT-, SPS-and UPS-schemes do not require a strict condition on the space step, unlike the MC-based scheme (∆x = c∆t). In the following, the external field the particle is subject to, is assumed space-independent, and is given by

A(t) = A 0 sin exp(-α (t -T f ) 2 4 sin 2Cπt/T f . (37)
In [START_REF] Yang | Fully numerical relativistic calculations for diatomic molecules using the finite-element method[END_REF], A 0 denotes the maximal amplitude to the field, T f the pulse duration, C the number of cycles and α a real parameter characterizing the pulse envelope.

Test 1.

In this first test, we compare the MC-based scheme on a very fine mesh (N (mc) = 54799), with the pseudo-spectral methods on coarser meshes. In [START_REF] Yang | Fully numerical relativistic calculations for diatomic molecules using the finite-element method[END_REF], we select A 0 = -100, the number of cycles C = 6, T f = 5 × 10 2 /137 2 , and α = 10 7 . The field is represented in Fig. 2 ψ (1) (0, x), ψ (2) (0, x)

T = e ik 0 x-x 2 /10 , 0 T where the wavenumber k 0 = 4. The scalar potential is hence null and, the nuclear potential is of Coulomb type with charge Z = 1000, and softcore parameter ε = 10 -1 :

V nuc. (x) = - Z √ x 2 + ε .
That is we have V (t, x) = V nuc. (x). The other physical data are as follows (in atomic units): c = 137, m = 1 and the computational domain is given by (-a, a), with a = 20 and ∆t = 10 -1 /mc 2 = 10 -1 /137 2 . We report in The MC-based scheme for a fixed ∆t, requires ∆x = c∆t = 1/mc, that is N (mc) ≈ 2amc dof. Unlike the MC-scheme, the pseudo-spectral methods allowed for choosing space steps (much) larger than c∆t, then at a cheaper computational complexity, while keeping a good accuracy.

For instance, it is reported in Figs. 3,4, the bispinor components at final time, using the MC-based scheme with N (mc) = 54799 and the UPS-scheme with N (ups) = 5481, 1097, 275. We observe that with much less dof (1097 vs 54799), we still keep a good accuracy. Notice that the simple FFT-method is the most efficient method, but requires an operator splitting unlike the UPS-method.

Test 2.

In the following, we still consider a simple test with small scales in space (and time). The use of the MC-based is shown to be more competitive than pseudo-spectral methods (more specifically the so-called FFT-method), as it is required to select ∆x at least as small as c∆t for an accurate spatial resolution. The EM is given by [START_REF] Yang | Fully numerical relativistic calculations for diatomic molecules using the finite-element method[END_REF] where A 0 = -500, the number of cycles C = 6, T f = 5 × 10 2 /137 2 , and α = 10 7 . The initial wavefunction is given by ψ (1) (0, x), ψ (2) (0, x)

T = e ik 0 x-x 2 /2 , 0 where the wavenumber k 0 = 4. The scalar potential is hence null and, the nuclear potential is a combination of three nuclear potentials with Z = 500, and a softcore parameter is ε = 10 -6

V nuc. (x) = - 2Z (x -2.01) 2 + ε - Z (x -2) 2 + ε - 3Z (x -1.99) 2 + ε .
This potential is represented in Fig. 5. The other physical data are as follows (in atomic units): c = 137, m = 1, and the computational domain is given by (-5, 5), ∆t = 10 -2 /mc 2 = 1/137 2 . We compare the MCbased solution with N (mc) = 1.37 × 10 5 grid points and the FFT-method with respectively N (mc) = 1.37 × 10 5 and N (fft) = 1.37 × 10 4 . More specifically, in Fig. 6, we represent the first and second components of the MC-(N (mc) = 1.37 × 10 5 ) and FFT-solutions with N (fft) = 1.37 × 10 4 , close to the nuclei (Top-left), (Top-right), as well as the spatial spectrum of ψ 1 (T f , •) (i.e.

k, F x ψ 1 (T f , k) 2 in logscale). Although, the computational time is much smaller using the FFT-method (5.8 seconds versus 48.7 seconds for the MC-based scheme) the accuracy in space is insufficient to capture the very small spatial scales in V nuc. . In comparison, if we also take N (fft) = 1.37 × 10 5 for the FFT-method, the two solutions are comparable as observed in Figs. 7, as well as the spatial spectrum of ψ 1 (T f , •). However the MC-method is now shown to be more efficient: 48.7 seconds versus 68.5 seconds for the FFT-method.

(T f , •) in logscale: k, F x ψ 1 (T f , k) 2 .
Although simple, this test illustrates that the MC-based method is sequentially more competitive than the FFT-method for multiscale problems in space.

Test 3.

We again compare the performance of the MC-based scheme with the pseudo-spectral scheme. The EM is assumed space-independent and given by [START_REF] Yang | Fully numerical relativistic calculations for diatomic molecules using the finite-element method[END_REF], with A 0 = -100, the number of cycles C = 8, T f = 2 × 10 3 /137 2 , and α = 2 × 10 3 . The initial wavefunction is given by

ψ (1) (0, x), ψ (2) (0, x) T = e ik 0 x-x 2 /2 , 0 T
where the wavenumber k 0 = 5. The scalar potential is hence null and, the nuclear potential is of Coulomb type with charge Z = 100, and softcore parameter ε = 10 -1

V nuc. (x) = - Z √ x 2 + ε .
Therefore, one gets V (t, x) = V nuc. (x). The other physical data are as follows (in atomic units): c = 137, m = 1, and the computational domain is given by (-20, 20) with ∆t = 10 -1 /mc 2 = 1/137 2 . For both schemes, a second-order splitting (Strang) scheme is used to approximate the TDDE. When using the SPS-scheme, we take ∆x = c∆t ≈ 7.3 × 10 -4 (resp.

c √ c∆t/2 ≈ 1.46×10 -3 , c √ c∆t/2 ≈ 4.27×10 -3 , c √ c∆t ≈ 8.54×10 -3 , 2c √ c∆t ≈ 1.71×10 -2 , 5c √ c∆t ≈ 4.27 × 10 -2 , c 2 ∆t ≈ 10 -1
) the corresponding number of spatial grid-points is given by N (sps) = 54800 (resp. 27400, 9364, 4682, 2341, 937, 401). For the MC-based scheme, we take ∆x = c∆t ≈ 7.3 × 10 -3 and the number of spatial grid-points is fixed and given by N (mc) = 54800. Although the MC-based scheme is much simpler to implement, its computational cost is naturally much higher whenever N (mc) N (sps) . We first report in Fig. 8 the bispinor components at final time T f by using the MC-based method with N (mc) grid-points and the SPS-method with N (sps) = 937 and 401 gridpoints. In Table 2, we report the CPU-time to solve the TDDE using the MC-based method on N (mc) = 54800 grid-points (∆x = c∆t = 1/mc), and the SPS-method using N (sps) = 54800

(∆x = c∆t), N (sps) = 9364 (∆x = √ c∆t/2 = 1/2 √ cm), N (sps) = 4682 (∆x = c √ c∆t = 1/m √ c), N (sps) = 2341 (∆x = 2c √ c∆t = 1/m √ c), N (sps) = 937 (∆x = 5 √ c∆t = 5/ √ cm)
, and N (sps) = 401 (∆x = c 2 ∆t = 1/m). Recall again that the CFL= 1 condition constraining the MC-method, no more occurs with the SPS-method which allows for choosing space steps much larger than c∆t.

Test 4.

The following test is devoted to the analysis of the scalability of the MC-based and FFTmethods, still in one-dimension. The domain we consider is (-64, 64). The initial density (ρ 0 (•) = 4 i=1 |ψ 4 (0, •)| 2 ) and the interaction potential (with 3 nuclei) are represented in Figs. 9 (Top-left) and (Top-right). The nuclei are respectively located in z = -10, 0, and 10, and their respective charge is Z 1 = 30, Z 2 = 60 and Z 3 = 70. Softcore constants are used to regularize the potential singularities, and are respectively equal to 5.0, 0.791 and 0.071. The time step is fixed to 1.78158 × 10 -6 and the number of degrees of freedom is N = 524288. The final physical time is T f = 1.78158 × 10 -4 . We also report the density at final time in the global domain in Figs. 9 (Bottom-left) and (Bottom-right). We compare the scalability of the MC-based and FFT-method for 2 2p processors with p = 0, • • • , 3. The MC-based scheme is parallelized by domain decomposition thanks to the algorithm described in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], and shortly recalled in Section 4. The efficiency of the parallelization for this method comes from i) the hyperbolicity, ii) the linearity of the Dirac equation and iii) the use of alternate direction splitting. The FFT-method which is used here, is based on the standard, and celebrated parallel version of fftw. In Fig. 10, we report the CPU-times for both the FFTand MC-based methods. More specifically, we provide i) the efficiency T 1 /pT p (where T p denotes the CPU time for solving the equation on p processors with p = 1, 4, 16, 64), where p denotes the number of processors in Fig. 10 (Top-left), and ii) the CPU-time in logscale as a function of the number of processors (Top-right).

We previously noticed on low dimensional problems, that for a fixed number of dof, the MC-based method was more efficient than the FFT-method. It is still observed on this higher dimensional problem that the MC-based method has a much better speed-up than the FFT-method, in particular when increasing the number of processors with a fixed number of dof.

Unlike the MC-based method, the FFT-method permits for fixed ∆t to select much larger spatial steps while keeping a very good accuracy in space. To illustrate this fact, we report on Fig. 10 (Bottom), and for 1 processor the CPU-time with the FFT-method on the same problem as above, but with a total number of dof: 524288/2 2p , with p = 0, • • • , 3. The MCbased method maintains an excellent efficiency up to 64 processors, unlike the FFT-method for which a strong discrepancy is observed beyond 16 processors.

Multi-dimensional experiments

In this section, we present large scale experiments in 2-d and 3-d in order to compare the efficiency of the MC-based and the pseudo-spectral FFT-methods. The latter is the most efficient and simple pseudo-spectral method. Due to the dimension of the problem, the methods are implemented in parallel following the approach described in Section 4. The MCbased algorithm is parallelized by using a simple domain decomposition strategy: the whole domain is divided into subdomains and the mesh data in each of these subdomains is managed to a different process. The communication between the subdomains which is required in the computation of the solution in the elements close to the subdomain boundaries is performed by using the Message Passing Interface (MPI) library.

Test 1.

This test is dedicated to the two-dimensional evolution in the (x, z)-plane of a wavepacket with wavenumbers k x = k z = 5, angular moment j z = 0.5, and its interaction with 3 atoms of equal charge Z 1 = Z 2 = Z 3 = 20, with softcore constant equal to 5. The nuclei are respectively located in (10, -10), (0, 0) and (-10, 10). We report the initial density ρ 0 (•, •) = 10 -4 , and the final physical time T f = 0.912169 (corresponding to 4000 iterations). The number of degrees of freedom is N x × N z = 4096 × 4096. We report in Fig. 12, the overall density at final time. The main purpose of this test is to illustrate the efficiency of the MCbased method. We denote by p x (resp. p z ) the number of processors in the x-direction (resp. z-direction), and by p the total number of processors with p = p x × p z . For the MC-based method, we respectively present some efficiency results in the case where i) p = p x × p z = p z (that is parallelism only in the z-direction), and in the case ii), where p x = p z = √ p. In both cases, we will use up to p = 256 processors. The FFT-method is parallelized in the (x, z)plane, by decomposing the domain only in layers in the z-direction (p = p z ). The main issue regarding the FFT-method is the use of parallel one-dimensional FFTs in the z-direction. Its poor efficiency when the number of dof is relatively "small" deteriorates the overall efficiency of the method despite a perfect scaling in the x-direction (as the one-dimensional FFTs in the x-direction are computed independently on different processors). As N z = 4096 is too small for efficient parallelization of the one-dimensional FFT, we only report the efficiency of the FFT-method for 1 to 4 processors. Even with this small number of processors, we indeed observe already a strong discrepancy of the efficiency. In comparison, the MC-based method shows an excellent efficiency (even higher than 1 due to cache-effects) up to 64 processors for both decompositions i) z-direction (p = p z ) and in x, z (p x = p z = √ p).

We observe that with both types of decomposition and parallelism, the efficiency is relatively identical. However, it is interesting to notice that at 256 processors, the second approach (p x = p z = √ p) has a better efficiency than the first one (p = p z ). In the following, we will only use a decomposition by layers in the z-direction.

Test 2.

In the following test, we still propose a two-dimensional simulation in the (x, z)-plane, in order to illustrate that due to the parallel one-dimensional FFT in the z-direction, the relative loss of efficiency of the method, can actually be reduced if N z is taken larger. A wavepacket with wavenumbers k x = k z = 5, angular moment j z = 0.5, is interacting with 2 atoms of equal charge Z 1 = Z 2 = 50, and with softcore constants equal to 1. The nuclei are respectively located in (-1, -10), and (-1, 10). We report the initial density and potential in Fig. 13 (Left) and (Right).

The overall computational domain is (-4, 4) × (-64, 64), the time step is ∆t = 5.70106 × 10 -5 , and the physical times T = 0.05521 and T f = 0.1142 (corresponding to 1000 and 2000 iterations). The number of degrees of freedom is N x × N z = 1024 × 16384 (resp. N x × N z = 256 × 4096). We report in Fig. 14, the overall density at final time. The MCbased method is parallelized in the z-direction. We observe in Fig. 15 that the efficiency of the FFT-method is indeed slightly improved when a larger N z is selected (4096 vs 16384). This is due to an improvement of the efficiency of the parallel one-dimensional FFT, for larger numbers of dof. As mentioned several times along this paper, although at fixed number of dof (that is fixed space step), the FFT-method is much less efficient (sequentially and in The magnetic field is defined by A = (A x , A y , A z ), with A x = A y = 0 and A z is represented at final time T f in the (x, y)-plane at z = 0, and in Fig. 17 (Left), and in the (x, z)plane and at y = 0 in Fig. 17 (Right). The initial data is a wavepacket centered at (-2.5, 0, 0) with wavenumber k x = 10, k y = -10, k x = -10, and the corresponding density ρ 0 (•) = 4 i=1 |ψ i (0, •)| 2 is represented in the (x, y)-plane at z = 0 in Fig. 16 (Left), and in the (x, z)-plane at y = 0 in Fig. 16 (Right). We select the time step equal to ∆t = ∆x/c, that is ∆t ≈ 9.1 × 10 -4 . We then report in Fig. 18 the density ρ(T f , x, y, 0) in the (x, y)-plan at z = 0 at final time T f (Left), as well as ρ(T f , x, 0, z) in the (x, z)-plan at y = 0 at T f (Right). We denote by p x (resp. p y , p z ) the number of processors in the x-direction (resp. y-,zdirection), and by p the total number of processors with p = p x × p y × p z . For the MC-based method, we respectively present efficiency results in the case i) p = p x × p y × p z = p z (that is parallelism only in the z-direction) with p = p z = 1, 4, 16, 64, 256, and ii) p x , p y and p z are greater than 1 (parallelism in x, y, z-directions): p x = p y = p z = 1 (p = 1), p x = p y = p z = 2 (p = 8), p x = p y = p z = 4 (p = 64) and p x = p y = 4, p z = 16 (p = 256). We report in Fig. 19, the efficiency and CPU-time as a function for the MC-and FFT-methods. We notice that the MC-based method is again the most efficient, with in particular, a better efficiency with parallelism in the z-direction. However there is a noticeable deterioration of the efficiency of the MC-based method beyond 16 processors, and unlike the 2d-case. This is due to the heavier load of data to exchange between processors in 3-d.

Conclusion

In this paper we have proposed a performance comparison of simple and efficient Dirac equation solvers: a real space method based on the characteristic equation (MC-based method), also referred as a Quantum Lattice Boltzmann method, with simple pseudo-spectral methods (PS-methods). We have established that for a given number of degrees of freedom, the MC-based method is much more efficient sequentially and in parallel than PS-methods. The price to pay is a strict condition imposed on the spatial discretization step. However and except for multiscale problems in space, the PS-methods do not impose such constraints, and allow then for taking much larger spatial discretization steps (then reducing the overall computational complexity), while keeping a good accuracy. In conclusion, the MC-based method should be preferred for problems involving very small spatial scales, such as those involving heavy ions, ultrashort laser-atom interaction, while the PS-methods are more adapted for other the physical configurations.

Acknowledgments. The authors would like Dr F. Fillion-Gourdeau (INRS) for helpful discussions.

Appendix A. One-dimensional numerical schemes

Classical FFT-scheme. This scheme necessarily requires the splitting of the Dirac Hamiltonian. From time t n to t n+1 the schemes reads as follows.

1. Source term integration, and 0 j N (ps) -1

ψ (1),n 1 h,j = exp -i∆t(mc 2 + V nuc.,j )/2 ψ (1),n h,j , ψ (2),n 1 h,j = exp -i∆t(-mc 2 + V nuc.,j )/2 ψ (2),n h,j .
2. We then solve the potential-free and mass-free Dirac Hamiltonian. h,k 1 e -iξp(x k 1 +ax) e iξp(x+ax j ) .

• We set φ n 1 h = φ (1),n 1 h , φ (2)
• Then, we get

ψ (1),n 2 h = φ (1),n 2 h + φ (2),n 2 h √ 2 , ψ (2)
,n 2 h = φ

(1),n 2 h -φ

(2),n 2 h √ 2 .

3. Finally, we integrate one more time the source term (Strang-splitting), for 0 j N (ps) -1 ψ (1),n+1 h,j = exp -i∆t(mc 2 + V nuc.,j )/2 ψ (1),n 2 h,j , ψ

(2),n+1 h,j = exp -i∆t(-mc 2 + V nuc.,j )/2 ψ (2),n 2 h,j . Split Pseudo-spectral (SPS) scheme. From time t n to t n+1 the scheme reads 1. Source term integration, and 0 j N (ps) -1 ψ (1),n 1 h,j = exp -i∆t(mc 2 + V nuc.,j )/2 ψ (1),n h,j , ψ

(2),n 1 h,j = exp -i∆t(-mc 2 + V nuc.,j )/2 ψ (2),n h,j . k 1 e -iξp(x k 1 +a) e iξp(x j -a) .

We set φ n

1 h = φ (1),n 1 h , φ (2)
• Then, for all 0 j N (ps) -1

φ (1),n * 1 h,j = φ (1),n 1 h,j -c∆t[[∂ x ]] φ (1),n 1 h,j + i∆teA n 1 φ (1),n 1 h,j , φ (2),n * 1 h,j = φ (2),n h,j + c∆t[[∂ x ]] φ (2),n 1 h,j -i∆teA n 1 φ
(2),n 1 h,j .

• Followed, for all 0 j N (ps) -1, by -φ

φ (1),n 2 h,j = φ (1),n h,j - c∆t 2 [[∂ x ]] φ (
(2),n 2 h √ 2

3. Finally, we integrate one more time the source term (Strang-splitting), for 0 j N (ps) -1 ψ (1),n+1 h,j = exp -i∆t(mc 2 + V nuc.,j )/2 ψ (1),n 2 h,j , ψ

(2),n+1 h,j = exp -i∆t(-mc 2 + V nuc.,j )/2 ψ (2),n 2 h,j .

Unsplit Pseudo-spectral (UPS) scheme. This scheme corresponds to the unsplit version of the SPS-scheme. From time t n to t n+1 the schemes reads, as follows. We set φ n 1 h = φ φ n k 1 e -iξp(x k 1 +a) e iξp(x j -a) .

• Then for all 0 j N (ps) -1 (1),n h,j + V nuc.,j φ

(2),n h,j .

• Followed, for all 0 j N (ps) -1, by φ MC-based scheme. From time t n to t n+1 the schemes reads 1. Source term integration and 0 j N (ps) -1 ψ

(1),n 1 h,j = exp -i∆t(mc 2 + V nuc.,j )/2 ψ

(1),n h,j , ψ

(2),n 1 h,j = exp -i∆t(-mc 2 + V nuc.,j )/2 ψ (2),n h,j .

2. We set φ n 1 h = φ

(1),n 1 h

, φ

(2),n 1 h T , where

φ (1),n 1 h = ψ (1),n 1 h + ψ (2),n 1 2 √ 2 , φ (2)
,n 1 h = ψ

(1),n 1 h -ψ

(2),n 1 h √ 2

• For any 1 j N (mc) -2 φ

(1),n 2 h,j = φ

(1),n 1 h,j-1 + i∆teA n 1 φ

(1),n 1 h,j , φ

(2),n 1 h,j+1 -i∆teA n 1 φ

(2),n 1 h,j .

• Then, we get

ψ (1),n 2 1 = φ (1),n 2 h + φ (2),n 2 h √ 2 , ψ (2)
,n 2 h = φ

(1),n 2 h -φ

(2),n 2 h √ 2 .

3. Finally, we integrate one more time the source term (Strang splitting) for 0 j N (mc) -1 ψ (1),n+1 h,j = exp -i∆t(mc 2 + V nuc.,j )/2 ψ (1),n 2 h,j , ψ

(2),n+1 h,j = exp -i∆t(-mc 2 + V nuc.,j )/2 ψ (2),n 2 h,j .

Figure 1 :

 1 Figure 1: Parallelization principle of the FFT in the z-direction.

  , as function of time in [0, T f ]. The initial wavefunction is given by

Figure 2 :

 2 Figure 2: Electric potentials A(t).

Figure 3 :

 3 Figure 3: Comparison: bispinor first component at final time T f , computed with MC-based scheme with N (mc) = 54799 gridpoints, and SPS-method (Left) and UPS-method (Right) with N (ps) = 275, 1097, 5481 gridpoints.

Figure 5 :

 5 Figure 5: Interaction potential V nuc. .

8 FFTFigure 6 :

 86 Figure 6: Comparison: bispinor first (Top-left) and second (Top-right) components at final time T f , computed with MC-based scheme with N (mc) = 1.37 × 10 5 gridpoints, and FFT-method with N (fft) = 1.37 × 10 4 gridpoints. (Bottom) Spatial spectrum of ψ 1 (T f , •) in logscale: k, F x ψ 1 (T f , k)

8 FFTFigure 7 :

 87 Figure 7: Comparison: bispinor first (Top-left) and second (Top-right) components at final time T f , computed with MC-based scheme with N (mc) = 1.37 × 10 5 gridpoints, and FFT-method with N (fft) = 1.37 × 10 5 gridpoints. (Right) Spatial spectrum of ψ 1 (T f , •) in logscale: k, F x ψ 1 (T f , k) 2 .

Figure 8 :

 8 Figure 8: Comparison: bispinor components at final time T f , computed with MC-based scheme with N (mc) = 54800 gridpoints, and SPS-method with N (sps) = 401 and 937 gridpoints.

Figure 9 :

 9 Figure 9: (Top-left) Initial spinor density. (Top-right) Interaction potential (Bottom-left) Spinor density at final time. (Bottom-right) Spinor density at final time: zoom in the central nucleus.

4

 4 i=1 |ψ i (0, •, •)| 2 , and potential in Fig. 11 (Left) and (Right).The overall computational domain is (-64, 64) 2 , the time step ∆t = 1/mc 2 = 2.28042 ×

Figure 10 :

 10 Figure 10: (Top-left) Efficiency T /pT p , (Top-right) CPU-time as a function of the number of processors (p = 1, 4, 16, 64), (Bottom) CPU time for FFT-method with N z = 524288/2 2p , with p = 0, • • • , 3, and MC-based method with 524288 dof.

Figure 11 :

 11 Figure 11: (Left) Density of the initial spinor. (Right) Interaction potential.

Figure 13 :

 13 Figure 13: (Left) Density of the initial spinor. (Right) Interaction potential.

Figure 14 :

 14 Figure 14: Density at time: (Left) T = 0.05521 and (Right) T f = 0.1142.

Figure 15 :

 15 Figure 15: (Left) CPU-time as function of the number of processors. (Right) Efficiency T 1 /pT p for the MC-based and FFT methods: p = 1, 2, 4, 16, 64 processors.

Figure 16 :

 16 Figure 16: Density at initial time.(Left)x, y, ρ 0 (x, y) , for (x, y) ∈ (-8, 8) 2 . (Right) x, z, ρ 0 (x, 0, z) , for (x, z) ∈ (-8, 8) × (-512, 512) .

Figure 17 :

 17 Figure 17: Electric potential. (Left)x, y, A z (T f , x, y, 0) , for (x, y) ∈ (-8, 8) 2 .(Right) x, z, A z (T f , x, 0, z) , for (x, z) ∈ (-8, 8) × (-512, 512) .

Figure 18 :

 18 Figure 18: Density at final time.(Left)x, y, ρ(x, y, 0, T f ) , for (x, y) ∈ (-8, 8) 2 . (Right) x, z, ρ(x, 0, z, T f ) , for (x, z) ∈ (-8, 8) × (-512, 512)

Figure 19 :

 19 Figure 19: (Left) Efficiency T 1 /pT p for the MC-and FFT-methods. (Right) CPU-time as function of number of processors (p = 1, 2, 4, 8, 16, 64, 256 processors).

NN

  (fft) /2-1 p=-N (fft) /2 e -icξp∆t N (fft) -1 k 1 =0 φ (j),n 1 h,k 1 e -iξp(x k 1 +ax) e iξp(x+ax j ) (fft) /2-1 p=-N (fft) /2 e icξp∆t N (fft) -1 k 1 =0 φ (j),n 1

• 1 N

 1 For any 1 j N (ps) -2[[∂ x ]] φ n 1 h = (ps) N (ps) /2-1 p=-N (ps) /2 iξ p N (ps) -1 k 1 =0φ n 1

• 2 [

 2 For any 1 j N (ps) -

( 2

 2 -c∆t[[∂ x ]] φ (1),n h,j + i∆t eA n φ (1),n h,j -V nuc.,j φ (1),n h,j -mc 2 φ + c∆t[[∂ x ]] φ (2),n h,j -i∆t eA n φ (2),n h,j + mc 2 φ

  a splitting of order at least 2, should be implemented, which reads

	ψ n+1 h	= P x	∆t 4	P y	∆t 2	P x	∆t 4	P y	∆t 2	P z (∆t)P y	∆t 2	P x	∆t 4	P y	∆t 2	P x	∆t 4	ψ n h .
	Basically, compared to

Table 1

 1

	# dof CPU-time MC CPU-time SPS CPU-time UPS CPU-time FFT
	54799	19.46	-	-	-
	5481	-	2.86	3.87	2.15
	1097	-	0.63	0.81	0.48
	275	-	0.13	0.19	0.10

, the CPU-time for solving the above problem, using the MC-based, FFT-, SPS-, UPS-schemes, and with a variable number of dof.

Table 1 :

 1 Computational time in second: MC-, SPS-, UPS-, FFT-schemes.

Table 2 :

 2 Computational time in second: MC-and SPS-schemes (MC-scheme with 54800 points) and SPSsolution.

	# dof CPU-time MC CPU-time SPS
	54800	39.52	87.63
	27400	-	42.1
	9364	-	13.62
	4682	-	6.71
	2341	-	3.35
	937	-	1.34
	401	-	0.58

parallel) than the MC-based method, it allows to select much larger space steps (less dof for given spatial domain, for fixed imposed ∆t) but still keeping a good accuracy, except for multiscale problems. In the latter case, a very fine resolution in space is indeed required, making the MC-method more attractive. This remark is also illustrated in Fig. 15, where it is observed that on a coarse mesh (N x = 256, N z = 4096) and still taking ∆t = 1/mc 2 , the FFT-method is more efficient than the MC-method when N x = 1024, N z = 16384, at least for a small number of processors.

Test 3.

In this last test, we consider the evolution of a wavepacket subject to an external dynamic and static fields in 3d. We again compare the performance of the MC-based and FFT-methods. The computational domain is (-8, 8) × (-8, 8) × (-512, 512), and the time