1 Introduction

During the development of modern and future electronic devices, frequently quantum mechanical effects in carrier transport processes have to be accounted for. In cases such as these, the Wigner formalism presents a convenient formulation in phase space, thereby allowing many classical notions and concepts to be reused [7, 8].

In our deterministic method, the critical discretization of the diffusion term is done through the utilization of an integral formulation of the Wigner equation [2, 12]. The developed method describes the time evolution of the solution as a superposition of fundamental wave packages.

In cases where physical quantities vary over many orders of magnitude in the phase space, the deterministic methods are the only possible approach, because stochastic methods lack the necessary precision [6].

2 The deterministic method

The basis for the following considerations is the Wigner equation in one dimension

$$\begin{aligned} \frac{\partial f_\mathsf {w}(x,k,t)}{\partial t} - \frac{\hbar \, k}{m^*} \; \frac{\partial f_\mathsf {w}(x,k,t)}{\partial x} = \int \mathrm {d}k^{\prime } \, V_\mathrm {w}(x,k-k^{\prime }) \; f_\mathsf {w}(x,k^{\prime },t). \end{aligned}$$
(1)

It describes the evolution of the Wigner function \(f_\mathsf {w}(x,k,t)\) under the influence of the Wigner potential \(V_\mathrm {w}\). This Wigner potential is derived from the electrostatic potential by

$$\begin{aligned} V_\mathrm {w}(x,k)=\frac{1}{i \hbar }\frac{1}{L_\mathrm {coh}} \int \limits _{-L_\mathrm {coh}/2}^{L_\mathrm {coh}/2} e^{-i 2 k s} \left[ V(x+s) + V(x-s)\right] \mathrm {d} s, \end{aligned}$$
(2)

with \(L_\mathrm {coh}\) the coherence length. A detailed explanation of the aspects of the Wigner formalism can be found, for instance, in [11]. In contrast to classical probability distributions, the distribution in the Wigner picture may also take negative values, therefore named quasi-distribution functions.

The promoted deterministic approach uses the integral formulation of the Wigner equation. The integral form [10, 13] is obtained by considering the characteristics of the Liouville operator on the left-hand side of (1), which are the Newtonian trajectories x(., t) initialized with \(x^{\prime }\), \(k^{\prime }\), \(t^{\prime }\)

$$\begin{aligned} x(x^{\prime },k^{\prime },t^{\prime },t) = x^{\prime } + \frac{\hbar }{m^*} k^{\prime } (t-t^{\prime }). \end{aligned}$$
(3)

This approach has already been used for the development of stochastic solvers, which rely on the corresponding Neumann series of the integral equation [1, 4, 9]. The existence of the partial time derivative guarantees the uniqueness of the solution to the studied evolution problem [3].

The discretization of the quantities in phase space (x, k) and time t is performed in the simulation domain on an equidistant ortho grid

$$\begin{aligned} x=n \varDelta x, \ k=m \varDelta k, \ t=l \varDelta t, \quad \forall \; n \in [0 \dots N],\ m \in [-M \dots M],\ l \in [0 \dots L]. \end{aligned}$$
(4)

The developed deterministic method describes the evolution of the discrete Wigner function \(f_\mathsf {w}(\upnu ,\upmu ,\uptau )\) at a certain point \((\upnu \varDelta x, \upmu \varDelta k)\) in phase space at time \(\uptau \varDelta t\) by

$$\begin{aligned} f_\mathsf {w}(\upnu ,\upmu ,\uptau ) = \sum _n \sum _m {f}_0(n,m) \cdot {q}_{\upnu ,\upmu ,\uptau }(n,m), \end{aligned}$$
(5)

which can be interpreted as the superposition of fundamental wave packages \({q}_{\upnu ,\upmu ,\uptau }\) weighted by the initial distribution \({f}_0\) of the Wigner function. For each point \((\upnu ,\upmu )\) in discrete phase space, an initial distribution \({q}_{\upnu ,\upmu ,0}\) at \(\tau =0\) defined by

$$\begin{aligned} {q}_{\upnu ,\upmu ,0}(n,m)=\delta _{\upnu ,n} \,\delta _{\upmu ,m} \end{aligned}$$
(6)

evolves into a single fundamental solution \({q}_{\upnu ,\upmu ,\uptau }\)

$$\begin{aligned}&{q}_{\upnu ,\upmu ,\uptau }(n^{\prime },m^{\prime }) = e^{-\sum \limits _{j=0}^1 \gamma (N^{\prime }(0,j))\varDelta t} {q}_{\upnu ,\upmu ,\uptau -1}(N^{\prime }(0,1),m^{\prime })\nonumber \\&\quad +\,\varDelta t \sum _{l=0}^{1} \sum _m \; \Gamma (N^{\prime }(0,l),m,m^{\prime }) \; e^{-\sum \limits _{j=0}^{l}\gamma (N^{\prime }(0,j))\varDelta t} {q}_{\upnu ,\upmu ,\uptau -l}(N^{\prime }(0,l),m) \;\omega _l \end{aligned}$$
(7)

in the simulation domain. Here, \(N^{\prime }(l^{\prime },l)\) describes the discrete version of the trajectories (3)

$$\begin{aligned} N^{\prime }(l^{\prime },l) = n^{\prime } + int \Bigl ( \frac{\hbar }{m^*} \frac{\varDelta k \varDelta t}{\varDelta x} m^{\prime } (l-l^{\prime })\Bigr ). \end{aligned}$$
(8)

The fact that the discrete version of the trajectory presumes integral numbers is accounted for by the int operator in (8) and will be discussed later.

The introduced functions \(\Gamma\) and \(\gamma\) depend on the discrete Wigner potential

$$\begin{aligned} \Gamma (x,m,m^{\prime })\!=\!V_\mathrm {w}^+(x,m\!-\!m^{\prime }) + V_\mathrm {w}^-(x,m^{\prime }\!-\!m) \!+\!\gamma (x)\,\delta _{m,m^{\prime }} \end{aligned}$$
(9)
$$\begin{aligned} \gamma (x)=\sum _m V_\mathrm {w}^+(x,m), \end{aligned}$$
(10)

with

$$\begin{aligned} V_\mathrm {w}^+(x,m) = \left\{ \begin{array}{ll} V_\mathrm {w}(x,m)&{} \text {if} \, V_\mathrm {w}(x,m)>0,\\ 0&{}\text {else,} \end{array} \right. \end{aligned}$$
(11)
$$\begin{aligned} V_\mathrm {w}^-(x,m) = \left\{ \begin{array}{ll} -V_\mathrm {w}(x,m)&{} \text {if} \, V_\mathrm {w}(x,m)<0,\\ 0&{}\text {else.} \end{array} \right. \end{aligned}$$
(12)

2.1 Computation

Equation (7) shows a mixture of terms \({q}_{\upnu ,\upmu }\) in \(\uptau\), which are the unknown values of the actual time step, and in \(\uptau -1\), which are the already calculated values of the previous time step. Reordering the unknowns to the left side delivers an equation system for each unknown distribution \({q}_{\upnu ,\upmu ,\uptau }\)

$$\begin{aligned}&{f}_{\upnu ,\upmu ,\uptau } = \mathbf {f}_0^T \cdot \mathbf {q}_{\upnu ,\upmu ,\uptau } \end{aligned}$$
(13)
$$\begin{aligned}&\mathbf {A}_{\upnu ,\upmu } \cdot \mathbf {q}_{\upnu ,\upmu ,\uptau } = \mathbf {B}_{\upnu ,\upmu } \cdot \mathbf {q}_{\upnu ,\upmu ,\uptau -1}. \end{aligned}$$
(14)

The vector \(\mathbf {q}_{\upnu ,\upmu ,\uptau }\) represents the unknown distribution of a single fundamental solution for \((\upnu ,\upmu )\) with its based equation system of complexity \(2\,M\!\times \!N\). For each combination \((\upnu ,\upmu )\) in the simulation domain, the equation system has to be solved consecutively for each time step \(\uptau\).

The weighting factors \(\omega _l\) originate from the numerical time integration in the continuous representation of (7). Unfortunately, these factors appear in both the left- and right-hand side of the equation system (14). According to the trapezoidal rule, factors of \(\omega _l=0.5\) are used.

In cases where the Wigner potential stays stationary over time, \(\mathbf {q}_{\upnu ,\upmu ,\uptau }\) stays stationary, too. Consequentially, the equation system may be calculated only once for the first time step. Subsequent values of \({f}_{\upnu ,\upmu ,\uptau }\) can be achieved by reinsertion of previously calculated distributions of \(f_\mathsf {w}\) in place of \({f}_0\) in (13), which reduces calculation times dramatically.

Equation (7) describes the computation of time step \(\uptau\) from its preceding time step. However, variants of calculating the fundamental solutions \(\mathbf {q}_{\upnu ,\upmu ,\uptau }\) from multiple previous time steps, or even from the initial distribution, seem feasible. Actually, the rank of the resulting equation systems will be boosted rapidly. As these equation systems have to be solved for each \((\upnu ,\upmu )\) , the applicability of that method will be problematic shortly for a few number of implied time steps.

2.2 Interpolation

For common device dimensions, the examination of the trajectory (8) detects areas of low velocities as problematic. By a simple truncation to integer values, the deterministic trajectories stop evolving, especially for small values of m'. Apart from a simple spatial resolution enhancement, which mainly impacts the simulation times, the movement could be adapted by accounting for the accumulation in positional error.

However, a resulting step-wise moving wavefront may result in increasing amplifying oscillations in the solution, especially near corners in the potential distribution.

In this case, an interpolation of \(N^{\prime }\) between the grid points and, as a consequence, to the involved values of \({q}_{\upnu ,\upmu ,\uptau }\) delivers best results. Attention has to be paid to the fact that information is lost in this connection and edges may broaden during this interpolation.

3 Application to the time evolution of wave packages

The Wigner methodology is often validated by the application to wave packages, as it describes the dependencies of the quantities, not only in spatial direction but also in k-space, in a comfortable way. Also the developed method is well suited to examine wave packages [5].

3.1 A single Gaussian wave package

The complex wave function of a minimum uncertainty wave packet [7] around the spatial position \(x_0\) is described as

$$\begin{aligned} \uppsi _{x_0}(x) = A\, e^{-\frac{(x-x_0)^2}{2\sigma ^2}} e^{i k_0 x}, \end{aligned}$$
(15)

with \(\sigma\) the standard deviation and \(k_0\) the wave number of the package. Via the factor A, the probability function is normalized to unity

$$\begin{aligned} \int \limits _{-\infty }^{\infty } \uppsi ^*(x) \, \uppsi (x) \,\mathrm {d}x = 1. \end{aligned}$$
(16)

The corresponding Wigner function allows a simple analytic presentation if the coherence length is let to infinity, giving rise to a continuous k. Accordingly, the factor \(1/L_\mathrm {coh}\) is replaced by \(1/2\pi\)

$$\begin{aligned} f_\mathsf {w}(x,k) = \frac{1}{2\pi }\int \limits _{-\infty }^{\infty } \uppsi ^{*}(x+\frac{y}{2}) \;\uppsi (x-\frac{y}{2}) \;e^{i k y} \,\mathrm {d} y. \end{aligned}$$
(17)

For the Wigner function of the one-dimensional Gaussian wave package in phase space (x, k), this yields

$$\begin{aligned} f_\mathsf {w}(x,k) = B \, e^{-\frac{(x-x_0)^2}{\sigma ^2}} e^{-(k-k_0)^2 \sigma ^2}. \end{aligned}$$
(18)

Here, the factor B normalizes the Wigner function over x and k, satisfying

$$\begin{aligned} \int \limits _{-\infty }^{\infty } \int \limits _{-\infty }^{\infty } f_\mathsf {w}(x,k) \,\mathrm {d}k \,\mathrm {d}x= 1. \end{aligned}$$
(19)

3.2 Two superposed states

In the following, a superposition of the wave functions of two wave packages will be examined. Two Gaussian packages are located at \(-x_0\) and \(+x_0\), respectively, which propagate in the same direction with same speed. The resulting normalized wave function can be evaluated as the sum of both packages

$$\begin{aligned} \uppsi (x) = \frac{1}{\sqrt{2}} \left[ \uppsi _{x_0} (x) + \uppsi _{-x_0} (x) \right] . \end{aligned}$$
(20)

Accordingly, the phase space representation delivers

$$\begin{aligned} f_\mathsf {w}(x,k) = C \Bigl [ e^{-\frac{(x-x_0)^2}{\sigma ^2}} e^{-(k-k_0)^2 \sigma ^2} +e^{-\frac{(x+x_0)^2}{\sigma ^2}} e^{-(k-k_0)^2 \sigma ^2} + 2\,e^{-\frac{x^2}{\sigma ^2}} e^{-(k-k_0)^2 \sigma ^2} \cos 2 \,x_0\, (k-k_0) \Bigr ]. \end{aligned}$$
(21)

In detail, the Wigner function consists of three parts. These are the two parts stemming from each separate package of the wave function plus an additional term in between from the mixed product of the separate wave packages in (17). The latter shows identical envelope to each single package and shows oscillatory behavior in k. Notably, the frequency depends directly on the distance of the packages \(x_0\) , while its amplitude is not affected with an amplitude two times of a separate package. A similar setup is discussed in [5], for instance.

Fig. 1
figure 1

Surface plot of the Wigner function of two superposed wave states. Two Gaussian packages at \(-x_0\) and \(+x_0\) (where \(x_0=60\) nm) plus the correlation state in the middle (at \(x=0\)) can be seen, which shows oscillatory behavior in k-direction. The Wigner function is drawn without the factor C to achieve a proper scaling.

Full size image

The according phase space representation is shown in Fig. 1. In the following, the Wigner functions and density distribution functions are drawn without the scaling factor C to achieve normalized functions in the range \([-2,2]\). The transformed system results in three contributions, which are the two contributions called the Gaussian states, which can be seen on the left- and right-hand sides, plus an additional term called the correlation state, which is shown in the middle. In this regime, in contrast to classical distribution functions, the Wigner function also takes negative values. This oscillatory behavior of the analytic Wigner function, Eq. (21), poses some numerical requirements for the discrete counterpart, which will be discussed later.

The density distribution function in x-space is calculated by integration or, in the discrete case, by summation over all k-values of the Wigner function

$$\begin{aligned} f(x,t) = \int f_\mathsf {w}(x,k,t) \,\mathrm {d}k. \end{aligned}$$
(22)

In the analytic Wigner function (21), the oscillations in k-space are canceled out. The resulting density function

$$\begin{aligned} f(x) = D \,\Bigl ( e^{-\frac{(x-x_0)^2}{\sigma ^2}} + e^{-\frac{(x+x_0)^2}{\sigma ^2}} + e^{-\frac{x^2}{\sigma ^2}} \cdot e^{-\frac{x_0^2}{\sigma ^2}} \Bigr ) \end{aligned}$$
(23)

shows contributions of three Gaussian packages, each with standard deviation \(\sigma\), at position \(x_0\) and \(-x_0\) and the third contribution in the middle at \(x=0\), where the latter one is weighted by the factor \(e^{-\frac{x_0^2}{\sigma ^2}}\), which is smaller by magnitudes than the other two Gaussians. The result of the discrete counterpart is shown in Fig. 2, where the oscillations are canceled by the summation over k.

Contrary to the previous calculation, in phase space the oscillatory correlated state must not be neglected as it shows contribution double the size than the other two Gaussians.

Fig. 2
figure 2

Density distribution function f(x) (without the factor D) of two superposed wave packages at the initial stage. The picture does not differ from the view of two classical wave packages, where only two separate wave packages at \(-x_0\) and \(+x_0\) (with \(x_0=60\) nm) could be seen. The correlation state nearly cancels out in the view of the density distribution

Full size image

In this way, the well-known analytical behavior of the minimum uncertainty Wigner function is used to validate our discrete numerical approach.

Figure 3 shows the situation in Wigner space when the packages traveled to the right. The density distributions of the Wigner function show a broadening in spatial direction. Additionally to the movement, the Wigner function shows a warping in dependence of the k-value. Contrary to the previous calculation, in phase space the oscillatory correlated state must not be neglected as it shows contribution double the size than the other two Gaussians.

Fig. 3
figure 3

Surface plot of the Wigner functions of the wave packages after time evolution. No barrier is used and the packages move undistorted. In phase space, a warping of the distributions is detected. This is caused by the fact that different regions in k move with different velocities.

Full size image

3.3 Introduction of a barrier

In this section, the behavior after passing a barrier will be examined. A barrier at the right side of the packages will be applied. One by one the packages will pass the barrier, where a fraction passes the barrier and another fraction is reflected. One might find the Gaussians are relatively smooth and wide, so the grid spacing may be noncritical and the coherence length may be chosen relaxed.

Contrarily, for the correlation state an appropriate estimate for the grid spacing has to be found. In solving the equation system, the fluctuations in k have to be resolved properly by the grid spacing \(\varDelta k\). Even though the Gaussian wave function is relatively smooth, the sampling in k-space may cause undersampling effects of the oscillation. This results in a nonphysical behavior of the reflected wave and peaks in the distribution function in the order of the entire distribution. An estimate for \(\varDelta k\) can be found by comparison with the period of the oscillation in (21)

$$\begin{aligned} 2 \,k\, x_0 = \frac{2 \pi }{T} k\quad \text {and} \quad T > 2 \varDelta k \qquad \Rightarrow \qquad \varDelta k < \frac{\pi }{2 x_0}, \end{aligned}$$
(24)

with a sampling two times the period in k. This gives an estimate of \(\varDelta k\) in dependence of the spatial distance of the initial packages.

As the coherence length \(L_\mathrm {coh}\) is related to \(\varDelta k\) by

$$\begin{aligned} \varDelta k = \frac{\pi }{L_\mathrm {coh}}, \end{aligned}$$
(25)

estimate (24) gives a minimum

$$\begin{aligned} L_\mathrm {coh}>2 x_0, \end{aligned}$$
(26)

which implies that the coherence length has to spawn the whole region of both packages. This seems also plausible in a physical point of view.

Figure 4 shows a comparison of simulations with and without correlation state after the right package reached a rectangular barrier.

Correlated packages with \(\sigma =10\,\mathrm {nm}\) at \(-x_0\) and \(x_0=40\,\mathrm {nm}\) are used. The rectangular barrier is located at \(x=60\,\mathrm {nm}\) with a height of \(0.3\,\mathrm {eV}\) and a width of \(5\,\mathrm {nm}\). Underlying simulations are performed on a simulation grid with N=500, M=90 which spawns a simulation domain of \(250\,\mathrm {nm} \times \pm 1\,\mathrm {nm}^{-1}\).

The introduced criterion delivers a maximal \(\varDelta k=0.039\,\mathrm {nm}^{-1}\) which corresponds to a minimal coherence length of \(L_\mathrm {coh}=80\,\mathrm {nm}\). The upper image shows a simulation with correlation state, and the lower image shows two wave packages calculated independently. In this case, the incident waves will be nearly completely reflected at the barrier. The right package was already reflected and changed its direction. The influence of the correlation state can be detected when the reflected right package reaches the left package. At this point, the correlation state nearly stopped and is about to change direction.

Fig. 4
figure 4

Comparison of two simulations after the right package was reflected by a barrier, located at \(x=60\,\mathrm{nm}\). The left package travels to the right, whereas the right package already changed its direction. The upper figure is calculated with a correlation state, whereas in the lower figure two separate wave packages are simulated. The correlation state nearly stopped moving at this moment and is on the point to change direction.

Full size image

3.4 Discussion

In the Wigner point of view, all three resulting distributions of the correlated wave packages can be simulated and studied independently, which shows the following observation:

Clearly, the left package did not reach the barrier and is not affected by it, whereas the right package already was reflected.

An interesting effect can be observed within the correlated package, although the envelope of the correlation state did not reach the barrier, it is already influenced by it. The movement nearly stopped and started turning its direction. This can be clearly explained by the fact that the correlation state also holds information of the right, already reflected package.

However, in the numerical method the originating behavior of the packages has to be accounted for by a proper selection of the grid spacing in k-direction or, reciprocally affected, an adequate coherence length.

4 Conclusion

In this paper, a deterministic method for solving the Wigner equation has been presented. Through its application to superposed wave packages, the adequate handling of the correlated state has been demonstrated, although strict limits on the resolution in k-space have to be satisfied. These bounds may represent a challenge for the solution, which has to be analyzed in further work.