Numerical constraints and non-spatial open boundary conditions for the Wigner equation

We discuss boundary value problems for the characteristic stationary von Neumann equation (stationary sigma equation) and the stationary Wigner equation in a single spatial dimension. The two equations are related by a Fourier transform in the non-spatial coordinate. In general, a solution to the characteristic equation does not produce a corresponding Wigner solution as the Fourier transform will not exist. Solution of the stationary Wigner equation on a shifted k-grid gives unphysical results. Results showing a negative differential resistance in IV-curves of resonant tunneling diodes using Frensley’s method are a numerical artefact from using upwinding on a coarse grid. We introduce the integro-differential sigma equation which avoids distributional parts at k=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0$$\end{document} in the Wigner transform. The Wigner equation for k=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0$$\end{document} represents an algebraic constraint needed to avoid poles in the solution at k=0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k=0$$\end{document}. We impose the inverse Fourier transform of the integrability constraint in the integro-differential sigma equation. After a cutoff, we find that this gives fully homogeneous boundary conditions in the non-spatial coordinate which is overdetermined. Employing an absorbing potential layer double homogeneous boundary conditions are naturally fulfilled. Simulation results for resonant tunneling diodes from solving the constrained sigma equation in the least squares sense with an absorbing potential reproduce results from the quantum transmitting boundary with high accuracy. We discuss the zero bias case where also good agreement is found. In conclusion, we argue that properly formulated open boundary conditions have to be imposed on non-spatial boundaries in the sigma equation both in the stationary and the transient case. When solving the Wigner equation, an absorbing potential layer has to be employed.


Introduction
Technology computer-aided design (TCAD) uses simulations to aid in semiconductor product development. Quantum transport models are indispensable to appropriately model modern devices and will further gain in relevance. The Wigner function method (WFM) [3,24,25] is based on a mathematical formulation of quantum mechanics which is formally close to a classical phase space description [5,26]. This allows for flexible mixed quantum-classical models and makes it an attractive approach for many applications where only parts of the system need to be modeled fully quantum mechanically. For example, the Wigner function approach can incorporate classical boundary conditions (BCs), semiclassical scattering models or classical force terms for the smooth part of the potential.
Quantum electron transport in modern semiconductor devices is described by a Wigner equation which is formally similar to the classical Liouville equation (also called Vlasov equation) [6]. Originally, the Wigner equation was used to simulate resonant tunneling diodes (RTDs) [4,21]. Now that the active region of transistors is getting closer and closer to the de Broglie wavelength, quantum mechanical effects already affect many currently manufactured devices. Several practical applications of the Wigner function method are demonstrated in the monograph [17].

Quantum distribution functions
Many key performance indicators of semiconductor devices are extracted from measured current-voltage characteristics. Those characteristics can be simulated by the numerical solution of stationary partial differential transport equations.
Stationary quantum transport is described by the Liouville-von Neumann equation for the density function (x, y) Here, V(x) is the potential energy. We will always assume constant mass m and only deal with a single spatial dimension in this paper. From the density function (x, y) , the Wigner function f(r, k) is defined as the result of two consecutive transformations described below.

Sigma function (characteristic function)
First, we introduce new coordinates for the quantum density Using these coordinates, the density matrix transforms into the sigma function and the stationary von Neumann equation transforms into the stationary sigma equation where the potential term U(r, s) is defined by In a single spatial dimension, Eq. (3) is the characteristic hyperbolic form of the stationary von Neumann Eq. (1). The sigma function has the symmetry property Its real part is an even function, the imaginary part is an odd function of s.
Particle density and local current are given as For practical simulation, a numerical cutoff parameter L in s-space has to be introduced which restricts s to an interval [−L, L] symmetric around s = 0 . This cutoff parameter is called the coherence length and is also used in the Wigner equation.

Wigner function
The Wigner function f(r, k) is then derived from the sigma function (r, s) via a Fourier transform in coordinate s With these conventions, a wave function of the form (x) = e ik 0 x has a corresponding sigma function (r, s) = e ik 0 s and a Wigner function f (r, k) = (k − k 0 ).
Applying the Fourier transform (7) to the sigma equation (3) gives the stationary Wigner equation Here, the Wigner potential V w (r, k) is defined as the Fourier transform of U(r, s) divided by iℏ The Wigner potential is a real function odd in variable k. For non-zero bias, it has a 1/k-pole at k = 0 . Particle density and local current are given as

Spatial boundary conditions
For open systems, classical inflow boundary conditions are imposed on the stationary Wigner equation (two-point boundary value problem) Here, f L and f R are prescribed distributions in the left and right reservoir depending on temperature and the doping in the electrodes.
The existence and uniqueness of solutions to the stationary Wigner Eq. (8) with inflow boundary conditions is a long-standing open problem even in a single spatial dimension [1,14,15].

Formal classical limit
For the classical limit of the sigma equation, we rescale the difference coordinate s and introduce a new coordinate = s∕ℏ . We label the function in rescaled coordinates as ̃(r, ) = (r, � ).
With p = ℏk we have k s = p and the Fourier transform of ̃(r, ) is f (r, p) . To perform the formal limit ℏ → 0 , we make the linear approximation The stationary sigma equation becomes in this formal limit where F(r) = −V � (r) is the force field. Equation (12) is the inverse Fourier transform of the stationary Vlasov equation

Numerical solution on a shifted k-grid
Frensley's original method [4] uses a shifted k-grid exluding the point k = 0 . Using a discrete Fourier transform, it is numerically much cheaper to solve for the solution in s-space where the equation is sparse in both coordinates, while the Wigner equation is not sparse in coordinate k. This has been undertaken in [11,12] and is discussed below. We use shooting methods to reduce the two-sided problem (inflow boundary conditions) to the one-sided case.

Integral form and Goursat problem
Integrating both sides of Eq.
Equation (14) is a two-dimensional integral equation of Volterra type. Existence and uniqueness of the solution can be proved [2,7,22]. We can arrive at a solution in full s-space by solving two Goursat problems separately in upper and lower half space where solutions are matched through the shared boundary condition (r) . In general, these solutions will not be L 2 -integrable and the Fourier transform does not exist as an ordinary function.
There is a striking discrepancy between Goursat boundary conditions and the boundary conditions for the Wigner equation. In the Wigner equation, no boundary conditions corresponding to (r) are specified (however, compare the discussion in Sect. 6). The Wigner and the sigma equation are formally equivalent via Fourier transform. But existence of the Fourier transform as an ordinary function is a nontrivial condition, and for this reason, the set of solutions to both equations is not in a simple bijective correspondence. This makes the mathematical analysis subtle.

Shooting methods
The Wigner function on a shifted k-grid (2N points) is related via a discrete Fourier transform to a discrete sigma function on a non-shifted s-grid ( 2N +1 points) plus antiperiodic boundary conditions in s-space [16].
The one-sided boundary value problem with spatial boundary conditions (s) for (r, s) given on an s-line (r 0 , s) and anti-periodic boundary conditions (r, L) = − (r, −L) can be reduced to the case of Goursat type boundary value problems by splitting into even and odd part and solving two separate equation systems in half space. The odd part fulfills o (r, 0) = 0 and the even part fulfills e (r, L) = 0 . Hence, the one-sided problem with anti-periodic boundary conditions is also well-posed. In the Wigner equation, usually two-sided inflow boundary conditions are given which can be imposed in the sigma equation via the discrete Fourier transform. Existence and uniqueness of the solution to the two-sided inflow problem has been proved in [1] for the discrete Wigner equation on a shifted k-grid. The same problem formulated in s-space is equivalent via the discrete Fourier transform and also well-posed.
The two-sided problem can be reduced to the solution of the one-sided problem using shooting methods (both for the sigma and the Wigner equation). In this way, no system matrix has to be stored, the method can easily be parallelized and numerical solution for very dense meshes is practical.

Failure of Frensley's method
The shooting method for the sigma equation with non-spatial anti-periodic boundary conditions was presented in [11]. No simulation results have been included in [11], as simulation results from the sigma equation did not reproduce the negative differential resistance expected in the simulation of IVcurves for resonant tunneling diodes and an implementation error could not be excluded.
The reason for the failure was not readily understood until further work was performed which is included in [10,12]. It turned out that simulation results which showed no negative differential resistance were actually numerically correct. In the case of the sigma equation, the mesh spacing (Δr , Δs) can be repeatedly refined until the solution does no longer change significantly and a numerical limit is reached. If the discrete Wigner equation is solved without upwinding using the same meshing parameters, we immediately closely reproduce the solution from the sigma equation showing no negative differential resistance in the IV-curve.
In contrast, IV-curves calculated using Frensley's method show a wished-for negative differential resistance. However, in Frensley's method, upwinding is used which introduces a huge discretization error in variable r. In the numerical limit, the solutions with and without upwinding have to agree. When using upwinding, the limit solution is only reproduced if the mesh is further refined in coordinate r as this reduces the discretization error. This was accomplished in [12] using shooting methods. With refinement in r, the solution with upwinding slowly converges towards the solution without upwinding. A very dense r-mesh has to be used. For illustration, the example from [12] (Fig. 2 in the original) is reproduced here in Fig. 2.
We conclude that the negative differential resistance shown in simulations of IV-curves for RTDs based on Frensley's method is nothing but a numerical artefact stemming from the use of upwinding on a coarse mesh. Furthermore, we noticed strong negative concentractions near k = 0 in simulation results. It is these observations which motivated this work.

Wigner equation at k = 0
In a single spatial dimension, Eq. (8) can be rewritten for k ≠ 0 as The latter form emphasizes the fact that the Wigner equation becomes singular at k = 0 , which complicates its analysis. We can derive two equations for k = 0 which are mathematically of different character: 1. Setting k = 0 in Eq. (8) gives the following integrability constraint In this degenerate case, the left hand side in Eq. (8) vanishes and we do not get a differential equation. This special case is called an algebraic constraint in [1,15] and is needed to avoid poles on the right hand side of Eq. (17). The integrability constraint has also a physical interpretation in the particle picture: The total potential inscattering rate at k = 0 must vanish in the steady state. 2. If constraint (18) is fulfilled, then we can take the limit k → 0 in Eq. (17). This gives the "transport" equation at k = 0 Here, f k (r, k) = f (r,k) k denotes the first order derivative of f(r, k) with respect to k. Alternatively, Eq. (19) can be derived by differentiating both sides of Eq. (8) with respect to k and then setting k = 0. Derivation of the equations at k = 0 assumes some regularity of f(r, k), respectively, its partial derivatives near k = 0 , symbolically written as Two equations at p = 0 can be derived also in the classical case. By setting p = 0 in (13), we get as the classical version of the quantum integrability constraint (18). Differentiating both sides of Eq. (13) with respect to p, we get the equation as the classical version of the transport Eq. (19).
The Wigner potential is odd in k, and in the stationary Wigner equation even and odd part only couple through inflow boundary conditions. The integrability constraint (18) only concerns the odd part of the Wigner function, while the transport equation at k = 0 (19) only concerns the even part. As the (inverse) Fourier transform conserves parity, this parity analysis transfers to the sigma equation.
In this work, we use the quantum transmitting boundary method (QTBM) [13] for comparison with the Wigner function method. The Schrödinger modes occurring in the QTBM solution represent unnormalizable states. Their Wigner transform can be calculated analytically for simple potentials (e.g., potential step and barrier structures in 1D) and in general contains distributional parts. When summed up according to the distribution of incoming modes, the Here, the zero bias case (see Sect. 6) is exceptional because for each QTBM mode (k)-singularities appear in its Wigner transform. Consequently, for zero bias, the Wigner transform of the full QBTM solution may contain distributional parts (k) . For non-zero bias, the Wigner transform of each QTBM mode is found to be smooth at k = 0 fulfilling both Eq. (18) and Eq. (19) in the appropriate sense.

Constrained equation
Solution of the Wigner equation on a shifted k-grid does not produce physically meaningful results as discussed in Sect. 3. Using a shifted k-grid, we can derive a discrete form of the transport equation at k = 0 (19) by subtracting the discrete equations for Δk∕2 and −Δk∕2 . But those solutions violate the integrability constraint (18).
A natural attempt to repair this is to include k = 0 in the mesh which incorporates the integrability constraint. Using a suitable discrete inverse Fourier transform [16], this corresponds to periodic boundary conditions in s-space. Unfortunately, the sigma equation with periodic boundary conditions (and hence also the corresponding discrete    26) gives what we call the integro-differential sigma equation where sgn denotes the sign function. In s-space there is a priori no unique antiderivative, the antiderivative is undefined by an additive constant. The derivative s turns into multiplication with ik in k-space. This operation has a canonical inverse operation in k-space, namely division by ik , because we do not allow for distributional parts (k) in the Wigner function. This is tacitly assumed in Eq. (26).
By splitting the integral on the right hand side of Eq. (27) we see that this expression picks out a specific antiderivative H(s) of a function h(s), in particular the arithmetic mean of upper and lower antiderivative.
Restricting the integro-differential sigma Eq. The only reasonable choice for the integration constant is to set on a s-domain symmetric around s = 0 , i.e., fully homogeneous boundary conditions. Fully homogeneous boundary conditions have been motivated by inverse Fourier transform of the integrability constraint (18) and the transport equation at k = 0 (19) in [10]. When solving the Wigner equation in k-space on a shifted k-grid, one has to incorporate the

Zero bias case
It has to be pointed out that the zero bias case is exceptional. At zero bias, independent of the energy, the Wigner function resulting from a single Schrödinger mode (as used in the QTBM method) has a singularity at k = 0 which stems from the long range correlation between the reflected and the transmitted part of the mode.
All the Wigner functions from the stationary Schrödinger modes contain distributional parts at k = 0 which are not smoothed out when modes are summed up. If distributional parts at k = 0 are allowed in the solution, then multiple solutions may actually exist.
A related kind of critique has been raised in the papers [18] and [23]. The authors claim a "failure of conventional boundary conditions" and suggest that the inflow boundary value problem might have multiple solutions. We note that those papers fail to identify the more fundamental cause of numerical problems (proper non-spatial open boundary conditions need to be employed, see Sect. 8).
Specific examples discussed by [18] use a symmetric and periodic potential, hence zero bias. In the zero bias case, the Wigner potential V w (r, k) is not singular at k = 0 and we may allow for (k)-terms and also for 1/k-poles in the solution f(r, k). Using the ansatz in the Wigner Eq. (8), we get the Wigner equation for h(r, k) Here, we have defined Ṽ w = m � V w . In the ansatz (32) the function c(r) can be freely chosen and acts like a boundary condition on the line k = 0 (i.e., an analogue of Goursat type boundary data). The integrability constraint (18) will in general not be fulfilled by the solution if we allow for poles 1/k.
In the zero bias case, the solution to the stationary Wigner equation is non-unique. Here, the right choice of a functional space could help to filter out the relevant solution for the boundary value problem. Our numerical strategy is to assume an equilibrium distribution function without distributional parts or poles at k = 0 , and in this work, we always impose the constrained sigma equation.

Simulation results
The constrained Wigner equation has first been introduced in [12]. Solving the constrained sigma equation in a least squares sense gives reasonable but not completely satisfactory results [10]. The natural way to incorporate double homogeneous boundary conditions is by employing an absorbing potential layer as recently introduced in [20]. This has the effect of damping down the solution near the non-spatial boundary. The imaginary potential is chosen as an even function in s and we still get a real Wigner potential which can be employed when solving the Wigner equation. The absorbing potential couples even and odd parts of the equation inside the simulation domain which otherwise only couple through inflow boundary conditions. In [8], a spectral method for the constrained sigma equation with an absorbing potential layer has been presented. This amounts to a discrete solution of the constrained Wigner equation in k-space. Memory requirements are huge and a special direct linear solver had to be implemented for the solution.
Below we present simulation results from solving the constrained sigma equation (double homogeneous boundary conditions) with an absorbing potential layer using finite differences. The corresponding discrete equation is sparse in both variables and memory requirements are much lower than if solving in k-space. The aim in this work is to investigate whether the constrained equation reproduces the solution from the QTBM in the numerical limit (fine mesh).
We use an orthogrid symmetric around s = 0 including s = 0 and a stencil for the unit square.
Given on the boundary, we calculate f on the boundary via the discrete Fourier transform and then impose inflow boundary conditions exactly as in Frensley's method on a shifted k-grid. As we are solving an overdetermined system, there is a freedom to choose equations which are imposed as strict constraints and which equations are only fulfilled in a least squares sense. In our over-determined method, the inflow boundary conditions are imposed as strict constraints.
The overdetermined system is then solved in a least squares sense with constraints exactly fulfilled. Mathematically this represents a quadratic programming problem with linear equality constraints. For the solution, we (33) use the method described in [8] which avoids the explicit building of normal equations. The numerical method is evaluated by simulating a GaAs/AlGaAs RTD and comparison with the quantum transmitting boundary method which we take as the reference model. The RTD structure used in the simulation is depicted in Fig. 3, it consists of a 4.5 nm-wide GaAs quantum well and 2.8 nm-wide AlGaAs barrier layers. The barrier height is 0.27 eV.
In order to get a good fit to QTBM a coherence length of at least 120 nm and a mesh spacing Δr < 0.1 nm should be employed. Here, we consider an even larger coherence length of 170 nm and a mesh size of N r = 6000, N s = 5601 with Δs = 2Δr . The length of the electrodes used in the simulation is 85 nm (half the coherence length). The coherence length and the mesh size have been chosen in order to study the method in the limit of a large coherence length and a fine mesh.
An absorbing layer is employed at upper and lower s-boundary. We use an imaginary quadratic potential with maximum value at the boundary corresponding to a bias of 2 V. Each absorbing layer comprises 10% of half the s-domain.
Simulation results for particle density are depicted in Fig. 4. In general, the fit with QTBM is very good. For zero bias, we see a small deviation in the density at the center of the well but also in this case the fit is still good.
In order to demonstrate the numerical capability of the method, we simulate up to a high bias which is beyond the usual operating conditions of the device. In Fig. 5, we see a weak second resonance in the IV-characteristics for both QTBM and sigma function method. The fit between both methods is excellent.

Conclusion
Frensley's method [4] uses a shifted k-grid which ignores the integrability constraint (18). Results showing a negative differential resistance when simulating IV-curves of resonant tunneling diodes are an artefact of numerical diffusion (upwinding) using a too coarse mesh. The method breaks down in the limit of a fine mesh [10].
In order to avoid a singularity at k = 0 , the integrability constraint has to be incorporated when solving the stationary Wigner equation. The integro-differential sigma equation (27) avoids distributional parts in the Wigner transform at k = 0 . We impose the inverse Fourier transform of the integrability constraint in the integro-differential sigma equation. After a cutoff, double-homogeneous boundary conditions are derived.
To get results closely reproducing QTBM, an absorbing potential layer has to be employed. When using an absorbing potential, double homogeneous boundary conditions are naturally fulfilled and results from QTBM are reproduced with high accuracy. In the zero bias case, a solution to the stationary Wigner equation may contain distributional parts at k = 0 . Exemplary simulation results still show a good fit for the density n(x) between QTBM and our method.
We want to stress that it is not necessary to actually solve an overdetermined system (neither for the sigma nor for the Wigner equation). Results presented in [20] demonstrate that combining an absorbing potential with antiperiodic boundary conditions gives good results closely reproducing results from QTBM. We have also experimented with combining an absorbing potential and periodic boundary conditions which again gives good results in numerical simulations but will not be further discussed here.
The overall insight from these studies is that proper open boundary conditions should be employed on non-spatial boundaries when solving the stationary sigma equation (characteristic Neumann equation). Analysis of the two equations for k = 0 [10] has been a detour to arrive at this conclusion. An alternative formulation (not based on an absorbing potential and not overdetermined) of non-spatial open boundary conditions for the sigma equation has been suggested in [9].
We point out that the considerations above concern the stationary case as well as the transient case. Transient simulation of RTDs using an absorbing potential is demonstrated in [19]. In transient simulations, the absorbing potential avoids artificial reflections on non-spatial boundaries. The requirement to impose properly formulated open boundary conditions is in line with simulation