Numerical methods for the Wigner equation with unbounded potential

Numerical methods for the Wigner equation with unbounded potential

Zhenzhu Chen222LMAM and School of Mathematical Sciences, Peking University, Beijing 100871, China.,    Yunfeng Xiong222LMAM and School of Mathematical Sciences, Peking University, Beijing 100871, China.,    Sihong Shao222LMAM and School of Mathematical Sciences, Peking University, Beijing 100871, China. 111To whom correspondence should be addressed. Email:
July 5, 2019

Unbounded potentials are always utilized to strictly confine quantum dynamics and generate bound or stationary states due to the existence of quantum tunneling. However, the existed accurate Wigner solvers are often designed for either localized potentials or those of the polynomial type. This paper attempts to solve the time-dependent Wigner equation in the presence of a general class of unbounded potentials by exploiting two equivalent forms of the pseudo-differential operator: integral form and series form (i.e., the Moyal expansion). The unbounded parts at infinities are approximated or modeled by polynomials and then a remaining localized potential dominates the central area. The fact that the Moyal expansion reduces to a finite series for polynomial potentials is fully utilized. Using a spectral collocation discretization which conserves both mass and energy, several typical quantum systems are simulated with a high accuracy and reliable estimation of macroscopically measurable quantities is thus obtained.

AMS subject classifications: 81Q05; 65M70; 81S30; 45K05; 82C10

Keywords: Wigner equation; Moyal expansion; spectral method; quantum dynamics; unbounded potential; uncertainty principle; double-well; Pöschl-Teller potential; anharmonic oscillator

1 Introduction

Unbounded potentials are ubiquitous in quantum mechanics, especially in simulating quantum tunneling phenomena ranging from various branches of physics and chemistry. As a typical example among them, the double-well potentials with two minima separated by a barrier have been widely used in understanding the transition of quantum states [1, 2] as well as in modeling the potential energy surface of small molecules like the ammonia and the methane [3, 4]. The often-used double-well potentials can be simply characterized by unbounded polynomial potentials plus bounded localized potentials, and related quantum observables are calculated through either the Schrödinger equation [3] or the path integral formalism [5]. In this work, we turn to adopt the Wigner function [6], a quasi-probability distribution, to investigate the quantum tunneling effects, because it grants us a natural description of quantum observables in a statistical form due to the Weyl correspondence [7, 8].

However, solving the Wigner equation that describes the time evolution of the Wigner function in the phase space is usually a tough task because of the difficulties in tackling the nonlocal and highly oscillating pseudo-differential operator. The situation becomes worse when the unbounded potentials are taken into account. The existed deterministic solvers, including the finite difference schemes [9, 10] and the spectral collocation methods [11, 12, 13], always require the potentials to decay fast and vanish at infinities, i.e., the localized potentials, since the Wigner kernel is evaluated by the Poisson summation formula, and thus is not applicable for the unbounded potentials. For the potentials of the polynomial type, on the other hand, as an equivalent series form of the pseudo-differential operator, the Moyal expansion reduces to a finite series [14, 15] and the resulting equation can be solved by either the spectral method [16] or the Hermite expansion [17]. In this work, we attempt to combine the advantages of the above two to evolve the Wigner quantum dynamics in the presence of unbounded potentials.

Figure 1: An asymmetric double-well potential [5]. It can be decomposed into a polynomial part and a localized part .

Specifically, we focus on the potentials that are sufficiently smooth, with the asymptotic behaviors at infinities governed or approximated by polynomials. Such kind of potentials covers the double-well potentials as mentioned above and is commonly used to fit the observed data (known as the polynomial regression model). For simplicity, we assume that the unbounded can be split into a polynomial potential and a localized one that decays at infinities, e.g., as displayed in Fig. 1. In this manner, two equivalent forms: the integral form and the Moyal expansion, can be employed to deal with and separately and the resulting equation can be solved by different techniques. In this paper, we choose the spectral collocation method because it is able to accurately resolve both the linear differential operators and the Fourier integrals for sufficiently smooth potentials. Meanwhile, discussions on the conservation of both mass and energy, as well as on their numerical counterparts, are performed. In some situations, the operator splitting technique can be further introduced to improve the performance.

The proposed scheme allows us to study the dynamics of many interesting quantum systems, like the Pöschl-Teller potential [18, 19] and the double-well systems [17, 20]. Several macroscopically measurable quantities, such as the differences in energy levels, the quantum tunneling rate, and the autocorrelation function can also be obtained with a satisfactory accuracy. Moreover, the uncertainty principle will be shown directly in phase space by simulating the superposition of a harmonic oscillator perturbed by a sixth-order anharmonic one.

The remaining is organized as follows. Section 2 gives a brief introduction to the Wigner equation with two equivalent forms. In Section 3, the treatment of unbounded potentials as well as the numerical scheme for the resulting equation is presented. We will prove that the proposed scheme can maintain the mass and total energy, and the demonstration of its performance is left for Section 4. Concluding remarks and further discussions are collected in the Section 5.

2 Quantum mechanics in phase space

As a classical mathematical representation of quantum mechanics in phase space, the Wigner function allows a direct connection with the classical picture and its dynamics equation reduces to the classical Vlasov equation as the reduced Planck constant vanishes. In this section, we will sketch the Wigner formalism for quantum mechanics and the exposition is restricted to one-dimensional one-body situation for simplicity. The Wigner function is defined in the phase space for the position and the wavenumber through the Weyl-Wigner transform of the density matrix as follows


If there are stationary states corresponding to energy with , then the dynamics is given by and thus the superposed state of these states is . In consequence, the Wigner function for such superposed state reads


Given a quantum operator at the instant , the expectation value can be calculated by averaging the corresponding Weyl symbol with the Wigner function:


Therefore we can easily deduce from Eqs. (2.2) and (2.3) that can be decomposed into components, the frequencies of which must be proportional to for all . That is, it is readily to obtain the differences in energy levels of the quantum system in question only via a direct spectrum analysis of .

Starting from the Schrödinger equation, it can be shown that the Wigner function follows the following dynamics, i.e., the time-dependent Wigner equation,


where is the mass and is nothing but from Eq. (2.1) through an inverse Fourier transform. Here is the so-called nonlocal pseudo-differential operator containing the quantum information and has different but equivalent expressions as follows.

  • For , the pseudo-differential operator is characterized by a convolution one:


    where is the so-called Wigner potential or Wigner kernel.

  • For , performing the Taylor series for at yields


    and substituting the above expression (2.9) into Eq. (2.5) leads to the Moyal expansion


    Here we adopt the compact notations: and for .

  • For with and , representing a general class of (unbounded) potentials which can be decomposed into two parts (e.g., see Fig. 1). Owing to the linearity of the Fourier transform, the pseudo-differential operator can be rewritten as a linear combination of Eqs. (2.7) and (2.10). More importantly, if is indeed a polynomial or can be approximated at infinities by a polynomial, denoted by , then we only need to consider a much simpler expression:


    where is the Wigner kernel corresponding to the localized potential via Eq. (2.8), give the coefficients from the polynomial potential via Eq. (2.11), and denotes the degree. A key observation in Eq. (2.12) is the Moyal expansion reduces to a finite series which can be readily resolved by standard numerical techniques.

A detailed comparison of the above three expressions (2.7), (2.10), and (2.12) for the pseudo-differential operator shows: Eq. (2.7) fails to hold when is unbounded, while Eq. (2.10) involves infinite terms even when with a compact support. By contrast, Eq. (2.12), albeit in a somewhat complicated form, only requires the asymptotic behavior of at infinities is governed by polynomials. It deals with the unbounded part with a finite series of linear differential term, and captures the fine structure in the central area through a twisted convolution. A simple example is the following unbounded rational-fraction-type potential


which converges to as . It can be easily verified that, the expression (2.7) breaks down because the Wigner kernel is not well defined any more in the classical sense, and the infinite series in the Moyal expansion (2.10) sticks in there and thus is difficult to handle with; on the contrary, the expression (2.12) gets rid of those problems by taking


In summary, the combined expression (2.12) serves as the starting point of this work for investigating the Wigner quantum dynamics in the presence of unbounded potentials.

Before proceeding, we would like to mention two conservation laws that are always used to guide the design of numerical methods. One is the mass conservation stated by the continuity equation


where is the particle density and the current density [13]. The other is the energy conservation


where is the quantum Hamiltonian operator.

3 Numerical methods

Now we turn to seek a numerical approximation to the Wigner equation, where the nonlocal term has the form as in Eq. (2.12). The convolution term poses the first challenge since it involves double integrations. In general, a simple nullification of the distribution outside a sufficiently large -domain is usually adopted [12, 13]. For a sufficiently large -domain , the truncated version of the Wigner equation is




denotes the discretized Wigner kernel for the localized potential and with being the spacing step in -space. Such approximation stems from the Poisson summation formula:


Here we assume that decays and thus ignore the periodic images. A necessary and sufficient condition for the truncated Wigner equation (3.1) to conserve the mass has been given in [13] and reads


where represents the length of -domain. In -space, the popular quantum transitive boundary condition will be adopted hereafter as did in [12].

The spectral (element) collocation methods have been demonstrated in [12, 13] to resolve the oscillations of the Wigner function and thus will be utilized in this work to discretize the truncated Wigner equation (3.1). In particular, a Fourier spectral collocation scheme is adopted in -space and a collocation spectral element method with Gauss-Lobbato points in -space. An explicit fourth-order Runge-Kutta discretization [21] is then employed for the time marching as did in [12]. Actually, the plane wave expansion in -space is a natural choice since both the linear differential operator and the convolution term can be accurately approximated in a compact form. Furthermore, we are able to prove the numerical conservation of both mass and energy for the resulting full discretization.

3.1 The Spectral collocation method

The uniform collocation points in are with and then the plane wave expansion in -space reads




and the coefficients are further determined by a collocation spectral element method with Gauss-Lobbato points in -space for easy implementation of boundary conditions. Let be the computational domain in -space and we divide it into non-overlapping elements as with , and . For simplicity we adopt a uniform mesh in which the number of collocation points keeps the same for all , denoted by , and use the Gauss-Lobbato points in each element. Then the spectral expansion for the coefficients in is



with and . That is, over and can be approximated by


Consequently, the partial derivative of with respect to can be directly obtained as



In a similar way, the partial derivatives of with respect to have a much simpler expression:


With the help of the orthogonal relation of the Fourier basis:


the truncated convolution term can be calculated analytically as follows


where the coefficients are determined by :


Finally, we obtain the following semi-discretizated scheme for the truncated Wigner equation (3.1)


and the fast Fourier transform (FFT) can be used to accelerate the computation.

3.2 Conservation laws

An explicit fourth-order Runge-Kutta method is used to evolve the semi-discretizated scheme (3.14) as we did in [12]. Below we will show that the resulting full discretization scheme conserves both the mass and the energy. To this end, it suffices to consider the one-step forward Euler method with the time step and the resulting fully discretizated scheme is


where and denote the numerical solutions of and at time , respectively.

To illustrate the numerical conservation laws, we need to consider the inner product in the computational domain


and the numerical current density

Proposition 1.

The numerical scheme (3.15) conserves the mass, i.e.,


provided that the total inflow and outflow are in balance, say,


Through integration by parts, it is easy to verify

due to the Eq. (3.19) as well as the periodic condition in -space

So we only need to show .

Splitting the summation with respect to into two parts, one for and the other for , it leads to


where we have applied in order Eq. (3.11) and the fact that for any according to Eq. (3.13). ∎

The numerical energy conservation, however, requires some additional conditions, as stated below.

Proposition 2.

The numerical scheme (3.15) conserves the energy, i.e.,

where , provided that

(a) ;

(b) , ;

(c) , .


We intend to prove the following relations:


Using the integration by parts leads to

where both conditions (b) and (c) are applied to eliminate the boundary terms, and then we arrive at the Liouville theorem

which is nothing but the first relation of Eq. (3.21).

For , by the Taylor theorem, we have

Substituting it into yields

implying that, for the remaining two relations of Eq. (3.21), we only need to verify


both of which must vanish through the integration by parts due to for and the condition (c).

3.3 Splitting treatment

Finally, we would like to mention that the splitting treatments of unbounded potentials used in Eq. (2.12) are very useful when some resulting subproblems allow exact solutions. Moreover, within the framework of the splitting schemes, different methods could be used to tackle and separately. For example, we may solve the truncated Wigner equation (3.1) in a splitting manner


where the subproblem (B) has explicit solutions for the quadratic potential [22].

Next, we consider the conservation laws of the splitting methods. It only needs to consider the simplest Lie-Trotter scheme, with the same spectral collocation method adopted for both subproblems, namely,


and thus we arrive at

Proposition 3.

The splitting scheme (3.24) conserves the mass, i.e.,

provided that

(a) ;

(b) , for .


Comparing the scheme (3.24) with (3.15) and according to Proposition 1, we only need to show

Actually, each term in the left-hand-side of the above equation must vanish.

As for the first term, it is a direct outcome of which can be readily verified through the integration by parts as well as using the periodicity of in -space.

From condition (b), a direct calculation shows

and then performing the integration by parts in -space for the second term yields

The proof is finished. ∎

Unfortunately, the numerical energy conservation fails to hold for the splitting treatment because


cannot be guaranteed for the scheme (3.24). This can be readily verified by taking, for instance, (see more details in Section 4.5).

4 Numerical experiments

In this section, several typical quantum systems are employed to test the performance of the proposed methods and the atomic units are used if not specified. We employ the -error and the -error to study the convergence rate of the spectral collocation method:


where is the numerical solution, and the reference solution which could be either the exact solution or the numerical solution on the finest grid mesh. In order to monitor the numerical conservation of mass and energy, the variations of total mass and energy are also chosen as the metrics,


As we did in [12, 13], all above metrics are evaluated by a simple rectangular rule over a uniform mesh.

4.1 The Pöschl-Teller potential

The Pöschl-Teller potential


serves as the first example for its energy levels and bound states with can be obtained analytically [19, 18].

For instance, when , we have

Figure 2: The Wigner function under the Pöschl-Teller potential. Left: The stationary state in Eq. (4.6); Right: The superposed state in Eq. (4.7).

The Wigner function corresponding to , as shown in the left plot of Fig. 2, reads

Figure 3: The Pöschl-Teller potential: The convergence rate with respect to (left) and (right). The spectral convergence in both -space and -space is observed for . All the errors are measured at the instant .

The Wigner function corresponding to the superposed state is