{'fraction_non_alphanumeric': 0.03735516227589785, 'fraction_numerical': 0.014152302991871793, 'mean_word_length': 4.654335071707953, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 6, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
The Wigner Function of Ground State and One-Dimensional Numerics

In this paper, the ground state Wigner function of a many-body system is explored theoretically and numerically. First, an eigenvalue problem for Wigner function is derived based on the energy operator of the system. The validity of finding the ground state through solving this eigenvalue problem is obtained by building a correspondence between its solution and the solution of stationary Schr\"odinger equation. Then, a numerical method is designed for solving proposed eigenvalue problem in one dimensional case, which can be briefly described by i) a simplified model is derived based on a quantum hydrodynamic model [Z. Cai et al, J. Math. Chem., 2013] to reduce the dimension of the problem, ii) an imaginary time propagation method is designed for solving the model, and numerical techniques such as solution reconstruction are proposed for the feasibility of the method. Results of several numerical experiments verify our method, in which the potential application of the method for large scale system is demonstrated by examples with density functional theory.


Introduction
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 twobody 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 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 = m = e = 1, where 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 −∞ to ∞ without extra explanation.

Wigner function and stationary Wigner equation
In this section, the Wigner formalism and the widely used stationary Wigner equation will be introduced briefly.
We define the density matrix where ψ 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 (x, p) in the 2D-dimensional phase space (x, p) ∈ R 2D is defined by applying Wigner-Weyl transform to the density matrix By taking integral of (3) with respect to 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 ψ(x), where ψ(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 p.
Following the method in [19], we can deduce the stationary Wigner equation from the time-independent Schrödinger equation where Θ[V] f is a non-local pseudo-differential operator defined by and the Wigner potential reads If the potential V ∈ C ω (R D ), (8) can be locally expressed by means of Taylor expansion where λ is a D-dimensional multi-index, λ! = D j=1 λ j !, x λ = D j=1 x λ j 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.
If F is a polynomial, then (14) can be further simplified as For A(x, p) = G(p), similarly, we have the corresponding eigenvalue problem and if G is polynomial, we have In [19], the general eigenvalue problem for an operatorÂ = A(x,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(x, p) = |p| 2 /2 + V(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 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 "stargenvalue 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.
Using the method of separation of variables we find Therefore the eigenfunction has the form To distinguish the Wigner function corresponding to pure state from other eigenfunctions in Remark 3.1, we have following relation: where δ i j is Dirac delta function.
Furthermore, the correspondence between Schrödinger eigenfunction and Wigner eigenfunction can be derived as follows. Proof. The derivation of the eigenvalue problem (19) and stationary Wigner equation (7) has clearly shown that the Wigner-Weyl transform of ψ satisfies these two equations.
On the contrary, let f (x, p) be an eigenfunction of H w with eigenvalue E, which satisfies (7). Applying the inverse Wigner-Weyl transform we have where A is the normalization constant. It can be verified by direct calculation and substitution of (7) that which means such construction yields a Schrödinger eigenfunction with the same energy.
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 Appendix B. Below we will focus on the Wigner formalism to develop our numerical schemes.

Hermite expansion of Wigner function
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 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 α is a D-dimensional multi-index. The basis function H α is defined as where 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]: Combining the last two relations we find 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 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 α (x) with |α| M, where M is an even positive number. Particularly, for two Wigner functions f a (x, p) and f b (x, p) corresponding to two orthogonal eigenstates ψ a and ψ b , respectively, following the property of Weyl transform, using integration by part we obtain (For detailed illustration, one can refer to [8]) where f (a) α (x) and f (b) β (x) are the coefficient functions of f a (x, p) and f b (x, p), respectively, and (For detailed derivation, see Appendix 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.

Stationary Wigner equation
For convenience, we take f α = 0 for any multi-index α with at least one negative component. Following the same method in [7], we get the series expansion fo the first term in (7) as Using (30), the pseudo-differential operator term Θ[V] f with expression (10) becomes Plugging (36) and (37) into (7) and equating the coefficient of each basis function to zero, we obtain This equation holds for every α with |α| being odd.

Eigenvalue problem
Using the recursion relation of Hermite polynomials, the first two terms in (19) can be expanded as With the help of the differential relation, the last term with local expression (21) can be written as Similar to the derivation of (38), the orthogonality of Hermite polynomials yields In the above equation, we choose α such that |α| is even. So far we have already derived the general eigenvalue problem and the constraint based on stationary Wigner equation for our model. In next section, we will restrict ourselves to one-dimensional case, and provide a feasible numerical framework on the strength of imaginary time propagation method.

One-dimensional numerics
Now we restrict ourselves to one-dimensional case, and take M = 2K, where K ∈ N. Taking α = 2k + 1 in (38) we find must be a finite expectation value, we have lim x→−∞ f 2k (x) = 0 for all k 0, integrating both sides of (42) from −∞ to x we obtain Thus f 2k+2 can be calculated by { f 2l } k l=0 . Once given the density ρ(x) = f 0 (x), we could construct { f 2k } ∞ k=1 recursively. Therefore in one-dimensional case, the Wigner function of pure state is determined by the density.
The corresponding eigenvalue problem is given by taking α = 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.

Imaginary time propagation method
Imaginary time propagation (ITP) method is a widely-used method for solving eigenvalue problems such as the time-independent Schrödinger equation. Its basic idea is to use the fact that for any Hermitian operator H and large t, we have exp(−tH)Ψ is approximately an eigenfunction of H associated with its smallest eigenvalue, given that Ψ is not orthogonal to the corresponding eigenspace. In order to prevent the eigenfunction from being too large/small, renormalization is applied during the time propagation. In the context of finding the ground state, we can illustrate this idea by considering the time-dependent Schrödinger equation where H s = −∇ 2 /2 + V is the Hamiltonian in the Schrödinger picture. We adopt Wick rotation of the time coordinate, i.e., t = −iτ, then (46) becomes with the formal solution Since the ground state eigenfunction ψ 0 has the lowest eigenvalue E 0 , given the initial state that is not orthogonal to ψ 0 , the normalized steady-state of (47) yields the ground state of the time-independent Schrödinger equation (2) as other exponentials decay more rapidly. To find the excited state, we only need to propagate several functions simultaneously with orthogonalization process after each step. Following the similar way to derive (38), with the help of (47), direct calculating the time derivative of Wigner function in the expression of (3) we obtain the equations for ITP. Recall the eigenvalue problem for the Wigner function (45), the ITP governing equations for coefficient functions are the ones by replacing the right-hand side of (45) with negative time derivative: where k = 0, 1, . . . , K.
It is desired to mention that the constraint of stationary Wigner equation makes the gap between the smallest eigenvalue and second-smallest eigenvalue larger, i.e., from (E 0 + E 1 )/2 − E 0 to E 1 − E 0 , which accelerates the convergence process. Although the infinite system satisfies the constraint of stationary Wigner equation if the initial condition satisfies such constraint. Due to the closure and numerical error arising during time evolution, we enforce this condition at each time step. The evolution of { f 2k } K k=0 is depicted by the following two-step procedure where (50a) is rewritten from (49). The projection step (50b) is to guarantee that (42) holds. In next subsection, we will propose a reconstruction method to deal with the underdetermined system (50a) and utilize the projection (50b).

Reconstruction
Notice that (50a) depicts a underdetermined system due to the appearance of f 2K+2 , to close the system we consider the reconstruction of f 2K+2 . As discussed at the beginning of this section, f 2K+2 is determined by { f 2k } K k=0 from stationary Wigner equation (42). On the other hand, stationary Wigner equation is a constraint to be met for ground state calculation. With appropriate boundary condition, we consider the reconstructing f 2K+2 based on (42). To achieve a system with finite unknowns, we adopt truncation on our domain, i.e., use [−a, a] instead of R, where a is a large positive number. Since lim x→±∞ f 2k = 0 as previous discussion, it is reasonable to use the approximation f 2K+2 (−a) = 0 as boundary condition for reconstruction.
It is deserved to mention that such reconstruction is not confined to f 2K+2 . Since the approximation f 2k+2 (−a) = 0 also works for 0 k < K, start from f 0 , we can construct f 2k for k = 1, 2, . . . , K in turn, and finally obtain coefficient functions { f 2k } K k=0 satisfying (50b). Therefore, such reconstruction is also adopted for the utilization of (50b).

Numerical discretization
Since the Crank-Nicolson scheme is unconditionally stable and second-order in time, it is adopted for temporal discretization as following approximation With regard to spatial discretization, we consider finite element method (FEM). The coefficient functions { f 2k } K k=0 are approximated by the linear combination of piecewise-polynomial basis functions {φ i }, where {φ i } is defined on a set of real space interpolation nodes {p i }. Denote coefficients by {ϕ k,i }, then we have: where N basis stands for the dimension of space V h = span{φ i , 1 i N basis }. And φ i is the i-th basis function that is typically chosen such that φ i (p j ) = δ i, j , δ i, j is Kronecker delta function. As a result, we have ϕ k,i = f 2k (p i ), for 1 i N basis and 0 k K. Let the superscript (n) denote the corresponded term at time τ = n∆τ, where ∆τ is the time step, and denote ϕ k = [ϕ k,1 , ϕ k,2 , . . . , ϕ k,N basis ] T . Introducing the approximation (51) to finite element discretization of the equation (49) within the subspace V h , we obtain where In the equations above, and Ω = [−a, a] is the truncated domain. The finite element discretization of (42) yields the equations for reconstruction and M (m) is defined as above. While the normalization condition of The main components of the flow chat are introduced as follows  A straightforward idea is to consider the strategy to restrict solution space. Based on our numerical experiment, although the convergence to ground state can be observed without such strategy, it can accelerate the convergence process. Furthermore, this strategy becomes necessary if the excited state is taken into consideration. Particularly, when we have a symmetric potential, it can be easily proved that the Schrödinger eigenfunction is either odd or even, which always leads to an even density. Therefore artificially taking f 0 to be even is a reasonable choice to restrict the solution space in such cases.

Numerical experiments
In this section, the numerical effectiveness of our method is verified by following four examples. The first example is a quantum harmonic oscillator, which is a fundamental example employed as a sanity check. Meanwhile, the influence of size of domain for the simulation is studied in this example, and the ability of calculating the excited states using our approach is also demonstrated by supplementing Algorithm 1 with additional orthogonalization process. Second, we consider a one-dimensional hydrogen system, where the solution contains singularity, and thus we need at least quadratic finite elements to produce a reliable ground state solution. The third and fourth examples are, respectively, a two-particle Hooke's atom system and its variation, where both systems are represented in the context of the density functional theory. Numerical results from these two examples show the potential of our method for simulating large scale systems.
In all examples, the computational domain is chosen as [−10, 10] with a uniform mesh unless otherwise specified. In each example, different mesh sizes h and truncation orders K are tested to show the convergence rate and the influence of truncation order. Error estimation is given by the infinity norm.

Harmonic oscillator
We consider the potential of harmonic oscillator For the sake of simplicity, we take ω = 1. Then the Wigner function of the n-th eigenstate is given by [18] where L n (x) is the n-th Laguerre polynomial. We first discuss the choices of truncation on domain. The numerical results obtained from different sizes of domain are demonstrated in Figure 6.1. It can be observed from Figure 6.1 that with a small domain a = 2.5, the desired numerical convergence cannot be observed. However, with the increment of the domain size, theoretical convergence rate of linear finite element method can be successfully obtained with a = 5, 10, 15. To balance numerical accuracy and efficiency, in our following simulations, a = 10 is always used.
Next, we would like to show the difference between the algorithms with/without the projection step (50b). It follows the discussion in the previous sections that introducing the constraint of stationary Wigner equation helps accelerate the convergence, as is also confirmed by our numerical experiments. It is illustrated in the Figure 6.2 that the error of the first excited state fails to decrease to 1.0 × 10 −10 even with sufficiently long time. Based on our numerical experiment, it is worth mentioning that the constraint of the stationary Wigner equation also contributes to suppressing some possible numerical or modelling errors accumulated during the simulation. In fact, divergent results have been observed in our experiments without such a constraint. Besides the ground state, it is noted that by adding an orthogonalization process in Algorithm 1, excited states of the system can be calculated simultaneously. Furthermore, the effect of truncation order K for the numerical solution is also studied in our numerical experiments. The numerical results are provided in Figure 6.3. The numerical observations of Figure 6.3 are summarized as follows: i) In the computation of all the states, secondorder convergence with respect to the mesh size is observed. ii) The calculation of ground state converges fast with respect to K; our numerical result for K = 5 almost coincides with that for K = 15. The convergence is slower for the excited state calculation, but the behavior still looks like the spectral convergence due to the smoothness of the solution. iii) Larger truncation order K is needed when higher excited state is calculated. This might be caused by the accumulation of numerical error from the orthogonalization process. In addition, since the 2nd excited state has higher energy, more terms in Hermite expansion are needed to depict its more oscillatory behavior.

One-dimensional hydrogen
The one-dimensional hydrogen system consists of an electron moving in the one-dimensional potential with boundary condition ρ(x) = 0. It is shown in [28] that the ground state density is It follows (3) that its Wigner function reads In our numerical experiment, it is found that linear finite element method fails to provide a convergent result, which can be explained as follows. Substituting f 0 = 2x 2 e −2|x| + ε 0 into the following numerical construction of f 2 we find Therefore in the element containing original point, we have the numerical error ε 0 /x 2 = O(1) which prevents the numerical convergence to the exact solution. To tackle this problem, we consider the quadratic finite element basis functions in this example, which yields theoretical error O(h) by (65). Numerical results are shown in Figure 6.2, where we present both the error of all coefficient functions and the error of density. The following observations can be made from the result: i) The error of density (orange curve) is better than that of all coefficient functions (blue curve), as the normalization condition acts on density. ii) It is observed that when the mesh size is large, the convergence rate is higher than the theoretical one. As the mesh is refined, the error gradually saturates due to the truncation of the infinite system, and eventually the overall error is dominated by the reconstruction of f 4 when the mesh size is sufficiently small. Note that for larger K, the higher-order derivatives of V turn out to be more singular, leading to the requirement of higher-order numerical methods to guarantee the convergence.

Hooke's atom
In this example we consider the Hooke's atom consisting of two electrons oscillating in the parabolic well. The two-particle Schrödinger equation is where the external potential of Hooke's atom is It follows the derivation in Appendix B that the ground state density is where C = (16 √ π + 10π) −1/2 , and erf(x) is the error function defined as Since erf(x) is odd by its definition, ρ(x) is an even function. This two-particle system can be transformed into the Kohn-Sham equation where The Kohn-Sham orbitals of (70) are given by With ψ = ψ 1 = ψ 2 and ε = ε 1 = ε 2 , and with the assumption that the functional derivative of V H + V XC vanishes at infinity, the Kohn-Sham potential can be exactly expressed as [24] V KS (x) = ε + 1 2 Since ρ(x) is even, V KS (x) is also an even function. The numerical results for K = 0, 1, 2 with our method are listed in following table. It can be observed in Table 6.1 that i) the numerical accuracy of the density is better than the high-order coefficients which is similar to the previous example. ii) For K = 0, the numerical convergence towards the theoretical result can already be obtained. However, with the refinement of mesh grids, the convergence order starts to decrease, since for small K, the reconstruction error dominates when mesh size is sufficiently small. iii) The above issue is remedied when K becomes larger. For K = 1, 2, the theoretical convergence order with respect to the grid size is obtained in all the simulations. Furthermore, comparable numerical accuracy can be observed from results with K = 1 and K = 2, which means that in this example, a small K can already provide sufficient numerical accuracy.

Contact-interacting Hooke's atom
Now we consider a more practical example by replacing the interacting function in (66) with delta function where V ext (x) is defined as (67). This two-particle system can be also transformed into the Kohn-Sham equation as (70) with (71). The energy of system is given by [13] where The local-density correlation energy functional is [29] where a = −1/24, b = −0.00436143, d = 0.252758 and e = 0.0174457. By these definitions, we have To solve (74), a self-consistent iteration is employed due to the nonlinearity [14]. To validate our method, the numerical result of the Schrödinger wave function is calculated as a reference. Both the numerical results of our method with K = 0, 1, 2 and the ones of Schrödinger wave function are listed in Table 6.2. Here both the Wigner function and the wave function are solved on different grids ranging from h = 0.2 to h = 6.25 × 10 −3 , and the "energy" columns list the ground state energy for all our simulations. For the "error" columns, we list the difference of the particle densities ρ w and ρ finest s , where ρ w stands for the Wigner function computed using different parameters h and K, and ρ finest s is obtained from the wave function computed on the finest grid. It is shown in Table 6.2 that i) for each fixed K, the convergence towards the Schrödinger results can be observed as the mesh is refined; ii) results of energies for two cases coincide with each other very well, especially when mesh size is sufficiently small (h ≤ 0.025). These observations again validate both the model and the numerical method proposed in this paper. In Table 6.3, the relative errors of numerical solution are demonstrated with a result obtained on the finest mesh. It can be found that the convergence rate is around the theoretical one of linear finite element method.

Conclusions
In this paper, a model of Wigner function of eigenstate is proposed, providing a theoretical foundation for ground state calculation in Wigner formalism. With a simplified model of [6], ITP method is adopted for the realization of one-dimensional ground state calculation. For validation purposes, we applied our method to the simulation of several benchmark systems. The aim of first two experiments is to show the capability of our method to handle excited state and singularity of potential. While in the last two examples, with the assistance of DFT, our approach successfully converges to the ground state. In particular, the consistency of results from our method with Schrödinger solution is observed in the last example, indicating the feasibility of our method in DFT regime.
Our ongoing work is to generalize the proposed method in this paper to three-dimensional case, in which the reconstruction of coefficient functions would be a nontrivial issue.
if all components of α are even; and C α = 0 if one of the component of α is odd.

Appendix B. Conversion between wave function and Wigner coefficient function
Given the Wigner function of pure state We have expanded it by series involving Hermite polynomials in Section 4. It follows (29) and (32) that we could define Conversely we have Thus the coefficient functions { f α } is equivalent to {h α } in the sense that one can be expressed by the linear combination of another. Using (B.1), h α could be calculated by the derivative of wave function: On the other hand, it follows (35) and (29) that Therefore with aid of (25), (27) and (28), we could recover the wave function from coefficient functions Taking Ψ(x 1 , x 2 ) = Ψ X (X)Ψ x (x), separating the variables we get − 1 4 ∂ 2 ∂X 2 + 1 4 X 2 Ψ X = E X Ψ X , and − ∂ 2 ∂x 2 + 1 16 The first equation is which is the one-dimensional Schrödinger equation with harmonic oscillator potential. So the ground state wave function is Ψ X = π −1/4 e −X 2 /2 , with energy 2E X = sgn(x)a j+1 x j + a 0 |x| = 0. (C.5) Thus a 0 = 0, matching the coefficients we get −2a 2 + sgn(x)a 1 = 0 (C.6) and a j+2 = sgn(x)a j+1 + ( j/2 + 1/4 − E x ) a j ( j + 1)( j + 2) for j 1. (C.7) For large j = 2k, the recursion formula approximately becomes a j+2 ≈ 1 2 j a j , with the approximate solution a j ≈ C(1/4) j/2 ( j/2)! = C(1/4) k k!