The Landau-Zener transition and the surface hopping method for the 2D Dirac equation for graphene
aa r X i v : . [ m a t h . A P ] M a y The Landau-Zener transition and the surface hoppingmethod for the 2D Dirac equation for graphene ∗ Ali Faraj † and Shi Jin ‡ Abstract
A Lagrangian surface hopping algorithm is implemented to study the two dimensionalmassless Dirac equation for Graphene with an electrostatic potential, in the semiclassicalregime. In this problem, the crossing of the energy levels of the system at Dirac pointsrequires a particular treatment in the algorithm in order to describe the quantum transition–characterized by the Landau-Zener probability– between different energy levels. We first derivethe Landau-Zener probability for the underlying problem, then incorporate it into the surfacehopping algorithm. We also show that different asymptotic models for this problem derived in[25] may give different transition probabilities. We conduct numerical experiments to comparethe solutions to the Dirac equation, the surface hopping algorithm, and the asymptotic modelsof [25].
Keywords:
Dirac equation; Wigner transform; Semiclassical model; Band crossing; Landau-Zener formula; Surface hopping algorithm; Spectral methods.
Subject classifications:
We are interested in the description of the transport of electrons in a single graphene layer. Thismaterial is a two-dimensional flat monolayer of carbon atoms which displays unusual and interestingelectronic properties arising from the bi-conically shaped Fermi surfaces near the Brillouin zonecorners (Dirac points). The electrons propagate as massless Dirac Fermions moving with the Fermivelocity v F , which is 300 times smaller than the speed of light v F ≈ c ≈ m.s − , and theirbehavior reproduces the physics of quantum electrodynamics but at much smaller energy scale.Although this model has been studied for a long time, see [3] for a bibliography, it has remainedtheoretical until the work of [27] where the graphene was produced for the first time. After theseresults, the interest of researchers on this material has shown a remarkable increase includingapplications in carbon-based electronic devices [21] and numerical simulations, see e.g. [9] andreferences therein.In this paper, we will consider a model of a two-dimensional Dirac equation ([1, 3, 26]) of agraphene sheet in the presence of an external potential. This model consists of a small parameter h directly related to the Planck constant. We are interested is the design of an efficient numericalmethod–the surface hopping method– for the graphene Dirac equation in the semiclassical regimewhere h <<
1. In this regime, the solution of the Dirac equation is highly oscillatory thus a huge ∗ Research supported by NSFC grant 91330203, NSF grant 1114546 and the NSF grant RNMS-1107291 “KI-Net” † Institute of Natural Sciences, Shanghai Jiao Tong University, Shanghai 200240 P.R. China. Current address:Grenoble INP, ESISAR, 26902 Valence Cedex 9, France. Email: [email protected] ‡ Department of Mathematics, Institute of Natural Sciences, MOE-LSEC and SHL-MAC, Shanghai Jiao TongUniversity, Shanghai 200240 P.R. China; and Department of Mathematics, University of Wisconsin, Madison, WI53706, USA. Email: [email protected] h →
0, to obtain the physical observables (such as density, flux, and energy) with a computationalcost far less than a direct quantum simulation. When the gap between different energy levels is oforder one (the so-called adiabatic case), the Wigner measure technique provides a simple descriptionof the motion: it can be well approximated by a fully diagonalized system, one classical Liouvilleequation for each energy level [10]. However, in the graphene Dirac equation, the energy levels cross at the Dirac points, where non-adiabatic transfers are observed and the particles can tunnelfrom one band to the other. The standard Wigner approach then needs to be revised to describethe non-adiabatic phenomena.One of the widely used approaches to simulate the non-adiabatic dynamics is the surface hop-ping method initially proposed by Tully and Preston [34] as an efficient computational method togo beyond the classical Born-Oppenheimer approximation. This method is widely used in chem-istry and molecular dynamics, see for examples [5, 30, 33, 34]. The basic idea is to combine theclassical transports of the system on individual potential energy surfaces with instantaneous tran-sitions from one energy surface to the other. The transition rates for this band-to-band hoppingare given by the well known Landau-Zener formula [35]. From the mathematical point of view, thefirst rigorous analysis of non-adiabatic transfer using Landau-Zener formula dates back to Hage-dorn [11]. More recently, the Wigner measure techniques for separated energy levels have beenextended in [7] and [8] to systems presenting band crossing. The proof is based on microlocalanalysis to justify the Landau-Zener formula. Departing from the results in [7] and [8], a rigoroussurface hopping algorithm in the Wigner picture was proposed in [20] for time-dependent two-levelSchr¨odinger systems with conically intersecting eigenvalues and implemented numerically in theLagrangian formulation in [19]. The corresponding Eulerian numerical scheme was proposed in[18] by formulating the hopping mechanism as an interface condition which is then built into thenumerical flux for solving the underlying Liouville equation for each energy level.In the present article we give a Lagrangian surface hopping algorithm for the graphene Diracequation similar to the algorithm in [19]. First, it is a classical result that the Wigner trans-form leads to two decoupled classical Liouville equations for each energy level, under the adiabaticassumption [10]. At the Dirac points where non-adiabatic transition occurs, we first derive theLandau-Zener formula, and then incorporate it into the surface hopping algorithm. We also showthat a reduced asymptotic model developed in [25] could give incorrect transition probability. Wethen compare through several numerical examples the solutions of the Dirac equation (solved bythe time-splitting spectral method), the surface hopping algorithm and the asymptotic models of[25]. Our numerical results show that, when there is no wave interference, the surface hoppingalgorithm indeed gives the correct non-adiabatic transition at Dirac points, with a much greatercomputational efficiency compared with the simulations based on solving directly the Dirac equa-tion.The article is organized as follows: in section 2 we give the graphene Dirac equation andits semiclassical limit via the Wigner transform in the adiabatic case. We give some examplesof potentials to show that non-adiabatic transition is indeed possible. In section 3, we derivethe Landau-Zener transition probability for the graphene Dirac equation. The surface hoppingalgorithm is given in section 4. In section 5, we study the asymptotic models introduced in [25]and show that a reduced model could give the incorrect transition probability. Numerical resultsare given in section 6 for comparisons of different models. For reader’s convenience we give thetime-splitting spectral method of the Dirac equation in Appendix A.2
Quantum transport in graphene and the Wigner measure
We consider the description of the transport of electrons in a single graphene layer in the presenceof an external potential. Following [1, 3, 26], it is modelled by the two dimensional Dirac equation: (cid:26) i ~ ∂ t ψ = [ − i ~ v F σ ∂ x − i ~ v F σ ∂ x + qV ] ψ, t ∈ IR, x ∈ IR ψ (0 , x ) = ψ I ( x ) , x ∈ IR (2.1)where ~ is the reduced Planck constant, v F is the Fermi velocity, q is the elementary charge of theelectron and V ( x ) ∈ IR is the electric potential. The Pauli matrices σ , σ are given by: σ = (cid:18) (cid:19) , σ = (cid:18) − ii (cid:19) . The initial wave function ψ I ( x ) ∈ C is normalized such that: Z IR | ψ I ( x ) | dx = 1 (2.2)and, using the mass conservation, the wave function ψ ( t, x ) ∈ C satisfies Z IR | ψ ( t, x ) | dx = 1 (2.3)where t, x = ( x , x ) are the time and space variables respectively. We consider the system (2.1) inthe semiclassical regime. For this purpose, we rewrite the equations such that there remains onlyone dimensionless parameter h . Proceeding as in [31], we change the variables to dimensionlessvariables as follows t → t/T , x → x/L (2.4)and define u I ( x ) = Lψ I ( Lx ) , u ( t, x ) = Lψ ( T t, Lx ) (2.5)where L is a reference length and T = L/v F . We remark that u I ( x ) and u ( t, x ) are chosen suchthat the change of variable preserves the normalization (2.2) and (2.3). Plugging (2.5) into (2.1),and dividing by the reference energy mv F where m is the effective mass of the electron, one getsthe dimensionless graphene Dirac equation: (cid:26) ih∂ t u = [ − ihσ ∂ x − ihσ ∂ x + U ] u, t ∈ IR, x ∈ IR u (0 , x ) = u I ( x ) , x ∈ IR (2.6)where h = ~ mv F L (2.7)is a small dimensionless parameter and U ( x ) is the dimensionless potential defined by U ( x ) = qmv F V ( Lx ) . We will consider the solution of (2.6) in the limit h → × P ( x, ξ ) = B ( ξ ) + U ( x ) I, ( x, ξ ) ∈ IR x × IR ξ , (2.8)where I is the 2 × B ( ξ ) is given by B ( ξ ) = ξ σ + ξ σ = (cid:18) ξ − iξ ξ + iξ (cid:19) , ξ ∈ IR , U ∈ C ∞ ( IR , IR ) is such that ∀ β ∈ IN there is a constant C β > (cid:12)(cid:12) ∂ βx U ( x ) (cid:12)(cid:12) ≤ C β , ∀ x ∈ IR . An easy computation shows that the matrix B ( ξ ) has eigenvalues ±| ξ | with corresponding or-thonormal set of eigenvectors given by χ ± ( ξ ) = 1 √ (cid:18) , ± ξ + iξ | ξ | (cid:19) T . (2.9)Moreover, it holds P ∈ S ( IR × IR ) × and the Dirac equation (2.6) can be rewritten as: (cid:26) ih∂ t u h = P ( x, hD ) u h , t ∈ IRu h (0) = u hI (2.10)where P ( x, hD ) is the Weyl operator defined for u ∈ C ∞ ( IR ) by the integral: P ( x, hD ) u ( x ) = 1(2 π ) Z IR ξ Z IR y P (cid:18) x + y , hξ (cid:19) u ( y ) e i ( x − y ) .ξ dξdy , (2.11)and D = − i∂ x . The operator P ( x, hD ) is essentially self-adjoint on the Hilbert space H = L ( IR ) and the domain of its self-adjoint extension is H ( IR ) . Therefore − ih P ( x, hD ) generates astrongly continuous group of unitary operators solution to (2.10).The eigenvalues λ ± ( x, ξ ) = U ( x ) ± | ξ | of the symbol P ( x, ξ ) satisfy λ + ( x, ξ ) = λ − ( x, ξ ) at thecrossing set { ξ = 0 } ⊂ IR x × IR ξ . The semiclassical limit away from the crossing set was performedin [10] for systems of the form (2.10) with more general symbols and initial data u hI subject toadditional conditions. In [10], the authors show that the semiclassical limit for the different bandscan be treated separately where, for each band, the related eigenvalue plays the role of a scalarclassical Hamiltonian. We recall now some basic notions of the Wigner analysis involved in the semiclassical limit. For u, v ∈ L ( IR ), the Wigner transform is defined by: w h ( u, v )( x, ξ ) = 1(2 π ) Z IR u (cid:16) x − h y (cid:17) v (cid:16) x + h y (cid:17) e iξ.y dy and for u ∈ H , the 2 × W h [ u ] = (cid:0) w h ( u i , u j ) (cid:1) ≤ i,j ≤ . We denote by w h [ u ] = tr W h [ u ] the scalar Wigner transform of u . For any bounded sequence f h in H , there is a subsequence of W h [ f h ] which converges in S ′ . Such a limit W is called a Wignermeasure associated to f h . If f h admits only one Wigner measure, we shall denote it by W [ f h ]and set w [ f h ] = tr W [ f h ].Another important object is the classical flow φ ± t ( x, ξ ) = ( x ± ( t ) , ξ ± ( t )) corresponding to theeigenvalues λ ± , i.e. the solution to: ( ddt x ± ( t ) = ± ξ ± ( t ) | ξ ± ( t ) | , x ± (0) = x ddt ξ ± ( t ) = − ∂ x U ( x ± ( t )) , ξ ± (0) = ξ , (2.12)where ( x, ξ ) ∈ IR x × IR ξ . Indeed, the decoupled semiclassical limit is valid on a set of the phasespace which is stable under the flow φ ± t and where no band crossing occurs. More precisely, ifthere exists an open subset Ω ⊂ IR x × IR ξ such that:Ω ∩ { ξ = 0 } = ∅ and φ ± t (Ω) ⊂ Ω , ∀ t ∈ IR u hI in (2.10) has a Wigner measure W I such that w I = tr W I satisfies w I | Ω c = 0then using the results in [10], it holds that w h [ u h ], the scalar Wigner transform of the solution u h to (2.10), converges to w ( t, x, ξ ) = w ( t, x, ξ ) + w − ( t, x, ξ ) (2.13)where w ± ( t, ., . ) is the scalar positive measure on IR x × IR ξ solving the classical Liouville equations: (cid:26) ∂ t w ± ± ξ | ξ | ∂ x w ± − ∂ x U.∂ ξ w ± = 0 , IR t × Ω w ± (0 , ., . ) = tr (cid:0) Π ± W I (cid:1) , Ω w ± ( t, Ω c ) = 0 , t ∈ IR . (2.14)For all ( x, ξ ) ∈ Ω, Π ± ( ξ ) in (2.14) denotes the projection on the eigenspace related to the eigenvalue λ ± ( x, ξ ). Using formula (2.9), it can be expressed explicitly as follows:Π ± ( ξ ) = 12 (cid:18) I ± | ξ | B ( ξ ) (cid:19) . (2.15)Moreover, the density n h ( t, x ) = | u h ( t, x ) | converges to n ( t, x ) = Z IR w ( t, x, ξ ) dξ . It follows that for initial data u hI such that the bands are separated initially, i.e. the support of w I does not intersect with the crossing set { ξ = 0 } , the measure w is described by w I if the supportof w I is stable under the characteristics solution to (2.12) (in that case, the characteristics startingfrom the support of w I do not reach { ξ = 0 } ). U = 0In the case of the trivial potential, the solutions to (2.12) are given by: (cid:0) x ± ( t ) , ξ ± ( t ) (cid:1) = (cid:18) x ± ξ | ξ | t, ξ (cid:19) (2.16)for ξ = 0 and we can take Ω = IR x × IR ξ \ { ξ = 0 } . Then the system (2.14) writes: ( ∂ t w ± ± ξ | ξ | ∂ x w ± = 0 , IR t × { ξ = 0 } w ± (0 , ., . ) = w I, ± , { ξ = 0 } w ± ( t, x,
0) = 0 , t ∈ IR, x ∈ IR (2.17)where w I, ± = tr (cid:0) Π ± W I (cid:1) , ξ = 0 . Using (2.16) in (2.17), we get that: w ± ( t, x, ξ ) = w I, ± (cid:18) x ∓ ξ | ξ | t, ξ (cid:19) , ξ = 0 . Therefore, the density n h ( t, x ) = | u h ( t, x ) | converges to n ( t, x ) = Z IR w I, − (cid:18) x + ξ | ξ | t, ξ (cid:19) dξ + Z IR w I, + (cid:18) x − ξ | ξ | t, ξ (cid:19) dξ . In the present case, if the support of the initial scalar Wigner measure w I does not contain thepoint { ξ = 0 } , then no hopping occurs: the bands do not communicate and at any time the measure w is described by the separate evolution of the level characteristics.5 .2.2 Case U = αx , α ∈ IR \ { } In that case, the solutions to (2.12) are such that: ξ ± ( t ) = ξ − ( αt,
0) (2.18)and we can take Ω = IR x × IR ξ \ { ξ = 0 } . Then, system (2.14) becomes: (cid:26) ∂ t w ± ± ξ | ξ | ∂ x w ± − α∂ ξ w ± = 0 , IR t × { ξ = 0 } ,w ± (0 , ., . ) = tr (cid:0) Π ± W I (cid:1) , { ξ = 0 } , w ± ( t, x, ( ξ , , t ∈ IR, x ∈ IR , ξ ∈ IR . (2.19)Unlike in section 2.2.1, non-adiabatic transfer may occur. Indeed, if the support of w I does notintersect with the crossing set { ξ = 0 } but contains points of the form ( x, ξ ) where ξ = ( ξ ,
0) and ξ = 0, then the characteristic curve ξ ± ( t ), given by (2.18), starting from ( ξ , t = ξ α and a band-to-band transition will take place.In these conditions, the system (2.19) does not describe correctly the asymptotics h → W h [ u h ].The goal of this paper is to develop efficient semiclassical methods to compute the non-adiabatictransition between different bands. To quantify the transition rates, we propose in section 3 aLandau-Zener formula which can be justified theoretically by using the results in [8]. Indeed, in[8] the following symbol was considered: P ( x, ξ ) = A ( ξ ) + U ( x ) I, ( x, ξ ) ∈ IR x × IR ξ , (2.20)where I is the 2 × A ( ξ ) = (cid:18) ξ ξ ξ − ξ (cid:19) . (2.21)Using the unitary equivalence B ( ξ ) = R ∗ A ( ξ ) R , (2.22)where R = 1 √ (cid:18) i − i (cid:19) , (2.23)it follows that the two symbols have the same eigenvalues and P can be reduced to P afterconjugation with the matrix R in the treatment of band crossing. Remark 2.1.
The functions w ± ( t, x, ξ ) in (2.13) are the diagonal terms of the Wigner measure W ( t, x, ξ ) := W [ u h ( t )]( x, ξ ). Indeed, it was shown in [10] that w ± is given by w ± = tr (cid:0) Π ± W (cid:1) Ω (2.24)and that W is diagonal in the sense that W = Π + W Π + + Π − W Π − on IR t × Ω. Moreover, aneasy computation leads to:Π ± W Π ± = tr (cid:0) Π ± W (cid:1) Π ± = w ± Π ± on IR t × Ω . (2.25) In this section we will study the Landau-Zener transition for the Hamiltonian P ( x, hD ).6 .1 Classical flow around the crossing set As remarked in section 2, non-adiabatic transfer happens only when the characteristics reach thecrossing set { ξ = 0 } . Proposition 3.1 below says that such characteristics will exist as soon as thepotential has points x ∈ IR such that ∂ x U ( x ) = 0. We refer to [8] for the proof. Proposition 3.1.
Consider x ∈ IR such that ∂ x U ( x ) = 0 , then there exist two unique curves s ( x ± ( s ) , ξ ± ( s )) which are continuous in s in a neighborhood of and C for s = 0 and suchthat: ( dds x ± ( s ) = ± ξ ± ( s ) | ξ ± ( s ) | , x ± (0) = x dds ξ ± ( s ) = − ∂ x U ( x ± ( s )) , ξ ± (0) = 0 . (3.1)To illustrate the fact that a spectral transfer for the symbol P happens at the crossing set,we will come back momentarily to the potential U = αx introduced in section 2.2.2. For such apotential, problem (3.1) writes: dds x ± ( s ) = ± ξ ± ( s ) | ξ ± ( s ) | , x ± (0) = x dds ξ ± ( s ) = − α (cid:18) (cid:19) , ξ ± (0) = 0 , and ∀ x ∈ IR , its unique solution is given by: x ± ( s ) = x ∓ sgn( α ) (cid:18) | s | (cid:19) , ξ ± ( s ) = − α (cid:18) s (cid:19) . Plugging this solution in the projectors defined in (2.15), we get for s < + ( ξ ± ( s )) = 12 (cid:18) α )sgn( α ) 1 (cid:19) , Π − ( ξ ± ( s )) = 12 (cid:18) − sgn( α ) − sgn( α ) 1 (cid:19) and for s > + ( ξ ± ( s )) = 12 (cid:18) − sgn( α ) − sgn( α ) 1 (cid:19) , Π − ( ξ ± ( s )) = 12 (cid:18) α )sgn( α ) 1 (cid:19) . It is easy to see that the projectors Π + ( ξ ± ( s )) and Π − ( ξ ± ( s )) interchange one with the other whenthe characteristics pass through the crossing set { ξ = 0 } . In this section, we give a heuristic argument, similar to the one in [19] and [20], to derive theLandau-Zener formula for the Hamiltonian P ( x, hD ).In general, the region for non-adiabatic transfer is not restricted to the crossing set { ξ = 0 } ,since quantum transition occurs as long as the two energy levels are sufficiently close, in the caseof avoided crossing [11]. As in [19] and [20], we define this region to contain the points where thedistance between the eigenvalues λ ± ( x, ξ ) of P ( x, ξ ) is minimal. In our case | λ + ( x, ξ ) − λ − ( x, ξ ) | =2 | ξ | and, when considered along the characteristics solution to (2.12), the necessary condition forminimal gap is: (cid:0) | ξ ( s ) | (cid:1) ′ = 0 ⇔ ξ ( s ) .∂ x U ( x ( s )) = 0and the hopping surface is chosen as the set: S = { ( x, ξ ) ∈ IR ; ξ.∂ x U ( x ) = 0 } . (3.2)The heuristics follows by inserting the characteristics in the trace-free part of the symbol to obtainthe system of ordinary differential equations: ihψ ′ ( s ) = B ( ξ ( s )) ψ ( s ) . R defined in (2.23) and using equation (2.22), one arrives at thefollowing new system: ihψ ′ ( s ) = A ( ξ ( s )) ψ ( s ) . Assume the particles defined by the trajectory ((2.12)) are near a point ( x ∗ , ξ ∗ ) ∈ S (due totranslation invariance we assume particles are at ( x ∗ , ξ ∗ ) initially), then the Taylor expansiongives: x ( s ) = x ∗ ± ξ ∗ | ξ ∗ | s + O ( s ) ξ ( s ) = ξ ∗ − ∂ x U ( x ∗ ) s + O ( s ) . Ignoring the O ( s ) terms, the system becomes: ihψ ′ ( s ) = (cid:18) ξ ∗ − ∂ x U ( x ∗ ) s ξ ∗ − ∂ x U ( x ∗ ) sξ ∗ − ∂ x U ( x ∗ ) s − ξ ∗ + ∂ x U ( x ∗ ) s (cid:19) ψ ( s ) . After conjugation with the rotation matrix: (cid:18) cos θ sin θ − sin θ cos θ (cid:19) where θ is the angle such that (cos 2 θ, sin 2 θ ) = ∂ x U ( x ∗ ) | ∂ x U ( x ∗ ) | we get: i h | ∂ x U ( x ∗ ) | ψ ′ ( s ) = − s ξ ∗ ∧ ∂ x U ( x ∗ ) | ∂ x U ( x ∗ ) | ξ ∗ ∧ ∂ x U ( x ∗ ) | ∂ x U ( x ∗ ) | s ! ψ ( s )where ξ ∧ ζ = ξ ζ − ξ ζ for ξ, ζ ∈ IR . After conjugation with the matrix (cid:18) (cid:19) it follows: i h | ∂ x U ( x ∗ ) | ψ ′ ( s ) = s ξ ∗ ∧ ∂ x U ( x ∗ ) | ∂ x U ( x ∗ ) | ξ ∗ ∧ ∂ x U ( x ∗ ) | ∂ x U ( x ∗ ) | − s ! ψ ( s ) . If we set ε = h | ∂ x U ( x ∗ ) | and η = ξ ∗ ∧ ∂ x U ( x ∗ ) | ∂ x U ( x ∗ ) | , the system becomes: iεψ ′ ( s ) = (cid:18) s ηη − s (cid:19) ψ ( s )which is the well known Landau-Zener problem (see [35]) for which the transition probability is: T = e − πε η . This allows us to propose the following non-adiabatic transition rate at the point ( x ∗ , ξ ∗ ): T ( x ∗ , ξ ∗ ) = e − πh ( ξ ∗∧ ∂xU ( x ∗ ) ) | ∂xU ( x ∗ ) | . (3.3)8 .3 About the rigorous justification of the Landau-Zener formula As already noticed in section 2, the symbol P involved in the Dirac equation (2.10) satisfies theidentity P ( x, ξ ) = RP ( x, ξ ) R ∗ where P and R are given respectively by (2.20) and (2.23). Therefore, if u h denotes the solutionto (2.10), the function v h = Ru h satisfies the equation: ih∂ t v h = P ( x, hD ) v h (3.4)where P ( x, hD ) is defined as in (2.11). A Landau-Zener formula was obtained rigorously in [8]for the two-scale Wigner measure ˜ ν of the function v h (see [8], [20] for the definition of two-scaleWigner measures). This result provides a rigorous proof of our Landau-Zener formula. Indeed, if ν I denotes the two-scale Wigner measure of u h , it follows that the two-scale Wigner measure of v h is ˜ ν = Rν I R ∗ . Then, the Landau-Zener formula for ν I can be deduced from the Landau-Zenerformula for ˜ ν . Remark 3.2.
The Landau-Zener formula obtained in [8] writes T = e − πη | ∂xU ( x ) | where η = δξ ∧ ∂ x U | ∂ x U | . When the direction δξ is equal to ξ √ h , this corresponds to the Landau-Zener formula (3.3) that weobtained heuristically in section 3.2. We give here a semiclassical Lagrangian algorithm for the evolution of the diagonal terms of theWigner matrix W h ( t, x, ξ ) of the solution u h ( t ) to (2.10). The algorithm is adopted from themethod proposed in [19] for time-dependent two-level Schr¨odinger systems with conically inter-secting eigenvalues.Define the level populations: P h ± ( t ) = || Π ± ( hD ) u h ( t ) || H (4.1)where for u ∈ H , Π ± ( hD ) u is the function defined via its Fourier transform Π ± ( hξ )ˆ u ( ξ ) andˆ u = F u is given by: F u ( ξ ) = 12 π Z IR u ( x ) e − ix.ξ dx (it is clear from (2.15) that Π ± ( hD ) u = Π ± ( D ) u ). With similar computations as in [19], we obtain: P h ± ( t ) = Z IR x Z IR ξ w h ± ( t, x, ξ ) dxdξ , (4.2)where w h ± ( t, x, ξ ) = tr (cid:0) Π ± ( ξ ) W h ( t, x, ξ ) (cid:1) . (4.3)As it appears from the following equationΠ ± W h Π ± = w h ± Π ± , which is obtained in a similar way as (2.25), w h ± are the diagonal terms of the Wigner matrix W h .Up to a small remainder, the function w h ± can be written in terms of the scalar Wigner transform9f the level projections of the solution u h ( t ). Indeed, using Lemma 2.3 in [10], it holds for all t ∈ IR : w h ± ( t ) = w h (cid:2) u h ± ( t ) (cid:3) + o (1) (4.4)in D ′ (cid:16) IR x × ( IR ξ \ { ξ = 0 } ) (cid:17) when h →
0, where u h ± ( t ) is the function with Fourier transformgiven by ˆ u h ± ( t, ξ ) = Π ± ( hξ )ˆ u h ( t, ξ )1 ξ =0 and ˆ u h ( t ) is the Fourier transform of u h ( t ). Comparing the relations (4.3) and (2.24), it followsthat, in the situation of section 2 without band-crossing, the partial differential equation in (2.14)satisfied by w ± can be used to approximate the time evolution of w h ± . In the case of band crossing,the idea is to use (2.14) for the time evolution of w h ± as long as the classical trajectories solutionto (2.12) are away from the hopping surface S defined by (3.2). When a trajectory reaches apoint ( x ∗ , ξ ∗ ) ∈ S a non-adiabatic transfer of weight occurs between w h + and w h − with transitionprobability T ( x ∗ , ξ ∗ ) given by (3.3). The algorithm
1. Initial sampling: in this step, an appropriate sampling of the function w hI, ± defined in (4.8)is chosen. Specifically, a set of sampling points { ( x k , ξ k , j k ) ∈ IR x × IR ξ × {− , + } ; k = 1 , ..., N } are chosen with associated weights w k ∈ IR given by: w k = w hI,j k ( x k , ξ k ) .
2. Hopping transport: away from the set S , each particle ( x k , ξ k , j k ) is transported by theassociated classical flow φ j k t solution to (2.12). In other words, for t ≥ x k ( t ) , ξ k ( t )) = φ j k t ( x k (0) , ξ k (0)) , w k ( t ) = w k (0) . (4.5)If there exists t ∗ > x k ( t ∗ ) , ξ k ( t ∗ )) =: ( x ∗ k , ξ ∗ k ) ∈ S , the weight is reduced usingthe transition rate T ∗ = T ( x ∗ k , ξ ∗ k ) , where T ( x, ξ ) is given by (3.3). Moreover, in order to describe completely the non-adiabatictransfer, a new particle with index l > N is created. Specifically, for t > t ∗ ( x k ( t ) , ξ k ( t )) = φ j k t ( x k (0) , ξ k (0)) , w k ( t ) = (1 − T ∗ ) w k ( t ∗ )and the new particle is created for t > t ∗ ( x l ( t ) , ξ l ( t )) = φ j l t − t ∗ ( x ∗ k , ξ ∗ k ) , j l = − j k , w l ( t ) = T ∗ w k ( t ∗ ) .
3. Final reconstruction: at the final time t f >
0, there are M ≥ N points { ( x k , ξ k , j k ) ∈ IR x × IR ξ × {− , + } ; k = 1 , ..., M } with associated weights w k which are approximations to w hj k ( t f , x k , ξ k ). Then, using equation(4.2), the surface hopping approximations P hsh, ± ( t f ) of the level populations P h ± ( t f ) are givenby: P hsh, ± ( t f ) = M X k =1 w k δ ± j k ω k (4.6)where δ ij is the Kronecker symbol related to i and j , and ω k is an appropriate quadratureweight. 10 emark 4.1. We note that this surface hopping algorithm is subject to some restrictions. Firstonly the dynamics of the diagonal components of the Wigner matrix away from the crossing set arewell approximated. Second, there are possible interferences which are not captured if no particulartreatment is performed. Indeed, if w h + ( t, x, ξ ) and w h − ( t, x, ξ ) arrive at the same time at some point( x, ξ ) close to the crossing set then a transfer of weight using only the Landau-Zener formula (3.3)might give an incorrect approximation of the dynamics. To avoid these interferences, we make theassumption that the initial data u hI in (2.10) satisfiesΠ − ( hξ )ˆ u hI ( ξ ) = 0 (4.7)where ˆ u hI is the Fourier transform of u hI . Then, it follows from (4.4) that the diagonal terms of theinitial Wigner matrix w hI, ± ( x, ξ ) = tr (cid:0) Π ± ( ξ ) W h (cid:2) u hI (cid:3) ( x, ξ ) (cid:1) (4.8)satisfy in D ′ (cid:16) IR x × ( IR ξ \ { ξ = 0 } ) (cid:17) w hI, + = w h (cid:2) u hI (cid:3) + o (1) , w hI, − = o (1) (4.9)when h →
0. Therefore, for the potential U ( x ) = αx given in section 2.2.2, the condition (4.7)insures that no interferences occur. Indeed, if w h + ( t, x, ξ ) reaches the crossing set { ξ = 0 } at sometime t ∗ >
0, the particle corresponding to w h − ( t, x, ξ ) immediately moves away as it appears fromthe equation ξ ± ( t ) = − α ( t − t ∗ ,
0) (4.10)for the momentum part of the characteristics solution to (2.12). Similarly, in the case where thepotential has no stationary points, i.e. ∂ x U ( x ) = 0, ∀ x ∈ IR , equation (4.10) is replaced by ξ ± ( t ) = − Z tt ∗ ∂ x U (cid:0) x ± ( s ) (cid:1) ds and the condition (4.7) insures that no interferences occur as long as | t − t ∗ | is small enough. Remark 4.2.
To deal with this interference, a possible solution might be to use a hybrid methodfor the Dirac equation (2.10). Such a method was introduced in [17] for the Schr¨odinger equationand mixes a Gaussian beam method or a Liouville equation in the region where no-interferencesoccur and a complete quantum solver in the region where the phase information of the wavefunction is required. In our case, the second region corresponds to the hopping region. Since theregion which involves the quantum solver can be chosen very small, a hybrid method allows anadapted treatment of the interferences without increasing dramatically the numerical cost even ifthe semiclassical parameter tends to 0. In addition, the resulting algorithm is supposed to workregardless of the conditions in Remark 4.1 on the potential and the initial data. Another possiblesolution is to keep the off-diagonal entries, which contain information about the non-adiabatictransition in the Wigner transform, and then derive the semiclassical models for the entire Wignermatrix, see [25, 4]. In the next section, we study such a model.
As it was mentioned, the surface hopping method breaks down where there are interferences. Apossible remedy for this is to derive improved models using the Wigner transform for the en-tire
Wigner matrix, in which the off-diagonal entries contain non-adiabatic transition information[25, 4]. In this section, we study such a model obtained in [25], and compare the Landau-Zenertransition rate of this model with the one we derived in section 3. To serve our purpose we willonly consider the model for the potential studied in 2.2.2, U = αx , α ∈ IR \ { } , for which theWigner matrix does not depend on the variable x . It was shown in [25] how this asymptotic model11an describe non-adiabatic transfer as a quantum correction of the decoupled system (2.19).When the Wigner matrix does not depend on the variable x , the system (2.19) reads (cid:26) ∂ t w ± ± ξ | ξ | ∂ x w ± − α∂ ξ w ± = 0 , IR t × { ξ = 0 } w ± (0 , ., . ) = tr (cid:0) Π ± W I (cid:1) , { ξ = 0 } ; w ± ( t, x, ( ξ , , t ∈ IR, x ∈ IR , ξ ∈ IR . (5.1)By taking into account the scaling (2.4)(2.5), the asymptotic model in [25] is the following correctionof (5.1): ( ∂ t w ± ± ξ | ξ | ∂ x w ± − α∂ ξ w ± = ± α ξ | ξ | Im (cid:16) ξ + iξ | ξ | w i (cid:17) , IR t × { ξ = 0 } ∂ t w i − i Λ( ξ ) w i − α∂ ξ w i = − i α ξ ( ξ − iξ ) | ξ | ( w + − w − ) (5.2)where Λ( ξ ) = − | ξ | h − αξ | ξ | . In the limit h →
0, the function w i ( t, x, ξ ) approximates the off-diagonal terms of the Wignermatrix of the solution u h to (2.10). As it appears from the last equation in (5.2), the function w i depends on w + and w − . As a consequence, w i provides a coupling term at the r.h.s. of the firsttwo equations in (5.2).Using the method of characteristics, w i can be expressed as an explicit function of the difference w + − w − . Inserting this solution in the first two equations in (5.2), the following approximateequations for w ± , in the limit h →
0, are obtained in [25]: ∂ t w ± ± ξ | ξ | ∂ x w ± − α∂ ξ w ± = ∓ τ ( ξ )( w + − w − ) , IR t × IR x × IR ξ (5.3)where ξ = 0 can be considered as a small parameter. In the domain (cid:26) ξ ∈ IR such that 0 < | ξ | ≤ √ αh and | ξ | ≤ αh | ξ | (cid:27) , the coefficient τ ( ξ ) is given by: τ ( ξ ) = α ξ | ξ | (cid:18) π αξ ) − arctan ξ ξ (cid:19) (5.4)where sgn( x ) denotes the sign of x . As it will be the case in section 6.2 for our surface hoppingalgorithm, the set { ξ = 0 } is the important region for the hopping. Indeed, the following lemma(proof left to the reader) shows that the set { ξ = 0 } plays the role of an interface where, in thelimit ξ →
0, the solution to (5.3) will have a discontinuity.
Lemma 5.1.
The function ξ τ ( ξ , ξ ) tends in D ′ ( IR ) to | α | βδ ξ =0 when ξ → , where β = π . (5.5)In the case α >
0, we show below that non-adiabatic transfer is possible using the model (5.3)and give the corresponding transmission matrix at the interface. Define ω ± as the set: ω ± = ( IR t × IR x × IR ξ ) ∩ {± ξ > } and consider an initial condition in the upper band and localized at the right of the interface. Inother words, the initial conditions for (5.3) are w ± (0 , ., . ) = w I, ± where w I, − = 0 and w I, + is a function in C ∞ ( IR x × IR ξ ) which is independent of ξ and hassupport in ω + . To perform the limit ξ →
0, the following assumptions are required on the solutionto (5.3): 12
1. For all ξ = 0, w ± ∈ C ( IR t × IR x × IR ξ ). A
2. There is a constant C which does not depend on ξ such that: | ∂ x w ± ( t, x , ξ ) | ≤ C, ∀ ( t, x , ξ ) ∈ ω + ∪ ω − .A
3. There is a function w ± which is continuous on ω + and ω − such that w ± tends to w ± when ξ → ω + and ω − . A
4. For all t ∗ > x ∗ ∈ IR , the limitlim ( t,x ,ξ ) → ( t ∗ ,x ∗ , t,x ,ξ ) ∈ ω + w ± ( t, x , ξ ) (resp. lim ( t,x ,ξ ) → ( t ∗ ,x ∗ , t,x ,ξ ) ∈ ω − w ± ( t, x , ξ ))exists and is denoted w r ± ( t ∗ , x ∗ ) (resp. w l ± ( t ∗ , x ∗ )).The characteristics related to (5.3) are the classical trajectories ϕ ± t ( x I , ξ I ) = ( x ± ( t ) , ξ ± ( t )) solutionto: ( ddt x ± ( t ) = ± ξ ± ( t ) √ ξ ± ( t ) + ξ , x ± (0) = x Iddt ξ ± ( t ) = − α, ξ ± (0) = ξ I . (5.6)Due to the localization of the support of the initial data w I, + , only positive initial momenta ξ I > { ξ = 0 } at the time t ∗ = ξ I α . Moreover, due to the condition w I, − = 0, the behavior at the interface isdescribed by the characteristics corresponding to the upper band only. To be more precise, bydifferentiating the map s w ± ( s, x +1 ( s ) , ξ +1 ( s )), one gets by using (5.3): dds (cid:18) w + ( s, x +1 ( s ) , ξ +1 ( s )) w − ( s, x +1 ( s ) , ξ +1 ( s )) (cid:19) = − τ ( ξ + ( s )) M (cid:18) w + ( s, x +1 ( s ) , ξ +1 ( s )) w − ( s, x +1 ( s ) , ξ +1 ( s )) (cid:19) + F ( s ) (5.7)where ξ ± ( s ) = ( ξ ± ( s ) , ξ ) , M = (cid:18) − − (cid:19) and F ( s ) = 2 ξ +1 ( s ) | ξ + ( s ) | (cid:18) ∂ x w − ( s, x +1 ( s ) , ξ +1 ( s )) (cid:19) . Applying the Duhamel formula to (5.7), it follows: (cid:18) w + ( t, x +1 ( t ) , ξ +1 ( t )) w − ( t, x +1 ( t ) , ξ +1 ( t )) (cid:19) = e − R t τ ( ξ + ( s )) dsM (cid:18) w I, + ( x I , ξ I ) w I, − ( x I , ξ I ) (cid:19) + Z t e − R ts τ ( ξ + ( µ )) dµM F ( s ) ds . (5.8)Using (5.4), a direct computation shows that there is a constant C which does not depend on ξ such that: ∀ s, t ∈ IR, | Z ts τ ( ξ ± ( µ )) dµ | ≤ C . (5.9)Moreover it holds lim ξ → Z t τ ( ξ ± ( s )) ds = (cid:26) < t < t ∗ β if t > t ∗ . (5.10)In addition, since the lower band is not occupied initially and most of the band-to-band transferoccurs at the time t ∗ , it holds: ∀ < t < t ∗ , lim ξ → ∂ x w − ( t, x +1 ( t ) , ξ +1 ( t )) = 0 . (5.11)In order to simplify the presentation, the proof of (5.11) is postponed until the end of the presentdiscussion. Consider first the case t < t ∗ in (5.8). Using assumption A Z t e − R ts τ ( ξ + ( µ )) dµM F ( s ) ds −→ ξ → . (5.12)Combining assumption A ξ → (cid:18) w ( t, x I + t, ξ I − αt ) w − ( t, x I + t, ξ I − αt ) (cid:19) = (cid:18) w I, + ( x I , ξ I ) w I, − ( x I , ξ I ) (cid:19) . (5.13)To obtain (5.13), we have used the limit below:lim ξ → ( x +1 ( t ) , ξ +1 ( t )) = (cid:26) ( x I + t, ξ I − αt ) if 0 < t < t ∗ ( x I + 2 t ∗ − t, ξ I − αt ) if t > t ∗ which follows from the explicit formula of the classical trajectories solution to (5.6). Then, bypassing to the limit t → t ∗ in (5.13), it follows from assumption A (cid:18) w r + ( t ∗ , x ∗ ) w r − ( t ∗ , x ∗ ) (cid:19) = (cid:18) w I, + ( x I , ξ I ) w I, − ( x I , ξ I ) (cid:19) (5.14)where x ∗ = x I + t ∗ . The case t > t ∗ follows the same line with the difference that (5.12) has to bereplaced by: Z t e − R ts τ ( ξ + ( µ )) dµM F ( s ) ds = Z t ∗ e − R ts τ ( ξ + ( µ )) dµM F ( s ) ds + Z tt ∗ e − R ts τ ( ξ + ( µ )) dµM F ( s ) ds . (5.15)The arguments used to deduce (5.12) can be applied here to show that the first integral at the r.h.sof equation (5.15) tends to 0 when ξ →
0. Moreover, using assumption A C ( t − t ∗ ), which tends to0 when t → t ∗ (here C is a constant which does not depend on ξ ). Then, by taking successivelythe limit ξ → t → t ∗ in (5.8), we obtain: (cid:18) w l + ( t ∗ , x ∗ ) w l − ( t ∗ , x ∗ ) (cid:19) = e − βM (cid:18) w I, + ( x I , ξ I ) w I, − ( x I , ξ I ) (cid:19) . (5.16)Putting together (5.14) and (5.16), we get (cid:18) w l + ( t ∗ , x ∗ ) w l − ( t ∗ , x ∗ ) (cid:19) = (cid:18) − T TT − T (cid:19) (cid:18) w r + ( t ∗ , x ∗ ) w r − ( t ∗ , x ∗ ) (cid:19) (5.17)where the transition probability T is given by: T = 1 − e − β β is defined in (5.5). The system (5.17) has the form of the solution of the wellknown Landau-Zener problem (see [35]). Since x I and ξ I > t ∗ > x ∗ ∈ IR .We can now give the proof of (5.11). Differentiating (5.3) with respect to x gives the followingPDE for ∂ x w − : ∂ t ( ∂ x w − ) − ξ | ξ | ∂ x ( ∂ x w − ) − α∂ ξ ( ∂ x w − ) = τ ( ξ )( ∂ x w + − ( ∂ x w − )) . (5.19)Consider ϕ − s (˜ x I , ξ I ) = ( x − ( s ) , ξ − ( s )), the lower band classical trajectory solution to (5.6) with aninitial condition (˜ x I , ξ I ) such that ϕ − t (˜ x I , ξ I ) = ϕ + t ( x I , ξ I ). Then it holds ∂ x w − ( t, x − ( t ) , ξ − ( t )) = ∂ x w − ( t, x +1 ( t ) , ξ +1 ( t ))14nd it is enough to show that ∂ x w − ( t, x − ( t ) , ξ − ( t )) tends to 0 when ξ tends to 0. Now, if onedifferentiates the map s ∂ x w − ( s, x − ( s ) , ξ − ( s )), one gets using (5.19): dds (cid:0) ∂ x w − ( s, x − ( s ) , ξ − ( s )) (cid:1) = τ ( ξ − ( s )) (cid:0) ∂ x w + ( s, x − ( s ) , ξ − ( s )) − ∂ x w − ( s, x − ( s ) , ξ − ( s )) (cid:1) . By integrating the previous equation with the Duhamel formula, it follows: ∂ x w − ( t, x − ( t ) , ξ − ( t )) = Z t e − R ts τ ( ξ − ( µ )) dµ τ ( ξ − ( s )) ∂ x w + ( s, x − ( s ) , ξ − ( s )) ds (5.20)where we have used that w I, − = 0. Using assumption A
2, together with the equations (5.9) and(5.10), we can conclude that ∀ < t < t ∗ , the integral at the r.h.s. of (5.20) tends to 0 when ξ → Remark 5.2.
Although the system (5.17) has the correct form, the formula (5.18) for the transitionprobability is different from the Landau-Zener transition probability (3.3) obtained in section 3.2and which writes for our particular potential: T = e − πξ h | α | . (5.21)It will be verified numerically in section 6.3 that the band transmission corresponding to theeffective model (5.3) is given by (5.18) whereas the band transmission corresponding to the model(5.2) is given by the correct transition probability (5.21). This signifies that the approximationmade in [25] to obtain (5.3) from (5.2) is not correct. Remark 5.3.
If the coefficient τ is replaced by ξ , (5.3) becomes a hyperbolic relaxation system(see [14], [28]) and the solution to (5.3) satisfies w + = w − when ξ →
0. Since R w + + R w − = 1this leads to: Z w + = Z w − = 12 . (5.22)We will observe numerically in Figure 9 that (5.22) is true when the time is big enough. In this section, the results provided by the surface hopping algorithm are compared with thereference level populations given by (4.1) where the solution u h ( t ) is computed numerically usingan accurate method to solve the Dirac equation (2.10). In particular, it is verified numericallyin section 6.1 that the spectral method is more accurate than the finite difference method. Theimage of the operator Π ± ( hD ) is computed by using discrete Fourier transform (DFT) and Fouriermultiplication. A comparison of the surface hopping algorithm with the models (5.2) and (5.3) isgiven in section 6.3. We suppose that the initial data u hI is such that its Fourier transform ˆ u hI satisfies the relation:ˆ u hI ( ξ ) = ˆ f h ( ξ ) χ + ( hξ ) (6.1)where χ + ( ξ ) is defined in (2.9) and ˆ f h ( ξ ) is the Fourier transform of some function f h ∈ L ( IR ).Such an initial data satisfies the non-interference condition (4.7). We remark that, like for the levelpopulation (4.1), χ + ( hξ ) can be replaced by χ + ( ξ ) in (6.1).If f h is bounded in L ( IR ), using Lemma 2.3 in [10] again, we obtain the following approxi-mation for the scalar Wigner transform of u hI : w h (cid:2) u hI (cid:3) = W h (cid:2) f h (cid:3) + o (1) (6.2)15n D ′ (cid:16) IR x × ( IR ξ \ { ξ = 0 } ) (cid:17) when h →
0, where W h (cid:2) f h (cid:3) = w h ( f h , f h ). By plugging (6.2) in(4.9), we obtain the following asymptotics for the initial value of the diagonal terms of the Wignermatrix: w hI, + = W h (cid:2) f h (cid:3) + o (1) , w hI, − = o (1) (6.3)in D ′ (cid:16) IR x × ( IR ξ \ { ξ = 0 } ) (cid:17) when h →
0. In the present section, the initial upper level functionis an h -scaled Gaussian wave packet: f h ( x ) = ( πh ) − e − | x − xh | h + i ξ . ( x − xh h with center x h ∈ IR , momentum ξ ∈ IR and norm k f h k L ( IR ) = 1. Its h -scaled Fourier transformand its Wigner transform can be computed explicitly. Indeed, F h f h ( ξ ) = ( πh ) − e − | ξ − ξ | h − i xh .ξh where ∀ u ∈ L ( IR ) , F h u ( ξ ) = h − F u ( h − ξ ). Moreover, W h (cid:2) f h (cid:3) ( x, ξ ) = ( πh ) − e − | x − xh | h − | ξ − ξ | h . In this section, the parameter h is equal to its physical value given by (2.7). The effective mass ofthe electron is given by m = 0 . m e where m e is the mass of the electron, the Fermi velocity is taken as v F = 10 m.s − and, havingin mind the simulation of devices of size equal to hundreds of nanometers as in [25, 12], we willtake the reference length equal to L = 500 nm . Then, the numerical value of the parameter h corresponding to formula (2.7) is: h = 3 . × − which is small enough for the problem to be considered in the semiclassical regime. The simulationdomain Ω is equal to Ω = [ − √ h, √ h ] × [ − √ h, √ h ] . In the next section, we solve the Dirac equation (2.10) for different choice of the potential U byusing the time-splitting spectral method (TSSM) presented in Appendix A. The potential is equal to U = v x ≥ which corresponds to a Klein step. For this potential, wecompare the TSSM and the Finite difference time domain method (FDTD) in [12]. The centerand the momentum of the initial Gaussian wave packet f h are chosen as x h = ( − √ h, , ξ = 12 (1 , . The height of the Klein step and the simulation time are respectively v = 2 | ξ | and t f = 13 √ h . For the TSSM, the number of grid points are given by N = 126 , N = 66and for the FDTD by N = 850 , N = 426 . t = ∆ x √ , where ∆ x and ∆ x are the mesh sizes in the x and x directions respectively. In the presentsection and in the following section, the supscript h will be omitted for the solution u h ( t, x ) to(2.10) and for the initial condition u hI ( x ). Using the discrete Fourier transform (DFT) (A.5), resp.its inverse, to approximate the Fourier transform, resp. the inverse Fourier transform, we get thefollowing approximation of the projectors Π ± ( hD ) u ( t n , x j ):Π ± ( hD ) u nj = 1 N N X k ∈K Π ± ( hξ k ) d ( u n ) k e iξ k . ( x j − a ) , j ∈ J (6.4)where d ( u n ) k is the Fourier coefficient defined in (A.5) and the discretization is the same as inAppendix A. Then, using formula (4.1), the approximation P hdir, ± ( t n ) of the level populations P h ± ( t n ) of the Dirac equation is given by: P hdir, ± ( t n ) = k Π ± ( hD ) u n k (6.5)where Π ± ( hD ) u nj is defined by (6.4) and for u = ( u j ) j ∈J , u j ∈ C , we have k u k = X j ∈J | u j | ∆ x ∆ x . (6.6)Similarly, the initial condition u I ( x j ) defined by (6.1) is approximated by the formula:( u I ) j = 2 π ( b − a )( b − a ) X k ∈K ˆ f h ( ξ k ) χ + ( hξ k ) e iξ k .x j , j ∈ J . We remark that in the previous formula, we used the exact value of the Fourier transform ˆ f h ( ξ k )instead of the DFT [ ( f h ) k . The initial data for the TSSM with the parameters described above isrepresented in Figure 1.In Figure 2, we depict the evolution with respect to the time t n of the level populations P hdir, ± ( t n ) provided by the TSSM and by the FDTD. The curve with title Upper level, resp.Lower level, refers to the plus sign, resp. minus sign, in (6.5) and the curve with title Total refersto the total population k u n k . Initially, the charge is carried completely by the upper level, thennon-adiabatic transfer occurs at time t = 0 .
28 and the charge is almost all transferred to the lowerlevel. Moreover, we remark that the TSSM is more accurate than the FDTD. Indeed, the firstmethod conserves the total mass (A.8) whereas for the second method the total mass decreasesat the hopping time (it was shown in [12] that the quantity conserved by the FDTD is not themass but a related functional). To reduce this mass loss, the number of spatial points has to bechosen big enough which increases the CPU time of the method. In particular, the CPU timescorresponding to the simulations of Figure 2 are the following: 493 . s for the FDTD, 1 . s for the TSSM.We remark in Figure 3 that the transition occurs when the charge reaches the Klein step whichis a classically forbidden region for the Upper level. The band transition to the Lower level makesthe Klein step an allowed region for the particles, which can tunnel in the region x ≥ U = αx , α > U = αx , α >
0. The Dirac equation (2.10) is solved usingthe TSSM. The center and the momentum of the initial Gaussian wave packet f h are taken to be x h = ( − √ h, , ξ = (1 , . | u ( x ) | x x | u ( x ) | -0.45 -0.4 -0.35 -0.3 -0.25 -0.2 -0.15-0.2 0 0.2-8-6-4-2 0 2 4 6 8 R e ( u ( x )) x x R e ( u ( x )) Figure 1: Initial condition for the TSSM: | u I ( x ) | (top) and, in a smaller region, Re (( u I ) ( x ))(bottom). 18 P opu l a t i on Time Total, TSSMUpper level, TSSMLower level, TSSMTotal, FDTDUpper level, FDTDLower level, FDTD
Figure 2: Time evolution of the level populations with a Klein step potential: comparison of theTSSM and the FDTD. 19 =0-0.2 0 0.2 x t=0.2150 t=0.2867 t=0.3583 t=0.5733-0.4 0 0.4x -0.2 0 0.2 x -0.4 0 0.4x -0.4 0 0.4x -0.4 0 0.4x -0.4 0 0.4x Figure 3: For different times and with respect to the space variable: representation of the positiondensity of the projection of the wave function u solution to (2.10) where U = v x ≥ and theinitial condition is given by (6.1). The upper half corresponds to | Π + ( hD ) u nj | and the lower halfto | Π − ( hD ) u nj | where the projectors are computed using (6.4). The solution u is computed withthe TSSM. 20he level position density of the solution u ( t, x ) to (2.10) is supported around ( x ± ( t ) , ξ ± ( t )),solution to: ( ddt x ± ( t ) = ± ξ ± ( t ) | ξ ± ( t ) | , x ± (0) = x h ddt ξ ± ( t ) = − α (1 , , ξ ± (0) = ξ . The solution of the above problem can be computed explicitly. It is given by: x ± ( t ) = x ∗ ∓ ( | t − t ∗ | , , ξ ± ( t ) = ξ − α ( t,
0) (6.7)where t ∗ = ( ξ ) α is the time such that ξ + ( t ∗ ) = 0 and the non-adiabatic transfer occurs. The point x ∗ = x + ( t ∗ ) = x h + ( t ∗ ,
0) is the point where the hopping occurs. The coefficient α is chosen suchthat the potential at the hopping point is equal to U ( x ∗ ) = v where v = | ξ | . This leads to: α = v − ( ξ ) ( x h ) . The simulation stops at time: t f = 13 √ h and the number of space points are given by N = 160 , N = 80 . The time step size is equal to ∆ t = ∆ x √ . The initial data ( u I ) j , the Lower level and Upper level populations P hdir, ± ( t n ) and the Totalpopulation are computed as explained in section 6.1.1. The time evolution of the level populationsis represented in Figure 4. The numerical hopping is observed around the time t = 0 .
39 which isaccurate enough compared to the predicted value t ∗ = 0 . x ∗ = (9 . × − , x ∗ . As for the Kleinstep potential, the band transition to the lower level allows the particles to tunnel in the region x ≥ x ∗ . However, for the potential considered here, we can observe that the part remaining onthe upper level is reflected. This could have been predicted from equation (6.7). Indeed, for t ≥ t ∗ ,the classical flow corresponding to the upper level (plus sign) moves to the left with respect to x ∗ whereas the classical flow corresponding to the lower level (minus sign) moves to the right withrespect to x ∗ . In the present section, the potential U is equal to U = αx where α = 15. The initial condition isas described in section 6.1. The center and the momentum of the initial Gaussian wave packet f h are taken as x h = ( − √ h, , ξ = (1 , . (6.8)The diagonal terms w h ± ( t, x, ξ ) of the Wigner matrix defined by (4.3) are computed using thesurface hopping algorithm presented in section 4. Then, for different values of the parameter h and of the simulation time t f , the surface hopping level populations P hsh, ± ( t f ) are computed from w h ± ( t f , ., . ) by using (4.6). The results provided by the surface hopping algorithm are compared tothe level populations P hdir, ± ( t f ) of the Dirac equation computed as explained in section 6.1.1.For the level populations P hdir, ± ( t f ), the Dirac equation (2.10) is solved using the TSSM in thesimulation domain Ω = [ − √ h, √ h ] × [ − √ h, √ h ] . P opu l a t i on Time TotalUpper levelLower level
Figure 4: Time evolution of the level populations for the potential U = αx , α > t f as follows: ∆ t = t f / . (6.9)The number of space points are given by: for h = 10 − and h = 10 − N = 250 , N = 126 , for h = 10 − N = 300 , N = 150 , for h = 10 − N = 850 , N = 426 . For the level populations P hsh, ± ( t f ), the initial sampling is chosen as follows. We consider a uniform J × J discretization of the domain of the x variable:[ a , b ] × [ a , b ]where a = a = x h − √ h and b = b = x h + 5 √ h and a uniform K × K discretization of the ξ variable domain [ c , d ] × [ c , d ]where c = c = ξ − √ h and d = d = ξ + 5 √ h . The grid points x j , j = 1 , . . . , J and ξ k , k = 1 , . . . , K are ordered such that | f h ( x ) | ≥ . . . ≥ | f h ( x J ) | , |F h f h ( ξ ) | ≥ . . . ≥ |F h f h ( ξ K ) | . =0-0.2 0 0.2 x t=0.1807 t=0.3976 t=0.5422 t=0.7591-0.4 0 0.4x -0.2 0 0.2 x -0.4 0 0.4x -0.4 0 0.4x -0.4 0 0.4x -0.4 0 0.4x Figure 5: For different times and with respect to the space variable: representation of the positiondensity of the projection of the wave function u solution to (2.10) where U = αx and the initialcondition is given by (6.1). The upper half corresponds to | Π + ( hD ) u nj | and the lower half to | Π − ( hD ) u nj | where the projectors are computed using (6.4). The solution u is computed with theTSSM. 23hen, we determine the minimal numbers N x of points of the x variable and N ξ of the ξ variablesuch that N x X j =1 | f h ( x j ) | ∆ x ≥ − tol x , N ξ X k =1 |F h f h ( ξ k ) | ∆ ξ ≥ − tol ξ where ∆ x = b − a J = b − a J , ∆ ξ = d − c K = d − c K and tol x , tol ξ are well chosen tolerances. Thephase space points ( x j , ξ k ), j = 1 , . . . , N x , k = 1 , . . . , N ξ are ordered such that: W h (cid:2) f h (cid:3) ( x , ξ ) ≥ . . . ≥ W h (cid:2) f h (cid:3) ( x N x N ξ , ξ N x N ξ )and we determine the minimal integer N such that: N X k =1 W h (cid:2) f h (cid:3) ( x k , ξ k )∆ x ∆ ξ ≥ − tolwhere tol is a well chosen tolerance. This gives the first step of the algorithm of section 4. Indeed,the initial sampling is the set: { ( x k , ξ k , +) ∈ IR x × IR ξ × {− , + } ; k = 1 , ..., N } where, using equation (6.3), the associated weights w k ∈ IR can be approximated by: w k = W h (cid:2) f h (cid:3) ( x k , ξ k ) . In the present case, the hopping surface S defined by (3.2) is equal to: S = { ( x, ξ ) ∈ IR ; ξ = 0 } and the second order derivative of the function s
7→ | ξ ± ( s ) | , defined by the characteristics solutionto (2.12), is equal to: (cid:0) | ξ ± ( s ) | (cid:1) ′′ = 2 α > . Therefore, for the potential considered here, the points of extremal gap are all minimas.Moreover, the classical flow ( x ± ( t ) , ξ ± ( t )) solution to (2.12) is given by x ± ( t ) = x ± Z t ξ ± ( s ) | ξ ± ( s ) | ds, ξ ± ( t ) = ξ − αt (1 , . (6.10)In equation (6.10), the formula of the momentum is explicit and the hopping transport step canbe simplified. Indeed, for k = 1 , . . . , N , if 0 < ( ξ k ) α < t f , the trajectory ( x k ( t ) , ξ k ( t )), 0 < t < t f defined by (4.5) will pass through a point ( x ∗ k , ξ ∗ k ) ∈ S at the time t ∗ = ( ξ k ) α and non-adiabatictransfer occurs. For such a k , the weight is changed such that for t > t ∗ w k ( t ) = (1 − T ∗ ) w k ( t ∗ )and a new particle is created on the lower band with index l > N . For the new particle, theassociated weight is such that for t > t ∗ w l ( t ) = T ∗ w k ( t ∗ ) . In the above equations, the transition rate T ∗ is equal to T ∗ = T ( x ∗ k , ξ ∗ k )where T ( x, ξ ) is given by (5.21). The surface hopping level populations defined by (4.6) are thengiven by: P hsh, + ( t f ) = N X k =1 w k ( t f )∆ x ∆ ξ , P hsh, − ( t f ) = M X l = N +1 w l ( t f )∆ x ∆ ξ . P hdir, + ( t f ) P hsh, + ( t f ) P hdir, − ( t f ) P hsh, − ( t f ) CPU dir ( s ) CPU sh ( s )10 − . × − . × − . . . . − . × − . × − . . . . − . × − . × − . . . . − . × − . × − . . . . t f = 0 .
13 and for different values of the semiclassical parameter h : level popula-tions obtained by the Dirac solver and the surface hopping method.For all the tests, the size of the sampling grids are equal to: J = K = 16 . The sampling tolerances are taken as:tol = 10 − , tol x = tol ξ = 10 − × tol . For such parameters, the number of particles obtained numerically for the initial sampling is equalto N = 6981.For the simulation time t f = 0 .
13 and for different values of the semiclassical parameter h , thelevel populations obtained by the two methods are listed in Table 1. We remark that the referencelower level population P hdir, − ( t f ) increases when h increases. This is due to the fact that for largervalues of h , the transition process is slower and, at the time t f , the post-transition relaxationobserved in Figure 7, has no yet happened. We notice that, for h smaller or equal to 10 − , thesurface hopping level populations P hsh, ± ( t f ) are almost constant with respect to h . This can beexplained by the fact that, for the potential considered, the transition rate depends only on thevariable ξ / √ h which does not depend on h when considered on the sampling points. The columnCPU dir, resp. CPU sh, denotes the CPU time required by the Dirac solver, resp. the surfacehopping method, to complete the simulations. The time required to choose the initial sampling isnot taken into account in CPU sh. It appears that the surface hopping method is very interesting.Indeed, when h decreases, the numerical cost of the Dirac solver increases whereas the surfacehopping method provides very accurate results for an almost constant and much smaller CPUtime.We depict in Figure 6 the absolute error of the level population corresponding to the upperlevel: E h + ( t f ) := | P hdir, + ( t f ) − P hsh, + ( t f ) | (6.11)with respect to the semiclassical parameter h . We remark that the surface hopping level populationsconverge numerically when h → h = 10 − , the time evolution of the level populations provided by the Dirac solver andthe surface hopping algorithm are depicted in Figure 7. The curve with title Upper level dir,resp. Lower level dir, refers to P hdir, + ( t ), resp. P hdir, − ( t ). The same is true when dir is replacedby sh. The total population of the Dirac equation, curve titled Total dir, is obtained as explainedin section 6.1.1 and the surface hopping total population, curve titled Total sh, corresponds to P hsh, + ( t ) + P hsh, − ( t ). We observe that the numerically obtained surface hopping total populationis conserved. For a smaller CPU time, the time evolution of the level populations provided bythe surface hopping algorithm agrees well with the one provided by the Dirac solver. Indeed, forthe two methods, the charge is initially carried completely by the upper level, then non-adiabatictransfer occurs at time t = ( ξ ) α = = 0 . { ξ = 0 } starting from the momentum ξ of the initial Gaussian wave packet) and thegreat majority of the charge is transferred to the lower level. The CPU time required to completethe simulation is 14 . s for the Dirac solver and 5 . s for the surface hopping method. Forthe second method, the CPU time is bigger than the time indicated in Table 1. This is due tothe fact that the surface hopping curves in Figure 7 are obtained by repeating the surface hoppingalgorithm described in the present section for the sequence of times t f = 0 . n/ ≤ n ≤ Le v e l popu l a t i on ab s o l u t e e rr o r h Figure 6: At time t f = 0 .
13, logarithmic plot of the absolute error (6.11) of the level populationcorresponding to the upper level while varying the semiclassical parameter h = 10 − p , p = 1 , , , P opu l a t i on Time Total dirUpper level dirLower level dirTotal shUpper level shLower level sh
Figure 7: For h = 10 − and U = αx , α = 15: time evolution of the level populations provided bythe Dirac solver and the surface hopping algorithm.26 .3 Simulations using the models of [25] In this section we will compare numerically the models (5.2) and (5.3) with the two dimensionalversion of our surface hopping algorithm. The comparison will be performed with h = 10 − , t f = 0 .
13 and with the same potential as is section 6.2. The initial condition is the 2D version ofthe initial condition in section 6.2. In other words: w h + | t =0 = W h [ f h ] , w h − | t =0 = 0is replaced by w + (0 , x , ξ ) = ( πh ) − e − ( x − ( xh h − ( ξ − ( ξ h , w − (0 , x , ξ ) = 0 . (6.12)The functions w ± and w i depend on ξ . However, since ξ plays the role of a parameter, thisdependence is not written on the l.h.s. of the equations in (6.12). The center and the momentumof the initial Gaussian wave packet are given by (6.8). The x -independent solutions to the surface hopping algorithm presented in section 4 are obtainedby replacing the system (2.14) with its two dimensional version: ∂ t w ± ± ξ | ξ | ∂ x w ± − α∂ ξ w ± = 0 , ξ = 0 . (6.13)More precisely, equation (6.13) is used for the time evolution of w ± as long as the classical trajec-tories ϕ ± t ( x , ξ ) solution to (5.6) are away from the hopping surface { ξ = 0 } . Then, the hoppingtransport is described as follows.Consider an initial set of sampling points { ( x ,k , ξ ,k , +) ∈ IR x × IR ξ × {− , + } ; k = 1 , ..., N } with associated weights w k ∈ IR given by: w k = w + (0 , x ,k , ξ ,k ) . For t ≥ x ,k ( t ) , ξ ,k ( t )) = ϕ + t ( x ,k , ξ ,k ) , w k ( t ) = w k where, using (5.6), we have ξ ,k ( t ) = ξ ,k − αt . Then, for k = 1 , . . . , N , if 0 < ξ ,k α < t f , the classicaltrajectory is such that ξ ,k ( t ∗ ) = 0 at the time t ∗ = ξ ,k α and non-adiabatic transfer occurs. Forsuch a k , the weight is changed such that for t > t ∗ w k ( t ) = (1 − T ∗ ) w k ( t ∗ )and a new particle is created on the lower band with index l > N . For the new particle, theassociated weight is such that for t > t ∗ w l ( t ) = T ∗ w k ( t ∗ ) . In the above equations, the transition rate T ∗ is equal to: T ∗ = e − πξ h | α | . The surface hopping level populations at the final time t f are then given by: P hsh, + ( t f ) = N X k =1 w k ( t f )∆ x ∆ ξ , P hsh, − ( t f ) = M X l = N +1 w l ( t f )∆ x ∆ ξ (6.14)where ∆ x and ∆ ξ are the mesh sizes in the x and ξ directions respectively.27 .3.2 Comparison with the models in [25] For equations (5.2) and (5.3), the simulation domain is:( x , ξ ) ∈ h − √ h, √ h i × h ( ξ ) − αt f − √ h, ( ξ ) + 5 √ h i . They are solved with periodic boundary conditions on a uniform grid with 500 points in the x -direction and 500 points in the ξ -direction. The time step size is chosen such that the conditionof stability of the upwind method is satisfied [22]. More precisely, we take:∆ t = (cid:18) x + α ∆ ξ (cid:19) − . Equations (5.2) and (5.3) are solved using a time-splitting method where the free equation (withoutsource term) is solved using a dimensional splitting: the free problem is splitted in two one-dimensional problems and each one dimensional problem is solved using the one-dimensional secondorder upwind MUSCL scheme (see[22]). The source term is integrated in time using a RK2 methodfor equation (5.2) (in order to preserve the second order accuracy) and exactly for equation (5.3).For equation (5.2), we take: w i | t =0 = 0in addition to the initial condition (6.12). We remark that, since the MUSCL method is writtenfor real valued functions, the third equation in (5.2) has to be splitted in two equations, one forthe real part of w i and one for its imaginary part.The level populations provided by the asymptotic model are defined by: P ham, ± ( t ) = Z IR w ± ( t, x , ξ ) dx dξ where w ± is the solution to (5.2) and the level populations provided by the effective model aredefined by: P hem, ± ( t ) = Z IR w ± ( t, x , ξ ) dx dξ where w ± is the solution to (5.3). For the level populations P hsh, ± ( t f ) defined in (6.14), the initialsampling is chosen as follows. We consider a uniform J × K discretization of the domain of the( x , ξ ) variables: [ a, b ] × [ c, d ]where a = ( x h ) − √ h , b = ( x h ) + 5 √ h and c = ( ξ ) − √ h , d = ( ξ ) + 5 √ h . The grid points( x ,j , ξ ,k ), j = 1 , . . . , J , k = 1 , . . . , K are ordered such that: w + (0 , x , , ξ , ) ≥ . . . ≥ w + (0 , x ,J × K , ξ ,J × K )and we determine the minimal integer N such that: N X k =1 w + (0 , x ,k , ξ ,k )∆ x ∆ ξ ≥ − tolwhere ∆ x = b − aJ , ∆ ξ = d − cK and tol is a well chosen tolerance. Then, the initial sampling is theset: { ( x ,k , ξ ,k , +) ∈ IR x × IR ξ × {− , + } ; k = 1 , ..., N } . For all the tests in this section, the number of grid points for the two dimensional surface hoppingalgorithm is: J = 100 , K = 10028 P opu l a t i on Time Total amUpper level amLower level amTotal shUpper level shLower level sh
Figure 8: For h = 10 − , ξ = 10 − and U = αx , α = 15: time evolution of the level popula-tions provided by the asymptotic model (5.2) and the two dimensional surface hopping algorithmpresented in section 6.3.1.and the tolerance is tol = 10 − . We depict in Figure 8 the time evolution of the level populations provided by the asymptoticmodel (5.2) and the two dimensional surface hopping algorithm. The curve with title Upper levelam, resp. Lower level am, refers to P ham, + ( t ), resp. P ham, − ( t ). The total population, curve titledTotal am, corresponds to P ham, + ( t ) + P ham, − ( t ). The same is true when am is replaced by sh.The surface hopping curves are obtained by repeating the surface hopping algorithm described insection 6.3.1 for the sequence of times t f = 0 . n/ ≤ n ≤ t = 0 . T ∗ is replaced by(5.18)(5.5). The curve with title Upper level em, resp. Lower level em, refers to P hem, + ( t ), resp. P hem, − ( t ). The total population, curve titled Total em, corresponds to P hem, + ( t ) + P hem, − ( t ). Thesame is true when em is replaced by sh. As in Figure 8, hopping occurs at time t = 0 . P hem, − ( t f ), thepopulation provided by the effective model (5.3) on the lower level and at the final time, for 3129 P opu l a t i on Time Total emUpper level emLower level emTotal shUpper level shLower level sh
Figure 9: For h = 10 − , ξ = 10 − and U = αx , α = 15: time evolution of the level populationsprovided by the effective model (5.3) and the two dimensional surface hopping algorithm presentedin section 6.3.1 where the transition probability T ∗ is replaced by (5.18)(5.5).30 T β Trans effTrans th
Figure 10: For h = 10 − , ξ = 10 − , t f = 0 .
13 and U = αx , α = 15: numerical verification of theformula (5.18) for the transition probability corresponding to the model (5.3).different values of the constant β distributed uniformly on the interval [0 , β givenby (5.5) is made arbitrary in (5.18) by replacing the coefficient τ appearing in (5.3) by ˜ τ = βπ / τ .The curve with title Trans th is the representation of the coefficient T given by (5.18) for thesame values of β . We remark that the two curves are very close which validates the limit ξ → Appendix A Time-splitting spectral method (TSSM)
The solution u ( t, x ) to (2.6) is computed on the domain:Ω = [ a , b ] × [ a , b ]using a time-splitting spectral method as in [13]. For r = 1 ,
2, we choose the spatial mesh size∆ x r = b r − a r N r in the r direction for a given integer N r . Define a uniform grid x j = a + ( j ∆ x , j ∆ x ) , j ∈ J (A.1)where J = { j = ( j , j ) ∈ N | ≤ j < N , ≤ j < N } and a = ( a , a ). For a given time step size ∆ t , let u nj denote the numerical approximation of u ( t n , x j ) at the time t n = n ∆ t , n ≥
0. Then, u n +1 j is computed from u nj by decomposing theproblem (2.6) in the two sub-problems ih∂ t u = [ − ihσ ∂ x − ihσ ∂ x ] u (A.2)and ih∂ t u = U u . (A.3)31he free Dirac equation (A.2) is solved using a spectral method in space and exact time integration,whereas equation (A.3) can be integrated exactly on [ t n , t n +1 ]. To discretize (A.2), we introducethe trigonometric interpolant of u :˜ u ( t, x ) = 1 N N X k ∈K d u ( t ) k e iξ k . ( x − a ) (A.4)where ˜ u ( t, x j ) = u ( t, x j ). In equation (A.4), we have K = { k = ( k , k ) ∈ Z | − N ≤ k < N , − N ≤ k < N } and ξ k = 2 π (cid:18) k b − a , k b − a (cid:19) . For a given function f : IR → C , ˆ f k is the discrete Fourier transform (DFT) of f defined byˆ f k = X j ∈J f j e − iξ k . ( x j − a ) (A.5)where f j = f ( x j ). When f = ( f j ) j ∈J is a sequence with f j ∈ C , the DFT of f is the sequenceˆ f k defined by equation (A.5). Inserting (A.4) into (A.2), and using the orthogonality of the set (cid:8) e iξ k . ( x − a ) , k ∈ K (cid:9) with respect to the scalar product of L (Ω), one gets the following system ofODE ddt d u ( t ) k = − iB ( ξ k ) d u ( t ) k . For any t ∈ IR , the exact solution to the previous equation is given by: d u ( t ) k = M ( t − t , ξ k ) [ u ( t ) k where M ( δ, ξ ) = e − iδB ( ξ ) . Applying an inverse discrete Fourier transform, one obtains the following expression for the ap-proximation u ( t ) j of u ( t, x j ): u ( t ) j = 1 N N X k ∈K M ( t − t , ξ k ) [ u ( t ) k e iξ k . ( x j − a ) . We remark that using the eigenvectors (2.9), the matrix B ( ξ ) can be diagonalized as follows: B ( ξ ) = P ( ξ ) D ( ξ ) P ( ξ ) ∗ where D ( ξ ) = diag ( | ξ | , −| ξ | ) , P ( ξ ) = ( χ + ( ξ ) , χ − ( ξ ))and therefore M ( δ, ξ ) = P ( ξ ) e − iδD ( ξ ) P ( ξ ) ∗ = (cid:18) cos( δ | ξ | ) − i sin( δ | ξ | )( ξ − iξ ) / | ξ |− i sin( δ | ξ | )( ξ + iξ ) / | ξ | cos( δ | ξ | ) (cid:19) . (A.6)The Strang splitting is the second order method constructed as follows: solve the first subproblem(A.2) over only a half time step of length ∆ t . Then, we use the result as data for a full time stepon the second subproblem (A.3) and finally take another half time step on (A.2). This leads tothe following method: u ∗ j = 1 N N X k ∈K M (∆ t/ , ξ k ) d ( u n ) k e iξ k . ( x j − a ) u ∗∗ j = e − ih U j ∆ t u ∗ j (A.7) u n +1 j = 1 N N X k ∈K M (∆ t/ , ξ k ) [ ( u ∗∗ ) k e iξ k . ( x j − a ) . emark A.1. It follows directly from the unitarity of the DFT and of the matrix M ( δ, ξ ) givenby (A.6) that the TSSM (A.7) conserves the discrete total charge, i.e.: k u n k = k u k , ∀ n ≥ k . k is defined in (6.6). Moreover, when the potential is equal to a constant U = U , the Dirac equation (2.6) admits the following plane wave solutions: u ( t, x ) = χ ± ( k ) e ih ( k.x − ( U ±| k | ) t ) , which are integrated exactly by the TSSM (A.7) if k ∈ IR satisfies kh = ξ k ′ for some k ′ ∈ K . Inthe case of the Schr¨odinger equation, an analysis of the stability of the TSSM was performed in[6] for initial conditions which are close to plane waves. Remark A.2.
In the applications, the solution of the Dirac equation (2.6) moves away fromthe simulation box. Therefore, in addition to the periodic boundary conditions provided by thespectral method, we use absorbing boundary layers at the edge of the simulation box, see e.g. [29].
Acknowledgements
A.F. acknowledges Clotilde Fermanian Kammerer for discussion about sec-tion 3.3.
References [1] C. W. Beenakker, Colloquium: Andreev reflection and Klein tunneling in graphene, Rev. Mod.Phys. 80 (2008) 1337.[2] D. Brinkman, C. Heitzinger, P. A. Markowich, A convergent 2D finite-difference scheme forthe Dirac-Poisson system with magnetic potential and the simulation of graphene, J. Comput.Phys. 257, A (2014) 318-332.[3] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, A. K. Geim, The electronicproperties of graphene, Rev. Mod. Phys. 81, 1 (2009) 109-162.[4] L. Chai, S. Jin, Q. Li, O. Morandi, A Multiband Semiclassical Model for Surface HoppingQuantum Dynamics, SIAM Multiscale Model. Simul. 13, 1 (2015) 205-230.[5] K. Drukker, Basics of surface hopping in mixed quantum/classical simulations, J. Comp. Phys.153 (1999) 225-272.[6] E. Faou, L. Gauckler and C. Lubich, Plane Wave Stability of the split-step Fourier methodfor the nonlinear Schr¨odinger equation, Forum Math. Sigma 2, e5 (2014) 45 pp.[7] C. Fermanian Kammerer, P. G´erard, A Landau-Zener formula for non-degenerated involutivecodimension three crossings, Ann. Henri Poincar´e 4 (2003) 513-552.[8] C. Fermanian Kammerer, P. G´erard, Mesures semi-classiques et croisements de mode, Bull.S.M.F 130, 1 (2002) 123-168.[9] G. Fiori, G. Iannaccone, Simulation of Graphene Nanoribbon Field-Effect Transistors, IEEEElectron. Device Lett. 28, 8 (2007) 760-762.[10] P. G´erard, P. A. Markowich, N. J. Mauser, F. Poupaud, Homogenization Limits and WignerTransforms, Comm. Pure Appl. Math. 50, 4 (1997) 323-379.[11] G. A. Hagedorn, Proof of the Landau-Zener formula in an adiabatic limit with small eigenvaluegaps, Commun. Math. Phys. 136 (1991) 433-449.3312] R. Hammer, W. P¨otz, A. Arnold, Single-cone real-space finite difference scheme for the time-dependent Dirac equation, J. Comput. Phys. 265 (2014) 50-70.[13] Z. Huang, S. Jin, P. A. Markowich, C. Sparber, C. Zheng, A time-splitting spectral schemefor the Maxwell-Dirac system, J. Comput. Phys. 208 (2005) 761-789.[14] S. Jin, Runge-Kutta methods for hyperbolic conservation laws with stiff relaxation terms, J.Comput. Phys. 122 (1995) 51-67.[15] S. Jin, P. A. Markowich, C. Sparber, Mathematical and computational methods for semiclas-sical Schr¨odinger equations, Acta Numerica 20 (2011) 121-209.[16] S. Jin, K. Novak, A coherent semiclassical transport model for pure-state quantum scattering,Comm. Math. Sci. 8 (2010) 253-275.[17] S. Jin, P. Qi, A Hybrid Schr¨odinger/Gaussian Beam Solver for Quantum Barriers and SurfaceHopping, Kinetic and Related Models 4 (2011) 1097-1120.[18] S. Jin, P. Qi, Z. Zhang, An Eulerian surface hopping method for the Schr¨odinger equationwith conical crossings, Multiscale Model. Simul. 9, 1 (2011) 258-281.[19] C. Lasser, T. Swart, S. Teufel, A rigorous surface hopping algorithm for conical crossings,Commun. Math. Sc. 5, 4 (2007) 789-814.[20] C. Lasser, S. Teufel, Propagation through conical crossings: an asymptotic semigroup, Comm.Pure Appl. Math. 58, 9 (2005) 1188-1230.[21] M. C. Lemme, T. J. Echtermeyer, M. Baus, H. Kurz, A Graphene Field-Effect Device, IEEEElectron. Device Lett. 28, 4 (2007) 282-284.[22] R. J. LeVeque, Finite Volume Methods for Hyperbolic Problems, Cambridge University Press,Cambridge, 2002.[23] O. Morandi, Multiband Wigner-function formalism applied to the Zener band transition in asemiconductor, Physical Review B 80, 2 (2009) 024301-024312.[24] O. Morandi, F. Sch¨urrer, Wigner model for Klein tunneling in graphene, Communications inApplied and Industrial Mathematics 2, 1 (2011).[25] O. Morandi, F. Sch¨urrer, Wigner model for quantum transport in graphene, J. Phys. A: Math.Theor. 44 (2011) 265301.[26] D. S. Novikov, Elastic scattering theory and transport in graphene, Phys. Rev. B 76 (2007)245435.[27] K. S. Novoselov, A. K. Geim, S. V. Morozov, D. Jiang, Y. Zhang, S. V. Dubonos, I. V.Gregorieva, A. A. Firsov, Electric Field Effect in Atomically Thin Carbon Films, Science 306,5696 (2004) 666-669.[28] L. Pareschi, Central Differencing Based Numerical Schemes for Hyperbolic Conservation Lawswith Relaxation Terms, SIAM J. Numer. Anal. 39, 4 (2001) 1395-1417.[29] U. W. Rathe, C. H. Keitel, M. Protopapas, P. L. Knight, Intense laser-atom dynamics withthe two-dimensional Dirac equation, J. Phys. B 30, 15 (1997) L531-L539.[30] D. Sholl and J. Tully, A generalized surface hopping method, J. Chem. Phys. 109 (1998)7702-7710.[31] C. Sparber, P. A. Markowich, Semiclassical asymptotics for the Maxwell-Dirac system, J.Math. Phys. 44, 10 (2003) 4555-4572. 3432] B. Thaller, The Dirac Equation, Springer, New York, 1992.[33] J. Tully, Molecular dynamics with electronic transitions, J. Chem. Phys. 93 (1990) 1061-1071.[34] J. Tully, R. Preston, Trajectory Surface Hopping Approach to Nonadiabatic Molecular Colli-sions: The Reaction of H + with D2