Numerical methods for the Wigner equation with unbounded potential
Abstract
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 timedependent Wigner equation in the presence of a general class of unbounded potentials by exploiting two equivalent forms of the pseudodifferential 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; doublewell; PöschlTeller 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 doublewell 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 oftenused doublewell 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 quasiprobability 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 pseudodifferential 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 pseudodifferential 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.
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 doublewell 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öschlTeller potential [18, 19] and the doublewell 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 sixthorder 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 onedimensional onebody situation for simplicity. The Wigner function is defined in the phase space for the position and the wavenumber through the WeylWigner transform of the density matrix as follows
(2.1) 
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
(2.2) 
Given a quantum operator at the instant , the expectation value can be calculated by averaging the corresponding Weyl symbol with the Wigner function:
(2.3) 
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 timedependent Wigner equation,
(2.4)  
(2.5)  
(2.6) 
where is the mass and is nothing but from Eq. (2.1) through an inverse Fourier transform. Here is the socalled nonlocal pseudodifferential operator containing the quantum information and has different but equivalent expressions as follows.

For , the pseudodifferential operator is characterized by a convolution one:
(2.7) (2.8) where is the socalled Wigner potential or Wigner kernel.

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 pseudodifferential 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:
(2.12) 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 pseudodifferential 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 rationalfractiontype potential
(2.13) 
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
(2.14) 
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
(2.15) 
where is the particle density and the current density [13]. The other is the energy conservation
(2.16) 
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
(3.1) 
where
(3.2) 
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:
(3.3) 
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
(3.4) 
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 GaussLobbato points in space. An explicit fourthorder RungeKutta 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
(3.5) 
where
(3.6) 
and the coefficients are further determined by a collocation spectral element method with GaussLobbato points in space for easy implementation of boundary conditions. Let be the computational domain in space and we divide it into nonoverlapping 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 GaussLobbato points in each element. Then the spectral expansion for the coefficients in is
(3.7) 
where
with and . That is, over and can be approximated by
(3.8) 
Consequently, the partial derivative of with respect to can be directly obtained as
(3.9) 
with
In a similar way, the partial derivatives of with respect to have a much simpler expression:
(3.10) 
With the help of the orthogonal relation of the Fourier basis:
(3.11) 
the truncated convolution term can be calculated analytically as follows
(3.12) 
where the coefficients are determined by :
(3.13) 
Finally, we obtain the following semidiscretizated scheme for the truncated Wigner equation (3.1)
(3.14) 
and the fast Fourier transform (FFT) can be used to accelerate the computation.
3.2 Conservation laws
An explicit fourthorder RungeKutta method is used to evolve the semidiscretizated 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 onestep forward Euler method with the time step and the resulting fully discretizated scheme is
(3.15) 
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
(3.16) 
and the numerical current density
(3.17) 
Proposition 1.
The numerical scheme (3.15) conserves the mass, i.e.,
(3.18) 
provided that the total inflow and outflow are in balance, say,
(3.19) 
Proof.
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 .
The numerical energy conservation, however, requires some additional conditions, as stated below.
Proposition 2.
(a) ;
(b) , ;
(c) , .
Proof.
We intend to prove the following relations:
(3.21) 
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
and
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
(3.22) 
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 LieTrotter scheme, with the same spectral collocation method adopted for both subproblems, namely,
(3.23) 
and thus we arrive at
(3.24) 
Proposition 3.
(a) ;
(b) , for .
Proof.
Comparing the scheme (3.24) with (3.15) and according to Proposition 1, we only need to show
Actually, each term in the lefthandside 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. ∎
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:
(4.1)  
(4.2) 
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,
(4.3)  
(4.4) 
As we did in [12, 13], all above metrics are evaluated by a simple rectangular rule over a uniform mesh.
4.1 The PöschlTeller potential
The PöschlTeller potential
(4.5) 
serves as the first example for its energy levels and bound states with can be obtained analytically [19, 18].
For instance, when , we have
The Wigner function corresponding to , as shown in the left plot of Fig. 2, reads
(4.6) 
The Wigner function corresponding to the superposed state is