A WENO-solver combined with adaptive momentum discretization for the Wigner transport equation and its application to resonant tunneling diodes.

We present a novel numerical scheme for the deterministic solution of the Wigner transport equation, especially suited to deal with situations in which strong quantum effects are present. The unique feature of the algorithm is the expansion of the Wigner function in local basis functions, similar to finite element or finite volume methods. This procedure yields a discretization of the pseudo-differential operator that conserves the particle density on arbitrarily chosen grids. The high flexibility in refining the grid spacing together with the weighted essentially non-oscillatory (WENO) scheme for the advection term allows for an accurate and well-resolved simulation of the phase space dynamics. A resonant tunneling diode is considered as test case and a detailed convergence study is given by comparing the results to a non-equilibrium Green's functions calculation. The impact of the considered domain size and of the grid spacing is analyzed. The obtained convergence of the results towards a quasi-exact agreement of the steady state Wigner and Green's functions computations demonstrates the accuracy of the scheme, as well as the high flexibility to adjust to different physical situations.

Introduction

Nowadays, well-established approaches for the simulation of charge transport in semiconductor devices are based on drift-diffusion and hydrodynamic models. Such macroscopic descriptions enable a rather quick computation and therefore, an industrial application. But, owing to the continuous progress made in device fabrication, microscopic descriptions of charge transport become of prime importance. One approach of this kind is to make use of the semi-classical Boltzmann transport equation (BTE) [1,2]. Simulations based on the BTE have attained great success in many different applications, but fail, as soon as quantum mechanical effects dominate the device behavior. A prototypical example for such a device is a resonant tunneling diode (RTD) [3,4]. In this case, a description of the charge carriers as localized particles becomes invalid and a fully quantum mechanical treatment is needed instead. In order to allow for a realistic device simulation, it is essential to employ methods which enable the description of open quantum systems where elastic and inelastic scattering mechanisms are present in addition. Two of the most prevalent approaches of this kind are the non-equilibrium Green's functions (NEGF) [3,5] and the Wigner transport equation (WTE) [4,6,7]. In this work, we will focus especially on a deterministic solution of the latter one. The Wigner function shares many analogies with its classical counterpart, but with the essential difference, that the function is in general not positive definite. Since its marginal distributions are proper probability distributions and since measurable quantities are calculated in the same manner as in the case of the Boltzmann phase space distribution, the Wigner function is commonly termed quasi-distribution as well.

The application of the WTE to the simulation of charge transport has been investigated by many groups in the past decades. Following the pioneering work of Frensley [8], extensions of the basic formalism by including a spatially dependent effective mass [9], by coupling the WTE to Poisson's or Schrödinger's equation [10,11], and by including scattering mechanisms have proven to enable a realistic device simulation [7,12-15]. The solution strategies employed for the WTE are twofold: On the one hand, Monte Carlo schemes have been devised and on the other hand, various discretizations of the WTE for a deterministic solution exist. The Monte Carlo (MC) approaches are either based on the concept of particle affinity [16] or on the concept of particle sign [17]. In particular, the second method has shown to be successful in tackling multi-dimensional [18] and many-body [19] quantum problems. In general, a close agreement with NEGF computations could be achieved and for an introduction to the MC-method we refer to [4,7]. In the case of a deterministic solution of the WTE, the treatment of the pseudo-differential operator, which accounts for the non-local action of the potential, is much more involved than that of terms in classical transport equations. The most common strategy is to make use of the fact that the eigenfunctions of the pseudo-differential operator are plane waves, so that a Fourier transformation allows for an efficient evaluation. The superior scaling of the fast Fourier transform [20] is also made use of in the spectral methods developed in [21] or [22]. The second term in the WTE coincides in the parabolic band approximation with the advection term (also called diffusion term) in the BTE. Well-developed discretization schemes from classical transport theory can thus be employed. In [11,12,23] various higher order finite difference stencils have been applied together with a convergence analysis. The adaptive, high-order scheme in [22] makes use of Gauss–Lobatto collocation points for the advection term.

In this work, we employ a weighted essentially non-oscillatory (WENO) [24] scheme for the advection term, well known from applications within the BTE [2,25,26]. The scheme allows for a highly accurate estimation of the flux in smooth regions of the underlying function, and at the same time, avoids the creation of spurious oscillations at discontinuities or regions with large gradients. Regarding the pseudo-differential operator, even though the fast Fourier transform allows for an efficient evaluation in order steps, certain restrictions are associated with it. On the one hand, the points for the momentum grid have to be chosen equidistantly and on the other hand, in order to be consistent with the continuity equation, the momentum grid and the spatial grid are interrelated and the discretizations of each one cannot be chosen independently. From our point of view, these restrictions may cause difficulties to resolve the Wigner function properly especially in situations where strong quantum effects are present, such as in the case of resonant tunneling devices. We thus propose a novel numerical scheme which allows for a highly flexible choice of the grid in order to adjust to different physical situations. The scheme follows similar ideas as in finite element or also finite volume methods, namely to approximate the Wigner function locally by piecewise polynomials. It is ensured by construction that the discretized WTE is consistent with the continuity equation without further conditions, so that the particle density is conserved for arbitrarily chosen grids. The local approximation enables one to change the resolution in different regions of the phase space by orders of magnitude. This allows us to properly resolve all details of the steady state Wigner function of a RTD, consisting of a rather smooth shape on which very short-scaled and large-valued oscillations are superimposed in certain regions of the phase space.

As a test case for the method we consider the simplest model of a RTD, given by a one-dimensional, static potential with a homogeneous doping profile and a constant effective mass. In addition, scattering mechanisms are neglected and fully coherent transport is treated instead. Despite of these simplifications the considered test case is representative and furthermore, the treatment of the fully coherent regime is most challenging for a discretization of the WTE. It is natural to expect that the oscillations in the Wigner function are damped as soon as scattering is included. Also due to the recent work outlined in [27,28], which questions the applicability of the inflow/outflow boundary conditions for the WTE, a further and more detailed analysis of the coherent transport regime is of interest.

The paper is organized as follows: In Section 2 we recap basic aspects of the Wigner function formalism before discussing the numerical method in Section 3 in detail. After performing an averaging of the WTE, we briefly introduce an operator splitting scheme, and thereafter, focus on the particular discretization of the pseudo-differential operator and of the advection term separately. In Section 4 we present results for the simulation of a RTD. At first, a detailed convergence study is given together with comparing the current–voltage characteristics with a NEGF computation, before presenting well-resolved phase space plots of the steady state Wigner function and results for a transient response simulation. Concluding remarks are given in Section 5.

Wigner transport equation

A possible quantum mechanical description of an open system of electrons where decoherence phenomena such as electron–phonon scattering are present, is given in terms of the von Neumann equation for the density operator . Frensley [8] pointed out that it is in some cases problematic to impose appropriate boundary conditions for the density operator. By making use of the Wigner function it is possible to circumvent this disadvantage [1,6-8,29-31]. The approach is based on transforming the density matrix in position space into a quasi-distribution function in phase space, with and when considering the special case of one dimension. It is done upon employing a Fourier transformation with respect to the difference coordinate ξ: known as the Wigner–Weyl transform, which is invertible [30]. The Wigner function enables an analogous phase space description as in the case of a classical distribution function with the same moments when integrating with respect to momentum p and position x. In particular, the particle density is given as the marginal distribution and the current density is related to the first moment by with the electron charge and denoting the effective mass. However, does not share all properties of a probability distribution since it may take on negative values in some regions of the phase space and is thus termed quasi-distribution function (for an interpretation of these negative-valued regions see for instance [32]). Due to the analogy to a classical distribution function, the same boundary conditions are applicable which clearly distinguish the incoming from the outgoing part of the distribution , see Section 3.4.

The dynamics of the density operator is governed by the von Neumann equation [30] where is the commutator with the Hamiltonian. Similar to Eq. (1), it is possible to map any Hilbert space operator onto a c-number phase space function by the Weyl–Wigner correspondence [30]. For Hermitian operators such as , the Wigner map results in a real-valued phase space function. For the special case of a Hamiltonian of the form , i.e. without products of non-commuting operators, the corresponding phase space function is simply obtained by replacing the operators by their eigenvalues, [29]. In place of the non-commutative product of operators, the so-called Moyal star product is encountered [29] where an arrow to the left (right) indicates that the derivative acts only on functions to the left-hand (right-hand) side. A product of operators is transformed as [29] Therefore, upon applying a Wigner transformation to the von Neumann equation (4), the phase space analogue, Moyal's equation, is found [29]: The Moyal bracket is defined by [29,33]

For the Hamiltonian considered in this work, Moyal's equation simplifies further. In particular we choose with the kinetic energy term in the effective mass approximation , where is the energy of the conduction band edge, and with an arbitrary electrostatic potential . Due to this approximation, the kinetic part of the Moyal bracket reduces to and is thus simply given by an advection term [34]. When rewriting the potential energy term in Eq. (7) in form of an equivalent integral expression, we arrive at the so-called Wigner transport equation (WTE) [1,6,35] with the pseudo-differential operator [1] also known as the Wigner kernel, and the Fourier transformed Wigner function defined as When comparing the WTE to the semi-classical Boltzmann transport equation [1], it is apparent that the kinetic part is represented in the same way by an advection term but the electrostatic potential enters the WTE in a more complicated, non-local way through .

The pseudo-differential operator acts in the Fourier transformed space as a simple multiplication of and the multiplicator which is called the symbol of the pseudo-differential operator [1].

One fundamental property of the WTE is that the continuity equation can be retrieved. A short calculation reveals that the contribution of vanishes when integrating it with respect to p: When making use of the relations of to and , Eqs. (2) and (3), respectively, and integrating the whole WTE (10) with respect to p, the continuity equation is readily obtained:

Numerical method

Averaging the Wigner transport equation

In the following we average the Wigner transport equation over grid cells defined by with In order to construct a conservative method, we also need the half-bounded intervals with The grid for the x variable is chosen to be equidistant but all of the following derivations can be directly adopted to a non-equidistant x grid. No constraints are assumed for the grid spacing of the p variable, which enables a highly flexible and adaptable grid for different physical situations. In particular, the cell boundaries for the x and p grids are defined by The central points of the grid cells are labeled by integer indices and given by

The particular approximations used for the Wigner function will be discussed later in Sections 3.3 and 3.4. At this point we only demand that for . This restriction to a compact support of may introduce finite size errors in actual computations but can always be controlled by increasing the domain size, as demonstrated in Section 4.2. The cell averages of the Wigner function over the grid cells are given by To arrive at a conservative scheme, we now average the whole Wigner transport equation (10) over grid cells. For the interior grid cells , with interior referring to the p variable, we arrive at with , where the left-sided term was identified as the time derivative of the cell average. To obtain the governing equations for (), we consider the corresponding finite intervals () together with the semi-infinite p intervals (): Since is assumed to vanish outside of , the left-hand side of Eq. (24) reduces to whereas we need to consider for the right-hand side the integration over the full semi-infinite p interval. On the whole, this results in the following expression for the boundary terms : and for the case we obtain:

The equations stated so far suffice to show that the numerical scheme developed in this work preserves the particle density, independent of the particular discretization of the momentum grid. From the knowledge of the cell averages, the averaged particle density in the interval is calculated byWhen inserting the expressions for the time derivatives of the cell averages, Eqs. (23), (26) and (27), into the time derivative of Eq. (28), we arrive at where the property Eq. (14) was made use of. As can be seen, the action of the pseudo-differential operator does not influence the particle density and the expression in the last line corresponds to the net flux into the interval . This tells us that the averaging procedure outlined here preserves the particle density and is consistent with the continuity equation.

Operator splitting

In the following two sections the numerical scheme is formulated in terms of the cell averages . This is done by choosing for a certain approximation for each operator, which is uniquely determined by the set of cell averages. When introducing the vector notation the averaged WTE, Eqs. (23), (26) and (27), is written in abstract notation as where the operators and correspond to the ones resulting from the advection term and from the pseudo-differential operator, respectively. For the case that and are not explicitly time-dependent, the formal solution is given by By employing a so-called Strang splitting [36,37], the time evolution due to the distinct operators can be treated separately, one after the other: The error encountered in each time step Δt is proportional to the commutator and to , such that the overall accuracy is of second order in Δt. The advantage of employing an operator splitting is that appropriate time stepping methods can be used for the individual operators in the three sub-time steps in Eq. (33).

Discretization of the pseudo-differential operator

In the following we focus on the sub-problem of discretizing the pseudo-differential operator, starting from Eqs. (23), (26) and (27) but leaving out the advection term. For each interior grid cell, the integration with respect to p in Eq. (23) is taken over a bounded interval, thus conceptually not difficult and results in where the definition of the pseudo-differential operator, Eq. (11), was used. To perform the integration over the semi-infinite p intervals in Eqs. (26) and (27) we recall the Fourier transform of the Heaviside step function, given by [38] with PV denoting the principal value, so that we obtain for the case of the lower boundary since . An analogous result is found for the upper boundary , in which is replaced by and an overall minus sign occurs.

After this general considerations, we restrict ourselves to Wigner functions which are representable as piecewise polynomials with respect to p and x. In particular we choose In principle one is not limited in the order γ of the polynomials, but in this work we consider at most a piecewise constant approximation with respect to x and first-order polynomials with respect to p, given by The slopes are determined from the cell averages by central finite differences with respect to p at the interior grid points and by one-sided finite differences at the boundaries . Since a direct numerical evaluation of the oscillatory integrals involved in is problematic, we also choose for a certain basis representation. To be specific, piecewise linear polynomials are used as well: with Here, a continuous form of is chosen but it should be mentioned that this is not mandatory. The grid points for the potential are selected to be a subset of the grid points for the Wigner function, with . Furthermore, two additional grid points outside the x interval for are introduced, and with and , which we will let go to ±∞ in the final equations to model the situation of semi-infinite leads under bias. Schematic drawings of the approximate forms chosen for and are depicted in Fig. 1.

Fig. 1

Schematic drawings of the piecewise linear approximations chosen for the Wigner function in (a) and for the potential in (b).

Inserting the particular approximations for and into Eqs. (34) and (36) results in with and the Fourier transform on the interval given by The integrations in Eq. (42) are carried out analytically, as far as possible at least, since sine and cosine integrals appear. In practice, standard routines are used to evaluate the trigonometric integrals numerically. After lengthy but conceptually not difficult calculations, for details see Appendix A, one arrives at the following set of equations: with the matrix elements given by and The abbreviations stand for whereby the expressions contain the sine and cosine integrals, see Appendix A. The intervals are given by with The corresponding expression for is readily obtained from Eq. (48) by the replacements , and leaving the same. For an actual implementation it is important to know the limits with k representing some wavenumber, and and labeling two distances. They are needed for the diagonal terms where occurs, for more details see Appendix A.

After having stated the governing equations, we now focus on the time stepping. The set of Eqs. (45) together with the chosen linear dependence of the slopes on the cell averages may be rewritten as The formal solution to advance from time step to is given by the system for each spatial index . In order to evaluate the matrix exponential numerically, we make use of the scaling and squaring method [39]. For steady state or transient response simulations, as considered in this work, the extra computational cost for evaluating the matrix exponential is negligible. Fully time dependent problems on the other hand require different approaches, such as a Runge–Kutta scheme for instance.

Discretization of the advection term

In the present section we focus on the advection term and consider the following sub-problem of the averaged WTE, Eqs. (23), (26) and (27): with , . This can be rewritten as with the flux at the two spatial boundaries of grid cell given by The advection term is acting via only on the x coordinate of , so that it is convenient to approximate in this term as a piecewise constant function with respect to p. With respect to x some higher-order approximation will be chosen. The assumption that is constant on the interval greatly simplifies Eq. (57) and we arrive at Eqs. (56) decouple now with respect to m and the problem reduces to solving a set of one-dimensional advection equations. To express the flux in terms of the cell averages , one has to rely on approximate schemes. High-order methods without slope or flux limiters [34,37] were found to be problematic, in the coherent transport regime at least, due to the creation of spurious oscillations. We obtained good results with a MC limiter (monotonized central-difference) scheme [40], but found the convergence rate of the method to be rather low, especially in regions with strong variations of and thus of as well. Since the problem of strong spatial variations of the distribution function is common in device simulations due to steep doping profiles, see e.g. [2], well developed schemes exist that can cope with such cases, known as weighted essentially non-oscillatory (WENO) methods. In this procedure, the numerical fluxes are obtained as the convex sum of a certain number of approximations on different stencils, where the corresponding weights are adjusted to the smoothness of the distribution function. The first such methods were developed in [24]. In this work we will focus on the WENO5 scheme presented in [2] or [25]. The method and especially its stability properties were extensively analyzed in [41].

In the WENO5 scheme, three different stencils are used to approximate the flux such that (Eq. (58)) at time is given by For the case the three single fluxes are determined by The single weights are normalized by the equation with where the values of are fixed by , and . The small quantity prevents the denominator to vanish. The smoothness indicators are determined by the following equations The fluxes and smoothness indicators for are constructed in an analogous way.

In order to treat the boundaries appropriately, ghost cells [34] are introduced which specify values for at the positions and . Boundary conditions for ohmic contacts are applied as described in [2]. For the case of the left-sided contact this results in for and with labeling the one-dimensional Fermi–Dirac distribution of the left contact, see [8,42,43]. The values for the ghost cells on the right-sided reservoir are determined in an analogous way.

To perform the time step, a third-order Runge Kutta method termed SSP(3,3) is employed, motivated by the work of Wang [41]. In this method three sub-steps are needed to advance one complete time step and the abbreviation SSP stands for the strong stability preserving property. When written in terms of the fluxes, the three sub-steps of the SSP(3,3) algorithm consist of where the individual fluxes , and are determined by the WENO5 scheme, Eqs. (59)–(63), out of the set of cell averages , and , respectively.

Results for the simulation of a resonant tunneling diode

Device under consideration and methodology

As a particular test case, we consider a resonant tunneling diode (RTD) consisting of an AlGaAs–GaAs-heterostructure, as depicted schematically in Fig. 2(a). The AlGaAs–GaAs-RTD is a standard problem in the field of device simulations to benchmark numerical methods for quantum transport calculations, see for instance [8,23,44,45]. Here, we restrict ourselves to the simplest case of a homogeneous doping profile, without aiming for a self-consistent solution by coupling the WTE to Poisson's equation. We assume that the bias voltage causes a linear potential drop in the region of the double barrier and that stays constant within the reservoirs, see Fig. 2(b). For the specific composition Al0.3Ga0.7As, the energy band offset has a magnitude of [23]. The length of the reservoirs included in the x domain for the Wigner function is labeled by and varied in different simulations. Furthermore, expressed in terms of the lattice constant of GaAs, we choose for the barriers a width of and a slope of length , as well as for the size of the well. A spatial dependence of the effective mass is neglected and the value for GaAs is used on the whole x domain. The electron distributions and the chemical potentials , inside the contacts are fixed by choosing a donor density of and a temperature of [8,44]. From the knowledge of the chemical potentials and temperatures for each reservoir, the cell averages of the ghost cells are set according to Eq. (64).

Fig. 2

Illustration of the considered RTD consisting of an AlGaAs–GaAs-heterostructure in (a) and the corresponding potential qV(x) in (b). The plot of qV(x) depicts the case of a rather short reservoir length of and a bias voltage of . The chemical potentials and of each reservoir are indicated.

The validity of semi-classical boundary conditions for the WTE as introduced in [8] is a topic under vivid debate, especially after recent works which address the non-uniqueness and the symmetry properties of the Wigner function [27,28,46]. The numerical test cases presented therein are for symmetric potentials for which we cannot provide reliable, i.e. well-resolved, results due to the presence of singular terms in the steady state Wigner functions, see Section 4.3. Other recent studies demonstrate the convergence of the WTE calculations upon increasing the size of the simulation domain [44] as well as possible improvements by adapting the boundary distribution to the physical state of the active device region [47]. Despite their approximate nature we employ inflow/outflow boundary conditions here as well and demonstrate that accurate and physically valid results can be achieved for sufficiently large values of . Due to the problematics with singular terms we present simulations only for non-zero bias voltages .

In order to evaluate the accuracy of the numerical method developed in this work, we compare the steady state curves and particle densities against a non-equilibrium Green's functions (NEGF) calculation. Details on the NEGF method may be found for instance in [3,5,42,45,48-50]. The NEGF simulations were performed with a grid spacing , for which we know that the calculated quantities are well converged.

For the WTE simulations, we choose an equidistant grid for the x variable with different values for the spacing Δx and the reservoir length . For the p grid we make extensive use of the possibility to apply a non-equidistant spacing, in order to resolve all the oscillation patterns of well enough. We specify a certain which determines the p domain by , as well as a maximum spacing for the outermost region of the p grid. All the other of the interior subdivisions are then expressed as a fixed fraction of . Thus, a refinement of by a certain factor causes a refinement of the whole p grid by the same factor. You are referred to Appendix B for some more details on the choice of the p grid.

For the calculation of steady state properties we employ the common strategy to evolve in time until the stationary state is reached. But, it is important to point out that for certain parameter sets of the x and p grid, a non-smooth convergence with a sudden build-up of error was observed. The problem arises due to the fact that the steady state distribution may exhibit heavily oscillating regions in phase space, especially for situations where tunneling and coherence phenomena are prominent. This is exactly the case for RTDs in the resonant tunneling regime. When simulating such a situation with a too coarse grained grid, one encounters that the quantities of interest initially converge towards the true steady state values, but large errors arise as soon as the oscillations in become too short scaled to be resolved by the chosen grid. It is, therefore, mandatory to carefully inspect the time evolution of and, if possible, to check the robustness of the results upon a refinement of the grid spacing. In all simulations performed in the course of this work, the spacing of the p grid was the crucial factor for a smooth convergent behavior and the spacing of the x grid was comparatively unproblematic.

Convergence with respect to the domain size and grid spacing

At first we investigate the impact of the considered domain size, i.e. the size of the p interval and the reservoir length on the calculated quantities. We set the spacing of the x grid to and for the p grid we choose . This combination of the x and p spacings is suited to properly resolve all of the oscillation patterns appearing in the steady state Wigner functions, at least for the presented sizes of the x domain (note that an increase of to larger values than considered here may also require a decrease of ).

In Fig. 3, the calculated curves for simulations with different values of the reservoir length , as well as with different values of are presented. A minimum value of at least is chosen. For smaller values of one observes that the cell averages of the Wigner function at take on significant values. This indicates that the underlying assumption (see Section 3.1) that vanishes outside of the interval , is a too severe restriction for the physical situation, for the chosen value of . As one can see from Fig. 3, the reservoir length has great impact on the accuracy of the calculated curves. Already a value of is large compared to the double barrier size of (see also Fig. 2), but the influence of the boundaries on the overall device behavior is prominent. Increasing the reservoir length to (not plotted) lessens this influence and for a fairly good agreement between the Wigner function and the NEGF calculations is observed in the whole voltage range . Taking the minor differences in the results for and into account, we conclude that a size of the domain determined by and is a reasonable choice for the considered physical problem. A convergence study for the particle density does not give additional insight and is not presented here. A similar dependence on as for the case of the -curves was found, whereby the magnitude of the errors in is much smaller and already results for agree fairly well with the NEGF reference.

Fig. 3

Comparison of the simulated j(V) curves for different values of the reservoir length and different extensions of the p domain, specified by , and . The results depicted on the left are for and those on the right are for . The relative differences of the currents to the NEGF reference are presented in the plots on the bottom. For all of the simulations the grid spacings are chosen to Δx = a and , so that (left plots) and (right plots), and .

We now focus on the question of the convergence of the calculated quantities with respect to the grid spacing. In particular we reduce Δx and consider one fixed but adjusted value of . Once the Wigner function is already well-resolved with respect to the momentum variable, a further decrease of the p spacing has only minor influence on the calculated quantities.

In Fig. 4, the results for three different spacings , and are compared. For the particle densities , the relative errors for the two voltages and are shown in Fig. 4(c), (d). For , only the region of negative differential resistance (NDR) is plotted due to its particular importance for device applications. It is obvious to see a monotonic decrease of the error in the current density with decreasing values of Δx and a very good agreement between the NEGF reference and the WTE calculations for the case , with a relative error below one percent. This clearly demonstrates the ability of the developed algorithm to produce very accurate results and furthermore, the convergence of the calculated quantities when refining the grid parameters.

Fig. 4

Plot (a) depicts the j(V) curve for the NDR-region, (b) the relative difference of the j(V) values to the NEGF solution and (c) and (d) the relative difference in particle density for and , respectively. The results are obtained by using and three different spacings . The size of the (p,x) domain is the same in all three simulations and determined by and . The number of grid points is and .

Wigner functions for the steady state

After these preliminary convergence studies we now focus on more physical aspects and present steady state Wigner functions for which we know that the results are well-converged. Fig. 3 displays the characteristic S-shaped curve of a RTD including a negative differential resistance (NDR) [7]. The peak in at is taken on when resonant tunneling is most prominent. The resonant tunneling decreases in the NDR-region until one enters the regime of conventional tunneling at the valley, , and for higher voltages.

To see how the two physically dissimilar regimes manifest itself in phase space, Fig. 5 displays steady state Wigner functions for the voltages and . In the conventional tunneling regime (Fig. 5(d)), the Wigner function is negative-valued only in small regions of the phase space and takes on a rather smooth shape, close to a classical Boltzmann distribution. In the resonant tunneling regime (Fig. 5(a)–(c)), on the other hand, the emergence of heavily oscillating regions in phase space is apparent. On the one hand, rather long-scaled oscillations are present throughout large parts of the phase space and on the other hand, a very sharp and large-valued stripe of oscillations builds up around . The detailed phase space dynamics are involved but one can see in the reflected part of the distribution on the left-hand side that momenta of approximately are extenuated, corresponding to the part for which a resonant tunneling process is accessible. On the right-hand side, one can see the outgoing, accelerated beam at . The sharp stripe of oscillations at is closely related to the resonant tunneling process itself and to the coherence between the reflected and transmitted part at and , respectively. At least this is what we conclude after analyzing various simulations and performing analytical calculations for simple tunneling situations. When solving for the Wigner function of a single, unbiased tunneling barrier and for the case of incoming plane waves, singular terms of the form and together with oscillating prefactors with respect to x emerge (see also [51]). A superposition of plane waves with different wavenumbers causes a change of the form . Such a behavior is consistent with . The question may arise why the simulations reveal a sharply peaked and heavily oscillating behavior of but no singular terms. To our understanding this is a consequence of the non-zero bias voltage combined with a continuous potential. Tests for the present RTD and a bias of were performed, with the result that sharp oscillations arose in the central region at which could not be properly resolved upon decreasing the p spacing. We thus believe that it is in general advisable to consider only non-zero bias voltages and preferably situations with a continuous . Due to the presence of the singular terms and for the case of unbiased barriers, we cannot provide reliable simulations with the present algorithm to the situations discussed in [27,28]. We can only report on the observation that unphysical solutions were encountered indeed, but solely in the case of under-resolved momentum grids.

Fig. 5

Illustrations of the steady state solutions of f(p,x,t) obtained by using the parameters , , and for two different bias voltages. Depicted in (a), (b), (c) are the solutions for and in (d) the one for . The number of grid points for the two cases are , and , for and , respectively.

From the knowledge of the detailed shape of the steady state Wigner functions in phase space it is possible to understand the particular dependence of the -curves on the reservoir length . As in the case of classical transport theory, the application of the inflow/outflow boundary conditions is only an appropriate approximation if, in the vicinity of the boundaries, the gradients of with respect to x are negligible. It corresponds to the requirement that the action of the pseudo-differential operator close to the boundaries and thereafter, of course, is not significant. From Fig. 3 it is obvious to see that the low-bias regime of , where resonant tunneling is prominent, converges much slower with respect to than the high-bias regime of , where conventional tunneling dominates. For the case of the conventional tunneling regime, one can note from Fig. 5(d) that exhibits only comparatively weak oscillations outside the region of the double barrier. In contrast, the results in Fig. 5(a) are characterized by strong and long-ranged oscillations of in the resonant tunneling regime. Clearly, the gradients of with respect to x at position are much smaller in the former case than in the latter one.

Transient response simulation

The Wigner function formalism gives direct access to time-resolved quantities and enables one to address time-dependent situations. In the following we present results for a simple example, namely the large-signal transient response of a RTD [8]. Of particular interest for this are the two dissimilar situations at and , corresponding to the peak and to the valley in the curves, respectively. For the simulations termed peak-to-valley we start from the steady state Wigner function for the peak voltage and consider an abrupt switching in bias voltage at , and for the case valley-to-peak the same is done in the opposite way.

To examine the dynamics in detail, Figs. 6–8 depict the time evolution of the current density as well as of the Wigner function. For the case of the peak-to-valley simulation, it is apparent to see from Fig. 6(a) an initial rise of the current density in the region of the double barrier, which then propagates to the contact on the right-hand side. This corresponds to the part of the electron distribution which occupied the well state at and is then accelerated to higher momenta by the increased bias voltage, see also Fig. 7 at and . Overall, the relaxation of to the new stationary value is much faster on the upwind side (left-hand side) of the barrier than on the downwind side. Fig. 6(b) depicts the time evolution of for the case of switching from valley to peak. Similar to the previous case, the time scale for the change in current density is shorter on the upwind side than on the downwind one. The abrupt rise in current density on the upwind side can be accounted to the reduction in reflection of electrons on the left-hand side of the barriers, as a result of the resonant tunneling of some of the electrons through the well state. In phase space, this effect is visible as an extenuation of the Wigner distribution on the upwind side at the corresponding momenta, see Fig. 8 at and .

Fig. 6

Time evolution of the current density j(x,t) for the transient response simulations, in (a) for peak-to-valley and in (b) for valley-to-peak. The parameters for the simulations are Δx = a, , and , i.e. and . In both cases the same p grid is used, namely the one optimized for the peak voltage.

Fig. 7

Time evolution of f(p,x,t) in phase space for the peak-to-valley case. For the simulation parameters see Fig. 6.

Fig. 8

Time evolution of f(p,x,t) in phase space for the valley-to-peak case. For the simulation parameters see Fig. 6.

The time-resolved phase space plots of in Fig. 8 allow one to investigate the detailed dynamics of the advent of the sharp stripe of oscillations at . In the beginning, the oscillations arise in the center at , rather parallel to the p axis and with a long-scaled modulation, before spreading out to larger values of x together with becoming more short-scaled. In the course of this evolution the oscillations turn more parallel to the x axis in sort of a shear movement. As a result, a very short-scaled modulation of with respect to p is formed, which requires an extremely fine p grid to resolve the steady state Wigner function properly (). Upon increasing , the required resolution with respect to p increases further. We believe that this effect is most prominent for the ballistic case as considered here, and that one can expect a damping of the oscillations as soon as dissipative mechanisms, such as electron–phonon scattering are included.

Conclusion and outlook

In this work, a novel numerical scheme for the deterministic solution of the Wigner transport equation has been presented. The central aspect of the method is to allow for a highly flexible and adaptive choice of the simulation domain. A detailed study of convergence is given by comparing the WTE computations to a reference solution which we obtained with the NEGF method. The results clearly demonstrate the convergence of the algorithm and a quasi-exact agreement between the WTE and the NEGF calculations could be obtained. As indicated already in other studies [44], the application of the inflow/outflow boundary conditions may yield physically valid and correct results, provided a sufficiently large distance between the active device region and the boundaries of the simulation domain is included. For large enough computational domains, a monotonic convergence of the current–voltage characteristics towards the NEGF reference could be achieved upon refining the grid spacing, with a minimal relative error below one percent. Subsequent to the convergence study, steady state Wigner functions for the RTD operating in the resonant and in the conventional tunneling regime have been compared and discussed. The two physically dissimilar regimes manifest itself in phase space by drastically different Wigner functions. In the resonant tunneling case, a sharp and large-valued stripe of oscillations has been discovered, accounting for the coherent superposition of reflected and transmitted states. In the final part of this work, a transient response simulation has been considered together with investigating the dynamics of the current density and the Wigner function.

Possible extensions of the method and further investigations are manifold. Since RTDs are meant to operate in the THz regime, fully time-dependent calculations are of interest. As investigated already thoroughly by other groups, a self-consistent solution of the Wigner–Poisson system including scattering mechanisms is important for a realistic description of transport at room temperature and could be also combined, of course, with the approach presented here.