A higher-order accurate operator splitting spectral method for the Wigner-Poisson system

An accurate description of 2-D quantum transport in a double-gate metal oxide semiconductor filed effect transistor (dgMOSFET) requires a high-resolution solver to a coupled system of the 4-D Wigner equation and 2-D Poisson equation. In this paper, we propose an operator splitting spectral method to evolve such Wigner-Poisson system in 4-D phase space with high accuracy. After an operator splitting of the Wigner equation, the resulting two sub-equations can be solved analytically with spectral approximation in phase space. Meanwhile, we adopt a Chebyshev spectral method to solve the Poisson equation. Spectral convergence in phase space and a fourth-order accuracy in time are both numerically verified. Finally, we apply the proposed solver into simulating dgMOSFET, develop the steady states from long-time simulations and obtain numerically converged current-voltage (I-V) curves.


Introduction
In the last two decades, the Wigner function approach [1,2] has provided a powerful tool for studying quantum effect in various electronic devices, such as the resonant tunneling diodes (RTDs) [3] and the metal oxide semiconductor filed effect transistors (MOSFETs) [4]. A coupled system of the Wigner equation and the Poisson equation is usually adopted for taking the space charge effects into account. Finite difference methods were often used to obtain numerical solutions of the Wigner equation [5,6] as well as of the Wigner-Poisson (WP) system [7,8], and several spectral methods were also tried [3,9,10]. In order to accurately capture 2-D quantum transport in a doublegate MOSFET (dgMOSFET), the WP system in 4-D phase space is required to be integrated with high resolution. However, all above-mentioned numerical methods were implemented in 2-D phase space, and highly accurate deterministic numerical methods for the WP system in 4-D phase space are very few up to now. This paper is intended to fill this gap by exploiting a recently developed operator splitting spectral method for the 4-D Wigner equation in quantum double-slit interference [11]. Specifically, we will take advantage of the operator splitting spectral method to solve the 4-D Wigner equation, in which the semi-discrete models resulted from spectral expansion in phase space for the sub-equations have analytical solutions, and continue to use a Chebyshev spectral method to solve the 2-D Poisson equation.
Detailed benchmark tests are performed with the Gaussian barrier scattering in 2-D and 4-D phase space, and demonstrate that the proposed operator splitting spectral method indeed has a spectral accuracy in phase space and a fourth-order accuracy in time. We also show that the electric field induced by the space charge has a great effect on the rate of quantum tunneling. After calibration, we apply our high-resolution solver into simulating RTD and dgMOSFET. Numerical experiments show that the steady states can be well developed from long-time simulations and the corresponding current-voltage (I-V) curves are numerically converged as the number of collocation points increases.
The remainder of this paper is organized as follows. Section 2 briefs the WP system. Section 3 presents the operator splitting spectral method. Section 4 conducts benchmark tests with the Gaussian barrier scattering. Simulations and discussions of RTD and dgMOSFET are given in Sections 5 and 6, respectively. The paper is concluded in Section 7 with a few remarks.

The Wigner-Poisson system
The Wigner function f (x, k, t) living in 2d-D phase space: (x, k) ∈ R 2d with position x ∈ R d and wavevector k ∈ R d , obeys the following Wigner equation [1] where d gives the dimension of position space, t denotes the time, is the reduced Planck constant, m is the mass, and Θ V [f ] is the so-called nonlocal pseudo-differential operator containing all the quantum information: Here V (x, t) gives the external potential, and can be rewritten into V (x, t) = V b (x) + V e (x, t) when taking the space charge effects into account, where V b (x) denotes the conduction band potential and V e (x, t) the effective electric potential. Actually, V e (x, t) can be determined by a Poisson equation with the electron density as its source term: Figure 1: A 16 nm dgMOSFET structure [15]. The gate length L G , equivalent gate oxide thickness EOT and silicon channel thickness T si are 6 nm, 1 nm and 3 nm, respectively. The source and drain doping is 5 × 10 19 cm −3 . The transistor is assumed to be wide, i.e., the y-direction is treated as infinite long.
where q e denotes the positive electron charge, (x) is the dielectric constant, N d (x) is the doping density and n(x, t) denotes the density of electrons given by And, the current density J (x, t) can be further calculated by In this work, we focus on developing a high-resolution solver for the WP system in 4-D phase space (i.e., d = 2) and let x = (x, z), k = (k x , k z ). In particular, our target is to simulate the dgMOSFET (as shown in Fig.1) and the working-equations read where we have chosen the commonly used inflow boundary conditions for the Wigner equation [5,12], and mixed boundary conditions for the Poisson equation: the Dirichlet boundary at the source/drain in x-direction and the gates in z-direction plus the Neumann boundary at the Oxide/Air interfaces in z-direction, the computational domain for position is Ω x = [x l , x r ] × [z l , z r ], V gu and V gl give the upper and lower gate voltage, respectively, V ds is the source/drain bias potential, and N + D refers to the ionized donor doping concentration.

Numerical methods
Considering the decay property of the Wigner function when |k| → +∞, a simple nullification outside a sufficiently large k-domain Ω k is usually adopted [11,13,14], thus we are in fact using a truncated pseudo-differential operator Θ T V [f ] in k-space as follows where Ω k = [k x,min , k x,max ] × [k z,min , k z,max ] and y x,µ = µ∆y x , y z,ν = ν∆y z with ∆y x , ∆y z being the spacing, which satisfy ∆y i = 2π/L k i with L k i = k i,max − k i,min for i = x, z in this paper.

Solving the 4-D Wigner equation
The operator splitting spectral method developed in [11] for simulating the quantum double-slit interference is employed here for solving the 4-D Wigner equation. A brief description is given below and the interested readers are referred to [11] for more details.
An s-stage exponential operator splitting method for the Wigner equation given in Eq.
(1) reads where f n (x, k) := f (x, k, t n ) denotes the exact solution at time t n := n∆t and s j=1 e a j ∆tA e b j ∆tB f n (x, k) gives the corresponding numerical solution. Here A, B are the convection operator and pseudo-differential operator, which correspond to two sub-equations of the Wigner equation, respectively: We adopt the advective approach to march the sub-equation (A) in Eq. (10) strictly along the characteristic lines as follows and the Chebyshev expansion of the Wigner function with respect to x is used to obtain function values at shifted points. Motivated by the intrinsic nature of Fourier transformation contained in the pseudodifferential term Eq. (8), we use a Fourier spectral method to solve the sub-equation (B) in Eq. (10). The interpolation operator I k,N reads where ψ ν i (k i ) = e 2πiν i (k i −k i,min )/L k i are the Fourier basis functions and N i the number of collocation points in k i -space for i = x, z. Substituting the interpolation function I k,N f (x, k, t) into the pseudo-differential term Eq. (8) also yields spectral approximation Accordingly, the orthogonal relation of the Fourier basis functions implies ∂ ∂t a νx,νz (x, t) = c νx,νz (x)a νx,νz (x, t), the solution of which has the following explicit form a n+1 νx,νz (x) = e cν x,νz (x)∆t a n νx,νz (x).
To match with the spectral accuracy in phase space, we adopt a fourth-order splitting scheme with s = 4 in Eq. (9): In subsequent numerical experiments, we usually choose N x = N z := N for convenience.

Solving the 2-D Poisson equation
The Chebyshev expansion in x direction continues to be used to solve the Poisson with Dirichlet boundary conditions: where the function V gl (x) (resp. V gu (x)) reduces to the lower (resp. upper) gate voltage V gl at the gate, and vanishes otherwise. In order to achieve the spectral convergence, we use a cubic polynomial to smoothly connect, for example, V gl to 0.
can be approximated by the truncated Chebyshev series as follows where φ n give the Chebyshev polynomials of the first kind.
For simplicity, we suppose that the numbers of collocation points in x and z directions are even and the same, denotes by M . It can be readily verified that the expansion coefficients satisfy the following relationships where c 0 = 2 and c n = 1 for n ≥ 1. The collocation equations for {a nm } that follow from Eqs. (13)- (14) are then n odd a nm = 0, n even m odd We define the column vectors X i , B i for i = 0, 1, . . . , M by and let P be the (M + 1) × (M + 1) matrix as shown in Eq. (31) in Appendix. Consequently, Eqs. (15)- (17) can be rewritten into The solution process and related specific solution form are detailed in Appendix. Although the size of the coefficient matrix in Eq. (18) is (M + 1) 2 × (M + 1) 2 , only the calculation of sub-matrices of order (M + 1) × (M + 1) is involved rather than directly inverting the original matrix. Therefore, the proposed Chebyshev spectral method for the 2-D Poisson equation is not only highly accurate but also efficient. It should be pointed out that we select the Dirichlet boundary condition above just for an example and the proposed numerical solver for the Poisson equation is able to deal with all kinds of boundary conditions.
In summary, we evolve the WP system (7) in 4-D phase space as follows Step I. Calculate the potential V e (x, t = 0) with initial density n(x, t = 0) via the 2-D Poisson equation (Eq. (4)) by using Chebyshev spectral methods; Step II. Using the obtained potential V e (x, t = 0) to solve the time-dependent 4-D Wigner equation with operator splitting spectral method to obtain f (x, k, ∆t) and then to calculate the density n(x, ∆t) via Eq. (5); Step III. Calculate the potential V e (x, ∆t) with the density n(x, ∆t), repeat Step I and Step II until to the final time t f .

Calibration
In this section, we first would like to verify the convergence rate and efficiency of the proposed solver. The L 2 -error ε 2 (t) and L ∞ -error ε ∞ (t): are employed to study the convergence rate in terms of the number of collocation points and the time step, where Ω = Ω k × Ω x gives the computational domain in 4-D phase space, f num and f ref denote the numerical solution and reference solution, respectively. To conveniently visualize the 4-D Wigner function, we plot the reduced 2-D Wigner function [13] in this paper as follows In addition, the Chebyshev collocation points in x-direction for the 4-D Wigner equation and 2-D Poisson keep the same, which may avoid additional interpolations when calculating the electron density n(x, t) in Eq. (4).
the reference solution of which is V (x, z) = (x 2 + z 2 )e x 2 +z 2 and the right terms are Table 1 gives the calculation time (the second column) and L ∞ -error (the third column) under the different number of collocation points (the first column). We set M = 8, 16, 32, 64 and 128. When M is equal to 64 and 128, the L ∞ -error has reached 10 −14 , but the calculation time only takes less than 0.1 second. The calculation time in the Table 1 is the serial time with 1 CPU (Intel Core TM i7-8550U CPU @ 1.80GHz). That is, the proposed Chebyshev spectral method for the 2-D Poisson equation is not only highly accurate but also efficient. The right plot in Table 1 clearly shows the spectral convergence with respect to M .  To further validate the overall performance of the operator splitting spectral method for the WP system, we simulate the Gaussian barrier scattering of the Gaussian wave packet (GWP) [11,13] to investigate its convergence rate. We first make tests in 2-D phase space: , Ω x = [−25 nm, 25 nm], and adopt the initial GWP as where x 0 is the center, a the minimum position spread and k 0 the initial wavenumber. The Gaussian barrier reads with ω = 1 nm, x b = 0 and H = 2.3 eV. The other parameters are: x 0 = −10 nm, k 0 = 1.4 nm −1 , a = √ 2 nm, = 1 eV · fs, the effective mass m e = 1 eV · fs 2 · nm −2 and the final time t f = 20 fs. The Poisson equation satisfies the Dirichlet boundary condition with bias potential V 0 = 0.5 eV, the dielectric constant = 10 Fm −1 and the doping density N d (x) = 0.
In order to study the convergence rate with respect to N (resp. M), we fix M = 400 (resp. N = 128) and ∆t = 0.01 fs. As shown in the left and middle plots of Fig. 2, the proposed splitting spectral method shows the spectral convergence with respect to both N and M . The right plot of Fig. 2 further displays the fourth-order convergence rate with respect to ∆t on a fixed mesh (N, M ) = (128, 400).
Next, we would like to use such Gaussian barrier scattering to study the effect of the space charge on quantum tunneling. The tunneling rate P r (t) [14] P r (t) = and = 10 Fm −1 , the tunneling rate is almost proportional to the bias potential V 0 and higher than the value 0.0353 indicated by the red line, which is the rate for the case without coupling the Poisson equation, as shown in the left plot of Fig. 3. The middle plot of Fig. 3 gives the negative correlation between the tunneling rate and the doping density N d when the bias voltage V 0 = 0 and = 10 Fm −1 . The tunneling rate is much smaller than that for the case without coupling the Poisson equation when the doping density gets larger than a certain value (about 0.005). The right plot of Fig. 3 also shows the negative correlation between the tunneling rate and the dielectric constant , but the rate is always larger than that for the case without coupling the Poisson equation when fixed V 0 = 0, N d = 0. We further compare the Wigner functions at instants t = 8, 13, 18 fs for only the Wigner equation (left) with those for the WP system (right) in Fig. 4 when setting V 0 = 0, N d = 0, = 4 Fm −1 . We are able to clearly see there that it is much easier for GWP to pass through the barrier when the Poisson equation accounting for the space charge effects is coupled.  It is clearly shown that the space charge effects helps GWP with its tunneling through the barrier.
Now we will calibrate the proposed solver in 4-D phase space still with the Gaussian barrier scattering. We choose the Gaussian barrier as and the initial GWP as where x 0 , z 0 are the center, k 0 x/z is the initial wavenumber and σ x/z is the minimum position spread. We set the parameters to be Ω x = [−20 nm, 20 nm]×[−20 nm, 20 nm], , and x 0 = −6 nm, z 0 = 6 nm. And we still choose = 1 eV · fs, m = 1 eV · fs 2 · nm −2 , = 1 Fm −1 , N d (x) = 0 and V g = 0.5 V. The numerical results are displayed in Fig. 5, where the left (resp. right) plot shows clearly the spectral convergence with respect to N (resp. M ) while fixing M = 200 (resp. N = 100) and ∆t = 0.01 fs. Moreover, we show the reduced Wigner functions of the WP system in 4-D phase space in Fig. 6. It clearly shows that GWP crosses the barrier even when its average kinetic energy (0.72 eV) is lower than the barrier height (1.3 eV) and the Wigner functions obviously have negative values.

Resonant tunneling diode
As a classical 1-D hetero-structure device with negative differential resistance, RTD exploits resonant tunneling through double barriers as its basic mechanism. Fig. 7 gives a typical type of RTD in which two thin layers (gray) are sandwiched by another three layers (white) to form two energy barriers and one quantum well [8]. In this work, we use constant effective mass m = 0.067m 0 with m 0 being the electron mass in vacuum and set the length of the device to 40 nm which means the computational domain in x-space is Ω x = [0 nm, 40 nm]. The barrier region is set to 3nm, the length of the quantum well is 4 nm and the length of the contact is 10 nm. The doping profile in both contacts is depicted as the Fig. 7, where the n-parts are doped with a concentration 4.446 × 10 17 cm −3 and the i-part is doped intrinsically.  The initial and boundary conditions are both taken to be fixed, and given by the equilibrium Fermi-Dirac distribution: where T is the temperature, k B is the Boltzmann constant and µ L , µ R are the Fermi levels at the left and right contacts, respectively. The parameters are set as: Ω k = [− 5 3 π nm −1 , 5 3 π nm −1 ], m 0 = 9.10956 × 10 −31 kg, = 1.0546 × 10 −34 J · s , ε = 13.1ε 0 , ε 0 = 8.85 × 10 −12 Fm −1 , k B T = 2.5852 × 10 −2 eV, q e = 1.602 × 10 −19 C and µ L = µ R = 0.01 eV. In order to get rid of possible Gibbs oscillation, a cubic interpolation (smoothing) is used over a unit near the discontinuities in the band potential V b (x) and the doping density N D (x).
We are mostly interested in the formation of steady states of RTD, which correspond formally to the limit as t → +∞. Once the steady state is attained, the current of RTD should not appreciably vary with time any longer. To this end, we regard the numerical solution to be the steady state only when the difference in L ∞ -norm of the electron density given in Eq. (5) between two successive time steps is less than 10 −5 . Here, the time evolution is performed with a step of 0.02 fs up to the final time t f = 500 fs, at which the Wigner function has reached a steady state.  circles). That is, high resolution plays a key role in producing an accurate I-V curve, which constitutes the main reason for us to develop an efficient WP solver with high accuracy. Second, the I-V curves in Fig. 8 show that an incoming distribution of electrons given in Eqs. (26) and (27) can still generate a current flowing through the device even though under low bias like V 0 = 0.1 V. Simultaneously, from the converged I-V curve with M = 300 (see the blue curve with asterisks in Fig. 8), we are able to observe there that, the current under V 0 = 0 is zero, but it reaches a peak under V 0 = 0.2 V and has a valley around V 0 = 0.4 V. At even higher bias potentials, electrons surmounting the double barriers again increase the current. However, the left plot of Fig. 9 shows that the height of barrier decreases with the increase of the bias V 0 . The current peak can be reached around V 0 = 0.2 V because the resonant level in the double well is aligned with the energy of the injected electrons, which is manifested by the central peak of the electron density inside the well (see the red curve in the right plot of Fig. 9).
Finally, in Fig. 10, we show the steady Wigner functions at the final time, and find out that there are few electrons cross the barrier when the bias is 0, which explains that the current is almost zero under V 0 = 0. Meanwhile, it is obvious that electrons pass though the double barrier and partially reside inside the well under V 0 = 0.2 V, thereby verifying the resonance again.
6 Double-gate MOSFET Fig. 1 cartoons typical structure of dgMOSFET [15,16]. The width of the device is assumed to be large, and the potential is invariant along y-direction. The silicon layer is sandwiched by two symmetric oxide layers. Source and drain are doped heavily. In this work, the size parameters of the dgMOSFET device are set as follows: the gate length L G , equivalent gate oxide thickness EOT and silicon channel thickness T si are 6 nm, 1 nm and 3 nm, respectively. The highly-doped source and drain access regions are 16  10 −12 Fm −1 , the temperature T = 300 K, the doping density N d = 5 × 10 19 cm −3 in the highly-doped regions.
The electrons in the real source/drain contacts are in equilibrium characterized by a Fermi level E F 1 /E F 2 : where E F 1 = E F 2 = 0.0307 eV and the total energy ε(x) of the electron is And, the initial distribution function in z-direction vanishes in the two oxide layers and stays constant in the semiconductor layer (see Fig. 11).  We still investigate the I-V curves but now V contains the source/drain bias potential V ds and gate voltage V g . Fig. 12 shows the current J of steady states against V ds and V g , where the y-coordinate is the integral of the current density over the contact area. When V g = 0 V is fixed, the left plot of Fig. 12 gives the I-V curves against V ds . We find that the current increases with V ds . Similarly, we plot the I-V curves against V g in the right plot of Fig. 12 for fixed V ds = 0.5 V and observe there that the current also increases as V g increases. The trends of these two I-V curves are consistent with the results in [15]. We further plot the potential along the channel in Fig. 13 to explain the characteristics of the I-V curves. The left plot of Fig. 13 shows that a larger source/drain bias V 0 help electrons to pass through, so the current is higher. The right plot of Fig. 13 displays that the increasing gate voltages V g make the the potential well shallow and eventually disappear, so the more easily the electrons pass through, the higher the current becomes.  and electric potential V e (x, z) (second row) under V ds = 0, V g = 0 (first column), V ds = 0.5 V, V g = 0 V (second column) and V ds = 0.5 V, V g = 0.5 V (third column).
In order to give more details on the steady states, we also plot in Fig. 14 the reduced Wigner function f (x, k x ) (first row) and electric potential V e (x, z) (second row) under different V ds and V g . The first row of Fig. 14 shows that electrons flow more easily from the left to the right as the bias V ds increases. And at the same time, increasing the gate voltage V g is also beneficial to the flow of electrons. The second row of Fig. 14 displays that the doping forms a barrier when V g = 0, so that the intermediate channel forms a potential well, but the well depth decreases as V g increases.

Conclusion
In this paper, we made the first attempt to solve the Winger-Poisson system in 4-D phase space with high accuracy, and succeeded to develop steady states and to obtain numerically converged I-V curves from reliable long-time simulations. We believe that the proposed high-resolution solver may provide more reference solutions to benchmark the stochastic algorithms which have recently attracted a lot of attention due to its simplicity as well as its satisfactory scaling on parallel high-performance machines [17][18][19][20].