A Deterministic Solution of the Wigner Transport Equation with Infinite Correlation Length

We propose a new formulation of the Wigner transport equation with infinite correlation length. Since the maximum correlation length is not limited to a finite value, there is no uncertainty in the simulation results owing to the finite integral range of the nonlocal potential term. For general and efficient simulation, the WTE is solved self-consistently with the Poisson equation through the finite volume method and the fully coupled Newton-Raphson scheme. Through this, we implemented a quantum transport steady state and transient simulator with excellent convergence.


I. Introduction
Semiclassical models based on Boltzmann transport equation (BTE) have been widely used for semic onductor device simulation [1]- [4].However, as devices continue to scale down, it is important to cons ider the quantum mechanical effects such as tunneling and quantization [5]- [7].Therefore, semiclassical models are no longer valid at gate lengths of 10 nm or less, and quantum transport models are require d for accurate prediction of device characteristics prediction.
For quantum transport simulation, non-equilibrium Green's function method (NEGF), Wigner transport equation (WTE), and master equation approach are mainly used.Although these models successfully pr edict the characteristics of nanoscale devices, the long computational time and lack of versatility of the simulator are still a critical issue.
In NEGF formalism, it is hard to consider the microscopic scattering mechanism because it requires the inversion of matrix of huge rank because the self-energy terms are generally nonlocal function [11].
Electron-phonon scattering can be efficiently calculated through local approximation, but large computa tional cost is required to include other scattering mechanisms.Also, although this method is well defin ed in steady state, it is inconvenient for transient simulation yet which is very important in device cha racterization [12].
Recently, Pratik B. Vyas et al reported a simulation of the dissipative quantum transport through the Pauli master equation (PME) [13].They show successful simulation results in an ultra-thin body double -gate FET based on the quantum transmitting boundary method (QTBM).This is an attractive model fo r efficiently handling the scattering mechanism, but it can be applied only when the perturbation is we ak and the device length is sufficiently short [13].
As an alternative to the above two methods, we focus on the WTE for the simulation of quantum tr ansport in this work [14]- [16].Transient simulation and dissipative transport simulation are possible bas ed on the WTE, and there have been several studies on this direction [17]- [22].However, more proper formulation of WTE for practical simulations should be introduced.There is a nonlocal potential term i n WTE, which includes an infinite range of integration (infinite correlation length): where χ is real space index, k is momentum space index and U is potential energy.To solve this num erically, previous studies have used finite correlation length, resulting in several problems.W. R. Frensl ey noted that if the computational domain is confined to a finite region, Wigner-Weyl transformation c annot be unitary, because some of the information will be lost.A. S. Costolandski et al confirmed that different simulation results were obtained depending on the correlation length, and mentioned that the appropriate correlation length is different depending on the device structure and there is no simple phy sics-based rule to determine it.As such, the simulation based on the finite correlation length has a pro blem in that there is uncertainty in the simulation result depending on the correlation length and may not be physically consistent with the density operator.
Recently, L. Schulz et al already proposed complex absorbing potential formalism to consider the un bounded computational domain.However, we present an alternative and much simpler method in this p aper.We derive a new formulation with an infinite correlation length by assuming an ideal semi-infinit e reservoir in the contact region.Through reconstruction of the nonlocal potential term, an equivalent e quation with a finite integral range is derived.Since our new formula considers the integral range of n onlocal potential terms up to infinity, it can solve the problem of uncertainty of simulation results acco rding to the integral range.To solve the proposed equation numerically, we use the unwind scheme, fi nite volume method (FVM), and backward Euler method.Through these, quantum transport steady-state and transient simulation with excellent convergence are successfully implemented.By applying our simu lator to resonant tunneling diode (RTD), it was confirmed that reliable results are obtained by showing the plateau region and transient oscillation in unstable bias.

II. Theory
The WTE can be expressed as follow: where the first term of LHS is scattering integral, t is time, ℏ is Dirac's constant, and ε is energy lev el.In general, to solve the WTE numerically, the integral range of the nonlocal potential term (Eq.( 1)) is limited to a finite range.
To consider the infinite correlation length, we first assume an ideal contact condition in which the re servoir is semi-infinitely long and has a constant potential energy.Such boundary conditions are comm only used in quantum transport simulation.In a one-dimensional space, the contact is located at χ=0 a nd χ=L, and the potential energy at the right contact is higher than the left one by Uex.To reconstruct the equation, the nonlocal potential term is divided into the sum of the two terms as follows: where u is a unit step function.In this representation, we just add and subtract the product of Vex and the unit step function to the nonlocal potential term.The integrand function of the first term in Eq. ( 3) becomes an odd function for ζ, and the second term can be calculated through the Fourier transfor m relational expression of the unit step function as follow: ..

 ⎯⎯ → +
(5) Thus, Eq. ( 3) can be rewritten as: When χ is between 0 and L, the integrand of the first term is always 0 if ζ is greater than 2L.There fore, the integration range can be reduced to [0, 2L]: Since the integral range is finitely limited through reformulation of nonlocal potential terms, it is possi ble to solve WTE with infinite correlation length numerically.

III. Simulation methods
To solve our new formulation numerically, we use the finite volume method (box integration metho d).In the steady state, WTE can be expressed as , , ( , ) 0, where W is quantum evolution term and C is collisional term.To apply the upwind scheme, the formu la can be divided into two cases according to the direction of the group velocity: where v is group velocity and the + sign represents when the group velocity is positive and negative, respectively.For easy box integration, the equation is transformed as follows using partial differentiatio n: In x space with a uniform mesh size, the box integration at node xi can be obtained by integrating E q. (13) from xi-1/2 to xi+1/2: 0.5 0.5 0.5 .
To calculate the first term of Eq. ( 14), we need to know the distribution function at xi and xj.The si mplest way to do this is to use the average value of two adjacent nodes.However, in this work we u se the Quadratic upstream interpolation for convective kinematics (QUICK) scheme for high numerical accuracy.The value at the cell face can be calculated as follows through the QUICK scheme: If the calculation range is outside the boundary, it is assumed to have the same value as in the bound ary.Dirichlet boundary conditions apply only to the left if the group rate is greater than 0, and only t o the right if it is less than 0 as upwind method.
The second term of Eq. ( 14) vanishes if the group velocity is constant.However, if the group velocity changes at any point for reasons such as partial varying effective mass, the second term should be cal culated.If the group velocity at xi abruptly changes from A to B, the second term can be written as 0.5 where a is the Dirac delta function which is the derivative of the unit step function.For example, ass uming a parabolic band, if the effective mass at xi changes from m1 to m2, Eq. ( 17) becomes 1 ( , ).
Before describing the quantum evolution term, mesh spacing in k-space should be considered.When E q. ( 13) is integrated over k-space, the equation becomes a continuity equation for charge density.In or der to satisfy the charge conservation, the integral of the quantum evolution term must also be 0. If a uniform mesh size is used, the integral can be expressed discretely as , , 0.
The above equation holds when the Fourier completeness relation is satisfied, and the mesh size at that time can be expressed as follow: The above equation holds when the Fourier completeness relation is satisfied, and the mesh size at that time can be expressed as follow: , where Nk is the number of meshes in k-space.
For calculating the quantum evolution term, we first need to calculate the nonlocal potential term in Eq (10).To calculate the first term of RHS, we assume that the potential changes linearly between me shes as shown in Fig. 2. In this way, the nonlocal potential is calculated through direct integration rath er than discrete integration in order to accurately account for changes in the sin function by position.Since a linear potential is assumed, the equation can be expressed in the form of a product of a sin f unction and a linear function, so that analytical calculations are easily possible.Eventually, the quantum evolution term is calculated discretely as follow: ( , ) ( , ) ( , ) 2 To simply include the scattering effect we use the relaxation time approximation: , ( , ) 1 ( , ) ' ( , ') ' ( , ') Where t is the relaxation time and f is the local equilibrium distribution function.The equilibrium distr ibution function uses the solution when the applied bias is zero.

Electrostatic potential Vp is obtained through the Poisson equation as
where Nd is the doping concentration and n(x) is the electron density obtained from the Wigner functi on: The potential energy v(x) used in the nonlocal potential term can be calculated as follow: ( ) where Uc is the band structure function which considers the band offset considering the barriers and w ells.
Because an iterative solver is required because WTE is a nonlinear system that needs to be solved t ogether with the Poisson equation, we simulate with two methods: the Gummel method [25], which is mainly used as a decoupled scheme, and the Newton Raphson method, which is a fully coupled schem e.
In the case of transient simulation, the calculation method is the same as that of steady-state simulat ion, and it is calculated using the backward Euler (implicit) method.Compared to the forward Euler (e xplicit) method, a much larger time step size can be used and more stable simulation is possible.

IV. Simulation results
For the verification of our new model, we simulate a GaAs-AlGaAs-GaAs resonant tunneling diode a s an example.The emitter and collector are 35 nm and have a doping concentration of 2*1018 /cm3, and the barriers and well are not doped.3nm barriers, 4nm well, and 7 nm spacer are used.The band offset at the barrier is 0.3 eV and the temperature is 77 K.The relaxation time approximation is use d to consider the scattering effects, and the relaxation time is 525 fs.Fig. 1. shows the convergence of the simulation.At both low and high bias, the coupled method sh ows much better convergence.The simulation does not even converge with the decoupled method at 0. 26 V.Of course, it is possible to enforce convergence by making the bias step smaller or multiplying the damping parameter, but the simulation efficiency deteriorates.Because the semiclassical assumption is used in the Gummel method, the convergence is poor in devices with strong quantum effect, but the Newton Raphson method shows quadratic convergence at all bias steps because it calculates an exact response matrix.The plateau region between 0.25 V and 0.32 V can be confirmed through the I-V characteristics in Fig. 2. A plateau appears in the case of a positive bias sweep.This result is obtained through the stea dy-state simulation method, so there is a possibility that an incorrect solution was found during the iter ation process, so transient simulation is performed under the same conditions.Fig. 3(a).shows the transient current characteristics when the bias is raised in the form of a step fu nction by 0.01V by positive sweep.In all other biases except 0.27V, the current converges to the same value as the steady-state simulation result after a sufficient time as at 0.28 V in Fig. 3(b).However, the current oscillates at 0.27 V even after a long period of time.Specifically, it can be confirmed that only bias with negative resistance among the plateau region has unstable transient current characteristic s.
Fig. 4 shows the band diagram at the plateau and the electron density at several biases.Fig. 4(a) sh ows the band diagram in steady-state when forward bias sweep is performed.The solid line from 0.27 V to 0.32 V is the plateau region, and we can know that the characteristic is clearly different from t he normal operation region, which is the dotted line.At the plateau, we can see that band banding oc curs in the emitter region in front of the first barrier.Therefore, the emitter region also shows characte ristics like another quantum well, and a quantized state exists.And when this state is similar to the re sonance energy level between the double barriers, a new current path is formed as shown schematically in Fig. 4(a).Therefore, as shown in Fig 4(b), not only the current in the plateau region but also the electron density in the quantum well does not drop significantly compared to the peak current.
Since our simulator clearly shows the unique characteristics of RTD such as plateau and this is cons istent with previous studies [18]- [22], our new formulation is reliable and is expected to be used as a versatile quantum transport simulator.

V. Conclusion
We derive a novel representation of nonlocal potential terms with infinite correlation length with the assumption of ideal contact.Through this, more accurate simulation is possible without uncertainty of t he WTE solution due to the finite correlation length.Not only steady-state simulation but also transient simulation are possible, and since the Newton-Raphson method is used, the accurate linear response of the equation can be calculated, and thus small signal or noise analysis will also be readily possible.Now that the numerical properties and simulation results in one dimension have been confirmed, it is e xpected that it will be possible to extend to a three-dimensional simulation through a mode space appr oach, and it is expected that more realistic device simulation will be possible.

Fig. 1 .
Fig. 1.Convergence of steady-state simulation according to iteration method.Coupled scheme (red line) shows better convergence than decoupled scheme (black line).

Fig. 2 .
Fig. 2. Current-voltage characteristics obtained through steady-state simulation.0.01V is used as the bias step, and bistability is shown in the case of forward bias sweep (black line) and backward bias sweep (red line).

Fig. 3 .
Fig. 3. Transient current characteristics up to (a) 2000 fs and (b) 15 ps.At all biases except 0.27 V, the current converges to a steady-state solution.(b) At 0.27 V, the current oscillates even after a long time, and the average value of the oscillation current is almost the same as that of the steady-state solution.

Fig. 4 .
Fig. 4. (a) Conduction band diagram and (b) electron density in steady-state when forward bias sweep is performed.