Wigner ensemble Monte Carlo simulation without splitting error of a GaAs resonant tunneling diode

A Monte Carlo technique for the solution of the Wigner transport equation has been developed, based on the generation and annihilation of signed particles (Nedjalkov et al. in Phys Rev B 70:115319, 2004). A stochastic algorithm without time discretization error has been recently introduced (Muscato and Wagner in Kinet Relat Models 12(1):59–77, 2019). Its derivation is based on the theory of piecewise deterministic Markov processes. Numerical experiments are performed in the case of a GaAs resonant tunneling diode. Convergence of the 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
Modelling of electronic transport in nanometer systems requires a theory that describes open, quantum-statistical systems which cannot be provided simply by the Schrödinger equation. Several formulations of quantum transport have been employed, such as those based on the density matrix, non-equilibrium Green's functions, and the Wigner function.
The Wigner function is a real-valued but not necessarily positive definite quasi-distribution, representing a quantum generalization of the Boltzmann distribution function. The Wigner function formalism is attractive as it provides the expression of quantum dynamics in a phase-space formulation, directly comparable with the classical analogue, more intuitive compared with the more abstract density matrix and Green's function approaches.
At the same time the Wigner equation can be augmented by a Boltzmann-like collision operator accounting for the process of decoherence.
The first finite-difference-based solver appears in the mid 1980s (see [6] for a review), and more efficient solvers have been developed nowadays [2,3,7,11,24,28]. We had to wait until the beginning of 2000, to have particle Monte Carlo (MC) solvers for the Wigner equation [20,25]. From that period until now, several papers have been published on this subject (see [27] for a review) and, recently, very interesting device simulations have been provided [4,19,22].
In the realm of the particle Monte Carlo methods, we have focused on the so-called signed Monte Carlo method [21]. Here, the Wigner potential is treated as a scattering source which determines the electron-potential interaction, and consequently new particles with different signs are stochastically added to the system. Recently this method has also been be understood in terms of the Markov jump process theory [12,17], producing a class of new stochastic algorithms. Algorithms that belong to this class are a standard time-splitting algorithm and a new no-splitting algorithm that avoids errors due to time-discretization [13,16,18]. The goal of this paper is twofold: i) to write, for the sake of clarity and dissemination purposes, a simplified set of rules for the splitting and no-splitting algorithms to help a potential user; ii) to perform, from a mathematical point of view, numerical experiments to compare the standard timesplitting algorithm with the the new no-splitting algorithm, for a simplified GaAS Resonant Tunneling Diode (RTD) structure.

The signed particle Monte Carlo method
The Wigner quasi-distribution function f w can be seen as a generalization of the semiclassical Boltzmann distribution function. If we separate the potential energy V tot felt by the electrons into a slowly-varying and rapidly-varying parts V tot = −e + V , the Boltzmann-Wigner transport equation (BWTE) has the form [23] where x ∈ ℝ d is the particle position, k ∈ ℝ d the wave vector (and ℏk the momentum), m * the particle effective mass, satisfying the Poisson equation where e is the absolute value of the electric charge, 0 the absolute dielectric constant, r the relative dielectric constant, N D , N A are the assigned doping profiles (due to donors and acceptors), and n the density The equation includes the Boltzmann scattering operator C(f w ) having the scattering rate (k) [9] and the quantum evolution term where V w is the Wigner potential The Wigner potential is a non-local pseudo-differential operator which is responsible of the quantum transport, is realvalued, and anti-symmetric with respect to k. Solving the Wigner equation, from the numerical point of view, is a quite difficult task. The main complication arising in the direct solution based on finite-difference scheme is the discretization of the diffusion term k ⋅ ∇ x f w due to the typically fast variation in the phase space. Particle-based MC techniques do not require the discretization of this term, but they need costly computational times.
According to the so-called signed particle Monte Carlo approach developed initially in [21], the quantum evolution term (4) looks like the gain term of a collisional operator in which the loss term is missing. But the Wigner potential (5) is not always positive and cannot be considered a scattering term. For this reason, it can be separated into a positive and negative parts V + w , V − w such that In this way, we can define an integrated scattering probability per unit time as and rewrite the quantum evolution term as the difference between gain and loss terms, i.e. Now the term w(k � , k) is interpreted as a new scattering rate which produces, from the old particle, a new pair of particles having weight u and −u . In conclusion, an initial parent particle (with sign) evolves on a free-flight trajectory and, according to a generation rate given by the function (x) , two new signed particles are generated in the same position having weight u and −u respectively. The momentum of the new particles is generated with probability V + w (x, k)∕ (x) . However this procedure suffers of an efficiency issue in particle generation, because usually is a rapidly oscillating function, and exponential growth of particle numbers. In order to contain the particle number, a cancellation procedure must be introduced.

Piecewise deterministic Markov process
A probabilistic model for the quantum evolution term of the BWTE is based on a particle system having time evolution of a piecewise deterministic Markov process [1]. Each particle is characterized by the state space where the first component is the weight u j ∈ {−1, +1} , the second component is the position vector, and the third component is the wave-vector. The time evolution of the particle system (10) is determined by a flow F(z j (t)) and a jump kernel Q(z j (t)) . Starting at a certain state z j (0) , the process performs a deterministic motion according to the flow and the random waiting time until the next jump satisfies 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 In the case of a "jump" in which the jth particle generates two new particles, added to the system such that the jump kernel takes the form [18] and the random waiting time, for the jth particle, satisfies where represents the generation probability. In order to assure finiteness in the integrals with respect to the wave vector, a cutoff parameter c > 0 must be introduced such that Eq. (7) is evaluated in Then the main result is that appropriate functionals of this process satisfy a weak form of the Wigner equation [17].

The splitting algorithm
The time evolution of the system (10) contains both a continuous component (movement in the position space) and a jump component (generation of new particles or phonon scattering). A splitting time step Δt is used in order to separate the transport and the jump processes. The total simulation time is divided into tiny time intervals Δt . Generation of new particles is dictated by Eq. (16). The probability that a particle will survive without generation during a ballistic flight of duration Δt is: Since the condition (20) 1 must be fulfilled during all the simulation time, then we have the following limitation for Δt which must be estimated at the beginning of the simulation from the bias condition, lattice temperature, device dimensions and other input parameters. Finally, this is the generation procedure: s0. For each jth particle s1. Let be r ∈ U[0, 1] do not generate anything, next particle GOTO s0. s2. Otherwise create a random vector k uniformly in (c) . s3. Let be r ∈ U[0, 1] do not generate anything, next particle GOTO s0 s4. Otherwise generate a couple of particle next particle GOTO s0. For the sake of clarity, a simplified set of rules is given in Table 1 for this splitting algorithm. Phonon scattering is tackled using standard MC techniques [9]. The cancellation is performed each time the total particle number N(t) is greater than the cancellation frequency N canc . Couples of particles which belong to the same phase-space cell with opposite affinity are canceled (see [17] for the details). The pros of this algorithm are that the flights progress synchronously in small increments, programming becomes much easier and a number of costly checking procedures, necessary for bookkeeping when dynamics calculations are done for the whole flight, may be avoided. The technique is also very effective for studying transient regimes and for vectorization, since the particles are naturally kept synchronous. The cons is that a discretization error depending on the time step is introduced, which could be minimized by decreasing Δt at the expense of a voracious CPU use.

The no-splitting algorithm
In this case we shall consider the majorant for the generation process (7) as well as one for the total phonon scattering rate and we have the total majorant The random waiting time (16) for all particles is: If ̂ does not depend on x, we have where r ∈ U[0, 1] and the random waiting time is completely determined. In the no-splitting algorithm the transport and the generation process are not separate. Let us indicate T the next observation time.   A set of rules for this algorithm is given in Table 2. The main advantage is that no discretization error due to the time step is introduced.

The resonant tunneling diode
Our goal is to perform numerical experiments with the splitting and no-splitting algorithms introduced in the previous section. We have considered the resonant tunneling diode structure introduced in [26], but with a simplified physics setup. The device simulated is a 1-D (in space) RTD of total length L = 150 nm , which consists of two barriers (with b = 3 nm, a = 0.3 eV) surrounding a quantum well (with b w = 5 nm). The barrier structure is centered in a 30 nm lightly doped ( N D = 10 16 cm −3 ) spacer region that is connected to 60 nm highly doped ( N + D = 10 18 cm −3 drift regions on either side). The quantum well region is symmetric with respect to the mid-point L/2, as shown in Fig. 1, where the distance between each barrier center and the mid-point is The device is considered to be entirely GaAs (with m * = 0.067 ), with scattering mechanisms due to polar optical phonons within a single Γ band [9]. The temperature is 300 K. In our simulation we do not use degenerate statistics in scattering rates and boundary conditions, as done in [26].
Let us suppose that the confining potential V depends only on the longitudinal dimension x. Then we can separate the vector = (x, y, z) in its x longitudinal component and the transversal one, as well as for = (k x , k y , k z ) i.e.
In this case, it is possible to consider a Wigner function of the type where the longitudinal variables come from the singledimensional definition, imposed by the fact that the potential depends only on x, while the transversal variables are introduced by the other parts of the Hamiltonian accounting for phonons. Since the Wigner potential (5), with d = 3 , can be written The vectors � ⟂ = (y � , z � ) and ⟂ = (k y , k z ) belong to the same plane, hence and finally we have  Fig. 1 The quantum well region The total barrier potential is one-dimensional and it is the sum of the single barrier potentials where After some manipulations we get with k → k x . The deltas in the formula (46) mean that the k y , k z of the particle do not change during the creation process. The majorant (13) in this case is and the creation probability (17) writes

Numerical results
We have used an uniform grid in (x, k)-space Other parameters used are: where c is the cutoff in Eq. (18), N ini the initial particle number, N canc the cancellation frequency, N r the repetition number. In the following we shall consider ohmic boundary conditions. The current at the contact is usually obtained by counting charges which are in the contact mesh. But the instantaneous current at the contact is also due to the induced current due to the electron flow, which can be evaluated using to the Ramo-Shockley theorem [5] and, in the 1D case, yields where N elec = elec × Area where elec is the total particle density and Area the cross section area. Then the current density is easily defined as J = I∕Area . Figure 2 shows the current versus the bias voltage, for some splitting time steps Δt and in the no-splitting case. The convergence of the splitting algorithm to the no-splitting one is clearly shown for Δt → 0 . Figure 3 shows the CPU time for both methods: a gain factor of 3 for the no-splitting case can be estimated for the smallest Δt = 0.03125 fs. Figures 4 and 5 show the potential energy V and density, respectively, obtained with the splitting algorithm, corresponding to the peak (on-resonance) and valley (off-resonance) of Fig. 2. The results presented have been obtained using a high performance cluster based on Intel Xeon E5-2620 V3 2.4 Ghz and 32 Gb RAM.

Conclusions
The Boltzmann-Wigner transport equation has been solved by using for the quantum evolution term the signed-particle Monte Carlo method, where a new pair of particles characterized by a sign is generated randomly and added to the system. This generation mechanism has been recently interpreted in terms of the Markov jump process, producing a class of new stochastic algorithms [17]. One of these algorithms has been implemented without time discretization error, and it has been applied to a resonant tunneling diode. The results are compared with the corresponding splitting algorithm, showing an excellent agreement as well as a low-cost computational effort. Future research will develop this MC methodology for the simulation of realistic devices, such as silicon nanowires according to the guidelines in [14,15] and new emerging materials such as graphene [8,10].