A stochastic algorithm without time discretization error for the Wigner equation

Stochastic particle methods for the numerical treatment of the Wigner equation are considered. The approximation properties of these methods depend on several numerical parameters. Such parameters are the number of particles, a time step (if transport and other processes are treated separately) and the grid size (used for the discretization of the position and the wave-vector). A stochastic algorithm without time discretization error is introduced. Its derivation is based on the theory of piecewise deterministic Markov processes. Numerical experiments are performed in a one-dimensional test case. Approximation properties with respect to the grid size and the number of particles are studied. Convergence of a time-splitting scheme to the no-splitting algorithm is demonstrated. The no-splitting algorithm is shown to be more efficient in terms of computational effort.


Introduction. The Wigner equation
was introduced in [18] as an auxiliary tool for specific quantum mechanical calculations. The solution f is real-valued, but not necessarily non-negative. It depends on the time t > 0 , the position x ∈ R d and the wave-vector k ∈ R d , where R d is the d-dimensional Euclidean space. Moreover, is Planck's constant (divided by 2π), m is mass, ∇ is the gradient and the central dot denotes the scalar product. The Wigner kernel V W is determined via the relation where V is potential energy and i denotes the imaginary unit. It is real-valued and anti-symmetric with respect to k . The solution f is related to the solution ψ of the 60 ORAZIO MUSCATO AND WOLFGANG WAGNER Schrödinger equation where ∆ is the Laplace operator. In particular, under some restrictions on ψ , the function f satisfies The initial condition is where f 0 is an integrable function. The Wigner equation turned out to be a convenient tool for modelling quantum effects in nanoelectronic devices, since it can be coupled easily to the scattering part of the semiconductor Boltzmann equation ( [6], [5], [12]). Comprehensive presentations of the field and, in particular, of numerical approaches to the Wigner-Boltzmann equation (including extensive lists of references) are given in [14] and [13]. Basic ingredients of the so-called Wigner Monte Carlo method are stochastic algorithms for solving the Wigner equation. The most common algorithms are using signed particles (weights ±1) and treat the integral with respect to the Wigner kernel via the generation of pairs of particles with opposite signs. Among the numerous recent studies related to this approach we mention [16] (comparison with a deterministic Wigner solver) and [2] (distributed-memory parallelization).
The purpose of this paper is to present an algorithm based on the "random cloud model" introduced in [17], which avoids any time discretization error. This approach via piecewise-deterministic Markov processes is efficient due to the introduction of a majorant for the Wigner kernel, which makes the corresponding rate function independent of the position. The idea with the majorant (leading to fictitious creation events) was used before in [10], where algorithms with splitting of transport and creation processes were studied. The paper is organized as follows. The algorithm is introduced in Section 2. Results of numerical experiments are provided in Section 3. Comments are given in Section 4.

Algorithm.
A "random cloud model" for the Wigner equation (1) was introduced in [17]. This is a class of piecewise deterministic Markov processes of the form Each process is a system of N (t) particles, which are characterized by a weight u j (t) , a position x j (t) and a wave-vector k j (t) . We use the notations for the states of single particles and for the states of the process.
The initial state of the processes (6) is generated according to the initial condition (5) of the Wigner equation. More precisely, the initial number of particles is deterministic, N (0) = N ini . For j = 1, . . . , N ini , the positions x j (0) and wave-vectors k j (0) are generated independently according to |f 0 | . The weights are obtained as u j (0) = signf 0 (x j (0), k j (0)) . Without loss of generality, we assume that |f 0 | is a probability density.
The time evolution of the processes (6) is determined by a flowF and a jump kernel Q . Starting at statez ∈ Z , the process performs a deterministic motion according toF . The random waiting time τ until the next jump satisfies where P denotes the probability measure and Then the process jumps into a new stateκ ∈ Z distributed according to The flow has the form (cf. (8)) so that particles move independently of each other. The single particle flow F is (cf. (7)) is the velocity corresponding to k . There is considerable freedom in choosing the jump kernel Q . The following choice is particularly well suited for numerical purposes. We introduce a majorant V W such that and consider jump kernels of the form where J(z, j, k) = z 1 , . . . , z N , z 1 (z, j, k), z 2 (z, j, k) and The majorantV W is a remaining degree of freedom. The kernel (13) without a majorant (or, withV W = |V W |) was considered in [17]. The kernel (13) with a majorant was used in [10], [7] in the context of a class of time splitting algorithms (transport and creation steps were separated).
According to (13), the waiting time parameter (10) takes the form so that the random waiting time satisfies (cf. (9)) When the majorantV W does not depend on x , then one obtains the following algorithm performing the evolution of the particle system on a time interval [0, T ] .
No-splitting algorithm. 0. The system time is denoted by t . The individual particle times t j , j = 1, . . . , N , indicate the moments of the last update of the corresponding position. Put t = 0 and t j = 0 , j = 1, . . . , N .
1. Generate the random waiting time τ until the next creation event, Put t := t + τ . 2. If t ≥ T , then update all positions (up to T ), Stop the algorithm. 3. Choose an index j ∈ {1, . . . , N } uniformly. Update the position and the individual time of the corresponding particle (up to t), 4. Generate a vectork according to the probability density 1 2γV 5. Check if the creation event is fictitious. With probability go to 1. 6. Create a pair of particles Put N := N + 2 and the individual times of these particles equal t .
7. If N ≥ N canc , then update all positions and individual times according to (16) and perform cancellation. This means that pairs of particles with similar positions and wave-vectors, but with opposite signs, are removed from the system. A detailed description of our cancellation procedure is given in [10, Sect. 2.1.2]. Go to 1.
Functionals. The above algorithm performs the time evolution of the stochastic process. Another part of the numerical procedure is taking measurements on the particle system. The connection between the particle system and the Wigner equation is given by the formula (cf. [17, Theorem 2.1]) where ϕ is an appropriate test function and E denotes mathematical expectation.
The expectation in (20) is approximated as the empirical mean over N rep independent runs of the algorithm. A functional of particular interest is the density (cf. (4)) In order to apply (20), the functional (21) is replaced by a spatially averaged version, with the test function where D ⊂ R d is a spatial cell, |D| denotes the volume and the symbol χ A denotes the indicator function of a set A .
Assumptions. The basic property (20) was established under the assumption (cf. According to (12), assumption (23) implies Example 1. For the one-dimensional rectangular barrier, the Wigner kernel has the form Assumption (24) is not fulfilled in Example 1, which is one of the common test cases. Therefore, we mention a version of the no-splitting algorithm, which works also in this case. We assume that the potential energy is integrable, which is fulfilled in Example 1. The modified algorithm is obtained by replacing the Wigner kernel V W by its truncation one obtainsγ so that assumption (23) is fulfilled.

Numerical experiments.
Here we perform numerical experiments with the no-splitting algorithm introduced in the previous section. First we introduce the test case and specify the various components of the algorithm. Then we study the discretization error due to the cancellation procedure. Finally we consider a timesplitting algorithm and illustrate its convergence with respect to the time step as well as its efficiency compared to the no-splitting algorithm.
3.1. Test case. We consider the one-dimensional Gaussian barrier The Wigner kernel has the form (cf. Remark 1 below) where The initial state is where x 0 , k 0 ∈ R and σ 0 > 0 . The corresponding initial state of the Schrödinger equation (3) is Algorithm. The generation of the initial set of particles is straightforward, since In particular, all initial particles are positive. We use the majorant (cf. (27)) which corresponds to the rateγ Note that condition (23) is fulfilled. According to (27) and (30), the steps (17)-(19) of the no-splitting algorithm take the form 4. Generate a vectork according to the probability density go to 1. 6. Create a pair of particles Put the individual times of these particles equal t and N := N + 2 . The approximation parameters of the no-splitting algorithm are related to the cancellation procedure (step 7). These parameters influence the accuracy and the frequency of the cancellations. The main parameter is the "cancellation grid". We use intervals [x min , x max ] and [k min , k max ] , which are divided, respectively, into N x and N k equal subcells. The first and last subcells are extended to infinity. Further parameters are the initial number of particles N ini and the cancellation bound N canc .
Functionals. We introduce several functionals, which illustrate the time evolution of the system. These quantities are the total average position (cf. (21)) the total average velocity (cf. (11)) and the total average kinetic energy where is the density with respect to the wave-vector. According to (29), one obtains the initial values The observation parameters of the no-splitting algorithm influence the measurements on the particle system. An "observation time step" ∆t obs is used for measuring the time evolution of the functionals. "Observation grids" are used for measuring functionals depending on the position or/and the wave-vector. They are determined by the parameters x obs min , x obs max , N obs where x ∈ [x obs j−1 , x obs j ) and j = 1, . . . , N obs x .

Cancellation error.
Here we study the approximation error in the no-splitting algorithm. This error is due to the cancellation of particles. We consider the test case (26), (28) with the parameters (as in [16]) (44) Figure 1 shows the position density (37), which is very close to the deterministic reference solution (39). Next we illustrate the dependence of the approximation error on the cancellation parameters. In these tests the measured position densities (37) are compared with the deterministic reference solution (39). Other functionals are compared with the corresponding results for the parameters (44), which provide a stochastic reference solution. Grid in the wave-vector space. For the parameters (44) there is ∆k = 0.05 nm −1 . Now we choose the parameters (44) with N k = 100 so that ∆k = 0.2 nm −1 . The error for this set of approximation parameters is illustrated using several functionals. Figure 2 shows the position density (37), which differs quite significantly from the deterministic reference solution (39). Figure 3 shows the wave-vector density (38) and the stochastic reference solution (parameters (44)). The error is quite significant, though the qualitative behaviour is described correctly. Figures 4-6 show the total average position, velocity and kinetic energy (32)-(34) as well as the corresponding stochastic reference solutions (parameters (44)). These quantities are measured at the cancellation times, with only one repetition. The figures show how the error accumulates with time. In particular Figure 6 shows that a coarse grid (i.e. N k = 100) has the effect to increase the particle energy.
Grid in the position space. For the parameters (44) there is ∆x = 0.15 nm . Now we choose the parameters (44) with N x = 100 so that ∆x = 0.6 nm . Figure 7 shows the position density (37), which in some region fluctuates around the deterministic reference solution (39).
Initial number of particles and cancellation bound. Finally we study the influence of the parameters N ini and N canc . Figure 8 shows the curves for the numbers of particles before and after cancellations versus the number of calls to cancellation for the parameters (44). The lower curve depends on N ini and on the cancellation grid, while the upper curve is determined by N canc . The corresponding curves for  20  358k  256  100 400 160k 480k  14  273k  218  400 100 160k 480k  15  301k  229  400 400 40k 480k  7  174k  151  400 400 160k 960k  9  402k  345  Table 1. Properties of the algorithm for various sets of cancellation parameters. The quantity "calls" denotes the number of calls to the cancellation procedure and N after is the average number of particles after the last cancellation.
other sets of parameters are not shown, but the relevant information is collected in Table 1.
The position density (37) was measured for the parameters (44), with N ini = 40000 and N canc = 960000 , respectively. We do not show the corresponding figures, since they are rather similar to Figure 1. In the first case there is a certain increase in fluctuations. Some other properties for these sets of parameters are provided in the last two lines of Table 1.
We conclude with some general remarks concerning the choice of the parameters N ini and N canc . Both parameters should not be below a certain level (depending on the cancellation grid). Decreasing the parameter N ini leads to larger fluctuations. Decreasing the parameter N canc makes cancellation inefficient and finally leads to a blow-up of the algorithm. Increasing these parameters leads to more effort.

Comparison with a time-splitting algorithm.
Here we compare the nosplitting algorithm with a time-splitting algorithm. We study convergence with respect to the time step as well as the effort in terms of the CPU time. We consider the test case (26), (28) with the parameters (41), (42). The observation parameters are given in (43). For reference purposes we provide an algorithm with time splitting (from [10]). In this algorithm a time step ∆t is used in order to separate the transport and the    Table 2. Properties of the time-splitting algorithm and the nosplitting algorithm. The quantities "err-max" and "err-aver" denote, respectively, the maximum and the average (over the cells) of the absolute differences between the measured density and the reference solution. The last column provides the measurements of the first, second and third cancellation time.
creation processes. The algorithm performs the evolution of the system (15) on the time interval [0, ∆t] according to the following steps: 1. Transport step The positions change according to The time-splitting algorithm is applied with the cancellation parameters (44) and various time steps. Figure 9 shows the position density (37), which approaches (for a decreasing time step) the deterministic reference solution (39). Further measurements are collected in Table 2, where the last line corresponds to the no-splitting algorithm. In particular, the data (third and fourth column) show that the error of the time-splitting algorithm quantitatively approaches the error of the no-splitting algorithm, when the time step decreases.
The last column of Table 2 provides the measurements of the first, second and third cancellation time. It is rather obvious that the convergence cannot be better than of the order ∆t , which is indeed observed. The value of the first cancellation time for the no-splitting algorithm is very close to the prediction 1.8925, obtained with (49) in the Remark 2. The prediction is so accurate, since the function γ is constant in most of the region, where the solution is concentrated. The function γ is shown in Figure 10. Note that γ(∞) ∼ 0.29 (×10 15 sec −1 ). Taking into account the average numbers of particles after the first and second call to cancellation, one might also make (less accurate) predictions for the second and third cancellation time.
Finally, the second column of Table 2 provides some insight into the important efficiency issue. The effort consists of three components: "transport and creation", "cancellation" and "measuring functionals". For the time-splitting algorithm, the first component is (roughly) inversely proportional to the time step. Thus, for small time step, this component is responsible for the major part of the effort. This gives the no-splitting algorithm a rather significant advantage.
Remark 3. Table 2 shows clearly that our time-splitting algorithm (cf. (45), (46)) convergences towards the reference solution (and, thus, to the results provided by the no-splitting algorithm) when time step decreases. But in [16] the authors report that "bigger time steps bring to higher accuracy solutions" and claim that "this can be explained in terms of the physical processes involved". This strange behaviour (error increasing with decreasing time step) is due to a wrong splitting rule. In fact, according to the algorithm in [16, p.173], each particle always creates some offspring during a time step. But the correct rule is that each particle creates offspring with a certain probability, as in (46).
Our algorithm has some more special features compared to other signed-particle algorithms for the Wigner equation (cf., e.g., [2]). In particular, a discretization of the state space is used only in the cancellation procedure, while transport and creation are treated with a continuous state space. An efficiency gain is obtained by the introduction of fictitious creation events via appropriate majorants and by using |V W | instead of V + W in the creation procedure (cf. (17), (18)). Some of these aspects, including the advantages of our cancellation procedure, were discussed in [10] for the corresponding time-splitting schemes.
The no-splitting algorithm is convenient for being combined with scattering processes in the context of the Wigner-Boltzmann equation. In this case a constant majorant has to be used (cf. (25)). In a similar way it should be possible to extend the algorithm towards the Wigner-Fokker-Planck equation (cf., e.g., [3], [1]). Here a Wiener process is added to the transport part, in analogy with the heated inelastic Boltzmann equation mentioned above. An interesting direction of further studies concerns more sophisticated cancellation procedures. A straightforward idea is to use a self-adjusting cancellation grid. This means that, given ∆x and ∆k , the parameters k min , k max , x min , x min are adapted to the (evolving in time) cloud of particles. Another opportunity would be using a non-uniform cancellation grid, i.e., taking a finer grid in regions with high densities or/and large gradients. This might be of particular importance in multi-dimensional situations.