The Wigner Function of Ground State and One-Dimensional Numerics

The study of the ground state of a many-body quantum system plays an important role in a variety of areas such as geometry optimization of molecules, photon absorption spectra of atoms and molecules, and linear response theory in the molecular dynamics.

The ground state can be obtained by solving the fundamental governing equation, i.e., Schrödinger equation, in quantum mechanics. However, due to the curse of dimensionality, direct numerical study of Schrödinger equation via classical mesh-based approaches is intractable even for small molecule such as methane. Hence, approximate solution of Schrödinger equation has been a long-standing research topic, in which many pioneer works have been done. For example, quantum Monte Carlo methods [2, 47] uses stochastic methods to evaluate integrals arising in the many-body problems. Quantum Monte Carlo method offers potential to describe directly many-body effect of a quantum system. However, the efficiency of the method suffers from its slow convergence in the simulations. By using a single Slater determinant as an ansatz for the many-body wavefunction, the main task of Hartree-Fock method[12] is to solve a set of equations derived from a variational method. With acceleration techniques for the self-consistent field iteration and quality solver for the generalized eigenvalue problem, the efficiency of Hartree-Fock method becomes acceptable for large scale system from practical problems, which makes the method very popular in the quantum computational chemistry community even nowadays. However, the lack of electron correlation would introduce large deviations from the experimental results, which limits the application of the Hartree-Fock method. Density functional theory[21, 30] combines the advantages from above two methods. Theoretically[20], it has been proved that three dimensional ground state electron density is a fundamental quantity in a given many-body system, and both the electron exchange and the electron correlation are described in derived Kohn-Sham model[40]. Numerically, the techniques developed for Hartree-Fock method can be borrowed for solving Kohn-Sham model. Even better, with the application of local basis functions for the wavefunction[4, 27], the numerical efficiency can potentially be further improved by using fast solvers for sparse system.

Besides the conventional Schrödinger wave function in Hilbert space, the Wigner phase-space quasi-distribution function [45] provides an equivalent approach to describe quantum object that bears a close analogy to classical mechanics [49]. Moreover, the intriguing mathematical structure of the Weyl-Wigner correspondence has also been employed in some advanced topics, such as the deformation quantization[50]. The Wigner formalism has been applied to a variety of situations ranging from atomic physics [43] to quantum electronic transport [33, 44] and many-body quantum systems[39]. Furthermore, both pure and mixed states of a quantum system can be handled in a unified approach by Wigner function[8]. All theoretical advantages motivate the research on developing models and numerical methods for finding the Wigner functions of a quantum system.

Different from the situation for Schrödinger equation that there have been lots of mature approximate models and numerical methods, more efforts are needed towards the Wigner functions. The first attempts to simulate quantum phenomena by Wigner function were [15, 16] for one-dimensional one-body case. Recently, several methods were designed for the simulation based on Wigner function, such as cell average spectral element method[41], moment method[26, 17], WENO-solver[11], Gaussian beam method[48], etc. While there were also various stochastic methods, e.g., signed particle Wigner Monte Carlo method[31, 32, 34] and path integral method[23, 22, 5]. In many-body situation, the Wigner based simulation was achieved by advective-spectral-mixed method[46], Monte Carlo method[35, 37, 38] and the method based on branching random walk[42]. It is known that in a dynamic study of a given system, an initial state of the system should be specified, which is the ground state of the system in most cases. It is noted that although there have been works mentioned above for dynamics of a given quantum system, the work towards the Wigner functions of the ground state is rare. To our best knowledge, only [36] proposed a feasible framework to handle both time-dependent and time-independent problem based on Monte Carlo method, while no result on the deterministic method for a many body quantum system can be found from the literature, even for the simplest one-dimensional two-body case.

In this paper, in the category of deterministic approach, with the aid of density functional theory, the ground state Wigner function is explored both theoretically and numerically for a given many-body quantum system. More specifically, we firstly derive an eigenvalue problem of energy operator for Wigner function based on the stationary Wigner equation. Then the correspondence between Wigner eigenfunction and Schrödinger eigenfunction in the sense of construction is deduced to guarantee the validity of calculation of the ground state Wigner function through solving this eigenvalue problem. Focusing on the one-dimensional case, a numerical method based on imaginary time propagation method [9, 25, 3] is designed for the solution of the proposed eigenvalue problem. A quantum hydrodynamic model proposed in [6] is simplified in our work to reduce the dimension of the problem, while a reconstruction method is proposed to resolve the well-posedness issue introduced by truncating the approximation. Four examples are tested to show the effectiveness of our method. The numerical convergence of the method can be observed clearly in all numerical experiments. Furthermore, the capability of the numerical method on calculating the excited states of the system, and on handling the case with singular potential, is also demonstrated successfully in the harmonic oscillator and the hydrogen examples, respectively. More importantly, the potential of our method for the ground state calculation of large-scale systems is also shown obviously in the last two examples with effective potentials in density functional theory. It is worth mentioning that compared with the Monte Carlo approach, our method is more robust in the sense that random initial guess of the Wigner function is adopted in our simulation. It is known that a general issue of Monte Carlo method is its slow convergence, while our method successfully demonstrates the theoretical convergence rate of linear finite element method. The method can be extended to higher-order cases in a natural way.

The rest of this paper is organized as follows: In Section 2 we briefly introduce Wigner function and stationary Wigner equation. The derivation of the eigenvalue problem and the correspondence between Wigner eigenfunction and Schrodinger eigenfunction are demonstrated in Section 3. In Section 4, the discretization along $\mathbf{p}$ direction is provided. One-dimensional numerics based on imaginary time propagation method is considered in Section 5. Four numerical examples are presented in Section 6. In Section 7 we conclude this paper and discuss the direction of future work. The conversion between wave function and the coefficient functions and related numerical benefit are illustrated in appendix A.

For convenience, we only consider the Hartree atomic units $\hbar=m=e=1$, where $\hbar$ is the reduced Planck constant, $m$ is the effective mass of electron, and $e$ is the positive electron charge. The range of appearing integrals is from $-\infty$ to $\infty$ without extra explanation.

In this section, the Wigner formalism and the widely used stationary Wigner equation will be introduced briefly.

We define the density matrix

where $\psi_{j}$ is the $j$-th eigenfunction of the time-independent Schrödinger equation

and $P_{j}$ is the probability to find the $j$-th eigenstate. The Wigner function $f(\mathbf{x},\mathbf{p})$ in the $2D$-dimensional phase space $(\mathbf{x},\mathbf{p})\in\mathbb{R}^{2D}$ is defined by applying Wigner-Weyl transform to the density matrix

By taking integral of (3) with respect to $\mathbf{p}$, we get the particle density

Following the basic property of the Weyl transform, the energy in the Wigner formalism can be expressed by

Specially, if the density matrix corresponds to a pure state $e^{iA}\psi(\mathbf{x})$, where $\psi(\mathbf{x})$ is a real-valued function, with the fact that Wigner function is also a real-valued function, we can derive that

Thus the Wigner function of any pure state with a constant phase factor is even with respect to $\mathbf{p}$.

Following the method in [19], we can deduce the stationary Wigner equation from the time-independent Schrödinger equation

where $\Theta[V]f$ is a non-local pseudo-differential operator defined by

and the Wigner potential reads

If the potential $V\in C^{\omega}({\mathbb{R}^{D}})$, (8) can be locally expressed by means of Taylor expansion

where $\mathbf{\lambda}$ is a $D$-dimensional multi-index, $\mathbf{\lambda}!=\prod_{j=1}^{D}\lambda_{j}!$, $\mathbf{x}^{\mathbf{\lambda}}=\prod_{j=1}^{D}x_{j}^{\lambda_{j}}$, and

We close this section by the following theorem, which illustrates the effect of stationary Wigner equation from the perspective of Fourier transform:

In order to find the Wigner function of ground state, the Wigner analogy of the eigenvalue problem with respect to energy operator is needed. Starting from the eigenvalue problem proposed in [19], Lemma 3.1 is deduced for deriving Wigner eigenvalue problem. Finally, to guarantee the validity of the ground state calculation based on Wigner function, the correspondence between Schrödinger eigenfunction and Wigner eigenfunction is established in Theorem 3.1.

In [19], the general eigenvalue problem for an operator $\hat{A}=A(\hat{\mathbf{x}},\hat{\mathbf{p}})$ has been formulated as

from which (14) and (16) can be immediately obtained. The alternative forms (15) and (17) can then be derived via integration by parts. Now we consider the Weyl transform of the energy operator $H(\mathbf{x},\mathbf{p})=|\mathbf{p}|^{2}/2+V(\mathbf{x})$. By Lemma 3.1, using (7) to cancel the imaginary part we can define the Wigner analogy of the eigenvalue problem corresponding to energy operator

where

If the potential $V$ is analytic on $\mathbb{R}^{D}$, the integral in (19) has local expression similar to (10)

The eigenvalue problem (19) and stationary Wigner equation (7) are exactly the real and imaginary parts of the “star-genvalue problem” introduced in [10], respectively. It is proved in [10] that the Wigner function satisfying star-genvalue problem corresponds to a pure state wave function.

Furthermore, the correspondence between Schrödinger eigenfunction and Wigner eigenfunction can be derived as follows.

The above theorem guarantees the correspondence between Schrödinger eigenfunction and Wigner eigenfunction in the sense of construction, which validates the ground state calculation based on Wigner function by solving the eigenvalue problem with the constraint of stationary Wigner equation. Furthermore, based on the construction in theorem 3.1, the conversion between Schrödinger eigenfunction and Wigner eigenfunction expressed by coefficient functions is established in B. Below we will focus on the Wigner formalism to develop our numerical schemes.

To deal with the global integral operators and reduce the dimension of problem, we consider the quantum hydrodynamic model proposed in [6]. It is noted that the hydrodynamic model in [6] is developed for time-dependent simulations, and is derived based on a series expansion of the Wigner function with respect to the momentum variable $\mathbf{p}$, where the basis functions vary for different positions and times. Here in the computation of ground states, we fix the bases and expand the Wigner function as

where $\mathbf{\alpha}$ is a $D$-dimensional multi-index. The basis function $\mathcal{H}_{\mathbf{\alpha}}$ is defined as

where $\operatorname{He}_{n}(x)$ is the $n$-degree Hermite polynomial

To deduce the equations of coefficient functions corresponding to stationary Wigner equation and the eigenvalue problem, we need the following useful properties of Hermite polynomial [1]:

1. 1.

   Orthogonality: $\int\operatorname{He}_{m}(x)\operatorname{He}_{n}(x)\exp(-x^{2}/2)dx=m!\sqrt{2\pi}\delta_{m,n}$;

2. 2.

   Recursion relation: $\operatorname{He}_{n+1}(x)=x\operatorname{He}_{n}(x)-n\operatorname{He}_{n-1}(x)$;

3. 3.

   Differential relation: $\operatorname{He}_{n}^{\prime}(x)=n\operatorname{He}_{n-1}(x)$.

Combining the last two relations we find

Therefore

A direct result of the orthogonality of Hermite polynomials is

Due to the fact that the Wigner function of any pure state is even with respect to $\mathbf{p}$, in the situation of ground state calculation, it is reasonable to assume

Finally, to avoid a system with infinite unknowns, we only consider $f_{\mathbf{\alpha}}(\mathbf{x})$ with $|\mathbf{\alpha}|\leqslant M$, where $M$ is an even positive number.

Particularly, for two Wigner functions $f_{a}(\mathbf{x},\mathbf{p})$ and $f_{b}(\mathbf{x},\mathbf{p})$ corresponding to two orthogonal eigenstates $\psi_{a}$ and $\psi_{b}$, respectively, following the property of Weyl transform, using integration by part we obtain (For detailed illustration, one can refer to [8])

where $f_{\mathbf{\alpha}}^{(a)}(\mathbf{x})$ and $f_{\mathbf{\beta}}^{(b)}(\mathbf{x})$ are the coefficient functions of $f_{a}(\mathbf{x},\mathbf{p})$ and $f_{b}(\mathbf{x},\mathbf{p})$, respectively, and (For detailed derivation, see A.)

Then we will use the properties mentioned above to derive governing equations of the coefficient functions corresponding to the stationary Wigner equation and the eigenvalue problem. The procedure follows the Petrov-Galerkin spectral method, which can also be regarded as taking moments on both sides of the equations.

Now we restrict ourselves to one-dimensional case, and take $M=2K$, where $K\in\mathbb{N}$. Taking $\alpha=2k+1$ in (38) we find

The corresponding eigenvalue problem is given by taking $\alpha=2k$ in (41)

It is noted that (45) depicts an underdetermined eigenvalue system of coupled functions, which is difficult to directly solve by commonly-used method such like block Schur factorization due to its coupled structure. Instead we consider solving it from the perspective of time propagation. In next subsection we will introduce imaginary time propagation method for the utilization of ground state calculation. It is worth to mention that the truncated system is underdetermined due to the appearance of $f_{2K+2}$, a reconstruction approach is proposed in Subsection 5.2 to achieve a well-posed system. Finally, the numerical detail is shown in the last subsection.