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

François Fillion-Gourdeau

Emmanuel Lorin
email: elorin@math.carleton.ca

Steve Mclean

Pseudospectral computational methods for the time-dependent Dirac equation in static curved spaces

Keywords: Dirac equation on curved space, pseudospectral approximation, strained graphene

Pseudospectral numerical schemes for solving the Dirac equation in general static curved space are derived using a pseudodifferential representation of the Dirac equation along with a simple Fourier-basis technique. Owing to the presence of non-constant coefficients in the curved space Dirac equation, convolution products usually appear when the Fourier transform is performed. The strategy based on pseudodifferential operators allows for efficient computations of these convolution products by using an ordinary fast Fourier transform algorithm. The resulting numerical methods are efficient and have spectral convergence. Simultaneously, wave absorption at the boundary can be easily derived using absorbing layers to cope with some potential negative effects of periodic conditions inherent to spectral methods. The numerical schemes are first tested on simple systems to verify the convergence and are then applied to the dynamics of charge carriers in strained graphene.

Introduction

The Dirac equation is one of the most important equations of Physics, giving a quantum description of relativistic spin-1/2 particles such as electrons and quarks [START_REF] Thaller | The Dirac equation[END_REF]. For this reason, it can be found in the theoretical description of many physical systems in nuclear physics, condensed matter physics, laser-matter interaction, cosmology, and many others. However, this partial differential equation is notoriously hard to solve, and one often has to resort to numerical methods for accurate and non-perturbative solutions. Motivated by quantum electrodynamics (heavy ion collision, pair production) [START_REF] Ackad | Calculation of electron-positron production in supercritical uranium-uranium collisions near the Coulomb barrier[END_REF][START_REF] Reinhardt | Theory of positron production in heavy-ion collisions[END_REF][START_REF] Gelis | Quark-antiquark production from classical fields in heavy-ion collisions: 1 + 1 dimensions[END_REF][START_REF] Wells | Asymptotic channels and gauge transformations of the time-dependent Dirac equation for extremely relativistic heavy-ion collisions[END_REF][START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF], by strong field physics (Schwinger's effect, intense laser-molecule interaction) [START_REF] Fillion-Gourdeau | Resonantly enhanced pair production in a simple diatomic model[END_REF][START_REF] Fillion-Gourdeau | Enhanced schwinger pair production in many-centre systems[END_REF][START_REF] Fillion-Gourdeau | Numerical computation of dynamical Schwinger-like pair production in graphene[END_REF][START_REF] Dunne | Catalysis of Schwinger vacuum pair production[END_REF][START_REF] Salamin | Relativistic high-power laser-matter interactions[END_REF], graphene modeling [START_REF] Fillion-Gourdeau | Time-domain quantum interference in graphene[END_REF][START_REF] Katsnelson | Chiral tunnelling and the Klein paradox in graphene[END_REF], tremendous efforts have been put these past two decades on the development of numerical methods for the computation of the time-dependent Dirac equation in flat Minkowski space. Real space methods such as Quantum Lattice Boltzman techniques [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][START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF][START_REF] Lorin | A simple and accurate mixed P 0 -Q 1 solver for the Maxwell-Dirac equations[END_REF], Galerkin methods [START_REF] Grant | Variational methods for Dirac wave equations[END_REF][START_REF] Fillion-Gourdeau | Galerkin method for unsplit 3-d Dirac equation using atomically/kinetically balanced B-spline basis[END_REF][START_REF] Ern | Éléments finis: théorie, applications, mise en oeuvre[END_REF], and pseudospectral or spectral methods [START_REF] Huang | A time-splitting spectral scheme for the Maxwell-Dirac system[END_REF][START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF][START_REF] Guo | Spectral and pseudospectral approximations using Hermite functions: Application to the Dirac equation[END_REF][START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF][START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Mocken | FFT-split-operator code for solving the Dirac equation in 2+1 dimensions[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF] were developed to efficiently and accurately solve the Dirac equation. Semi-classical regimes were considered in [START_REF] Wu | Gaussian beam methods for the Dirac equation in the semi-classical regime[END_REF] using Gaussian beams or Frozen Gaussian Approximations [START_REF] Chai | Frozen gaussian approximation for the Dirac equation in semi-classical regime[END_REF][START_REF] Swart | A mathematical justification for the Herman-Kluk propagator[END_REF], while the non-relativistic limit has been studied in several recent papers [START_REF] Bao | Numerical methods and comparison for the Dirac equation in the nonrelativistic limit regime[END_REF][START_REF] Bao | A uniformly accurate multiscale time integrator pseudospectral method for the Dirac equation in the nonrelativistic limit regime[END_REF]. The computational difficulties for solving the time-dependent Dirac equation include the fermion doubling problem, related to numerical dispersion [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF], and the Zitterbewegung [START_REF] Thaller | The Dirac equation[END_REF], resulting in highly oscillating solutions whose origin can be traced back to the presence of the mass term βmc 2 . Finally, drastic stability conditions can lead to numerical diffusion related to the finite wave propagation speed, the speed of light c.

Another numerical challenge for the Dirac equation shared by any other wave equation in real space solved on a truncated domain is the need of imposing special boundary conditions in order to avoid spurious wave reflections at the computational domain boundary. Therefore, the computational methods on truncated domains require non-reflective boundary conditions [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF][START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schroedinger equation[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Hammer | A dispersion and norm preserving finite difference scheme with transparent boundary conditions for the Dirac equation in (1+1)d[END_REF][START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF], absorbing or perfectly matched layers (PML) [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF][START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF][START_REF] Turkel | Absorbing PML boundary layers for wave-like equations[END_REF][START_REF] Zeng | The application of the perfectly matched layer in numerical modeling of wave propagation in poroelastic media[END_REF][START_REF] Tsynkov | Numerical solution of problems on unbounded domains. a review[END_REF], or the introduction of an artificial potential [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF]. On the other hand, Fourier-based methods applied on bounded domains naturally induce periodic boundary conditions, which can be problematic when dealing with delocalized wave functions. In this respect, the technique developed in [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with perfectly matched layers[END_REF][START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF] for the Dirac equation in flat space, where a spectral method is combined with PML, is an interesting alternative. One of the goals of this article is to extend this numerical scheme to the Dirac equation in curved space.

Recently, the Dirac equation in curved space-time has gained important interest in some applications such as condensed matter physics for describing the dynamics of charge carriers in deformed 2D Dirac materials [START_REF] Cortijo | Effects of topological defects and local curvature on the electronic properties of planar graphene[END_REF][START_REF] Cortijo | Electronic properties of curved graphene sheets[END_REF], as well as astrophysics for fermion tunneling in black holes [START_REF] Kerner | Fermions tunnelling from black holes[END_REF][START_REF] Criscienzo | Fermion tunneling from dynamical horizons[END_REF][START_REF] Li | Dirac particles tunneling from BTZ black hole[END_REF][START_REF] Chen | Hawking radiation of Dirac particles via tunnelling from rotating black holes in de Sitter spaces[END_REF]. In its discrete version, it has been considered from the lattice Boltzmann technique point of view [START_REF] Succi | Quantum lattice Boltzmann is a quantum walk[END_REF][START_REF] Flouris | Confining massless Dirac particles in two-dimensional curved space[END_REF][START_REF] Debus | Shifted Landau levels in curved graphene sheets[END_REF] and as continuous limits of quantum walks [START_REF] Di Molfetta | Quantum walks as massless Dirac fermions in curved space-time[END_REF][START_REF] Arrighi | Quantum walking in curved spacetime[END_REF][START_REF] Mallick | Simulating Dirac hamiltonian in curved space-time by split-step quantum walk[END_REF]. However, the literature on numerical methods for this equation is scarse. Therefore, this article is an attempt to fill this gap. In particular, we present a numerical method based on the pseudodifferential representation [START_REF] Taylor | Pseudodifferential Operators[END_REF] of the Dirac equation in curved space in combination with perfectly matched layers (PML). In the pseudodifferential representation, it is possible to efficiently use Fourier-based methods, even though the Dirac equation under consideration has non-constant spatial coefficients. A similar methodology was successfully developed in [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with perfectly matched layers[END_REF][START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF] for the Dirac equation in flat space and in [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods[END_REF] for the Gross-Pitaevskii equation. In these two cases, the pseudodifferential representation was used within some absorbing layers at the truncated domain boundary to implement wave absorption at the boundary using PML. Within the domain however, the scheme was an usual spectral numerical method. In curved space, as we will show below, it is beneficial to employ the pseudodifferential representation in the whole domain because convolution products do not appear explicitly. When combined with an implicit scheme for the time evolution, this allows for a Fourier-based method benefiting from spectral convergence and unconditional 2 -stability. In addition, the PML can be straightforwardly included, reducing the effects of the inherent periodic boundary conditions. This paper is organized as follows. In Section 2, we describe the Dirac equations under consideration. In Section 3, we present the type of absorbing (perfectly matched) layers for the Dirac equation. In Sections 4 and 5, we propose and analyze two numerical methods for approximating the Dirac equation in static curved space-times. Section 6 is dedicated to numerical experiments. We conclude in Section 7.

Dirac equation

This section is devoted to the presentation of the Dirac equations studied in this paper. We first recall the basics of the usual Dirac equation in flat space and then, we present its extension to curved space. Finally, we reformulate the latter in "Hamiltonian form", similar to the one in flat space but with space-dependent coefficients.

Dirac equation in flat space

The time-dependent Dirac equation in Cartesian coordinates reads [START_REF] Itzykson | Quantum Field Theory[END_REF]

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

H flat (t, x) = α • [-i∇ -eA(t, x)] + βm + I 4 V (t, x), (2)
where ψ(t, x) ∈ L 2 (R 3 ) ⊗ C 4 is the time t = x 0 and coordinate (x = (x 1 , x 2 , x 3 )) dependent four-spinor, 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 β, α = (α i ) i=1,2,3 are the Dirac matrices. In this work, the Dirac representation is used, where

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

σ 1 = 0 1 1 0 , σ 2 = 0 -i i 0 and σ 3 = 1 0 0 -1 , (4)
while I 2 is the 2 × 2 unit matrix. Note that natural units are used where = c = 1.

To simplify the notation further and to parallel the one used in the next section for the Dirac equation in curved space time, it is convenient to write the Hamiltonian as

H flat (t, x) = -iα • ∇ + βm + F flat (t, x), (5)
where the function F flat contains the contribution from the electromagnetic potential:

F flat (t, x) := -eα • A(t, x) + I 4 V (t, x). (6)

Dirac equation in curved space

In this section, the generalization of the Dirac equation to curved space is presented. Throughout, Einstein's notation is assumed with the following conventions: greek indices relate to general curved space (characterized by the general metric g µν (x), where x := (t, x) denotes a space-time point), uppercase latin indices relate to flat space (characterized by the Minkowski metric η AB = diag(1, -1, -1, -1)) while lowercase latin indices are summed over spatial coordinates only (g ij (x) for i, j = 1, 2, 3).

The extension of the Dirac equation to curved space follows by imposing general covariance under arbitrary coordinate transformations. In covariant notation and general background space, the Dirac equation takes the form [START_REF] Pollock | On the Dirac equation in curved space-time[END_REF]

iγ µ (x) [∂ µ + Ω µ (x) -ieA µ (x)] -m ψ(x) = 0, (7)
where A µ is the four vector electromagnetic potential. The generalized gamma matrices used in the Dirac equation define a Clifford algebra:

{γ µ (x), γ ν (x)} = 2g µν (x), (8)
where the notation {•, •} stands for the anticommutator and g µν is the metric characterizing the curved space. These matrices are a generalization of the Dirac matrices in flat space, which are not space-dependent and are related to the Minkowski metric as

{γ A , γ B } = 2η AB . (9)
The two sets of matrices are related via the tetrad formalism as

γ µ (x) = γ A e µ A (x), (10)
where e µ A (x) is the tetrad. The tetrads are used to link the metric in curved and flat spaces and thus, obey the property:

g µν (x) = e µ A (x)e ν B (x)η AB . (11)
The spinorial affine connection Ω µ (x) was introduced in the Dirac equation to preserve the covariance. It is given by

Ω µ (x) = - i 4 ω AB µ (x)σ AB , (12)
where σ AB = i[γ A , γ B ]/2 is the commutator of the "flat space" Dirac matrices while the spin connection is

ω AB µ (x) = e A ν (x) ∂ µ e νB (x) + Γ ν µσ (x)e σB (x) , (13)
where the Christoffel symbols Γ ν µσ (x) were introduced. It is also important to notice that in curved space, the usual 2 -norm is not preserved. Instead, denoting

ψ, ψ γ = ψ 2 γ = |g(x)|ψ † [γ 0 γ 0 (x)]ψd 3 x , (14)
where g(x) is the determinant of the metric, the norm ψ γ is preserved in time, see [START_REF] Leclerc | Hermiticity of the Dirac hamiltonian in gravitational field[END_REF][START_REF] Parker | One-electron atom as a probe of spacetime curvature[END_REF][START_REF] Parker | One-electron atom in curved space-time[END_REF]. To develop a numerical scheme, it is convenient to rewrite [START_REF] Fillion-Gourdeau | Enhanced schwinger pair production in many-centre systems[END_REF] in a form similar to (1) in flat space, i.e. in "Hamiltonian form". This is performed straightforwardly multiplying (7) by γ 0 (x) and using the anticommutation relation [START_REF] Fillion-Gourdeau | Numerical computation of dynamical Schwinger-like pair production in graphene[END_REF]. The Dirac equation in curved space can then be written as [START_REF] Parker | One-electron atom as a probe of spacetime curvature[END_REF] i∂

t ψ(t, x) = H(t, x)ψ(t, x), (15)
where H(t, x) is the Dirac Hamiltonian operator in curved space given by

H(t, x) = -i[g 00 (x)] -1 γ 0 (x)γ i (x) [∂ i + Ω i (x) -ieA i (x)] +[g 00 (x)] -1 γ 0 (x)m -I 4 [iΩ 0 (x) + eA 0 (x)] . (16)
Defining generalized Dirac matrices as

β(t, x) := [g 00 (x)] -1 γ 0 (x), (17)
α i (t, x) := [g 00 (x)] -1 γ 0 (x)γ i (x), (18)
the Dirac Hamiltonian becomes

H(t, x) = -iα(t, x) • [∇ + Ω(t, x) -ieA(t, x)] + β(t, x)m -I 4 [iΩ 0 (t, x) + eA 0 (t, x)] . (19)
At this point, it is convenient to simplify the notation further by defining a function F that allows for writing the Hamiltonian as

H(t, x) = -iα(t, x) • ∇ + β(t, x)m + F (t, x) , (20)
where

F (t, x) := -iα(t, x) • [Ω(t, x) -ieA(t, x)] -I 4 [iΩ 0 (t, x) + eA 0 (t, x)] .
This is the general form of the Dirac Hamiltonian in curved space-time. The functions β(t, x), α(t, x) and F (t, x) need to be determined a priori from the metric and/or from the electromagnetic field potential. Some explicit examples are presented below. It is well-known that the preceding Hamiltonian is not self-adjoint with respect to the covariant inner product when the metric is time-dependent [START_REF] Leclerc | Hermiticity of the Dirac hamiltonian in gravitational field[END_REF][START_REF] Parker | One-electron atom as a probe of spacetime curvature[END_REF][START_REF] Parker | One-electron atom in curved space-time[END_REF][START_REF] Huang | Hermiticity of the Dirac Hamiltonian in curved spacetime[END_REF], casting some doubts on the conservation of probability. This can be remedied by adding a new term in the Hamiltonian as [START_REF] Leclerc | Hermiticity of the Dirac hamiltonian in gravitational field[END_REF][START_REF] Huang | Hermiticity of the Dirac Hamiltonian in curved spacetime[END_REF]

H (t, x) = H(t, x) + i 2 ∂ t ln |g(x)|g 00 (x) . (21)
The new Hamiltonian H is now self-adjoint. The new term can be interpreted in many ways, such as in the pseudo-hermitian operator formalism [START_REF] Huang | Hermiticity of the Dirac Hamiltonian in curved spacetime[END_REF][START_REF] Gorbatenko | Solution of the problem of uniqueness and hermiticity of hamiltonians for Dirac particles in gravitational fields[END_REF][START_REF] Gorbatenko | Uniqueness and self-conjugacy of Dirac hamiltonians in arbitrary gravitational fields[END_REF] or as the time-dependence of the position eigenstates [START_REF] Huang | Hermiticity of the Dirac Hamiltonian in curved spacetime[END_REF]. However, the generality of these results are disputed by other authors, who proposed a different approach [START_REF] Arminjon | Basic quantum mechanics for three Dirac equations in a curved spacetime[END_REF]. In this work, we do not dwell into this controversy as it is outside the scope of this article. For the sake of simplicity, and throughout the rest of this article, we will assume the metric is time-independent, allowing us to write the Hamiltonian as

H(t, x) = -iα(x) • ∇ + β(x)m + F (t, x), (22)
F (t, x) = -iα(x) • [Ω(x) -ieA(t, x)] -I 4 [iΩ 0 (x) + eA 0 (t, x)] . ( 23
)
This Hamiltonian is now self-adjoint and is the starting point for the development of the numerical schemes. The main difference between the flat and curved space versions of the Dirac equation is twofold: 1) the Dirac matrices are space dependent functions α(x) and β(x), assumed here to be smooth, and 2) the function F (t, x) contains the contribution coming from the spin affine connection and the metric. As long as the latter are smooth enough, they do not lead to any particular computational issues. On the other hand, because the Dirac matrices are space dependent, it is obviously not possible to solve directly this equation with a Fourier-based method without avoiding convolution products. However, we can rewrite the equation in pseudodifferential form as follows in R 3

i∂ t ψ(t, x) = -iα(x) • F -1 x iξF x {ψ}(t, ξ) + β(x)m + F (t, x) ψ(t, x),
which will be the ground of the proposed methodology. In this last equation, F x {•}(t, ξ) is the Fourier-transform operator on spatial coordinates and ξ is the transform variable. In this formulation, the Dirac partial differential equation becomes an integral equation where the derivative is expressed through its spectral representation. Mathematically, the Dirac equation [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] along with the Hamiltonian (20) is a hermitian non-strictly hyperbolic system of equations [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]. In principle, other numerical methods such as finite volumes or Galerkin could be used. An example of a finite volume discretization in 1D and its interpretation as a lattice Boltzmann method and quantum walk can be found in [START_REF] Succi | Quantum lattice Boltzmann is a quantum walk[END_REF].

Absorbing Layers

From a practical point of view, the time-dependent Dirac equation is considered on a bounded truncated physical domain denoted by D Phy . The pseudospectral method used to solve the Dirac equation naturally induces periodic boundary conditions on a bounded domain. What follows, is a general strategy which simultaneously i) avoids/reduces artificial wave reflections at the domain boundary, ii) limits the transfer of the wave from one side to the opposite one by periodicity. To reach this goal, we add a layer D PML surrounding D Phy , and stretch the coordinates in all the directions. The overall computational domain is next defined by: D = D Phy ∪ D PML . We refer to [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for the construction of PMLs for quantum wave equations and more specifically to [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] for the derivation and analysis of PMLs for the Dirac equation. Here, we outline the main features of this technique which is detailed in [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with perfectly matched layers[END_REF]. The first step of the implementation of PMLs is the following change of variables [START_REF] Zheng | A perfectly matched layer approach to the nonlinear Schrödinger wave equation[END_REF] involving only the space variable:

x i = x i + e iθ x i L * i Σ(s)ds, (24)
where θ ∈ (0, π/2), i = 1, 2, 3 and Σ is a function to be determined. We then define

S i (x i ) := 1 + e iθ i Σ(x i ),
the function Σ i being given by

Σ i (x i ) = Σ(|x i | -L i ), L i * |x i | < L i , 0, |x i | < L i * ,
where L i * < L i , and such that {x i ∈ R :

L i * |x i | < L i }
is the absorbing layer. The partial derivatives are then transformed into

∂ i → 1 S i (x i ) ∂ i = 1 1 + e iθ i Σ i (x i ) ∂ i , (25)
with Σ vanishing while S i is equal to 1 in D Phy , respectively. On truncated domains, we will consider the transformation [START_REF] Beerwerth | Krylov subspace methods for the Dirac equation[END_REF], and the associated new Hamiltonian

H PML = -iα(x) • T (x) + β(x)m + F (t, x), (26)
where

T := ([S 1 (x 1 )] -1 ∂ 1 , [S 2 (x 2 )] -1 ∂ 2 , [S 3 (x 3 )] -1 ∂ 3 ) T .
Several types of functions can be selected. An exhaustive study of the absorbing functions Σ is proposed in [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods[END_REF] for Schrödinger equations. Here are some examples:

Type I: Σ 0 (x i + δ i ) 2 , Type II: Σ 0 (x i + δ i ) 3 , Type III: -Σ 0 /x i , Type IV: Σ 0 /(x i ) 2 , Type V: -Σ 0 /x i -Σ 0 /δ i , Type VI: Σ 0 /(x i ) 2 -Σ 0 /(δ i ) 2 ,
where

δ i := L i -L i * .
From the pseudospectral point of view, the space-dependence of the coefficients (S i ) -1 again prevents the direct application of the Fourier transform on the equation, even in the flat-space case. In the latter case, the same pseudodifferential operator representation will still allow for combining the efficiency of the pseudospectral method and the computation of the non-constant coefficient Dirac equation (see [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with perfectly matched layers[END_REF]).

Time-discretization

This section is devoted to the time discretization of the Dirac equation, in flat and curved spaces. The main tool which is often used is the operator splitting technique. The complications related to the appearance of spatial differential operators are relegated to the next section where the spatial discretization is discussed.

Time-discretization for Dirac equation in flat space

In order to solve the Dirac equation in flat space, a natural approach consists in splitting the equation as follows. Here, an order-2 Strang splitting [START_REF] Strang | On the construction and comparison of difference schemes[END_REF] is considered, but higher order operator splittings can naturally be used. Let us consider a time interval from t n to t n+1 and assuming ψ(t n , •) is given, the exact formal solution to the Dirac equation in flat space is

ψ(t n+1 , x) = T exp -i t n+1 tn H flat (s, x)ds ψ(t n , x), (27)
where T stands for the time-ordered exponential. The latter can be approximated to second order by the following symmetric decomposition [START_REF] Suzuki | General decomposition theory of ordered exponentials[END_REF][START_REF] Suzuki | Fractal decomposition of exponential operators with applications to manybody theories and Monte-Carlo simulations[END_REF]:

ψ(t n+1 , x) = e -i ∆t 2 [βm+Fflat(tn+1/2,x)] e -∆tα•∇ e -i ∆t 2 [βm+Fflat(tn+1/2,x)] ψ(t n , x) + O(∆t 3 ), ( 28
)
where three exponential operators have been introduced. The first and third exponential operators can be evaluated analytically using the fact that

e i[βG(t,x)+α•G(t,x)] = I 4 cos(|G|) + i [βG(t, x) + α • G(t, x)] |G| sin(|G|), (29)
for arbitrary functions G and G, where we defined

|G| = G 2 (t, x) + G 2 (t, x).
A common method to deal with the second differential operator is to use the Fourier transform F x , as follows (still using an order 2-splitting).

Step by step, denoting t n * = t n + ∆t, we have

     ψ(t n+1/2 , x) = e -i ∆t 2 [βm+Fflat(tn+1/2,x)] ψ(t n , x), t ∈ [t n , t n+1/2 ], ψ(t n * , x) = F -1 x e iα•ξ∆t ψ(t n+1/2 , ξ) , t ∈ [t n , t n * ], ψ(t n+1 , x) = e -i ∆t 2 [βm+Fflat(tn+1/2,x)] ψ(t n * , x), t ∈ [t n+1/2 , t n+1 ]. (30)
Then, the second step in Fourier space can also be evaluated via [START_REF] Wu | Gaussian beam methods for the Dirac equation in the semi-classical regime[END_REF]. After discretizing spatially, the second equation is commonly solved using the Fast Fourier Transform (FFT), resulting in an operator splitting pseudospectral scheme [START_REF] Bao | An efficient and stable numerical method for the Maxwell-Dirac system[END_REF][START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF][START_REF] Mocken | FFT-split-operator code for solving the Dirac equation in 2+1 dimensions[END_REF][START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF][START_REF] Mocken | Quantum dynamics of relativistic electrons[END_REF].

Time discretization for Dirac equation in curved space

Following the same procedure as in flat space, described in the last section, an operator splitting approach can be introduced in curved space. For a time interval from t n to t n+1 and assuming that ψ(t n , •) is given, the exact formal solution to the Dirac equation in curved space is

ψ(t n+1 , x) = T exp -i t n+1 tn H(s, x)ds ψ(t n , x), (31)
where T stands for the time-ordered exponential. Again, the latter can be approximated to third order accuracy by a symmetric decomposition [START_REF] Suzuki | General decomposition theory of ordered exponentials[END_REF][START_REF] Suzuki | Fractal decomposition of exponential operators with applications to manybody theories and Monte-Carlo simulations[END_REF]:

ψ(t n+1 , x) = e -i ∆t 2 [β(x)m+F(tn+1/2,x)] e -∆tα(x)•∇ e -i ∆t 2 [β(x)m+F(tn+1/2,x)] ψ(t n , x) +O(∆t 3 ). ( 32
)
This has the same form as [START_REF] Braun | Numerical approach to solve the time-dependent Dirac equation[END_REF], except for the second exponential, which now has a spacedependent Dirac matrix. The latter makes the direct use of the Fourier transform, as in the flat case, challenging because the efficient FFT cannot be used and the computational complexity would be O(N 2 ), where N is the number of lattice points. Our strategy is to approximate the exponential operator to simplify the problem and to be able to exploit the pseudodifferential form of the derivative operator. Two different types of approximation are introduced, leading to two different classes of numerical schemes:

1. Crank-Nicolson approximation: This is obtained by formally approximating the exponential operator in its lowest order unitary form (1/1 Padé's approximant of the exponential function):

e -∆tα(x)•∇ = I 4 -∆t 2 α(x) • ∇ I 4 + ∆t 2 α(x) • ∇ + O(∆t 3 ). ( 33
)
Then, the operator splitting can be implemented by the following sequence

         ψ(t n+1/2 , x) = e -i ∆t 2 [β(x)m+F(tn+1/2,x)] ψ(t n , x), t ∈ [t n , t n+1/2 ], I 4 + ∆t 2 α(x) • ∇ ψ(t n * , x) = I 4 - ∆t 2 α(x) • ∇ ψ(t n+1/2 , x), t ∈ [t n , t n * ], ψ(t n+1 , x) = e -i ∆t 2 [β(x)m+F(tn+1/2,x)] ψ(t n * , x), t ∈ [t n+1/2 , t n+1 ]. (34)
As shown below, this yields a semi-implicit numerical scheme. 2. Polynomial approximation: This is obtained by formally approximating the exponential operator e -∆tα(x)•∇ , by a polynomial of the form Np q=0

a q P q (-∆tα(x) • ∇), (35)
where (P q ) q=0,••• ,Np are some set of orthogonal polynomials (such as Taylor series, second order differencing or Chebychev polynomials, for example [START_REF] Tal-Ezer | An accurate and efficient scheme for propagating the time dependent Schrödinger equation[END_REF]) and (a q ) q=0,••• ,Np are the polynomial coefficients, fixed to have an accurate approximation of the exponential.

Then, the operator splitting can be implemented by the following sequence

           ψ(t n+1/2 , x) = e -i ∆t 2 [β(x)m+F(tn+1/2,x)] ψ(t n , x), t ∈ [t n , t n+1/2 ], ψ(t n * , x) = Np q=0 a q P q (-∆tα(x) • ∇)ψ(t n+1/2 , x), t ∈ [t n , t n * ], ψ(t n+1 , x) = e -i ∆t 2 [β(x)m+F(tn+1/2,x)] ψ(t n * , x), t ∈ [t n+1/2 , t n+1 ]. ( 36
)
As shown below, this yields explicit numerical schemes which do not require the solution to a linear system.

These two strategies are well-known for the Schrödinger equation, but are adapted here to the Dirac equation in curved space and combined to the pseudodifferential representation of the derivative. In the next sections, some methods will be given based on the discretization of Eqs. ( 34) and [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] with spectral accuracy in space. The Strang splitting was introduced mostly to be consistent with the traditional splitting strategy used in flat space and to isolate the part of the equation that requires a special treatment with pseudodifferential operators. In principle, this is not mandatory and unsplit schemes could also be used in combination with an approximation of the time-ordered exponential. Nevertheless, notice that the overall second-order of accuracy in time is still preserved by the operator splitting.

Space-discretization for the Dirac equation in curved space

In this section, the spatial discretization of the Dirac equation in curved space is described, based on the two different approaches described in the last section. Throughout, we assume that the Dirac equation is solved in curved space on a truncated domain

[-a 1 , a 1 ] × [-a 2 , a 2 ] × [-a 3 , a 3 ] R 3 .
The cases of flat space or the one with time-dependent coefficients are standard and not recalled here. We define two sets of grid-points in real and Fourier spaces by

D (x) N = x k := x k 1 ,k 2 ,k 3 = (x 1 k 1 , x 2 k 2 , x 3 k 3 ) k∈O (x) N , D (ξ)
N = ξ p := ξ p 1 ,p 2 ,p 3 = (ξ 1 p 1 , ξ 2 p 2 , ξ 3 p 3 ) p∈O (ξ) N ,
where N , k and p are multi-indices (N := (N 1 , N 2 , N 3 ), k = (k 1 , k 2 , k 3 ) and p = (p 1 , p 2 , p 3 ), respectively and N i ∈ 2N * ), and with

O (x) N = k ∈ N 3 / (k i = 0, • • • , N i -1) i=1,2,3 , O (ξ) N = p ∈ N 3 / p i = - N i 2 , • • • , N i 2 -1 i=1,2,3
.

The set D (x)

N defines a mesh with equidistant point positions in each dimension with sizes (for i = 1, 2, 3)

x i k i +1 -x i k i = h i = 2a i /N i .
One can deduce that the discrete wavenumbers in Fourier space are given by (for i = 1, 2, 3)

ξ i p i = p i π/a i .
The wave function ψ(t, x) is discretized spatially by a projection onto the spatial mesh while ψ(t, ξ) is discretized on the momentum mesh. Thus, we denote by ψ n k the approximate wavefunction at time t n and position x k and by ψ n p the wave function in momentum space at time t n and momentum ξ p . The discrete wave functions ψ n k and ψ n p are related by the discrete Fourier transform pair:

ψ n p = N -1 k=0 ψ n k e -iξ p •(x k +a) , ψ n k = 1 N N/2-1 p=-N/2 ψ n p e iξ p •(x k +a) ,
where a = (a 1 , a 2 , a 3 ). Armed with this notation, we can write the partial discrete Fourier coefficients in each dimension as, for i = 1, 2, 3:

ψ n k|k i →p i = F i (ψ n k ) = N i -1 k i =0 ψ n k e -iξ i p i (x i k i +a i ) ψ n k = F -1 i ( ψ n k ) = 1 N i N i /2-1 p i =-N i /2 ψ n k|k i →p i e iξ i p i (x i k i +a i ) ,
where the notation k|k i → p i means that the index k i in the set k is replaced by the index p i and where the partial discrete Fourier transform operator in the i coordinate is denoted by F i (•).

In order to approximate the partial derivatives, we use pseudospectral approximations of the pseudodifferential representation of the derivative operators. This leads to the following approximate first-order partial derivatives:

∂ i ψ(t n , x k ) ≈ [[∂ i ]]ψ n k := 1 N i N i /2-1 p i =-N i /2 iξ i p i ψ n k|k i →p i e iξ i p i (x i k i +a i ) . ( 37
)
This is the spectral representation of the derivative which, under standard assumptions on the smoothness of the wave function, has a spectral accuracy [START_REF] Canuto | Spectral methods. Scientific Computation[END_REF][START_REF] Gottlieb | Spectral methods for hyperbolic problems[END_REF]. Another representation of the spectral derivative can be found in Appendix B in terms of the differentiation matrix.

In the next sections, we will exploit these relations to obtain accurate numerical schemes. This approach not only allows to select the spatial steps as large as wanted, but it also preserves the very high spatial accuracy, the parallel computing structure and the scalability of the split method developed in [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]. In practice, we use the FFT to implement the Discrete Fourier Transform (DFT). This strategy is now combined with the two time discretization proposed in the last section to solve the Dirac equation in curved space.

5.1. Numerical scheme I: Crank-Nicolson scheme A Crank-Nicolson scheme is obtained by discretizing [START_REF] Antoine | A review of transparent and artificial boundary conditions techniques for linear and nonlinear Schroedinger equation[END_REF] and by using the spectral operator for computing the derivative. The main idea then consists in approximating F -1

x iξF x {ψ}(t, ξ) by using the discrete operator [[∇]] defined in [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. This yields

ψ n+1/2 k = e -i ∆t 2 β k m+F n+1/2 k ψ n k , t ∈ [t n , t n+1/2 ], G k ψ n * k = G k ψ n+1/2 k , t ∈ [t n , t n * ], ψ n+1 k = e -i ∆t 2 β k m+F n+1/2 k ψ n * k , t ∈ [t n+1/2 , t n+1 ], (38)
where we defined

α k = α(x k ), β k = β(x k ) and F n k = F (t n , x k ).
We also introduced the operators

G k := I 4 + ∆t 2 α k • [[∇]], (39)
G k := I 4 - ∆t 2 α k • [[∇]], (40)
for convenience. An explicit example of this procedure for a given 2D metric can be found in Appendix A.

Then, the numerical solution can be obtained by implementing the following algorithm. From time t n to t n+1 , and adding the PML T (see [START_REF] Bauke | Accelerating the Fourier split operator method via graphics processing units[END_REF]), the 3-steps scheme [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF] explicitly reads:

• Step 1. From t n to t n+1/2 such that t n+1/2 = t n + ∆t/2, with initial data ψ n k :

ψ n+1/2 k = e -i ∆t 2 β k m+F n+1/2 k ψ n k . (41)
To perform this step, an exact or approximate expression of the exponential operator is required. This can be achieved in different ways, depending on the metric chosen and the form of β and F . The first thing to note here is that both β and F are 4-by-4 matrices. One possibility is then to use one of the numerical techniques described in [START_REF] Moler | Nineteen dubious ways to compute the exponential of a matrix, twenty-five years later[END_REF] to compute the exponential of the matrix.

• Step 2. The second step, given by

G k ψ n * k = G k ψ n+1/2 k , (42)
can be written as a linear system of equations, by using the discrete pseudospectral representation of the derivative. Naively, one would construct the matrix G k explicitly as done in Appendix C, and solve the corresponding linear system. However, this procedure has the same computational complexity as the evaluation of convolution products. The computational efficiency of this step can be improved significantly by using a Krylov iteration solver (GMRES, conjugate gradient). This follows a technique developed before and we refer the reader to [START_REF] Antoine | Perfectly Matched Layer for computing the dynamics of nonlinear Schrödinger equations by pseudospectral methods[END_REF] for more details.

• Step 3. The third step is given by

ψ n+1 k = e -i ∆t 2 β k m+F n+1/2 k ψ n * k . ( 43
)
The same technique as in Step 1 can be used to evaluate the matrix exponential.

In flat space, it is well-known that the 2 -norm of the 4-spinor must be preserved, while in static curved space, the 2 γ -norm is preserved.

Proposition 5.1. Assume that (α i,n k ) i=1,2,3 are hermitian, and that F is a bounded function. The numerical scheme (41)-( 43) is unconditionally 2 -stable, and preserves the 2 -norm in flat space.

Proof. The proof is straighforward, as it mainly relies on i) the hermitivity of the Dirac matrices (α i k ) † = α i k , with i = 1, 2, 3, for steps [START_REF] Tsynkov | Numerical solution of problems on unbounded domains. a review[END_REF], [START_REF] Antoine | Computational performance of simple and efficient sequential and parallel Dirac equation solvers[END_REF], ii) the definition of G and G. Using the matrix representation of the derivative given in Appendix B, we start from the hermitian transpose of the differentiation matrices:

A i k i k i = - 1 N i N i /2-1 p i =-N i /2 iξ i p i e -iξ i p i (x i k i -x i k i ) . (44)
Setting p i = -p i and using the fact that ξ i

-p i = -ξ i p i
, we obtain

A i k i k i = - 1 N i N i /2-1 -p i =-N i /2 iξ i -p i e -iξ i -p i (x i k i -x i k i ) , (45)
= 1 N i N i /2-1 -p i =-N i /2 iξ i p i e iξ i p i (x i k i -x i k i ) . (46)
This implies that, by construction:

[[∇]] = -[[∇]],
and as a consequence

G † = I -(∆t/2) α • [[∇]] † = I + (∆t/2)α • [[∇]].
As the coefficients of the Dirac equation are bounded in space and time, we trivially conclude about the 2 -stability. In flat space, α is constant and F is purely real the 2 -norm is trivially preserved (Steps 1 to 3 are unitary).

Numerical II: polynomial scheme

Polynomial schemes are obtained using the same procedure as for the Crank-Nicolson method, i.e. by discretizing [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] and by using the spectral operator for computing the derivative. Again, we approximate F -1

x iξF x {ψ}(t, ξ) by using the discrete operator [[∇]] defined in [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. This yields

ψ n+1/2 k = e -i ∆t 2 β k m+F n+1/2 k ψ n k , t ∈ [t n , t n+1/2 ], ψ n * k = Np q=0 a q P q (-∆tα k • [[∇]])ψ n+1/2 k , t ∈ [t n , t n * ], ψ n+1 k = e -i ∆t 2 β k m+F n+1/2 k ψ n * k , t ∈ [t n+1/2 , t n+1 ]. (47)
The first and second steps are exactly the same as in the Crank-Nicolson scheme (see [START_REF] Pinaud | Absorbing layers for the Dirac equation[END_REF]) and thus are not discussed here. The second step, on the other hand, is different because it does not require a solution of a linear system. The main challenge is in computing powers of the operator [[∇]]. This is performed by using FFTs, where the number of FFT pairs is given by the order of the polynomial. Every scheme of this form is explicit and thus, is a priori at best conditionally stable. As a matter of fact, numerical experiments often show an instability of the numerical solution. However, the Crank-Nicolson scheme derived above naturally requires the solution to a large linear system at each time iteration. In order to improve the efficiency while keeping a reasonable accuracy, we propose a 2-steps polynomial scheme with directional splitting. In the following, in order to simplify the presentation, we will assume that α(x) is of the form

α(x) = a(x) • α, (48)
where a(x) := (a 1 (x), a 2 (x), a 3 (x)) T and where (a i ) i=1,••• ,3 are space dependent scalar functions. Using directional splitting and the form [START_REF] Li | Dirac particles tunneling from BTZ black hole[END_REF] for the matrices, the scheme ( 47) is slightly modified to

ψ n+1/2 k = e -i ∆t 2 β k m+F n+1/2 k ψ n k , t ∈ [t n , t n+1/2 ], ψ n * 1 k = Np q=0 a q P q (-∆ta 1 k α 1 [[∂ 1 ]])ψ n+1/2 k , t ∈ [t n , t n * ], ψ n * 2 k = Np q=0 a q P q (-∆ta 2 k α 2 [[∂ 2 ]])ψ n * 1 k , t ∈ [t n , t n * ], ψ n * 3 k = Np q=0 a q P q (-∆ta 3 k α 3 [[∂ 3 ]])ψ n * 2 k , t ∈ [t n , t n * ], ψ n+1 k = e -i ∆t 2 β k m+F n+1/2 k ψ n * 3 k , t ∈ [t n+1/2 , t n+1 ]. (49)
Let us remark that PMLs can easily be included by simply replacing a i (x) by a i (x)/S i (x i ), for i = 1, 2, 3 in the scheme below. Steps 1. and 5. are identical to the Numerical Scheme I (5.1). The principle of the second step of Numerical Scheme II consists in approximating the evolution of each direction by a Taylor expansion and by diagonalizing the Dirac matrix. We denote by Λ = diag(1, 1, -1, -1) and Π i the transition matrices, such that α i = Π i ΛΠ i, † , for i = 1, 2, 3. In addition, in the i-direction, we set φ n k = Π i, † ψ n k . Then, the time evolution is approximated as follows (with n * 0 = n + 1/2):

φ n * i k = φ n * i-1 k -∆ta i k Λ([[∂ i ]]φ n * i-1 ) k + O(∆t), = a i k F -1 i 1 -i∆tΛξ i F i (φ n * i-1 k ) + 1 -a i k φ n * i-1 k + O(∆t) . (50)
The first line of this equation is a polynomial Taylor scheme and is usually unstable while the second line is just a re-writing of the first line. Stability is recovered when the first term on the right-hand-side is approximated by an exponential, as 1 -i∆tΛξ i = e -i∆tΛξ i + O(∆t 2 ), which is accurate to second order in time. Then, the scheme reads

φ n * i k = a i k F -1 i e -i∆tΛξ i F i (φ n * i-1 k ) + 1 -a i k φ n * i-1 k . (51)
Finally, to recover the wave function, we set ψ

n * i k = Π i φ n * i k .
We proceed similarly in the other directions. As it is not totally obvious, we next prove the consistency of Numerical Scheme II. Proposition 5.2. The Numerical Scheme II is consistent with [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF] .

Proof. The analysis of the consistency only requires a focus on one of the steps 2-4, as the other steps are similar or standard. From e -i∆tξ i p Λ = I 4 -i∆tξ i p Λ -∆t 2 (ξ i p ) 2 Λ 2 + O(∆t 3 ) , we can write the first term of (51) as 51) is written as

Ξ n * i k := F -1 i Π i I 4 -i∆tξ i Λ -∆t 2 (ξ i ) 2 Λ 2 Π i, † F i (ψ n * i-1 k ) + O(∆t 3 ), = ψ n * i-1 k -∆tα i ∂ i ψ n * i-1 k + ∆t 2 I 4 ∂ 2 i ψ n * i-1 k + O(∆t 3 ), since Π i ΛΠ i, † = α i . Thus, (
ψ n * i k = a i (x k )Ξ n * i-1 k + 1 -a i (x k ) ψ n * i-1 k (52) = ψ n * i-1 k -a i (x k )∆tα i ∂ i ψ n * i-1 k + a i (x k )∆t 2 I 4 ∂ 2 i ψ n * i-1 k + O(∆t 3 ), ( 53
)
which is equivalent to

i ψ n * i k -ψ n * i-1 k ∆t = -ia i (x k )α i ∂ i ψ n * i-1 k + ia i (x k )∆tI 4 ∂ 2 i ψ n * i-1 k + O(∆t 2 ). ( 54
)
We have proven that is consistent at order 1 in time, with

i∂ t ψ(t, x) = -ia i (x)α i ∂ x ψ(t, x),
and this transport-like equation is equivalent to the 2 to 4 steps in the operator splitting.

We now provide a stability result.

Proposition 5.3. Let us assume that 0 sup x a i (x k ) C with i = 1, 2, 3. Then, the Numerical Scheme II approximating [START_REF] Fillion-Gourdeau | A split-step numerical method for the time-dependent Dirac equation in 3-d axisymmetric geometry[END_REF], which is assumed to be well-posed, is unconditionally 2 -stable.

Proof. Again the proof relies on Step 2. Let us focus on one of the direction and assume first that C 1. As Π i is unitary, we trivially get :

Ξ n k 2 2 = ψ n k 2 2
, where | • | 2 denotes the 2 -norm in C 4 on spinor components. Consequently, this yields

Ξ n k 2 := h 3 N k=0 Ξ n k 2 2 1/2 = ψ n k 2 .
Moreover, denoting R{z} (resp. I{z}) the real (resp. imaginary) part of z and using ( 52), we get

ψ n * i k 2 2 = [a i (x k )] 2 Ξ n * i-1 k 2 2 + 1 -a i (x k ) 2 ψ n * i-1 k 2 2 +2a i (x k ) 1 -a i (x k ) R{Ξ n * i-1 k }R{ψ n * i-1 k } + I{Ξ n * i-1 k }I{ψ n * i-1 k } , [a i (x k )] 2 Ξ n * i-1 k 2 2 + 1 -a i (x k ) 2 ψ n * i-1 k 2 2 +2 a i (x k ) 1 -a i (x k )| Ξ n * i-1 k 2 ψ n * i-1 k 2 , ( 55
)
where the Cauchy-Schwarz inequality was used to obtain [START_REF] Mallick | Simulating Dirac hamiltonian in curved space-time by split-step quantum walk[END_REF]. Expanding (55) and using Ξ

n * i-1 k 2 2 = ψ n * i-1 k 2 2
, we easily obtain the inequality:

ψ n * i k 2 2 ψ n * i-1 k 2 2
. We deduce that ψ

n * i k 2 ψ n * i-1 k
2 , for i = 1, 2, 3. Thus, we obtain:

ψ n * k 2 ψ n k 2 .
The stability analysis from Steps 1 and 5 is straightforward, hence leading to

ψ n+1 k 2
ψ n k 2 . This concludes the proof for C 1. For C > 1, we proceed similarly, and conclude with the Gronwall's inequality.

We can extend the above result to a second-order time scheme, by replacing the second step (52) by

ψ n * i k = a i (x k )Ξ n * i-1 k + 1 -a i (x k ) ψ n * i-1 k + a i (x k )∆t 2 I 4 [[∂ 2 i ]]Ξ n * i-1 k . ( 56
)
The second order in time is obtained thanks to the addition of the rightmost term to the scheme [START_REF] Debus | Shifted Landau levels in curved graphene sheets[END_REF]. The stability occurs from the implicitation of the correction/anti-diffusion term

a i (x k )∆t 2 I 4 [[∂ 2 i ]]ψ.

Computational complexity analysis

This paragraph is dedicated to the analysis of the computational complexity of the presented methods. In particular, we compare the complexity with i) the direct implicit method based on the direct application of the FFT on the equation, involving spatial convolution products, ii) the Crank-Nicolson scheme and the polynomial scheme. Each time iteration requires on a N -point grid the following operations:

• The implicit direct method requires O(N 2 + N ν ) operations, with ν > 1. The first term is due to the approximation of the convolution products by standard quadrature rules, and the second term comes from the numerical computation of the solution to the linear system.

• The Crank-Nicolson scheme needs O(N log N + N ν ) operations, for ν > 1. The first term is related to the approximation of the convolution products by FFT-method, while the second one comes from the solution to the linear system.

• The polynomial scheme implies O(N log N ) operations since the convolution products are computed by an FFT-method, the rest of the scheme being linear.

Numerical experiment for the curved-space Dirac equation

Several specific examples will be considered, which basically correspond to different metrics. The objective is to demonstrate how simple and efficient is the proposed methodology. We will start with simple one-dimensional tests, then will consider more elaborated twodimensional physical configurations. In 1-D and 2-D, the Dirac equation in curved space is slightly modified compared to 3-D. In particular, the dimension of Dirac matrices is 2-by-2, instead of 4-by-4. The numerical schemes described previously can be straightforwardly adapted to these cases.

Static spacetime

We consider the metric ds 2 = e 2Φ(x) dt 2 -e 2Ψ(x) dx 2 , such that Φ and Ψ are two spacedependent functions. This leads to the following one-dimensional Dirac equation [START_REF] Koke | Dirac equation in 2-dimensional curved spacetime, particle creation, and coupled waveguide arrays[END_REF]

i∂ t ψ = -ie Φ(x)-Ψ(x) σ x ∂ x + Φ (x) 2 ψ + e Φ σ z mψ .
Numerical Experiment 1. We propose a benchmark with Ψ(x) = e -10 -2 x 2 , Φ(x) = e -5×10 -3 x 2 , and for φ 0 (x) = exp(-x 2 /2 + ik 0 x), where k 0 = 5 and c = 1. The computational domain is [-5, 5], while the discretization parameters are set through: ∆t = 5 × 10 -4 and N 1 = 18027. We compare the proposed method, with a real-space method, at CFL = 0.99, which degenerates into the Quantum-Boltzmann method for flat space [START_REF] Fillion-Gourdeau | Numerical solution of the timedependent Dirac equation in coordinate space without fermion-doubling[END_REF]. The secondorder splitting implicit pseudospectral method (Numerical Scheme I) is implemented by using GMRES [START_REF] Saad | GMRES -A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems[END_REF]. We plot x → exp(Φ(x) -Ψ(x)), corresponding to a velocity field, in Fig. 1 (Left). We report the real and imaginary parts of the first component ψ 1 in Fig. 1 (Middle, Right) at time T = 0.5, corresponding to 1000 time iterations. Unlike the real space method for CFL = 1, the pseudospectral method is linearly stable. As an illustration, we compare in Fig. 2 (Left, Middle) the real and imaginary parts of ψ 1 on a coarse grid (h = 1.1 × 10 -2 , ∆t = 10 -2 ) and fine grid (h = 5.5 × 10 -4 , ∆t = 5 × 10 -4 ). Finally in Fig. 2 (Right), we report in logscale the 2 -norm of the error (

ψ h (T, •) -ψ ref,h (T, •) 2
) as a function of the space size h.

Numerical Experiments 2. We compare the pseudospectral method in flat and curved spaces. For the curved space case, we select Ψ(x) = cos(x/10)e -10 -2 x 2 , Φ(x) = e -10 -2 x 2 , with φ 0 (x) = exp(-x 2 /2 + ik 0 x), where k 0 = 5 and we again take c = 1. The computational domain is [-5, 5], with discretization parameters ∆t = 5×10 -4 and N 1 = 20001. The second order implicit splitting pseudospectral method is again solved by using GMRES [START_REF] Saad | GMRES -A Generalized Minimal Residual Algorithm for Solving Nonsymmetric Linear Systems[END_REF]. We plot x → exp(Φ(x) -Ψ(x)), corresponding to the velocity field in Fig. 3 (Left), and the real and imaginary parts of the first component ψ 1 in Fig. 3 (Middle, Right) at time T = 1, corresponding to 2000 time iterations. This illustrates the effect of the spatial curvature on -5 the solution to the Dirac equation.

Numerical Experiments 3. We now consider a two-dimensional Dirac equation in curved space, defined by the metric ds 2 = e 2Φ(x) dt 2 -e 2Ψ(x) dx 2 , such that Φ, Ψ are two spacedependent functions, with x = (x, y),

i∂ t ψ = -ie Φ(x)-Ψ(x) σ x ∂ x + ∂ x Φ(x) 2 + σ y ∂ y + ∂ y Φ(x) 2 ψ + e Φ(x) σ z mψ .
We assume that Φ(x) = e -10 -2 x 2 , Ψ(x) = e -5×10 -3 x 2 , with φ 0 (x) = e -x 2 /2+ik 0 •x , where k 0 = (5, 5) T and setting c = 1. The computational domain is [-5, 5] 2 . We report in Fig. 4 (Left) the velocity field on

D = [-5, 5] 2 , x → exp Φ(x) -Ψ(x) .
We here implement the Numerical scheme II. The numerical data are as follows: ∆t = 1.14 × 10 -4 , and N 1 = N z = 512. The initial data is a wavepacket

ψ 0 (x, z) = φ 1 (x, z), 0, 0, 0 T ,
where φ 1 (x, y) = e -(x 2 +y 2 )/2+5i(x-z) , which is plotted in Fig. 4 (Right). We report in Fig. 5 (resp. Fig. 6) the modulus of the first component (resp. real part of the first component) of the Dirac 4-spinor at times (in atomic units) t 1 = 0.57 × 10 -2 , t 2 = 1.14 × 10 -2 , t 3 = 2.28 × 10 -2 and t 4 = 4.56 × 10 -2 in the flat (Left) and curved (Right) spaces.

Massless Dirac particles in curved graphene

Charge carriers in graphene can be described theoretically by a 2-D Dirac equation in curved space-time, where the metric is related to the graphene sample deformation [START_REF] Cortijo | Electronic properties of curved graphene sheets[END_REF]. The dynamics of these charge carriers was recently studied numerically [START_REF] Flouris | Confining massless Dirac particles in two-dimensional curved space[END_REF][START_REF] Debus | Shifted Landau levels in curved graphene sheets[END_REF] with a lattice Boltzmann method. Using our numerical schemes, we now propose some tests with similar configurations. One-dimensional test: rippled graphene sample. This test is dedicated to the simulation of strained graphene, and a simple comparison with non-strained graphene (corresponding to flat space). More specifically, we consider a rippled graphene sheet parameterized by the coordinate transformation map [START_REF] Flouris | Confining massless Dirac particles in two-dimensional curved space[END_REF]. We set h(x) = a 0 cos(2πk 0 x/ ) and f (x) = h (x)

2 /2 = 2π 2 a 2 0 k 2 0 sin 2πk 0 x/ 2 / 2 , where a 0 and k 0 denote the amplitude and wave vector of the surface ripples, and is the length of the sheet. Moreover, in onedimension, the spatial part of the metric and the tetrad are given by

g(x) = 1 -f (x) 2 , e x x (x) = 1 1 -f (x)
.

The Dirac equation for modeling strained graphene with external electromagnetic potentials (A, V ) in 1-d then reads

∂ t ψ + σ x e x x ∂ x -iA x ψ = -iγ 0 (m -V )ψ . (57)
This equation can be solved using the numerical schemes developed in this article. It is easy to show that the following 2 γ -norm is conserved

ψ γ := 1 -f (x) |ψ(x)| 2 dx 1/2
. Indeed, we multiply ( 57) by (1 -f (x))ψ † , take the real part and integrate in space, and directly get:

d dt 1 -f (x) |ψ(x)| 2 dx = 0 .
Experiment 4. In this first experiment, we assume that the initial data is

ψ(0, x) = (1, i) T βe -βx 2 /2 / √ 4π .
Numerically, we take = 2, a 0 4 × 10 -1 , k 0 = 2, c 1, and = 5. Moreover, we fix A x (x) = V (x) = 5x. We plot in Fig. 7 (Right), the graph of x ∈ [-10, 10] → f (x). Numerically, we choose the discretization parameters to be ∆t = 10 -2 and h = 10 -2 . We report in Fig. 8 the density, defined by d

F (t, •) = |ψ 1 (t, •)| 2 + |ψ 2 (t, •)| 2 in
flat space, and the density d C (t, •) in curved space at different times t = 0.4, t = 0.8, t = 1.2 and t = 1.6. In Fig. 9 (Left), we report in logscale the 2 -and 2 γ -norms of the solution as a function of time iterations, in flat and curved space as well as illustrating the 2 -stability, and the 2 -norm conservation in flat space, and 2 γ -norm conservation in curved space. We also report the graph of convergence in Fig. 9 (Right.), with ∆t = 10 -5 , and T = 10 -1 (corresponding to 1000 time-iterations) and the computational domain is still [-10, 10]. We represent the 2 -norm error between a solution of reference ψ ref. (•, T ) (computed on a very fine mesh) and the approximate one ψ approx. (•, T ), computed with a mesh-size h of 1/2 i with i = 4, ..., 11.

Numerical Experiments 5.

A more severe test is performed with different physical data. The computational domain is given by [-5, 5] and the initial data is with β = 2. Numerically, we take ∆t = 10 -2 , and h = 10 -2 . Now, we take V (x) = 1/(|x|+1), A x = 10x 2 and a 0 = 0.4, = 10, k 0 = 5. The initial data, the potential, the functions G and f are reported in Fig. 10. We plot in Fig. 11, the density d F (t, •) in flat space, and the density d C (t, •) in curved space at different times t = 0.2, t = 0.4, t = 0.6 and t = 0.8. In Fig. 12, we report the 2 -norm of the solution as a function of time iterations, in flat and curved space, illustrating the 2 -stability (and 2 -norm preserving in flat space) of the proposed scheme. The straining on graphene is enhanced compared to the above setting.

ψ(0, x) = (1, i) T βe -βx 2 /2 / √ 4π ,
Numerical Experiments 6. The last experiment is dedicated to the use of PML in order to absorb the waves reaching the computational domain boundary. We compare the solution and its norm using and without using PML. The objective is to show that PML allow to circumvent the effect of periodic boundary conditions, by avoiding the propagation of nonphysical waves entering back into the physical domain from one side to the other. We refer to [START_REF] Antoine | A simple pseudospectral method for the computation of the time-dependent Dirac equation with perfectly matched layers[END_REF], for a detailed study of PML for the Dirac equation using the pseudospectral method presented in this paper. The objective in this example is not to construct the best PML possible (this would require a fine study of absorbing functions and their parameters), but rather to show how efficient these PML can be. In the following, we consider absorbing functions of type I (with σ 0 = 1 and θ = 0, see Section 3 for definitions), with a PML corresponding to 10% of the overall computational domain. The computational domain is given by [-4.5, 4.5] and the initial data is

ψ(0, x) = (1, i) T βe -βx 2 /2 / √ 4π ,
with β = 2. Numerically, we take ∆t = 10 -2 , and h = 10 -2 . Now, we take V (x) = A x (x) = 0, m = 0, and a 0 = 0.4, = 5, k 0 = 2. We report in Fig. 13, the density d F (t, •) in flat space, and the density d C (t, •) in curved space at different times t = 0.75, t = 1.5, t = 2.25 and t = 4, with (Left column) and without PML (Right column). In Fig. 14, we report the 2 -norm of the solution as a function of time iterations, in flat and curved space with and without PML. The latter illustrate the conservation of the norms, when using periodic boundary conditions without PML, and vice-et-versa.

This experiment shows that, although the proposed method naturally imposes perdiodic boundary conditions, its negative effects can be tackled thanks to Perfectly Matched Layers, which can very easily be implemented within the pseudospectral method.

Conclusion

We have derived and analyzed simple pseudospectral computational methods for solving the Dirac equation in curved space with perfectly matched layers at the computational domain boundary, and more generally for Dirac-like equations with non-constant coefficients. Interestingly, the proposed methods can easily be implemented from existing Fourier-based methods. Some numerical one-and two-dimensional experiments illustrating the properties of the numerical schemes were proposed. In a forthcoming paper, we will apply the developed methodology to an extensive study of strained graphene.

Appendix A. Example of space-discretization in the Crank-Nicolson scheme

As an example, let us consider a two-dimensional Dirac equation in curved space, defined by the metric ds 2 = e 2Φ(x) dt 2 -e 2Ψ(x) dx 2 , such that Φ and Ψ are two space-dependent functions. This leads to the following two-dimensional Dirac equation [START_REF] Koke | Dirac equation in 2-dimensional curved spacetime, particle creation, and coupled waveguide arrays[END_REF]

i∂ t ψ = -ie Φ(x)-Ψ(x) σ x ∂ x + ∂ x Φ(x) 2 + σ y ∂ y + ∂ y Φ(x) 2
ψ + e Φ(x) σ z mψ . We can rewrite this equation in the form ∂ t ψ = -e Φ(x)-Ψ(x) σ x ∂ x ψ -e Φ(x)-Ψ(x) σ y ∂ y ψ -e Φ(x)-Ψ(x) ∂ x Φ(x) 2 σ x + ∂ y Φ(x) 2 σ yie Φ(x) σ z m ψ .

We denote by ψ n k the approximate wavefunction at time t n and at x k . The numerical scheme reads:

ψ n+1 k -ψ n k ∆t = -e Φ(x k )-Ψ(x k ) σ x [[∂ x ]]ψ n+1 k -e Φ(x k )-Ψ(x k ) σ y [[∂ y ]]ψ n+1 k -e Φ(x k )-Ψ(x k ) ∂ x Φ(x k ) 2 σ x + ∂ y Φ(x k ) 2 σ y -ie Φ(x k ) σ z m ψ n+1 k .
Adding PML (see Section 3) to the numerical scheme formally leads to

ψ n+1 k -ψ n k ∆t = - e Φ(x k )-Ψ(x k ) S x (x k ) σ x [[∂ x ]]ψ n+1 k - e Φ(x k )-Ψ(x k ) S y (x k ) σ y [[∂ y ]]ψ n+1 k -e Φ(x k )-Ψ(x k ) ∂ x Φ(x k ) 2 σ x + ∂ y Φ(x k ) 2 σ y -ie Φ(x k ) σ z m ψ n+1 k .

Appendix B. Representation of the spectral derivative

Another representation of the spectral derivative can be obtained by taking the Fourier transform on the RHS of [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. One obtains

[[∂ i ]]ψ n k = 1 N i N i /2-1 p i =-N i /2 N i -1 k i =0 iξ i p i e iξ i p i (x i k i -x i k i ) ψ n k|k i →k i . (B.1)
This can be simplified further by noting that the sum on (p i ) i=1,2,3 can be performed explicitly. Therefore, the final result is that

[[∂ i ]]ψ n k = N i -1 k i =0 A i k i k i ψ n k|k i →k i , (B.2)
where the differentiation matrices are given by

A i k i k i = 1 N i N i /2-1 p i =-N i /2 iξ i p i e iξ i p i (x i k i -x i k i ) (B.3) = i π N i a i      2e iB i k i k i sin N i B i k i k i 2 e iB i k i k i -1 2 + i N i cos N i B i k i k i 2 e iB i k i k i -1      , (B.4)
where

B i k i k i := π a i (x k i -x k i ).
α i k S i k N i -1 k i =0 A i k i ,k i ψ n * k|k i →k i .
Re-arranging the sums and introducing Kronecker's symbol, the last expression can be written in the form of

G n+1/2 k ψ n * k = N -1 k =0 G n kk ψ n * k , (C.2)
where we define the matrix representation of G The RHS of Eq. (C.1) is simpler because it contains the initial data, which is a known quantity. Therefore, it is possible to use the spectral representation of the derivative in Fourier space as in Eq. [START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF]. This yields

G n+1/2 k ψ n+1/2 k := H n k , = ψ n+1/2 k - ∆t 2
i=1,2,3

α i k S i k N i N i /2-1 p i =-N i /2
iξ i p i ψ n k|k i →p i e iξ i p i (x i k i +a i ) .

With these results, the linear system has the form

G n ψ n * = H n , (C.4)
where G n is the matrix defined in (C.3) and the bold symbols represent vectors in real space, with components V = [V 0,0,0 , V 1,0,0 , • • • , V N 1 -1,N 2 -1,N 3 -1 ] T . Solving this linear system yields ψ n * k , the value of the wave function after the second steps of the Crank-Nicolson scheme. The main problem with this approach is the evaluation of the matrix G n , which is not efficient (the computational complexity is O(N 2 )), and the storing of the same matrix which can be problematic since it requires O(N 2 ) of memory.

Figure 1 :

 1 Figure 1: Experiment 1. (Left) Velocity field x → exp Φ(x) -Ψ(x) . (Middle) Real part of ψ 1 (•, T ) with real space and pseudospectral methods. (Right) Imaginary part of ψ 1 (•, T ) with real space and pseudospectral methods. The final time is T = 0.5.

Figure 2 :

 2 Figure 2: Numerical Experiments 1. (Left) Real part of ψ 1 (•, T ) on coarse and fine grids. (Middle) Imaginary part of ψ 1 (•, T ) on coarse and fine grids. (Right) 2 -norm error.

Figure 3 :

 3 Figure 3: Numerical Experiments 2. (Left) Velocity field x → exp Φ(x) -Ψ(x) . (Middle) Real part of ψ 1 (•, T ) for flat and curved spaces. (Right) Imaginary part of ψ 1 (•, T ) for flat and curved spaces. The final time is T = 1.

Figure 4 :

 4 Figure 4: Numerical Experiments 3. (Left) Velocity field. (Right) Initial wavefunction ψ 1 (0, •).

Figure 5 :

 5 Figure 5: Numerical Experiments 3. Modulus of the first component of the 4-spinor (Left) Curved space: at times t 1 = 0.57 × 10 -2 , t 2 = 1.14 × 10 -2 , t 3 = 2.28 × 10 -2 and t 4 = 4.56 × 10 -2 . (Right) Flat space: at times t 1 = 1.14 × 10 -2 , t 2 = 2.28 × 10 -2 and t 4 = 4.56 × 10 -2 .

Figure 6 :Figure 7 :

 67 Figure 6: Numerical Experiments 3. Modulus of the first real part of the component of the 4-spinor (Left) Curved space: at times t 1 = 0.57 × 10 -2 , t 2 = 1.14 × 10 -2 , t 3 = 2.28 × 10 -2 and t 4 = 4.56 × 10 -2 . (Right) Flat space: at times t 1 = 1.14 × 10 -2 , t 2 = 2.28 × 10 -2 and t 4 = 4.56 × 10 -2 .

Figure 8 :Figure 9 :

 89 Figure 8: Experiment 4. From Top-Left to Bottom-Right: density in flat and curved spaces at time t = 0.4, t = 0.8, t = 1.2 and t = 1.6.

Figure 10 :Figure 11 :

 1011 Figure 10: Numerical Experiments 5. (Left) Initial density and potential V . (Right) Graph of f .

Figure 13 :

 13 Figure 13: Numerical Experiments 6. Density in flat and curved spaces at time t = 0.75, t = 1.5, t = 2.25 and t = 4. (Left column) with PML. (Right column) without PML.

Figure 14 :

 14 Figure 14: Numerical Experiments 6. 2 -norm of the solution as a function of time iterations in flat and curved space, and 2 γ -norm of the solution as a function of time iterations in curved space. (Left) with PML. (Right) without PML.

  Appendix C. Explicit construction of the matrix in the Crank-Nicolson schemeThe second step, given explicitly by G as a linear system of equations, by using the discrete pseudospectral representation of the derivative. For the LHS of (C.1), we use the representation given in Eqs.[START_REF] Antoine | Absorbing boundary conditions for relativistic quantum mechanics equations[END_REF] and (B.2) to obtain G

  i ,k i δ kk |k i =k i . (C.3) Equation (C.2) is just the matrix-vector product with a matrix defined in (C.3).

  Figure 12: Numerical Experiments 5. 2 -norm of the solution as a function of time iterations in flat and curved space, and 2 γ -norm of the solution as a function of time iterations in curved space.

	0.9
	0.88
	0.86
	0.84
	0.82
	0.8
	0.78
	0	0.1	0.2	0.3	0.4	0.5	0.6	0.7	0.8

2, t = 0.4, t = 0.6 and t = 0.8.