Comparison of deterministic and stochastic methods for time-dependent Wigner simulations

Journal of Computational Physics 300 (2015) 167–185

Contents lists available at ScienceDirect

Journal of Computational Physics www.elsevier.com/locate/jcp

Comparison of deterministic and stochastic methods for time-dependent Wigner simulations Sihong Shao a , Jean Michel Sellier b,∗ a b

LMAM and School of Mathematical Sciences, Peking University, Beijing 100871, China IICT, Bulgarian Academy of Sciences, Acad. G. Bonchev str. 25A, 1113 Sofia, Bulgaria

a r t i c l e

i n f o

Article history: Received 20 March 2015 Received in revised form 2 June 2015 Accepted 1 August 2015 Available online 6 August 2015 Keywords: Quantum mechanics Wigner formalism Monte Carlo methods Spectral element methods Signed particles Time-dependent simulations

a b s t r a c t Recently a Monte Carlo method based on signed particles for time-dependent simulations of the Wigner equation has been proposed. While it has been thoroughly validated against physical benchmarks, no technical study about its numerical accuracy has been performed. To this end, this paper presents the first step towards the construction of firm mathematical foundations for the signed particle Wigner Monte Carlo method. An initial investigation is performed by means of comparisons with a cell average spectral element method, which is a highly accurate deterministic method and utilized to provide reference solutions. Several different numerical tests involving the time-dependent evolution of a quantum wave-packet are performed and discussed in deep details. In particular, this allows us to depict a set of crucial criteria for the signed particle Wigner Monte Carlo method to achieve a satisfactory accuracy. © 2015 Elsevier Inc. All rights reserved.

1. Introduction The Wigner formulation of quantum mechanics [1] represents an interesting and useful tool for both physicists and mathematicians as it is intuitive and describes quantum systems in terms of a (quasi) distribution function defined over a phase space. Its main advantage resides in the fact that both observables and states are represented by ordinary c-numbers in a non-commutative phase space and that this representation strongly mimics classical statistical mechanics [2,3]. As such, it defines physical systems in a language which is very close to the language of experimentalists and has found wide applications in many physical fields [4–6]. As a matter of fact, a whole branch of experimental physics exists, known as quantum tomography, which purpose is the reconstruction of the Wigner function from measurements [7,8]. Despite the great advantage of being intuitive, the resolution of the time-dependent Wigner equation remains an incredible challenge, from both analytical and numerical points of view, as it is a partial integro-differential equation where the unknown distribution function (i.e. the Wigner function) is defined over a 2 × D × N-dimensional phase space. Here D (= 1, 2, 3) is the dimension of space and N (= 1, 2, 3, · · · ) the number of bodies involved. Nowadays, numerical techniques exist to simulate the time-dependent Wigner equation and can be roughly categorized into deterministic and stochastic methods. Among the deterministic methods for time-dependent Wigner simulations, the first one is a first-order upwind scheme finite difference method (FDM) proposed by Frensley in simulating the quantum transport in a resonant tunneling diode [9,10] and after that several second-order FDMs have been used [11]. A detailed summary about FDMs for the

*

Corresponding author. E-mail addresses: [email protected] (S. Shao), [email protected], [email protected] (J.M. Sellier).

http://dx.doi.org/10.1016/j.jcp.2015.08.002 0021-9991/© 2015 Elsevier Inc. All rights reserved.

168

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

time-dependent Wigner equation can be found in Ref. [12]. However, it was reported that general FDMs are not very accurate for transient Wigner simulations [12] and what is worse for uniform FDMs, a constraint due to the requirement of electron conservation should be applied on the size of the computational domain in the wavenumber space [13]. To address these issues, several higher order methods have been introduced for transient Wigner simulations, e.g. a cell average spectral element method (SEM) [14], a hyperbolic moment method [15], a WENO-solver [16], etc. Besides of its high accuracy, it was also shown that the cell average SEM is able to maintain analytically the electron conservation for numerical solutions even with non-uniform adaptive meshes in the phase space [14]. Like all mesh-based numerical methods, the extension of these deterministic methods into multi-dimensional many-body quantum problems is not likely to be an easy task. Usually, researchers will use higher order deterministic methods to generate highly accurate benchmark solutions in the one-dimensional space, and recourse to stochastic particle (mesh-free) methods for multi-dimensional simulations. Recently a Monte Carlo method (MCM) based on the concept of signed particles [17] has been introduced [18,19]. This novel MCM is capable of simulating the Wigner equation for quantum single-body problems and has been further generalized to the important case of quantum many-body problems [20–22]. It should be noted that there exist other MCMs for time-dependent Wigner equation, e.g. [23,30,31]. In order to validate the signed particle Wigner MCM for time-dependent, multi-dimensional, single- and many-body problems, the usual strategy has been consisting of comparisons with well known physical systems and check if the method recovers the expected behavior or, in other words, if it mimics the physics involved in the process. For instance, when the single-body Wigner MCM was presented it was compared to the Schrödinger solution [19], and in the case of many-body problems it was shown to correctly rotate the entanglement of two particles in the phase space, and to correctly reconstruct a Fermi (or exchange-correlation) hole in the case of identical Fermions. From this perspective, one may claim that no meticulous study on the mathematical accuracy of the method has been done so far. This paper is the first step to correct this situation. This work performs a thorough comparison between the conservative cell average SEM [14] and the signed particle Wigner MCM [19] in the one-dimensional space for a single-body problem. The purpose of such comparison is twofold. On the one hand, we want to demonstrate the numerical accuracy of both methods. On the other hand, this allows us to depict criteria for the Wigner MCM in order to choose simulation parameters such as the time step and the frequency of annihilation processes. As a matter of fact, at the present time the choice of such parameters represents a weak point as it is a very complex procedure (practically by trials and errors) and can affect the accuracy of the solution (e.g. appearance of negative probabilities in the configuration space). In this paper, we report a set of criteria which gives clear directions on how to choose such parameters in order to achieve an acceptable accuracy. The rest of the paper is organized as follows. In the next section we describe the Wigner formulation of quantum mechanics. Then we comment on the discretization of the Wigner quasi-distribution function in a truncated phase space. A short description of the conservative cell average SEM and the signed particle Wigner MCM follows. Finally, a series of numerical experiments and comparisons is described and thoroughly analyzed, offering a set of criteria for the Wigner MCM in choosing the time step, the annihilation frequency and so on. 2. The Wigner function formalism In this section, we sketch the Wigner formulation of quantum mechanics. For the sake of clarity, we restrict the exposition to a single particle in a one-dimensional phase space. Namely, the Wigner function f (x, k, t ) for a one-dimensional quantum system resides in the phase space (x, k) ∈ R2 for the position x and the wavenumber k. It is defined by using the Weyl–Wigner transform of the density operator ρ (x, x , t ) as follows

+∞ f (x, k, t ) =

exp(−iky )ρ (x + y /2, x − y /2, t ) d y .

(1)

−∞

Starting from the Schrödinger equation, we can easily show that the Wigner function f (x, k, t ) satisfies the following dynamic equation (i.e. the time-dependent Wigner equation) [23]

∂ h¯ k ∂ f (x, k, t ) + f (x, k, t ) = Q V [ f ](x, k, t ), ∂t m ∂x

(2)

where m is the mass of the electron, h¯ is the reduced Planck constant, and the right hand side is the so-called Vlasov– Wigner pseudo-differential operator

Q V [ f ](x, k, t ) =

i h¯

+∞

exp(−iky )( V (x − y /2) − V (x + y /2)) f (x, y , t )d y .

−∞

(3)

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Here i is the imaginary unit, V (x) is the potential, and  f (x, y , t ) is just another notation for Eq. (1) through an inverse Fourier transform for the k-variable,

169

ρ (x + y /2, x − y /2, t ) from

+∞

1  f (x, y , t ) =

2π

exp(iky ) f (x, k, t )dk.

(4)

−∞

It has been shown that the L 2 -norm of the Wigner function f (x, k, t ) is time invariant for L 2 initial data and bounded potential, and so does  f (x, y , t ) due to Parseval’s identity [24]. Usually, we employ an equivalent formula for Q V as follows

+∞

V w (x, k − k) f (x, k , t )dk

Q V [ f ](x, k, t ) =

(5)

−∞

+∞

V w (x, k ) f (x, k + k , t )dk ,

=

(6)

−∞

where the non-local Wigner potential V w (x, k) is calculated from the physical potential V (x) by a Fourier transform as

V w (x, k) =

=

=

+∞

i 2π h¯

(7)

−∞

+∞

i 2π h¯

exp(−iky )( V (x + y /2) − V (x − y /2))d y

(8)

−∞

+∞

i

π h¯

exp(iky )( V (x − y /2) − V (x + y /2))d y

exp(−2iky )( V (x + y ) − V (x − y ))d y .

(9)

−∞

The Wigner function can be utilized to calculate the electron density n(x, t ) by

+∞ n(x, t ) =

f (x, k, t )dk,

(10)

−∞

and the current density j(x, t ) by

h¯

j(x, t ) =

m

+∞ kf (x, k, t )dk.

(11)

−∞

In addition, we denote the integration of the potential term with respect to k as

p(x, t ) =

+∞ Q V [ f ](x, k, t )dk.

(12)

−∞

Accordingly, integrating the Wigner equation (2) with respect to k from −∞ to +∞ leads to the so-called continuity equation,

∂ ∂ n(x, t ) + j(x, t ) = 0, ∂t ∂x

(13)

where we have used p(x, t ) vanishes (i.e. p(x, t ) ≡ 0) due to the anti-symmetry of V w (x, k) in k. This continuity equation corresponds to conservation of the zeroth moment, in other words, the mass conservation.

170

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

3. Cell average spectral element method An adaptive conservative cell average SEM for time-dependent Wigner simulations was proposed in Ref. [14]. For the sake of completeness, in this Section, we will briefly describe the main idea of the cell average SEM. More details about the method and numerical performance, especially on the adaptive part, are referred to Ref. [14]. The design of deterministic solvers for the time-dependent Wigner equation (2) starts from a necessary truncation of the infinite phase space. For the x-space, we just take the computational domain X = [xmin , xmax ] where the quantum system under consideration is placed. For the k-space, we adopt here a simple nullification of the distribution outside a sufficiently large k-domain K = [kmin , kmax ], and then the k-integration range in Eqs. (10)–(12) will be changed into K. That is, the size |K| = kmax − kmin needed depends on the decay of the Wigner distribution for large wavenumber k and affects the accuracy of the numerical solutions. Moreover, we need actually to consider  f (x, y , t ) for | y | ≤ 2Y for some large Y (i.e. truncating the y-space) because of the decay of  f (x, y , t ) in Eq. (4) as | y | → ∞. After implementing all those truncations and adopting the trapezoid rule for the y-integration in the Wigner potential, deterministic solvers try to integrate instead the following truncated (in both x and k) and semi-discrete (in y) Wigner equation

∂ h¯ k ∂ f (x, k, t ) + f (x, k, t ) ∂t m ∂x k max h Vw (x, k − k) f (x, k , t )dk , =

(x, k) ∈ X × K,

(14)

kmin

where h Vw (x, k)

=

Ny 2 y 

π h¯

N yy = Y ,

 





sin(2ky μ ) V x + y μ − V x − y μ



(15)

,

μ=1

y μ = μ y ,

μ = 1, · · · , N y .

(16)

Before presenting the SEM, an important issue on the mass conservation must be emphasized. It has been pointed out that conservation of mass is the most important property that a discretization of the Wigner equation should fulfill and only the numerical solution ensuring the mass conservation captures the correct physics and is reliable [13]. In order to keep the mass conservation (13), it is required that k max

k max

dk V w (x, k − k) f (x, k , t ) = 0.

dk kmin

(17)

kmin

To this end, Frensley suggested in [9,10] a sufficient condition

 y |K | = π ,

(18)

which guarantees that verified by noting that

h Vw (x, k)

k max

is not only odd but also periodic in k with a period |K|. This condition can be readily

sin[2μ y (k − k)]dk =

sin(μ y |K|)

μ y

sin(μ y (2k − kmin − kmax )).

(19)

kmin

We also see here that the condition (18) can be generalized into  y |K| = zπ with z being any positive integer and setting z = 1 leads to (18). Actually, the condition (18) results in a constraint on the size of the domain K. While using uniform FDMs to solve the Wigner equation, Eq. (18) says that the length of the computational domain K in the k-space should be doubled after halving the y-mesh size. If we use the same spacing in the x-space and the y-space, then we need four times as many grid points in the k-space and twice as many grid points in the x-space to double the resolution of f (x, k, t ). That is, for FDMs, it will result in huge cost for accurate numerical solutions in high-dimensional cases. But this is not the case for our cell average SEM as described below, which is able to keep discrete mass conservation analytically even with a non-uniform mesh. Let

K=

R  r =1

Kr =

N R  

I j,

Kr = [dr −1 , dr ],

d0 = kmin , d R = kmax ,

r =1 j =1

I j = [k j +1/2 , k j −1/2 ],

k j ∓1/2 = dr −1 +

d r − d r −1 2

 1 + cos

j−

1 2

∓

1 2



π N

 .

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

171

Here Kr are the non-overlapping elements in the k-space and I j are the cells in each element. Then we define local quantities corresponding to n(x, t ), j(x, t ) and p(x, t ) in each computational cell as follows



n j (x, t ) :=

f (x, k, t )dk, Ij

j j (x, t ) :=

h¯

(20)



m

kf (x, k, t )dk,

(21)

h Vw (x, k − k) f (x, k , t )dk dk.

(22)

Ij

  p j (x, t ) :=

Ij K

Accordingly, a local continuity equation for the cell I j is defined as

∂ ∂ n j (x, t ) + j j (x, t ) − p j (x, t ) = 0. ∂t ∂x

(23)

From Eqs. (20)–(22), we can see that both cell averages and point values of the Wigner function f (x, k, t ), expressed in terms of the Chebyshev polynomials of the k-variable over each cell, are involved in Eq. (23). To achieve both high accuracy and efficiency, we will take advantage of the fact [25,26] that, for any function u (k) in the space of the Chebyshev polynomials, there is an exact relation between its cell averages and point values. This relation allows us to determine the expansion coefficients of the numerical solution with a Chebyshev collocation method for local continuity equations, in an efficient manner with the help of fast sine and cosine transforms. That is, in order to determine the unknown expansion coefficients, we integrate the local continuity equation (23) where there are only spectral errors associated with the Chebyshev polynomial expansion because all the integrals in Eqs. (20)–(22) are calculated exactly. More importantly, it can be readily checked that the electron conservation for the numerical solution is maintained analytically even with non-uniform adaptive meshes in the phase space by using the additive property of the integral. After obtaining the above conservative cell average SEM in the k-space, we are able to solve the local continuity equation (23) to get the expansion coefficients by using a traditional (standard) collocation SEM with Gauss–Lobatto points in the x-space for easy implementation of the inflow boundary conditions [9] and fast cosine transforms [27]. For the time discretization, we employ explicit multi-step Runge–Kutta methods, e.g. the fourth-order Runge–Kutta scheme [28]. As demonstrated in Ref. [14], this cell average SEM is a high-order, conservative, stable and efficient numerical methods for transient Wigner simulations. 4. Signed particle Monte Carlo method In this section, we describe the Wigner MCM based on signed particles. For the sake of clarity, we restrict the exposition to a single particle in a one-dimensional configuration space. The interested reader can find further mathematical details in [19,21] for the single- and many-body Wigner equations respectively. For practical details on the implementation, the reader is invited to download the source code as well [29]. The signed particle Wigner MCM starts with the Lagrangian description of the time-dependent single-body Wigner equation (2),

d f (x(t ), k(t ), t ) dt

=

∂f ∂f +v = Q V [ f ], ∂t ∂x

x ∈ X,

(24)

with the velocity v = h¯ k/m, where attention is focused on a particle whose path x(t ) is obtained by integrating x˙ = v with x = x(0) at t = 0. By discretizing the wavenumber space k in multiples of the quantity k, one can express the semi-discrete version of the Wigner equation (24) as a Fredholm integral equation of the second kind

f (x, M , t ) − e−

t 0

∞ γ (x(ξ ))dξ f (x(0), M , 0) =

f (x , M  , t  )(x , M , M  )e−

t

t

0

dt 

+∞ 

+∞

dx ×

M  =−∞−∞

γ (x(ξ ))dξ θ(t − t  )δ(x − x(t  ))θ (x ), X

(25)

where

γ (x) =

+∞ 

+ VW (x, M ),

(26)

M =−∞

+ + (x, M , M  ) = V W (x, M − M  ) − V W (x, −( M − M  )) + γ (x)δ M , M  .

(27)

172

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Here the (shorter) notation f (x, M , t ) defines the quantity f (x, M k, t ), θ(t ) is the Heaviside step function, θX (x) is the + indicator function of the domain X , δ(x) is the Dirac delta function, δ M , M  is the Kronecker delta, and V W (x, M ) refers to the positive part of the product kV W (x, M k). The solution of the Fredholm integral equation of the second kind (25) can always be written as a (infinite) Liouville– Neumann series. Consequently, the average value A  of a macroscopic variable A = A (x, k) can also be written into a series. In particular, for Eq. (25), by writing down the first three terms of the series, one clearly sees how the expansion of A  branches, and the total value is given by the sum of all branches:

τ +∞ 2 +∞  

A s (τ ) = dt i dxi s =0

e−

e

−

ti 0

γ (xi (ξ ))dξ

⇑xi , M i ,0 t1 ti

−

τ A (x1 , M i ) δ(t i − τ ) + ⇑x1 =xi (t i )

γ (x1 (ξ ))dξ

t2 t1

τ ⇑x2 =x1 (t 1 )

γ (x2 (ξ ))dξ

A (x3 , M 2 ) δ(t 2 − τ )

dt 2 t1

,

⇑x3 =x2 (t 2 )

θX (x1 )(x1 , M 1 , M i ) ×

M 1 =−∞

A (x2 , M 1 ) δ(t 1 − τ ) +



⇑x2 , M 2 ,t 1

+∞ 

dt 1 ti

⇑x1 , M 1 ,t i

e

−∞

0

f ( x i , M i , 0) ×

M i =−∞

+∞ 

θX (x2 )(x2 , M 2 , M 1 )×

M 2 =−∞

(28)

where the coordinates of a new trajectory are denoted by the symbol ⇑. Therefore, a physical interpretation in terms of three appearing particles can be given, based on (x, M , M  ) in Eq. (27), along with the exponentials appearing in Eq. (28). A initial parent particle (with sign) evolves on a free-flight trajectory and, according to a generation rate given by the function γ (x) (depending on the actual position of the particle, and denoting the presence of a given potential), two new signed particles, one positive and one negative, are generated with one momentum state L = M − M  generated with a probability equal to + VW (x, L )

γ (x)

,

and, with the same probability, another momentum L  with L  = M  − M [19]. This interpretation can essentially be seen as a Monte Carlo method for the time-dependent simulation of the Wigner equation, and which involves only the evolution of Newtonian signed particles. One should note that the grow of (virtual) signed particles is exponential therefore showing the necessity for an annihilation process. Such process can be implemented by exploiting the fact that two particles with opposite signs and in the same phase-space cell have a total contribution to the Wigner function f = f (x, k, t ) equal to zero (as they cancel out due to their opposite signs). Finally, for the sake of clarity and duplication purposes, a simplified set of rules is given in Algorithm 1, which offers a relatively uncomplicated overview of the algorithm. We restrict ourselves to the case of a time-independent external potential (the extension to a time-dependent one being trivial). 5. Numerical experiments We now report a detailed description of the numerical experiments performed in the present study where length, time and energy are measured in nanometer (nm), femtosecond (fs) and electron volt (eV), respectively. The aim of this section is twofold. On one side we want to put the readers in the position to duplicate these simulations (and the signed particle Wigner MCM code can be actually downloaded [29]). On the other side, the details of the experiments are important as they give important hints on the physics involved. Comments and discussion of the results are given in the next section. We employ the relative L 2 error and the relative L ∞ error to study the accuracy. Let f ref (x, k, t ) denote the reference solution which could be the exact solution or the numerical solution on a relatively fine mesh, and f num (x, k, t ) the numerical solution. Then, the relative errors are written as



ε2 (t ) =

X ×K ( f

ε∞ (t ) =

2

X ×K ( f (x, k, t ))

dxdk

ref (x, k, t ))2 dxdk

max(x,k)∈X ×K  f (x, k, t ) max(x,k)∈X ×K f ref (x, k, t )

,

1

2

,

(29)

(30)

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

173

Algorithm 1: A simple diagram of the signed particle Wigner Monte Carlo method. The symbol N A refers to the frequency of annihilation processes happening during the simulation, namely, every N A time steps an annihilation is performed. Input : The initial Wigner function f (x, k, t = 0), the external potential V (x), the initial number of particles N P , the time step t, the final time T , the frequency of annihilation N A . Output: The Wigner function f (x, k, t = T ) and the number of particles N P . 1 Initialize signed particles according to f (x, k, t = 0) and N P ; 2 Calculate the total number of time steps N = T /t; h 3 Compute the Wigner potential V w in Eq. (15); 4 Compute the function γ (x) in Eq. (26); 5 for n = 1, 2, . . . , N do 6 foreach particle i do 7 Set the local time t local = 0; 8 while t local < t do 9 Compute the time interval δt in Eq. (33) at the end of which a creation event happens; 10 Evolve particle i for the free time interval δt; 11 Generate signed particles according to γ (x); 12 Update the total number of particles N P ; 13 Set t local ← t local + δt;

end while

14 15 16 17 18

end foreach if mod (n, N A ) = 0 then Perform the annihilation; Update the total number of particles N P ;

19

end if

20 end for

where  f (x, k, t ) = | f num (x, k, t ) − f ref (x, k, t )|, and the integrals above are evaluated using a simple rectangular rule over a uniform mesh,

xi = xmin + (i − 1/2)h x ,

k j = kmin + ( j − 1/2)hk ,

with i = 1, 2, · · · , 200, j = 1, 2, · · · , 400, where

hx =

xmax − xmin 200

hk =

,

kmax − kmin 400

.

In order to obtain a more complete view of the accuracy, we also measure the relative errors for physical quantities, e.g. the electron density n(x, t ), in a similar way. For the cell average SEM, we set the computational domain in the x-space X = [xmin , xmax ] = [0 nm, 60 nm], the spacing  y = 0.3 nm and thus the computational domain in the k-space K = [kmin , kmax ] = [−π /(2 y ), π /(2 y )] ≈ [−5.2360 nm−1 , 5.2360 nm−1 ] due to the constraint (18) from the mass conservation, the number of elements Q × R = 10 × 20 (meaning the domain X × K is divided into 10 × 20 non-overlapping elements), and adopt a uniform mesh (choosing the same ( M q,r , N q,r ) for all elements) with N q,r = M q,r ≡ N e for r = 1, 2, · · · , R, q = 1, 2, · · · , Q . For the signed particle MCM, we take the same computational domain X × K and set k = hk for convenience. The motion of a Gaussian wavepacket [12,14] is simulated to investigate the performance of the deterministic cell average SEM and the stochastic signed particle MCM for time-dependent Wigner equation (2). The Gaussian wavepacket in the free space is



 (x − x0 − v 0t )2 f (x, k, t ) = exp − π 2a2 (1 + β 2 t 2 )     β t (x − x0 − v 0 t ) 2 2 2 2 × exp −2a (1 + β t ) (k − k0 ) − , 2a2 (1 + β 2 t 2 ) 1

hk where x0 is the center at t = 0, a is the minimum position spread, v 0 = ¯m0 is the average velocity, and β = h¯ 2 k20

(31) h¯ . 2ma2

The kinetic

energy of such Gauss wavepacket is E 0 = 2m . Actually, the Gauss wavepacket (31) is the analytical solution to the Wigner equation without a Wigner potential [12]. In the following numerical simulations, we take a = 2.825 nm, m = 0.0665me with me being the effective mass of electron, x0 = 15 nm, k0 = 0.7 nm−1 , and the final time T = 22 fs. 5.1. Free Gaussian wavepacket The first experiment utilized in investigating the accuracy of both SEM and MCM, which involves a Gaussian wavepacket in the absence of any external potential, i.e. V (x) ≡ 0, ∀x ∈ X × K. This is an interesting case since the exact analytical

174

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Fig. 1. Free Gaussian wavepacket: Wigner functions at different time instants t = 0, 22 fs. The left-hand-side column is for the reference solution (i.e. the exact solution given in Eq. (31)), while the right-hand-side column for the numerical solution obtained from the stochastic MCM with the number of particles N P = 1.5 M. (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

solution (chosen to be the reference solution) exists. On the other hand, it is also an significant way to validate the signed particle MCM when no potential is applied and, therefore, when the function γ (x) is identically zero over the whole domain. In fact, in this case, the method reduces to simply drifting every signed particle of the initial ensemble defining the initial conditions and no particle creation is required in the process (see definition (26)). That is, the accuracy of MCM does not depend on the chosen time step t, but we still take the same t of SEM for comparison in MCM simulations. We run four MCM simulations with different initial number of particles N P = 0.1, 0.5, 1.5, 4.5 million (M) to see how the accuracy depends on N P . For SEM, we also run three simulations with different mesh size (i.e. the order of expansion) N e = 10, 20, 30 using the same time step t = 1/500 fs to see the convergence behavior with respect to N e . It can be easily checked that this time step satisfies the Courant–Friedrichs–Levy condition [14]. The simulation results for the first numerical experiment are shown in Figs. 1, 2, 3, and 4. The Wigner functions in the phase space at time 0 fs and 22 fs are plotted in Fig. 1. It is found there that the difference between the MCM solution with N P = 1.5 M and the reference one virtually could not be seen with the naked eye. However, as shown in Fig. 2, the relative L 2 (resp. L ∞ ) errors at the final time t = 22 fs are about 0.14 (resp. 0.12) for the Wigner function and about 0.004 (resp. 0.005) for the electron density. It is quite clear that the relative errors are comparatively smaller in the configuration space, which is in accordance with the fact that, when computing the electron density, stochastic noise is smoothed out in the numerical integration over the momentum space (i.e. positive and negative random contributions cancel out). The biggest absolute difference between the MCM solution and the reference one exists in the area where the peaks of the Gaussian wavepacket are located (see Fig. 4(b)). To a stochastic method, such accuracy is

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

175

Fig. 2. Free Gaussian wavepacket: Numerical errors in both Wigner function and electron density for the stochastic MCM with different number of particles N P = 0.1, 0.5, 1.5, 4.5 million (M). (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

very satisfactory for real applications. Actually, an accuracy improvement can be clearly observed in Fig. 2 as the number of particle N P increases, but the improvement is very limited when N P goes from 1.5 M to 4.5 M. Hereafter we will initially take N P = 1.5 M for MCM simulations. Relatively speaking, the cell average SEM provides much more higher accurate numerical solutions for time-dependent Wigner equation than the signed particle MCM. When the order of expansion is N e = 30, as shown in Fig. 3, both L 2 and L ∞ relative errors at the final time t = 22 fs have the magnitude of order 10−10 for both Wigner function and electron density. As expected and depicted in Fig. 4(a), the biggest absolute difference between the SEM solution and the reference one exists in the peak area, too. By increasing the order of expansion N e from 10 to 30 (i.e. refining the mesh), as shown in Fig. 3, we obtain the spectral convergence of L 2 and L ∞ relative errors in both Wigner function and electron density [14]. In the following study, we will take the SEM solution with N e = 30 as the reference one in case of nonexistence of an exact solution. 5.2. Gaussian potential barrier In the second numerical experiment, a Gaussian potential barrier is added to the spatial domain which reads

V (x) = H exp −

(x − x B )2 2w 2

,

(32)

176

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Fig. 3. Free Gaussian wavepacket: Numerical errors (in the common logarithm) in both Wigner function and electron density for the deterministic SEM with different order of expansion N e = 10, 20, 30. (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

where H = 0.3 eV, x B = 30 nm (i.e. the center of the computational domain), and w = 1 nm (the dispersion of the potential barrier). The scattering of the Gaussian wavepacket taken in the first experiment by this Gaussian barrier will be simulated. It should be noted that now we have no exact solution for such scattering process. Instead, we will take the numerical solution obtained from the cell average SEM with the mesh size N e = 30 as the reference solution. But up to now, we  still do not know how large the truncation length Y should be used. As already demonstrated in Ref. [14], the matrix χ (xi ;q,r , y μ ) 1≤i ≤ N ,1≤μ≤ N is sparse and its entries become zero for compact potentials or extremely small for fast e

y

decaying potentials as μ goes to infinity. So we can relate the number of nonzero elements to the size of the support by observing how this number behaves if increasing μ. Once this number keeps unchanged, we will know we have a large enough N y . For the potential barrier considered here, we find that the number of nonzero elements (a matrix entry less than 10−16 is regarded as zero) becomes unchanged if μ ≥ 127. Hence we set N y = 127 and thus Y = N y  y = 38.1 nm. After that, we run three SEM simulations with different mesh size N e = 10, 20, 30 using the same time step t = 1/500 fs. The reference Wigner function in the phase space at 7, 12, 17, 22 fs are plotted in Fig. 5, where we observe during the interaction, the Gaussian wavepacket is separated into two wavepackets moving away from each other; many small oscillations appear around k = 0; and the Wigner function takes negative values in some areas. Fig. 6 gives the error history curves and depicts a very similar convergence behavior to that in the free space shown in Fig. 3. That is, we still now have the spectral convergence of L 2 and L ∞ relative errors in both Wigner function and electron density, which also implies that the numerical solution obtained with the mesh size N e = 30 has really high accuracy and can be used as the reference solution.

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

177

Fig. 4. Free Gaussian wavepacket: Absolute difference  f (x, k, t ) at t = 22 fs between the reference solution and the numerical solution obtained from the deterministic SEM with the mesh size N e = 30 as well as from the stochastic MCM with the number of particles N P = 1.5 M. (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

We also found that the biggest absolute difference between the SEM solution with N e = 10 and the reference one exists in the area around k = 0 (see Fig. 8(a)), which corresponds to the small oscillations observed in Fig. 5. For the signed particle MCM, we take the same truncation for the y-integration in calculating the Wigner potential, set N P = 1.5 M and run seven simulations with different time step and annihilation frequency as follows 1 1 1 (t , N A ) = ( 125 , 125), ( 125 , 500), ( 250 , 125), 1 1 1 1 , 250), ( 250 , 500), ( 500 , 500), ( 500 , 1000), ( 250

where the symbol N A refers to the frequency of annihilation processes happening during the simulation (i.e. every N A time steps an annihilation is performed). Fig. 7 shows the relative errors for the distribution functions in the phase space (i.e. the Wigner function) and in the configuration space (i.e. the electron density). We can clearly observe there that: (a) the errors for t = 1/125 and 1/250 are comparable, but much more smaller than those for t = 1/500, which seems very counterintuitive (usually, the smaller t we take, the more accurate numerical results we will have); (b) the errors are comparatively smaller in the configuration space due to the cancellation of positive and negative random contributions, which already appeared in the free Gaussian wavepacket simulations (see Fig. 2). By close comparison of Fig. 7 and for convenience, we regard (t , N A ) = (1/125, 125) as the best case, (t , N A ) = (1/500, 500) the second-worst case, and (t , N A ) = (1/500, 1000) the worst case. In Fig. 5, we plot the time-dependent Wigner functions in the phase space obtained from the best MCM simulation. As the reader can see there, it is very difficult to spot any qualitative difference from the naked eyes and thus a detailed mathematical quantitative study is necessary. Actually, the biggest absolute difference between the best MCM solution and the reference one also exists in the area around k = 0 (see Fig. 8(b)), which corresponds to the small oscillations observed in Fig. 5. In order to get a feeling about how bad the MCM solutions could be, we plot in Fig. 9 the time-dependent electron density in the configuration space for the best case and the second-worst case against the reference solution. An obvious mismatch between the second-worst case and the reference one is quite clear there, on the contrary, the agreement between the best case and the reference one is excellent. In order to further figure out possible reasons for these counterintuitive accuracy tests, Figs. 10 and 11 focus on the number of (virtual) particles involved during the MCM simulations. The total number is shown in Fig. 10 for different (t , N A ), where the number is reported every 1 fs, and when an annihilation corresponds to the time of counting the particles, then total number after annihilation is reported. One clearly sees an interesting “bottom line” structure, which also shows how the behavior of the curves strongly depends on the parameters chosen. Fig. 11 reports the number of particles created at each generation (note that the colors are consistent with those of Fig. 10). For every initial generation (parents), we report there the total number of particles created in the next generation and we show every generation until no further particle is created (happening usually after no more than 3 generations). As one can easily see from the above figures, the behavior and accuracy of the signed particle Wigner MCM strongly depends on the choice of the two parameters t and N A . Therefore, so far, it has remained complicated to set up a reliable simulation without running several instances with different parameters until a qualitatively acceptable solution is found. In the following section, we discuss in more details the results obtained from these numerical experiments. Indeed, an attentive analysis allows us to clearly depict a set of criteria to efficiently and properly choose these initial parameters.

178

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Fig. 5. Gaussian potential barrier: Numerical Wigner functions at different time instants t = 7, 12, 17, 22 fs obtained from the deterministic SEM with the mesh size N e = 30 (left-hand-side column) as well as from the stochastic MCM with (t , N A ) = (1/125, 125) (right-hand-side column). (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

179

Fig. 6. Gaussian potential barrier: Numerical errors (in the common logarithm) in both Wigner function and electron density for the deterministic SEM with different mesh size N e = 10, 20. Here the numerical solution with the mesh size N e = 30 is chosen to be the reference solution. (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

6. Discussion on parameters for Wigner Monte Carlo method The aim of this section is a detailed investigation of the results described in the previous section in order to understand how to set up reliable and efficient Wigner MCM simulations of a given quantum system. As a first step we discuss the case of a Gaussian wavepacket in the free space, we then proceed with the inclusion of a (Gaussian) potential barrier in the middle of the domain and, finally, we conclude with a list of significant criteria which must be taken into account in order to achieve reliable signed particle Wigner MCM results. 6.1. Free space The free space case, i.e. a Gaussian wavepacket evolving in the phase space without the influence of any external potential, is the first step towards the comprehension of setting up a Wigner MCM calculation. As a matter of fact, when no potential is present (V (x) = 0 identically) then the function γ (x) = 0, ∀x ∈ X × K, therefore preventing the creation of any further particle. The method consists only in evolving the initial set of particles (as Newtonian field-less particles) defining the initial conditions of the system. This situation allows us to focus on the drift process only and see how this affects the overall accuracy of the method. Furthermore, the fact that the exact analytical solution is known in this case makes the accuracy study easier. From Fig. 1 it is practically impossible to find any difference between the results obtained from the exact solution (left-hand-side column) and our MCM solution (right-hand-side) but a closer look to their difference in absolute value,

180

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Fig. 7. Gaussian potential barrier: Numerical errors in both Wigner function and electron density for the stochastic MCM with different (t , N A ). The time steps t reported in the legends are measured in femtoseconds, while the symbol N A refers to the frequency of annihilation processes happening during the simulation (i.e. every N A time steps an annihilation is performed). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

shown in Fig. 4, reveals that the absence of differences is only illusory. In fact, not only the solutions obtained from the deterministic SEM and the stochastic MCM are different, they also compare differently with the exact solution (showing a better agreement in the deterministic case). Fig. 2 shows the situation in even more details by reporting the relative L 2 and L ∞ errors in comparing the deterministic and stochastic Wigner functions (phase space) and probability densities (configuration space). Although these differences exist, the MCM solution still remains a good solution. As a matter of fact, when interested in macroscopic predictions, i.e. in the configuration space, one sees that the errors are also lower than those in the phase space by at least 95%. The situation, of course, drastically changes with the introduction of an external potential in the simulation domain. 6.2. Potential barrier When a potential barrier is included in the spatial domain, the situation becomes physically and numerically different. In particular, from a purely Wigner Monte Carlo perspective, this corresponds to the general situation in which now γ (x) = 0 for (at least) some point of the domain X × K, thus introducing (in a direct way) the process of particle creation and (in a indirect way) the process of annihilation during simulation runs. These two processes can affect the accuracy of the solution and great attention should be given when setting up a new simulation. In more details, in this section we investigate the correlations between the choice of the time step t and the frequency of annihilations with the process of particle creations. We will see that the accuracy can be affected by these choices. Several comments on this topic are given below.

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

181

Fig. 8. Gaussian potential barrier: Absolute difference  f (x, k, t ) at t = 22 fs between the reference solution and the numerical solution obtained from the deterministic SEM with the mesh size N e = 10 as well as from the stochastic MCM with (t , N A ) = (1/125, 125). Here t is the time step and the symbol N A refers to the frequency of annihilation processes. (For interpretation of the colours in this figure, the reader is referred to the web version of this article.)

• The time step. First of all, a proper time step t for the signed particle Wigner MCM has to be chosen. We here assist to a rather counterintuitive process: bigger time steps bring to higher accuracy solutions. This can be explained in terms of the physical processes involved in the simulation. In fact, the presence of a non-zero potential V (x) corresponds to a non-zero γ (x) function, which can be interpreted as a creation rate for signed particles. Therefore, according to the position x of a (virtual) particle the quantity

γ (x(t )) dt is the probability of creating a pair of signed particles in the time interval dt, and one can choose a time δt at which a creation of particles happen by taking a random number r ∈ [0, 1], i.e.,

δt = −

log r

γ (x(t ))

.

(33)

Thus, in order to capture the dynamics of a quantum system the following condition must be valid

t ≥ δt .

(34)

This gives us a first criterion for the choice of t. If condition (34) is not valid, the creation of many signed particles will not happen as the time interval t is too small and, consequently, many of the creation events are not captured, making the MCM simulation unreliable. This becomes clear from an investigation of Figs. 9 and 7. Indeed, from Fig. 9 one sees how a too small t (star-magenta-curve) brings quantitatively and even qualitatively to the wrong solution, while a big enough time step (diamond-black-curve) brings to a very good agreement with the reference solution (continuous-blue-curve). Furthermore, in Fig. 7 one observes (e.g. for a fixed frequency of annihilations N A = 500) the tendency for the accuracy to improve for bigger t. Indeed, the errors in both phase and configuration spaces (and in both relative L 2 and L ∞ 1 1 norms) tends to decrease for the case t = 125 fs (o-blue-curve), and to drastically increase for the case t = 500 fs (∗-magenta-curve). The very same trend, especially in the phase-space comparisons, is observable for the case N A = 125 1 1 fs (diamond-black-curve) and t = 250 fs, in accordance with our above given explanation (condition (34)). when t = 125 • Particle creation and annihilation. Fig. 10 sheds light on further directions to take in our search of criteria for the choice of proper MCM simulation parameters. As matter of fact it shows us that the frequency of annihilation processes N A is as much important as the choice of the time step. It is well-known in the Wigner Monte Carlo community that if the number N A is too low the numerical accuracy of the method is affected since an artificial diffusion of the particles is introduced. Therefore, choices such as N A = 1, while seeming to recover the solution qualitatively, actually provides a definitely numerically non-accurate solution. At the same time, knowing that the creation of signed particles is a exponentially growing phenomenon, one wants to keep the frequency of the annihilation as high as possible. Therefore a good balance must be found between the creation and annihilation processes during a given simulation. This is a non-trivial aspect of the MC method which is clearly illustrated in Figs. 7 and 10. Indeed, from Fig. 7 one clearly observes the effects of the parameter N A over the accuracy of the method. In particular,

182

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Fig. 9. Gaussian potential barrier: Numerical electron density at different time instants t = 7, 12, 17, 22 fs obtained from the deterministic SEM with the mesh size N e = 30 (i.e. the reference solution) as well as from the stochastic MCM with different (t , N A ) = (1/125, 125), (1/500, 500). (For interpretation of the references to colour in this figure, the reader is referred to the web version of this article.)

the case t =

1 250

fs for N A = 125 (triangle-brown-curve), N A = 250 (o-green-curve), and N A = 500 (square-orange-

1 curve) is self-explanatory and confirms our claim. The very same situation can be observed for the case t = 500 fs when N A = 500 (star-magenta-curve) and N A = 1000 (triangle-violet-curve). Furthermore, a very interesting phenomenon seems to emerge from Fig. 10. The best accuracy seems to be obtained when the corresponding curve for the total number of particles involved during a simulation is smooth. This is, indeed, confirmed by the cases, for example, 1 1 t = 125 fs, N A = 125 which is very smooth (diamond-black-curve) and t = 125 fs, N A = 500 which is highly oscillating (∗-blue-curve). Another very interesting observable phenomenon is related to the quantity of quantum information stored into the system. It is easy to see how the annihilation process always forces the total number of particles to a certain smallest amount 1 fs, N A = 125 of particles (i.e. the “bottom line”). For clarity a good example is represented again by the cases t = 125 1 (diamond-black-curve) and t = 125 fs, N A = 500 (∗-blue-curve), it easy to see that, at every point where annihilation happens, the second curve coincides with the first curve. Somehow, it seems that a minimal amount of information is required to simulate a given quantum system and that the annihilation processes is limited to the only task of removing any kind of redundant information. But it also becomes clear that the removal of redundant particles introduces a further source of lost of accuracy. • Particle generations. Finally Fig. 11 provides a further criterion for the choice of the (concurrent) parameters t and N A . 1 First of all, by focusing on the best case (top left-black-curves), t = 125 fs, N A = 125, we note that: a big enough number of particles per each generation must be created to capture the dynamics of a quantum system (in particular we observe

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

183

Fig. 10. Gaussian potential barrier: Number of particles in function of time (every 1 fs). An interesting “bottom line” structure is clearly shown here, which gives the minimum number of particles required to capture the quantum information in the time-dependent simulations. The diamond (black) curve corresponding to (t , N A ) = (1/125, 125), the star (magenta) curve corresponding to (t , N A ) = (1/500, 500) and the back arrow (violet) curve corresponding to (t , N A ) = (1/500, 1000) are related to the best case, the second-worst case and the worst case, respectively (see Fig. 7). If an annihilation process happens at exactly a multiple of 1 fs, then the counting is done after the annihilation has been applied. The time steps t reported in the legends are measured in femtoseconds, while the symbol N A refers to the frequency of annihilation processes (i.e. every N A time steps an annihilation is performed). The solid (resp. dashed) curve in the inset depicts the ratio between the number of particles for (t , N A ) = (1/250, 250) (resp. (t , N A ) = (1/250, 125)) and that for (t , N A ) = (1/125, 125). It can be easily observed there that these three cases have almost the same number of particles during the simulation, which is also shown in the original plot. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

an exponential pattern in the first generation of every test performed); the second generation has to create a big enough number of signed particles but no exponential pattern must be created (and this is further confirmed by the case t = 1 fs, N A = 500 – orange curve – where the second generation is exponential and, indeed, the accuracy in this case is 250 worse); a third generation must be created as well in a big enough total number. Therefore, we have now a set of criteria to control the accuracy of the solutions provided by the MCM and involving

• the time step t, • the annihilation frequency N A , • the initial number of particles N P . First of all, condition (34) must hold. If this rule is not respected then the solution will never be accurate, whatever value the other parameters have. Secondly, in general, bigger time steps give better accuracy. Finally, it is important to note that a balance must be found between the particle creation process, which have to be created in order to reliably catch the physics involved, and the particle annihilation process, which tends to keep the minimal amount of information in the simulated system. These criteria are further developed in the next section. 7. Conclusions A typical deterministic method – the cell average spectral element method (SEM) and a typical stochastic method – the signed particle Monte Carlo method (MCM) have been presented for solving the time-dependent Wigner equation. We demonstrated that SEM is a high-order, conservative, stable and efficient method in one-dimensional spatial space and can be used to provide the reference solution. Precisely because we have high accurate SEM solutions as the reference, a detailed investigation of the performance of the MCM, especially the effect of several parameters (e.g. the time step, the annihilation frequency) on accuracy, has been be done in this work. Indeed, by close comparison and serious study, we can conclude with a set of criteria to make a proper choice of parameters for reliable signed particle Wigner MCM simulations.

• Although counterintuitive, the time step has to be big enough to capture the physics (creation of signed particles): smaller time steps usually mean worse accuracy (exactly the opposite of what happens in classical and deterministic methods).

184

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

Fig. 11. Gaussian potential barrier: Number of generated particles in function of time (every t). Every row refers to the particles created by the previous generation of particles, namely, for every (t , N A ) (i.e. color), the first row is generated by the parent particles, the second row is generated by the first generation and the third row is generated by the second generation. The meaning of the colors is consistent with the legends shown in Fig. 10. The total number of the third generation particles is: 126 for (a), 344 (b), 40 (c), 28 (d), 42 (e), 6 (f) and 12 (g). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

S. Shao, J.M. Sellier / Journal of Computational Physics 300 (2015) 167–185

185

• The creation of particles has to be a smooth process in time (no big oscillations have to be observed in the timedependent total number of particles curve, Fig. 10). Therefore a good balance between the time step t and the annihilation frequency N A must be found.

• At every time step, the process of creation must generate a big enough amount of signed particles in the first generation (exponential grow is allowed). The second generation must be big enough too (but non exponential grow must be present). A third generation must be created with still a big enough number of particles. • The annihilation process seems to kill the redundant information embedded in a quantum system, introducing a lost of accuracy as a consequence. We hope that these criteria offer a better way to set up a time-dependent Wigner simulation with the signed particle MCM, as it allows to remove the usual initial process of trial and error for setting parameters once for all. This is, of course, only a preliminary investigation and further mathematical analysis on the signed particle Wigner MCM are envisaged. Acknowledgements Sihong Shao acknowledges financial support from the National Natural Science Foundation of China (Nos. 11471025, 91330110, 11421101). Jean Michel Sellier acknowledges financial support from the European project EC AComIn (FP7-REGPOT-2012-2013-1). References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31]

E. Wigner, On the quantum corrections for thermodynamic equilibrium, Phys. Rev. 40 (1932) 749–759. C. Jacoboni, P. Bordone, The Wigner-function approach to non-equilibrium electron transport, Rep. Prog. Phys. 67 (2004) 1033–1071. N.C. Dias, J.N. Prata, Admissible states in quantum phase space, Ann. Phys. 313 (2004) 110–146. C. Zachos, Deformation quantization: quantum mechanics lives and works in phase-space, Int. J. Mod. Phys. A 17 (2002) 297–316. J. Hancock, M.A. Walton, B. Wynder, Quantum mechanics another way, Eur. J. Phys. 25 (2004) 525–534. H. Kosina, Wigner function approach to nano device simulation, Int. J. Comput. Sci. Eng. 2 (2006) 100–118. U. Leonhardt, Measuring the Quantum State of Light, Cambridge University Press, New York, 1997. D. Leibfried, T. Pfau, C. Monroe, Shadows and mirrors: reconstructing quantum states of atom motion, Phys. Today (April 1998) 22–28. W.R. Frensley, Wigner-function model of a resonant-tunneling semiconductor device, Phys. Rev. B 36 (1987) 1570–1580. W.R. Frensley, Boundary conditions for open quantum systems driven far from equilibrium, Rev. Mod. Phys. 62 (1990) 745–791. K.L. Jensen, F.A. Buot, The methodology of simulating particle trajectories through tunneling structures using a Wigner distribution approach, IEEE Trans. Electron Devices 38 (1991) 2337–2347. B.A. Biegel, Quantum electronic device simulation, PhD thesis, Stanford University, 1997. R. Kosik, Numerical challenges on the road to NanoTCAD, PhD thesis, Institute for Microelectronics, TU Vienna, 2004. S. Shao, T. Lu, W. Cai, Adaptive conservative cell average spectral element methods for transient Wigner equation in quantum transport, Commun. Comput. Phys. 9 (2011) 711–739. R. Li, T. Lu, Y. Wang, W. Yao, Numerical validation for high order hyperbolic moment system of Wigner equation, Commun. Comput. Phys. 15 (2014) 569–595. A. Dorda, F. Schürrer, A WENO-solver combined with adaptive momentum discretization for the Wigner transport equation and its application to resonant tunneling diodes, J. Comput. Phys. 284 (2015) 95–116. M. Nedjalkov, H. Kosina, S. Selberherr, C. Ringhofer, D.K. Ferry, Unified particle approach to Wigner–Boltzmann transport in small semiconductor devices, Phys. Rev. B 70 (2004) 115319. M. Nedjalkov, P. Schwaha, S. Selberherr, J.M. Sellier, D. Vasileska, Wigner quasi-particle attributes – an asymptotic perspective, Appl. Phys. Lett. 102 (2013) 163113. J.M. Sellier, M. Nedjalkov, I. Dimov, S. Selberherr, A benchmark study of the Wigner Monte-Carlo method, Monte Carlo Methods Appl. 20 (2014) 43–51. J.M. Sellier, I. Dimov, A Wigner Monte Carlo approach to density functional theory, J. Comput. Phys. 270 (2014) 265–277. J.M. Sellier, I. Dimov, The many-body Wigner Monte Carlo method for time-dependent ab-initio quantum simulations, J. Comput. Phys. 273 (2014) 589–597. J.M. Sellier, I. Dimov, On the simulation of indistinguishable fermions in the many-body Wigner formalism, J. Comput. Phys. 280 (2015) 287–294. D. Querlioz, P. Dollfus, The Wigner Monte Carlo Method for Nanoelectronic Devices: A Particle Description of Quantum Transport and Decoherence, Wiley-ISTE, London, 2010. P.A. Markowich, C.A. Ringhofer, C. Schmeiser, Semiconductor Equations, Springer-Verlag, Wien–New York, 1990. W. Cai, D. Gottlieb, C.W. Shu, Essentially nonoscillatory spectral Fourier methods for shock wave calculations, Math. Comput. 52 (1989) 389–410. W. Cai, D. Gottlieb, A. Harten, Cell averaging Chebyshev methods for hyperbolic problems, Comput. Math. Appl. 24 (1992) 37–49. J.P. Boyd, Chebyshev and Fourier Spectral Methods, 2nd edition, Dover, New York, 2001. S. Gottlieb, C.W. Shu, Total variation diminishing Runge–Kutta schemes, Math. Comput. 67 (1998) 73–85. J.M. Sellier, Nano-archimedes, accessed 01 Feb. 2015, www.nano-archimedes.com. M.S. Torres Jr., G. Tosi, J.M.A. Figueiredo, Simulation of the time evolution of the Wigner function with a first-principles Monte Carlo method, Phys. Rev. E 80 (2009) 036701. D. Querlioz, J. Saint-Martin, A. Bournel, P. Dollfus, Wigner Monte Carlo simulation of phonon-induced electron decoherence in semiconductor nanodevices, Phys. Rev. B 78 (2008) 165306.