Wigner Ensemble Monte Carlo Simulation Without Discretization Error of a GaAs Resonant Tunneling Diode


 A Monte Carlo technique for the solution of the Wigner transport equation has been developed during these years, based on the generation and annihilation of signed particles. A stochastic algorithm without time discretization error has been recently introduced. 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 from the Schrödinger equation. Several formulations of quantum transport have been employed practically, 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 Boltzmann distribution. The Wigner function formalism is attractive as it allows 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 Boltzmannlike collision operator accounting for the process of decoherence.
The first finite-difference based solver appears in the mid 1980s (see [5] for a review), and more efficient solvers have been developed nowadays [8], [20], [2], [24], [6], [18]. We have to wait until the beginning of 2000, to have particle Monte Carlo (MC) solvers for the Wigner equation [21,15]. From that period up to now, several papers have been published on this subject (see [23] for a review) and, recently, very interesting device simulations have been provided [3,14,17].
In the realm of the particle Monte Carlo methods, we have focused ourselves on the so called Signed Monte Carlo method [16]. 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,9,13], producing a class of new stochastic algorithms. In this paper a thorough validation of a stochastic algorithm without time discretization error has been introduced in the case of a GaAs Resonant Tunneling Diode.

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 and rapidly parts V tot = −eϕ+V , the Boltzmann-Wigner transport equation (BWTE) has the form [19] where x ∈ R d is the particle position, k ∈ R d the wave vector (and k the momentum), m * the particle 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 total scattering rate λ(k) [7] and the quantum evolution term where V w is the Wigner potential The Wigner potential is a non-local potential operator which is responsible of the quantum transport, is real-valued, 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 [16], 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 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 ma-jorantV w such that In the case of a "jump" in which the j-th particle generates two new particles, added to the system such that and the random waiting time, for the j-th 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. (17) is evaluated in Then the main result is that appropriate functionals of this process satisfy a weak form of the Wigner equation [12].

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: Ifγ 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 j-th particle s1. Let be r ∈ U [0, 1] do not generate anything, next particle GOTO s0. s2. Otherwise create a random vectork uniformly in B(c) . s3. Let be r ∈ U [0, 1] do not generate anything, next particle GOTO 2.0 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 [7]. The cancellation is performed each time the total particle number N (t) is greater than the cancellation frequency N canc . Couple of particles which belong to the same phase-space cell with opposite affinity are canceled (see [12] 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 (17) as well as one for the total phonon scattering rate and we have the total majorant Γ = Γ s +γ .
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. 4. Choose an index j ∈ {1, ...., N } uniformly. Move the j-th particle up to t solving Newton's equation of motion. If the particle is out of the boundary erase it and GOTO 2, otherwise update the individual time of the j-th particle into t j := t. The component u j does not change. The position and wave-vector are x j , k 0 = (k x0 , k y0 , k z0 ) respectively.

Let be r ∈ U [0, 1]
if r ≤γ Γ GOTO 6 ELSE GOTO 10 (31) 6. Create a vectork according to the probability density 1 2γV 7. Check if the creation event is fictitious. Let be r ∈ U [0, 1] The individual times of these particles equal t. Put N := N + 2 . 9. If the total number of particles exceeds a certain bound, then move all particles up to time t according to Newton's equations. Erase particles which are out of the boundary. Pairs of particles with similar positions and wave-vectors, but with opposite signs, are removed from the system. GOTO 2 10. check/perform phonon scattering, GOTO 2 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 nosplitting algorithms introduced in the previous section. We have considered the Resonant Tunneling Diode (RTD) structure introduced in [22], but with a simplified physics setup. The device simulated is a 1-D (in space) RTD of total lenght L = 150 nm, which consists of two barriers (with b = 3 nm, a =  Table 2 The no-splitting algorithm 1. Start with initial conditions, observation time T , the bias Voltage. Put t = 0 2. t = t + T , if t ≥ t f inal STOP otherwise continue 3. Drift/Generate/Scattering particles in T according to the DGS procedure 4. update b.c. 5. assign charge in the grid points and evaluate density 6. solve Poisson equation 7. goto 2 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 figure 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 [7]. The temperature is 300 K. Let us suppose that the confining potential V depends only on the longitudinal dimension x. Then we can separate the vector x = (x, y, z) in its x longitudinal component and the transversal one, as well as for k = (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 single-dimensional 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 vectors x ⊥ = (y , z ) and k ⊥ = (k y , k z ) belong to the same plane, hence and finally we have 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 iŝ and the creation probability (17) writeŝ

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 [4] and, in the 1D case, yields where N elec = ρ charge × Area where ρ charge is the total charge density and Area the cross section 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.

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 are 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 [12]. 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 algo, showing an excellent agreement as well as a low-cost computational effort. Future researchers will develop this MC methodology for the simulation of realistic devices, such as silicon nanowires according to the guidelines in [10,11].