Formulation of a phase space exponential operator for the Wigner transport equation accounting for the spatial variation of the effective mass

A novel numerical approximation technique for the Wigner transport equation including the spatial variation of the effective mass based on the formulation of an exponential operator within the phase space is derived. In addition, a different perspective for the discretization of the phase space is provided, which finally allows flexible discretization patterns. The formalism is presented by means of a simply structured resonant tunneling diode in the stationary and transient regime utilizing a conduction band Hamilton operator. In order to account for quantum effects within heterostructure devices adequately, the corresponding spatial variation of the effective mass is considered explicitly, which is mostly disregarded in conventional methods. The results are validated by a comparison with the results obtained from the nonequilibrium Green’s function approach within the stationary regime assuming the flatband case. Additionally, the proposed approach is utilized to perform a transient analysis of the resonant tunneling diode including the self-consistent Hartree–Fock potential.


Introduction
For the state-of-the-art design of nanoelectronic devices including quantum effects a holistic view including the stationary and transient case is needed. Investigating the literature, there are mainly two approaches being predominantly used for the numerical investigation of quantum transport, namely the nonequilibrium Green's function approach [1] and the Wigner formalism [2].
From the engineering point of view, the Wigner formalism is preferable [3]. On the one hand, the implementation of scattering mechanisms is straightforward due to the phase space representation [4]. On the other hand, the numerical analysis of transient effects with regard to high frequency characteristics is computationally less demanding.
In contrast to the Wigner formalism, for each time-dependent injected mode at a single energy a spatial, discretized grid has to be stored over the time within the NEGF framework [5,6].
For the numerical solution of the equation of motion within the Wigner formalism, namely the Wigner transport equation (WTE), a wide variety of deterministic numerical methods has been established beside the quantity of stochastic Monte-Carlo methods, e.g., [7].
When quantum electronic devices are analyzed with regard to their current-voltage characteristics, the conventional finite-difference methods are combined with the inflow boundary conditions [8] to solve the WTE. These finitedifference methods are either based on upwind difference schemes (UDS) [8][9][10][11][12][13] or hybrid difference schemes [14] (HDS), which combine upwind fluxes and central fluxes. Unfortunately, these methods have several drawbacks as discussed in the following.
Along with the UDS applied onto the diffusion operator, the forward and backward traveling waves are separated, e.g., [8,9], following the idea of characteristic lines as if no potential is present. To conclude, from the physical point of view, the approach suffers from inherent errors, when a potential distribution is considered, as mathematically derived in [15]. From the numerical point of view, the commonly utilized UDS can be interpreted as an operator splitting scheme regarding the diffusion operator, which tends to show poor convergence properties. The application can have far reaching consequences, such as nonphysical results [16,17].
Beyond, the quantum-statistical distribution function only interacts with the diffusion operator within the UDS framework, whereas the non-local kinetic drift operator should be accounted for in the same manner.
To overcome these limitations, an approach based on the formulation of an exponential operator (EO), namely the phase space operator has been proposed to solve the WTE [15].
Furthermore, the definition of inflow boundary conditions is not sufficient on a finite computational domain. Rather, open boundary conditions have to be included to avoid nonphysical reflections due to the finiteness of the computational domain [16,[18][19][20]. This aspect can be accounted for by the complex absorbing potential formalism [16,20]. The combination of the EO along with the complex absorbing potential results in an excellent agreement with the results obtained from the quantum transmitting boundary method (QTBM) as shown in [16].
Another point of criticism is the poor convergence [21] with regard to the convolution integrals and the kernels therein. This property is strongly related to the coupling of the discretized phase space variables, which are linked in an anti-proportional manner [9][10][11][12]. To overcome these limitations, an advanced finite-volume scheme [21] and sophisticated spectral element methods have been introduced in [22,23].
Nonetheless, all these approaches neglect the spatial variation of the effective mass or are based on UDS [11,12], which are, therefore, inherently combined with errors as discussed above. But, when heterostructure devices are investigated, the adequate inclusion of the effective mass is a prerequisite. Recently, an approach based on the formulation of a slowly varying envelope function for the inclusion of the spatial variation of the effective mass has been proposed [24]. Based on the assumption that the second-order derivative is negligible, an envelope function can be formulated. Unfortunately, under nonequilibrium conditions rapid changes of the Wigner function are expected, which the envelope function does not include.
In this contribution, the concept of the EO is generalized, which is originally utilized to solve the single band WTE for the case of a constant effective mass distribution [15,16,20]. Along with this generalization, second-order derivatives regarding the transport direction can be addressed. Here, the proper numerical treatment of second-order derivatives within the Wigner formalism is discussed. Indeed, when multi-band Hamiltonians are considered, a numerical approximation of second-order derivatives is mandatory [25]. The formalism is utilized to include the spatially varying effective mass, which represents by far the simplest system, in which the Wigner formalism considers a second-order derivative. Nonetheless, the formalism provides a sound basis for future investigations, since only submatrices have to be interchanged in the case of multi-band systems. Beyond this objective, a different perspective onto the discretization of the bounded phase space with regard to the continuity equation is provided leading to more flexible discretization patterns, as well.
For clarity, the fundamentals of the Wigner formalism for quantum well devices including the spatial variation of the effective mass distribution are briefly summarized in Sect. 2. In this section we also provide a different perspective onto the discretization of the bounded phase space. Then, the phase space exponential operator is derived in Sect. 3, on which the new discretization scheme for the WTE relies. The formalism is applied onto the test case of a resonant tunneling diode in Sect. 4. Furthermore, a comparison and discussion of the results with regard to the NEGF results is performed. In addition, the proposed approach is analyzed with regard to the transient regime. The contribution ends with a summary and a conclusion in Sect. 5.

Quantum transport in phase space
For the sake of clarity, the fundamentals of the WTE including the spatial variation of the effective mass are briefly summarized in the following.
To begin with, the Wigner function defined within the phase space with the position and the wave vector k represents the Wigner-Weyl transform of the density operator ̂ according to Due to the relation of the wave vector k with the momentum p according to k = ℏ −1 p , the wave vector is briefly referred to as (normalized) momentum in this contribution, too. For the derivation of the corresponding equation of motion of the Wigner function, the Wigner-Weyl transform is applied onto the Liouville von-Neumann equation. Assuming a spatially constant effective mass, the WTE reads which seems to be the standard equation when the Wigner function formalism is utilized [3,4,[8][9][10]. The integral ( kernels included in (2) will be discussed in the following after taking into account the spatial variation of the effective mass.
For the adequate investigation of heterostructure devices, the corresponding spatial variation of the effective mass has to be included [11,12,24]. In contrast to the standard case (2), the extension of the WTE regarding the spatial variation of the effective mass results in a far more complicated expression Before going into details of the definitions of integral kernels, it is noteworthy that this governing equation conceptually differs from its constant effective mass counterpart due to the second-order derivative with regard to .
In (3) the integral kernels within the convolution integrals accounting for the spatial variation of the effective mass are introduced, where m represents the effective mass. The integral kernel including the potential distribution is given by where V contains both, the conduction band potential as well as the self-consistent Hartree-Fock potential.
The additional potential term W( ) represents the complex absorbing potential [16], which is utilized to account for open boundary conditions with respect to the -direction. As a result, nonphysical reflections caused by the finiteness of the computational domain can be effectively suppressed as outlined in the following. The main idea behind W( ) is to enclose the computational domain by an artificial absorber region, which is non-reflective and attenuates the statistical density matrix towards a value close to zero near the boundaries of the computational domain leading to a square-integrable density matrix. For a detailed discussion of the complex absorbing potential formalism within the WTE we refer to [19].

Computational Domain and its relation to the continuity equation
For the discretization of the phase space an alternative ansatz is developed, which inherently decouples the and k-directions related to the limitations of the conventional UDS based approaches, e.g., [9][10][11]. The difference to the conventional strategy is mainly based on two facts. At first, the discretized values of and do not necessarily have to share the same discrete locations. Of course, the conventional case is included within the proposed approach, too. Secondly, the conservation of mass with regard to the continuity equation is realized by including = 0 within the computational domain. This constraint results in an odd discretization number with regard to the -direction as derived in the following.
For the development of the approach, the computational domain is bounded within the real space. Hereby, we introduce symmetrical intervals with regard to the -direction according to I = −L ∕2, +L ∕2 and the -direction according to I = −L ∕2, +L ∕2 , respectively.
Within the latter interval I , the statistical density matrix can be expanded into a plane wave basis where the expansion coefficients can be easily identified with the discretized Wigner function. As a consequence, a modified formulation of the Wigner function on the bounded computational domain is obtained To ensure a unique solution, the plane wave basis functions have to be orthogonal, i.e., holds, the values of k n are required to be k n = 2 L n + 0 .

3
The phase constant 0 is introduced to indicate the symmetrical distribution of the values of k around zero, as usual.
Hereby, singularity issues related to the bounded states k = 0 can be effectively avoided. From the insertion of (7) into (6) it can be concluded that Along with these versatile relations we are equipped to tackle the conservation of mass with regard to the continuity equation.
To derive the continuity equation, the WTE (3) has to be integrated over the momentum variable k∕2 . Exemplarily, the term including the potential distribution is investigated on the computational domain resulting in Noting that L = 2 ∕ k , the relation (10) can be exploited such that holds. Along with this identity (12), the expression (11) can be reformulated in a more meaningful result since B( , 0) = 0 holds along with the complex absorbing potential being W(0) = 0 as defined. To ensure that this identity and, therefore, the continuity equation hold at any level of approximation, the location = 0 has to be included within the computational domain. In contrast to the conventional approach, this result provides a new perspective onto the interrelationship of the variables.
The integration of the remaining terms within the WTE is carried out in the same manner. Therefore, these calculations are not carried out explicitly here. Nonetheless, the general procedure is discussed in the following.
With regard to the terms including the integral kernels m − ( , k − k � ) the integrals are zero valued due to the corresponding functions in real space being odd with respect to . As a consequence, these are zero-valued at the location = 0 and the terms do not contribute to the continuity equation. Utilizing (12), the remaining terms within the WTE including the kernels m + ( , k − k � ) result in exploiting the product rule on the right-hand side. These terms can be easily assigned to the carrier density as well as the current density per carrier which together describe the conservation of the zeroth moment of the Wigner function, namely the carrier density, which is also known as the conservation of mass.
To conclude, the Wigner formalism has been investigated on a bounded computational domain, which inherently leads to a discretization of the momentum variable k. Furthermore, it has been derived that = 0 has to be included in I so that the continuity equation is fulfilled.
The question left open is how to discretize the intervals I and I . Here, we propose an equidistant finite volume discretization. Since I is defined as a symmetrical interval, the number of discretization points N have to be odd to include = 0 . The corresponding discretization width is represented by . For the -direction, there is no such a constraint and the application of N discretization points results in N s = N − 1 slices as depicted in Fig. 1, which are Fig. 1 Schematic representation of the discretized computational domain within the phase space with the parameters defined in the text separated by an amount of . With regard to the discretized momentum k n as described by (9), there is in general no restriction. In practice, the number of discretization points N k is evenly distributed around k = 0 without containing the latter one due to singularity issues Furthermore, for numerical purposes N k ≤ N holds.
Investigating the discretized values of the momentum variable, it can be concluded that k = 2 ∕L = 2 ∕N holds. As aforementioned, the conventional case assuming = 2 is inherently included leading to the relationship k = ∕N . Nonetheless, the latter expression is not preferable since it introduces an artificial coupling between the momentum k and the real space variable .
It is noteworthy, that the Wigner function decays with increasing values of k depending on the actual dispersion relation E(k) of the material, the chemical potential as well as the externally applied forces. To conclude, there is an upper limit k max obtained from (17), in which the Wigner function is significant. From this upper limit, the corresponding real space length can be determined resulting in a flexible discretization pattern As a consequence, a larger number of discretization points lead to a denser mesh, without influencing the remaining parameters such as as well as k max , which is the case for the conventional approaches.
In addition, an adaptive algorithm can be established utilizing the dispersion relation E(k) and the externally applied forces to determine k max . For instance, when larger voltages are applied, a larger interval has to be accounted for the determination of the Wigner function.

Novel discretization Scheme for the Wigner transport equation
In this section, the novel discretization scheme for the WTE including a spatially varying effective mass is established. But, the conceptional significance is much more far-reaching since multi-band Hamilton operators such as kp-Hamiltonians exhibit the same formal structure. As a consequence, they can be handled on an equal footing, nonetheless, the investigation of multi-band kp-Hamiltonians is beyond the scope of this work.
To start with, the momentum variable k is discretized, whereas the remaining variables are kept continuous for the moment. Because the structure of each term on the right-hand side within the WTE (3) consists of a weighted convolution integral including an integral kernel as well as the Wigner function, the discretization procedure can be generalized. For this purpose, a prototype function U with the corresponding Fourier transform Ũ is introduced, which is conceptually reflecting each term within the WTE according to with the exponents and taking into account integer values and the integral kernel defined by representing the Fourier transform of the prototype function as aforementioned. The prototype function equally typifies either the effective mass functions m ± (4) or the complex drift operator (5). Along with this prototype function, all convolution integrals as apparent within the WTE can be conceptually discretized utilizing the midpoint rule, which is O( k 2 ) accurate. Carrying out the summation with regard to j ′ for each index j and introducing the vector containing the discretized values f ( , k j , t) , the expression (21) can be replaced by a matrixvector representation As a consequence, the vector contains N k elements and the matrices show a dimension of N k × N k . For the determination of the discretized integral kernel Ũ , the Fourier transform has to be implemented on the discretized computational domain. The application of the midpoint rule results in with the truncation error given by O( 2 ) . Along with the different perspective onto the discretization of the phase space, the discretized values l and later do not necessarily share the same grid points as is usually adopted, e.g., [9][10][11]. As aforementioned, this conventional case is inherently included in the proposed approach by assuming = 2 . As long as an analytical expression for the potential V is available, the corresponding values of U( , ) can be determined exactly in any case. Unfortunately, when the self-consistent Hartree-Fock potential is considered, this discretized potential V( i ) is only available at certain points i within the computational domain due to the coupling of the discretized WTE along with the discretized Poisson equation. As a consequence, the discretized values of the potential term U( i , l ) can be located between two different values V( n ) and V( n+1 ) , from which the value U( i , l ) is interpolated linearly. The extension to higher-order approximations is a topic for future investigations. Now, the analytical WTE (3), which represents a partial integro-differential equation, is transformed into a system of coupled partial differential equations, which is only dependent on t and by discretizing the momentum variable k and applying the generalized discretization procedure (21)- (23). In addition, the matrices [D 2 ( )] , [D 1 ( )] , and [D 0 ( )] are introduced with the corresponding matrix elements defined by which summarize the terms with regard to the different orders of derivatives of the Wigner function. Accordingly, the WTE can be rewritten into a system of coupled partial differential equations The phase space exponential operator is derived in the next section, which represents an analytical solution with regard to the spatial differential operator on the right-hand side of (25).

Derivation of the phase space exponential operator
For the derivation of the phase space exponential operator we exemplarily analyze the behavior of the i-th slice within the computational domain considering the stationary regime.
In doing so, appropriate approximants for the real-space derivatives with regard to the -direction can be derived, which is a prerequisite in order to account for an adequate inclusion of the coherent effects [15,19,24].  (26), (27) and utilizing (24), an auxiliary equation can be defined according to which represents a simple system of coupled differential equations of the first order. For this type of differential equation the solution is well known. Concerning the i-th slice within the interval ∈ i , i + of the computational domain, the formal solution of this differential equation (28) is given by Along with the assumption of the super-vector ( ) being continuous at the interfaces between the adjacent slices, the super-vector located at i + can be easily replaced by i+1 . The occurring integral within the exponential operator in (29) can be approximated utilizing standard numerical techniques as there are for instance the mid-point integration rule as well as the trapezoidal rule. For simplicity, the integral is approximated utilizing the mid-point rule. As a result, the expression (29) can be rewritten as a propagation algorithm between adjacent locations i and i+1 The exponential operator is defined as the phase space operator. Unfortunately, the analytical evaluation of this phase space operator is combined with high computational costs with regard to the memory requirement as well as the computational time. Thus, an adequate series expansion technique is desirable to resolve the evaluation of the phase space operator utilizing for instance Chebyshev polynomials or Taylor series expansion techniques. Here, for demonstration purposes, the phase space operator is approximated utilizing Padé approximants, which results in The approximation of the phase space operator has been carried out to the linear term in both, the nominator as well as the denominator. Additionally, the abbreviations The expression (31) can be rewritten in terms of the Wigner function and its corresponding derivatives, such that and holds, where are introduced for the sake of clarity. After a straightforward calculation, (32b) can be rewritten as exploiting the interrelationship (32a) on the right-hand side that takes into account the sum of the derivatives with regard to the Wigner function.
The left-hand side of (33) considers the difference of the derivatives of the Wigner function. For the determination of the left-hand side of expression (33) in terms of the discretized points of the Wigner function i a lengthy but conceptional straightforward calculation is carried out as presented in the following.
To begin with, the nearest neighbored cells, namely the i + 1-th slice and the i − 1-th slice have to be taken into account. By utilizing the interrelationship (32a), it can be easily concluded that holds. Taking the difference of these interrelationships, (34a) and (34b), results in Unfortunately, this expression (35) is not closed with regard to i+1 and i . As a consequence, the remaining derivatives i+2 and i−1 have to be re-expressed in terms of i+1 and i as well as the corresponding discretized values of the Wigner function. In order to determine these relationships, the phase space operator (30) has to be applied twice onto the super-vector so that the next nearest slices are connected with each other.
For the approximation of the resulting operator, the same Padé approximant is applied. When connecting the i-th and the i + 2-th slices by means of the phase space operator, the following relationship can be found so that i+2 can be easily found. In the same manner, when connecting the i − 1-th and the i + 1-th slices by the use of the phase space operator (30), the following relationship holds which leads to an expression for the term i−1 . Now, a closed expression for the left-hand side of (33) can be established by inserting (36) and (37) in (35) according to Following this concept, higher-order approximants can be easily found applying higher-order Padè approximants onto the phase space operator. Now, the approximant (38) is inserted into (33) and the whole term is multiplied by ] . Hereby, numerical singularity issues can be avoided for the case of a constant effective mass, since [D 2 ( i+ 1 2 )] = 0 is zero-valued function in this case as apparent from (24). Furthermore, when the spatial dependence of the effective mass distribution is neglected, the result is in coincidence with the conventional approximation of the phase space operator [16].
After taking the transient portion of the WTE into account, the discretized WTE finally results in To summarize, by comparing (25) to this results, the proper numerical approximation of the corresponding differential operators can be determined resulting in The same discretization scheme could have been obtained from a standard finite volume scheme in the same manner. But, along with the proposed scheme higher-order approximants can be derived straightforwardly. It is noteworthy that the conventional UDS approximation [11] cannot be derived from the quasi-analytical solution of the phase space exponential operator.
For the numerical approximation of the derivatives concerning the effective mass terms within the matrices [D 1 ( i+ 1 2 )] and [D 0 ( i+ 1 2 )] (24), the adjacent interface values are utilized to approximate the value of the vertex. Considering the element m − ( , K jj � = k j − k � j ) exemplarily, which represents a matrix entry included within the matrix [D 1 ] , the expression can be rewritten according to When the derivative is not related to the values of the adjacent interfaces, numerical artifacts arise. These numerical artifacts result in large spikes of the spatial current density located at the heterojunctions.
In addition, the proposed formalism reduces to the conventional discretization as presented in [16], when a constant effective mass is assumed, which had been found in an excellent agreement with the NEGF approach. Therefore, these characteristics are inherently included within the proposed approach.
Before going into the implementation details of the system matrix, we first discuss the boundary conditions. The boundary conditions can then be partially incorporated into the system matrix leading to a non-underdetermined system matrix.

Formulation of boundary conditions
For the formulation of boundary conditions the so-called inflow boundary condition scheme is utilized [8,9]. Following the concept of the inflow boundary condition scheme the values of the Wigner function at the boundaries of the computational domain are given by a quantum statistical distribution function.
As usual, the quantum statistical distribution function at the boundaries is assumed to be given by the Fermi-Dirac statistics [8] where k B is the Boltzmann constant, T is the ambient temperature, E(k) represents the parabolic energy dispersion, and being the Fermi level.
The inflow of particles at the location 1 is assigned to the positive values of the momentum k and the inflow of particles at the location N is assigned to the negative values of the momentum k, respectively [8]. Along with these definitions, the boundary conditions for the WTE read and Since the Fermi-Dirac statistics represent a semi-classical boundary condition, these are only valid infinitely far away from the quantum structure [26]. As a consequence, the application can lead to errors tending to pose problems [17,21,[26][27][28][29]. To address this aspect, device adaptive boundary conditions based on the solution of the Schrödinger equation have been presented in [29]. Nonetheless, an alternative method, which is self-contained within the Wigner framework would be preferable. Hence, alternative formulations of the boundary conditions remain a topic for the future, yet.
In practical situations, the limitations of the semi-classical inflow boundary conditions are avoided by adding sufficiently large contact lengths [19,21]. A material dependent study for simple potential structures with regard to the validity of the semi-classical inflow boundary condition is presented in [28]. For III-V semiconductors the lower limit of the contact length is given by 60 nm [28].
To address the character of the second-order differential with respect to , the Wigner function is assumed to be a constant function with regard to the -direction at the boundaries [12]. On this basis, it can be concluded that the derivative of the Wigner function is given by and This assumption is justified for the present case since a sufficiently large contact length is presumed as discussed above. In the same manner, when the contact length shrinks, for instance below the 60nm for III-V systems, these conditions (44a) and (44b) at the boundary would have been to be reformulated.

Assembling of the system matrix
For the sake of clarity, the discretized WTE (39) is rewritten summarizing the matrices with regard to the corresponding location of the discretized Wigner function according to As apparent from (45), when addressing the first slice i = 1 and the last slice i = N s , the Wigner function 0 and N +1 are taken into account, which are located outside of the computational domain. Exploiting the boundary conditions, (44a) and (44b), it can be easily concluded that 0 = 1 holds for the first slice i = 1 as well as N = N +1 holds for the last slice i = N s , respectively. Along with these points, the righthand side of expression (45) can be rewritten with regard to the first slice i = 1 resulting in introducing the abbreviation [̃1] according to 2 . (47) In the same manner, the last slice i = N s can be analyzed introducing the abbreviation [̃N s ] which is obtained by exploiting the fact, that N = N +1 holds.
As discussed in Sect. 2.1, the computational domain takes into account N s slices leading to N = N s + 1 discretization points located on the interfaces. Thus, when connecting all slices, a number of N s (×N k ) sub-equations is obtained with N (×N k ) unknown discretized interface values of the Wigner function i . Along with the incorporation of the conventional inflow boundary conditions, (43a) and (43b), the number of equations as well as unknown discretized values of the Wigner function are equal. For the inclusion of inflow boundary conditions within the system matrix, the abbreviation [1] being a N k ∕2 × N k ∕2 matrix is introduced, such that the system matrix results in Hereby, the matrix [ ] is introduced considering the spatial behavior of the Wigner function, which is now summarized in the vector (t) . With regard to the stationary solution of the WTE, the memory requirements can be estimated by means of this system matrix [ ] . All submatrices within this matrix have a dimension of N k × N k and are full-matrices. Fortunately, all matrices must be real-valued, so that no imaginary part has to be stored. An amount of N − 3 interfaces represent the interior nodes, which do not interact with the boundary ( i = 2 to N s − 1 ). Along with the four-point stencil due to the second-order derivative, a number of4 ⋅ (N − 3) ⋅ N 2 k elements within system matrix are non-zero valued. The von-Neumann boundary conditions have been incorporated in the system matrix at the outer interfaces, which, therefore, add only three full submatrices. Finally, to account for the inflow boundary conditions, 2 . (49) identity matrices of dimension N k ∕2 × N k ∕2 contribute N k ∕2 non-zero values on each site. To conclude, the worst case memory requirement regarding the non-zero values of the system matrix is given by assuming a 64-bit precision. For instance, a 500 × 500 grid results in approximately 4GByte non-zero values within the system matrix. When large contacts with constant characteristics are assumed, the matrices [ i ] are locally zero-valued, due to the function m − ( , k) being zero-valued, reducing the memory requirement.
The system matrix for the transient portion of the discretized Wigner function is carried out in the same manner leading to introducing the matrix [ ] . In contrast to the conventional approach, a transient reaction of the boundary conditions is inherently included within this discretization scheme as can be observed from (51). This property immediately reflects the non-local character of quantum mechanics. In this contribution, we assume that the quantum statistical distribution function at the boundaries is time-independent Nonetheless, the consideration of a time-dependent reaction of the system represents an interesting topic for future investigations. Finally, the discretized WTE can be rewritten with the abbreviations being apparent from the context (50) and (51).
For the determination of the Wigner function, the equation system (53) has to be solved either for the stationary case or the transient case taking into account the corresponding boundary conditions. With regard to the transient operator, standard methods such as explicit Euler or Runge-Kutta schemes can be utilized.

Evaluation of the proposed approach
For the assessment of the proposed approach a simple structured double barrier resonant tunneling diode based on the material system GaAs/Al x Ga 1−x As is investigated. The spatial variation of the effective mass is considered. Two different alloy contents according to x = 0.2 and x = 0.3 are investigated within the proposed approach. A larger alloy content does not only lead to an increase of the effective mass, it also leads to an increase of the values of potential within the barriers. This aspect affects the discretization applied onto the computational domain and is discussed here.
For the validation of the proposed approach, the results are compared to the results obtained from the NEGF approach [1] in the stationary and ballistic regime. Furthermore, the flatband case is considered because the coupling to the Poisson equation can result in additional numerical errors. In addition, the proposed approach is extended in a straightforward manner to the analysis of the transient case in contrast to the NEGF case. Therefore, this analysis is provided as well, exemplarily.

Resonant tunneling device
For the design of the resonant tunneling diode, we have basically utilized the length parameters according to [11], in which the spatial effective mass variation has been included for the first time being with regard to the numerical simulation of resonant tunneling diodes along with the WTE.
The structure of the resonant tunneling device as shown in Fig. 2 represents the basis for the numerical simulations. The length of the GaAs quantum well L w is given by 4.5nm , which is surrounded by L b = 2.8nm Al x Ga 1−x As barriers. On each site, the electrodes are in total L c + L d , with L c being a structure dependent quantity and L d = 17.5nm. In contrast to the structure in [11], longer contacts are introduced. As aforementioned, the semi-classical boundary conditions are only valid infinitely far away from the heterostructure device [19,21,26]. Beyond, the voltage is assumed be to be constant within the domains addressed by L c , whereas the voltage drop takes place within the remaining domain, which is addressed by L d , L b and L w .
The contact regions are highly doped by an amount of N d = 2 ⋅ 10 +18 cm −3 . The temperature is set to be T = 300K . Since the conduction band is analyzed, the potential energy at the -point of the material GaAs is utilized as the reference energy and set to zero. Consequently, the potential within the AlGaAs barriers is given by approximately 0.18eV for the alloy content x = 0.2 and 0.28eV for x = 0.3 , respectively. The effective masses are m GaAs = 0.063 ⋅ m 0 in the material system GaAs as well as m AlGaAs = 0.0796 ⋅ m 0 for x = 0.2 and m AlGaAs = 0.088 ⋅ m 0 for x = 0.3 in the two investigated AlGaAs barrier materials. The constant m 0 represents the rest mass of an electron. For all simulations, the length of the layer including the complex absorbing potential [19] with respect to the -direction has been set to 12.6nm. The values of the complex absorbing potential are defined to increase quadratically towards the boundaries taking a value of 1eV.

Investigation of the stationary regime
For the investigation of the resonant tunneling diode within the stationary case, the dynamic portion of (53) is omitted, such that the system of equations given by has to be solved with containing the boundary conditions given by the Fermi-Dirac statistics (43a) and (43b), respectively.
To start with, the numerical simulation results for both resonant tunneling diodes for the different alloy contents are investigated for a fixed discretization. Then, the discretization of the phase space is carried out utilizing different grids allowing a first study of the convergence. Hereby, the discretization width is set to be a, a/2, a/4 with a = 0.565 nm and an interval of [−1.5 nm −1 , +1.5 nm −1 ] utilizing N k = 100 up to N k = 600 discretization points are taken into account with regard to the momentum variable k.

Numerical simulation results
For the presentation of the results with regard to the alloy content x = 0.2 , a discretization width = a∕2 combined with a number of N k = 400 momentum discretization points is employed and the contact length is set to L c = 35 nm . At first, the carrier densities are investigated, which are obtained from the corresponding Wigner function utilizing (15). For the alloy content x = 0.2 , the carrier densities as well as the corresponding conduction band including different applied bias are depicted in Fig. 3. In Fig. 3a, the equilibrium case is depicted and a qualitatively good agreement with the reference solution stemming from the NEGF approach is observed. In a similar manner, the nonequilibrium results agree well, as there are an applied bias U = 0.2 V depicted  Fig. 3b and an applied bias U = 0.4 V, which is shown Fig. 3c, respectively. Nonetheless, for larger applied biases, the deviation between both methods slightly increases in regions, where due to a voltage drop an additional confinement region emerges, as can be seen from Fig. 3c.
For the determination of the current-voltage characteristic, a bias U ranging between 0 V and 0.4 V is applied in steps of U = 0.02 V. The spatial current density is calculated from the Wigner function for the corresponding bias U utilizing (16). In particular, the discretized current operator is a local function in , in contrast to the conventional methods [8][9][10][11][12], which are non-local in this variable.
Comparing the current density obtained from the reference method given by the NEGF approach along with the proposed approach, the curves correspond well with each other. The characteristic negative differential resistance is adequately reproduced within the proposed approach. Furthermore, the bias at approximately U = 0.14 V of the resonant peak coincides (Fig. 4).
Since the conservation of mass with regard to the continuity equation within the Wigner framework represents a critical quantity, which has to be adequately accounted for, the spatial dependency of the current density is shown in Fig. 5 for different applied biases U. From the results depicted in Fig. 5 it can be concluded that the conservation of mass is properly described within the proposed approach due to the current being a constant function of .
The case of the resonant tunneling diode with the alloy content according to x = 0.3 has to be considered more carefully due to the larger potential barrier. As a consequence, large contact regions have to be employed according to L c = 105nm in order to approximate the requirements of the application of the semi-classical Fermi-Dirac statistics.
For smaller contact lengths, nonphysical quantities can be obtained, i.e., negative carrier densities. In a similar manner, this aspect is observed and discussed in [21].
For the discretization of the phase space, the -direction is discretized utilizing a discretization width = a∕2 and the interval within the k-direction is discretized utilizing N k = 600 , discretization points.
The nonequilibrium carrier densities as well as the corresponding potential of the conduction band with regard to the alloy content x = 0.3 are shown in Fig. 6. The equilibrium results are not shown here, since these coincide very well with the results obtained from the reference NEGF approach, but do not provide additional insights. Considering an external bias according to U = 0.2 V, the carrier densities are depicted in Fig. 6a, whereas the carrier densities for an external bias U = 0.3 V are depicted in 6b. For the values of the contact length L c smaller than 105nm , negative carrier densities have been obtained along with the same discretization widths in and k. Indeed, this effect is observed to be prominent in the regions, where due to the applied bias quantum confinement occurs leading to an oscillating carrier density. Along with the proposed contact length, this nonphysical behavior related to the validity of the semi-classical Fermi-Dirac boundary conditions can be more or less effectively avoided. The implementation of a nonuniform grid would be preferable for this case and represents a topic for future investigations.
The current-voltage characteristic for the resonant tunneling diode considering an alloy content x = 0.3 is depicted in Fig. 7. When comparing the result obtained from the proposed approach to the results obtained from the NEGF approach, it can be concluded that the device behavior is predicted in a similar manner. The proposed approach exhibits the negative differential resistance occurring at the same resonance energy. Nonetheless, the values of the current beyond the peak resonance are quite above the reference values obtained from the NEGF, although the results seem to be converged. The investigation of this behavior and its possible relation to the contact length represents a topic for future works. The spatial dependency of the current is depicted in Fig. 8 for five different applied biases, as there are the equilibrium case as well as different nonequilibrium cases ranging from 0.1 V to 0.4 V spaced by 0.1 V. As can be observed from this figure, the steady state current is spatially constant as required by the continuity equation. To conclude, the continuity equation is preserved at a numerical level considering large discontinuities within in the effective mass distribution.

Analysis of the convergence behavior
Due to the numerical solution of two distinct approaches, namely the Schrödinger equation and the WTE, different errors arise. On the one hand, different numerical schemes are employed to discretize these equations, which lead to different numerical errors. On the other hand, different type of boundary conditions [21,26] as well as the finiteness of the computational domain [19] leading to an unavoidable finiteness of the phase space, result in different model errors. While a quantum statistical distribution function is defined within the Wigner formalism, open boundary conditions are adequately accounted for by the concept of self-energies within the NEGF framework.
As observed in the latter section, that for instance the carrier density takes nonphysical values when the contact length is too short within the Wigner formalism, cannot be observed within the NEGF formalism. To address this aspect, the convergence of the result by means of the relative errors is not only compared to the results obtained from the NEGF, rather the relative error is investigated with regard to a Wigner solution on a fine grid. Note that a detailed investigation of the convergence with respect to all the different sources of errors within the Wigner formalism is beyond the scope of this paper. Consequently, the results shown here, are not the last word, but should be more taken as the first word.
For the assessment of the convergence, the relative error is defined with regard to the carrier density nodes summarized in the vector according to where the reference solution ref is either taken from the NEGF approach or from the Wigner formalism on a fine mesh, as mentioned. For = a the convergence of the results is not guaranteed, since the resolution of the barriers seems to be to coarse. Therefore, the corresponding results are not considered. Furthermore, the decrease of the discretization width to a/4 does not provide significant further insights. Therefore, we restrict ourselves onto the analysis of the = a∕2 for the first. The relative error with regard for both resonant tunneling diodes in dependence of the number of discretization points N k is depicted in Fig. 9. Hereby, the reference solution is obtained from the WTE on the finest mesh ( N k = 600 and = a∕2 ). As can be observed from Fig. 9, the relative error decreases with increasing number of points. To conclude, the method seems to be stable under grid refinement. Nonetheless, it is noteworthy that there are certain grid configurations, which do not converge, occurring only in the nonequilibrium case. This aspect is strongly related to the Wigner function, which then represents a highly oscillating function that has to be adequately resolved, which is also denoted in [21]. Note, the relative error with regard to the current voltage characteristic considering the NEGF reference solution is in the order of approximately 6% for the alloy content according to x = 0.2 for all N k as well as about 45-30% percent, decreasing towards the finest mesh for the alloy content given by x = 0.3 . For both resonant tunneling diodes the relative error regarding the carrier densities is below 1%.
Along with the application of higher-order methods onto the discretization regarding the k-direction, as for instance the spectral element method [22,23], an improved convergence behavior is expected for the proposed approach being an interesting topic for future investigations.

Comparison to the constant effective mass case
From the solution of the Schrödinger equation, it is well known that the spatial variation of the effective mass has a deep impact on the current-voltage characteristic. Indeed, with respect to the class of intraband resonant tunneling diodes, in which the effective mass within the resonant double barrier structure usually increases, the increasing effective mass lowers the overall current density. Unfortunately, the conventional UDS-based approach [11] does not resolve this behavior. Along with the UDS-based approach, the increasing effective mass within the barriers results in an increasing overall current-density.
For that reason, the proposed approach is analyzed with regard to the constant effective mass case, too. In Fig. 10, the current densities for the RTD are depicted applying the case of a constant effective mass as well as the case of a spatially varying effective mass. For these investigations, the above defined structure with the alloy content x = 0.3 is utilized. The corresponding parameters of the computational domain are adopted.
As can be observed from Fig. 10, the overall values of the current density for the case of a constant effective mass (const m = 0.063 ⋅ m 0 ) are larger than the values of the current density for the case of a spatially varying effective mass. Indeed, the peak current density is approximately two times higher assuming a constant effective mass. Furthermore, a fairly well agreement along with the reference method obtained from the NEGF method is found for the case of a constant effective mass. In this case, the proposed approach reduces to the conventional phase space operator approach presented in [15], which has been found to be in an excellent agreement with the NEGF reference method when applying the complex absorbing potential formalism [19]. To conclude, these characteristics are inherently included within the proposed approach as well.
Nonetheless, we expect major improvements for the UDSbased approach by including the complex absorbing potential formalism, too, as has been conceptually shown in [19] for the constant effective mass case utilizing a second-order

Transient regime
For the investigation of the transient effects within the resonant tunneling diode, (53) is solved utilizing Heun's method for the approximation of the transient operator resulting in To account for the self-consistent Hartree-Fock potential, the transport equation is mutually solved along with the Poisson equation [8]. As a consequence, the potential is updated for each time step. Since the Heun method is an implicit method, a large time step with a value of t = 1 fs can be chosen. When conventional explicit methods such as Runge-Kutta schemes, explicit Euler, and multi-step methods are utilized the time step has to be several orders of magnitude smaller to guarantee a stable numerical scheme. To overcome these limitations with regard to the small timesteps of explicit schemes, exponential integrators seem to represent an attractive alternative. Nonetheless, the investigation of exponential integrators for the WTE is a topic for future analysis.
For the numerical presentation purposes, the resonant tunneling diode with the alloy content according to x = 0.2 is considered, due to the minor requirements with respect to the computational domain. Furthermore, the well region has been slightly doped with N w = 1 ⋅ 10 +17 cm −3 .
For the efficient numerical solution of the system (56) an iterative method, namely the bi-conjugate gradient stabilized method, is utilized [30]. In addition, the matrix on the lefthand side of (56) is exploited for the formulation of a preconditioner by applying the incomplete LU factorization [31] in order to effectively reduce the overall computation time.
The device is initially in the nonequilibrium case applying a bias U = 0.14V and abruptly driven towards 0.24V for t > 0 . The initial Wigner function at t = 0 is obtained by solving the stationary system of equations (54) combined with the Poisson equation applying the corresponding bias U = 0.14V , which is depicted in Fig. 11a.
After the system evolves 500fs in time, the Wigner function shows significant values at larger values of k as depicted in Fig. 11b, which arise due to the larger applied external voltage. As a consequence, due to this behavior a larger current can be expected stemming from the symmetry considerations with regard to the definition of the current within the Wigner framework (16). This behavior is also reflected in Fig. 12, where the transient evolution of the current density is depicted for different characteristic regions. These regions include the source, the drain and the quantum well region.As can be observed, the time-dependent current at the source and drain regions oscillate strongly in contrast to the current within the quantum well. Each contact, partially described by the Fermi-Dirac statistic, tries to establish an  Fig. 12 Transient evolution of the current density j(t) within different locations equilibrium situation with its surroundings. Consequently, the charge injection and depletion are prominent in these locations until a local equilibrium condition emerges.

Summary and conclusion
A novel scheme for the numerical solution of the Wigner transport equation including the spatial variation of the effective mass has been derived. The scheme is based on the Padé approximation of the quasi-analytical solution given by the phase space exponential operator. Beyond, a different perspective onto the discretization of the phase space is provided. The approach has been utilized to simulate to AlGaAs/ GaAs resonant tunneling diodes with different alloy contents. The results obtained have been compared to the results obtained from the nonequilibrium Green's function approach in the stationary regime. For the smaller alloy content fairly good agreement can be observed with regard to the current density as well as the carrier density. Regarding the larger alloy, the results coincide qualitative. The observed deviations have been discussed with regard to the boundary conditions. Furthermore, a brief investigation of the convergence is provided. From this investigation it could be concluded that the method seems to be stable under grid refinement.
Finally, the proposed approach is utilized to investigate the transient behavior considering the self-consistent potential. To put it in a nutshell, the proposed approach offers major improvement with regard to the numerical simulation of heterostructure devices along with the Wigner transport equation.