Quantum eigenstates from classical Gibbs distributions

We discuss how the language of wave functions (state vectors) and associated non-commuting Hermitian operators naturally emerges from classical mechanics by applying the inverse Wigner-Weyl transform to the phase space probability distribution and observables. In this language, the Schr\"odinger equation follows from the Liouville equation, with $\hbar$ now a free parameter. Classical stationary distributions can be represented as sums over stationary states with discrete (quantized) energies, where these states directly correspond to quantum eigenstates. We show that this correspondence is particularly pronounced for canonical Gibbs ensembles, where classical eigenstates satisfy an integral eigenvalue equation that reduces to the Schr\"odinger equation in a saddle-point approximation controlled by the inverse temperature. We illustrate this correspondence by showing that some paradigmatic examples such as tunneling, band structures, and quantum eigenstates in chaotic potentials can be reproduced to a surprising precision from a classical Gibbs ensemble, without any reference to quantum mechanics and with all parameters (including $\hbar$) on the order of unity. Interestingly, it is now classical mechanics which allows for apparent negative probabilities to occupy eigenstates, dual to the negative probabilities in Wigner's quasiprobability distribution. These negative probabilities are shown to disappear when allowing sufficient uncertainty in the classical distributions.


Introduction
Quantum mechanics and classical mechanics differ not just in the physics they describe, but also in the mathematical ways they are generally expressed. Following Lagrange and Hamilton, classical mechanics is usually expressed in terms of time-dependent phase space variables, such as either coordinates and momenta or wave amplitudes and phases, and the dynamics is determined by Hamilton's equations of motion. Quantum mechanics in its nonrelativistic formulation is generally expressed in the language of operators acting on states, where canonical variables are replaced by non-commuting operators. The time-dependence is then absorbed in either the states (Schrödinger representation) or the operators (Heisenberg representation). The non-commutativity of the operators leads to the fundamental quantum uncertainty relation set by Planck's constant . Classical uncertainties arise when considering e.g. an ensemble of particles, for which the equation of motion for the probability distribution is given by Liouville's equation. Similar ensembles of quantum particles are described by density matrices satisfying von Neuman's equation of motion.
It is possible to formally consider the limit → 0, where the fundamental quantum uncertainties become small and classical mechanics can be recovered (see e.g. Refs. [1][2][3][4]). One way of observing this emergence is through the Wigner-Weyl phase space language [5][6][7]. This effectively rewrites the quantum equations of motion in terms of classical phase space variables, which we will denote as x and p for concreteness 1 , and von Neumann's equation for the Wigner (quasiprobability) function at → 0 reduces to Liouville's equation for a classical probability distribution.
However, the reverse seems to be much less known: the classical Liouville equation can be rewritten in the language of state vectors (wave functions), and truncating third-and higher-order derivatives in the equations of motion leads to the Schrödinger equation. On the level of equations, the inverse Wigner transform of a probability distribution describing an ensemble of classical particles leads to a Hermitian quasi -density matrix. The Liouville equation now naturally returns von Neumann's equation of motion for this quasi-density matrix, with its eigenstates satisfying the time-dependent Schrödinger equation. This mapping and the subsequent emergence of the Schrödinger equation appears to have been (re)discovered at various occasions [8][9][10][11][12], where the original derivation seems to be first presented by S. Hayakawa, who in turn attributed it to S. Olbert [13]. Still, to the best of our knowledge no detailed discussion of its implications exist. One of the aims of the current paper is also to give a pedagogical motivation for the introduction of state vectors and Hermitian operators in quantum mechanics through a detailed overview of the emergent operator-state representation of classical mechanics.
However, and as should be expected, there are some crucial differences between this representation of classical mechanics and quantum mechanics. First, Planck's constant here appears as a scale that can be freely chosen and can be taken to be arbitrarily small. The proposed mapping to a quasi-density matrix and the operator representation is exact at all values of , but the resulting equations of motion only return the Schrödinger equation provided the initial classical distribution is sufficiently smooth on the scale of . In this sense, an initial condition corresponding to a single phase space point, i.e. fixed initial position and momentum, can never be properly described by the Schrödinger equation. Second, the weights arising in the classical formulation of the density matrix, which would correspond to probabilities in an actual density matrix, can be negative. Negative probabilities similarly appear in Wigner's quasiprobability function, and the negative probabilities here can be seen as dual to the ones in Wigner's function. As also discussed by Feynman, such negative (conditional) probabilities do not pose a problem as long as the probabilities of verifiable physical events remain positive [14]. These negative probabilities are in fact necessary if we want to avoid classical uncertainty relations, as also observed in Ref. [15]. Furthermore, in much the same way that the negative probabilities in the Wigner function disappear upon some averaging over phase space, we show that negative probabilities in classical mechanics disappear upon the introduction of sufficient uncertainty on the phase space distribution, in particular in the energy.
In the second part of the paper we focus on the operator-state representation of stationary distributions. We focus on two such distributions, important in the context of statistical mechanics: microcanonical and canonical. Remarkably, and this is perhaps the main finding of our paper, eigenstates of the canonical (classical) Gibbs ensemble are shown to satisfy an integral equation, which in turn reduces to the stationary Schrödinger equation in the saddlepoint approximation controlled by the inverse temperature β. This result has the advantage that the accuracy of the approximation is set by both and the temperature, where we show that a remarkably accurate correspondence between classical and quantum eigenstates can be obtained even for relatively large values of . We illustrate these results with several paradigmatic examples, where we reproduce nearly-exact quantum wave functions and energy levels from the classical Gibbs ensemble even in the absence of particularly small parameters.
As an initial illustration, in Fig. 1 we present the classical eigenstates for the canonical Gibbs ensemble for a single-particle system in a two-dimensional double-well potential. These are visually indistinguishable from the eigenstates obtained by solving the Schrödinger equation, where we have chosen to present the two lowest-energy states and two higher-excited states. The two lowest states correspond to the expected symmetric and antisymmetric combination of localized states, whereas the other states represent generic excited states. Figure 1: Illustration comparing classical eigenstates obtained from the Gibbs distribution (surface plots) and quantum eigenstates obtained solving the stationary Schrödinger equation (wire mesh) for a two-dimensional potential. Two lowestlying states and two excited states are shown for a two-dimensional potential, with an added shift in the vertical direction in order to visually separate different states. Potential (blue surface plot) corresponding to V (x, y) = ν 4 (1 − x 2 ) 2 + 1 2 mω 2 y 2 with m = ω = ν = 1, inverse temperature β = 0.1, and = = 0.1. See Section 5 for more details. This paper is organized as follows. In Section 2, we provide a short recap of relevant results from the Wigner-Weyl formalism, after which an overview of our results is presented in Section 3, discussing the classical Schrödinger equation, the representation of classical observables as operators, and arguing for the necessity of negative probabilities. The proposed mapping allows general classical dynamics to be mapped to an eigenvalue problem, without any approximation, as shown in Section 4. This also highlights the importance of stationary states, where the quasi-density matrices for stationary microcanonical and canonical distributions are calculated in Section 5. Analytic results for negative eigenvalues of the quasi-density matrix are presented for the linear potential and the harmonic oscillator, and various examples of eigenstates and eigenvalues for the Gibbs ensemble are shown in Section 6. These include the harmonic oscillator, the quartic potential, and the double-well and periodic potential, where it is shown how tunneling states can be obtained, before concluding with some examples of eigenstates of two-dimensional potentials. Section 7 is reserved for conclusions.

Wigner-Weyl formalism
For completeness, we first provide a short overview of the Wigner-Weyl formalism [5][6][7]. Readers familiar with this formalism can immediately skip to Sec. 3.
Given a quantum wave function ψ(x, t), the Wigner function is defined as More generally, for a density matrixρ(t), This function has the important property that it returns the correct marginal distributions for the canonical variables This clearly suggests that W could be interpreted as a joint probability distribution for x and p. However, while this function is real and normalized, it is not necessarily positive. Rather, it belongs to a class of quasiprobability distributions. The equation of motion for the Wigner function is highly similar to the classical Liouville equation, namely where H W ≡ H W (x, p) is the Weyl symbol of the quantum HamiltonianĤ, essentially a correctly-ordered classical Hamiltonian describing the particle, {·, ·} M is the Moyal bracket, and Λ is the symplectic operator, or loosely speaking simply the Poisson bracket operator: In the limit → 0 the sin-function can be linearized, sin( Λ/2) → Λ/2, and von Neumann's equation (5) reduces to the Liouville equation [5,6,16] such that in this limit the function W is conserved on continuous trajectories defined through the classical equations of motion In this formulation of quantum mechanics, the key differences with classical mechanics are (i) the non-positivity of the Wigner function, not allowing a straightforward interpretation of W as a probability distribution, and (ii) with the exception of harmonic systems, the impossibility of representing solutions of Eq. (5) through smooth characteristics x(t) and p(t) on which W is conserved in time. In other words, while the Wigner function itself smoothly changes in time, this change can not be represented by smooth phase space trajectories of the particle, necessitating so-called quantum jumps.

Operator-State representation of the classical Liouville equation
Before presenting the full derivation, we will present some of the results of the following section. Starting from a classical probability distribution P (x, p), a quasi-density matrix W can be defined as where we have introduced a dimensionful parameter that will play the role of . Switching to variables (x 1 , x 2 ) = (x + ξ/2, x − ξ/2), we can write which holds for any choice of parameter , such that both the states and weights have an implicit dependence on . The coefficients w α are real and satisfy α w α = 1 for normalized states ψ α (x). Classical averages over observables can be represented as Hermitian operators acting on the wave functions ψ α (x), returning the coordinate representation of x and p where operator averages are calculated in the usual way. This operator representation always holds, irrespective of the choice of and without any approximation. Making the timedependence explicit, in the limit → 0 (the precise scale determining this limit will be made clear in the derivation), the equation of motion for W can be rewritten to show that all weights w α are time-independent and the states ψ α (x) satisfy the Schrödinger equation Both in quantum and classical mechanics an important role is played by stationary, i.e. time-independent, distributions corresponding to stationary states. In classical mechanics such distributions are usually associated with statistical ensembles, with canonical and microcanonical ensembles predominant in systems with a conserved number of particles. Applying the expansion (10) to a stationary distribution W(x 1 , x 2 ) naturally leads to classical stationary states. A crucial result of this paper is that the stationary states corresponding to a canonical Gibbs ensemble at inverse temperature β have a striking resemblance to quantum stationary states. As shown in Section 5, the exact eigenvalue equation for such canonical eigenstates can be written as which is similarly shown to return the Schrödinger equation when evaluating the integral using a saddle-point approximation controlled by β, with the eigenvalues returning the expected Boltzmann weights. Interestingly, this close correspondence between quantum and classical states remains even in the absence of small parameters (see Sec. 6).
Of course, one immediate difference between quantum and classical mechanics in this language is that the parameter playing the role of is arbitrary. There are, however, other key differences highlighting the complementarity of the quantum and classical formalisms. In particule, the w α can now be interpreted as quasiprobabilities, and it is now classical mechanics that can lead to negative probabilities of occupying states, defined as eigenfunctions of W(x 1 , x 2 ). For the Gibbs ensemble in particular we find that at high temperatures the weights w α are all positive and coincide with the Boltzmann factors w α ∝ e −βEα , forming a broad distribution. In the opposite limit of vanishing temperatures β → ∞, where the differences between quantum and classical states are most pronounced, the distribution of this weights is oscillatory and broad, with w α ∝ (−1) α e −βEα , withβ ∝ 1/β. The distribution is maximally narrow at some specific temperature β −1 on the order of the ground state energy. It is precisely this negativity of probabilities that leads to the violation of the uncertainty principle in classical systems for small temperatures or, more generally, for narrow phase space probability distributions. Despite these subtleties, this "state language" of formulating classical mechanics is very useful as it provides both practical and conceptual tools for understanding the connections and differences between quantum stationary states and classical equilibrium distributions; quantum and classical chaos and integrability, entanglement and more. This connection also can help us to understand in what way some of the postulates of quantum mechanics are simply reformulations of classical results in the state language.

Classical Schrödinger equation.
We start this section by showing how the Liouville equation 2 for the classical probability distribution P (x, p, t) [16] can be rewritten entirely in the coordinate space. This derivation essentially follows that of Ref. [13], but because it is not widely known we will repeat it here for completeness. To simplify the derivation that follows, we consider a Hamiltonian describing a classical particle in a one-dimensional potential where m is the mass of the particle and V (x) is the potential in which it moves. The assumption of a one-dimensional/single-particle system will be unimportant in the following, and this derivation is readily generalized to more general Hamiltonians H(x, p) that are arbitrary analytic functions of the phase space variables. Substituting the Hamiltonian (15) into the Liouville equation we find It is convenient, by taking the Fourier transform with respect to momentum, to go from a representation of P in terms of position and momentum to a representation in terms of a pair of position coordinates In order for ξ to have the dimension of position, we also introduce a dimensionful parameter , which is kept arbitrary for the time being but will end up playing the role of . The reason for choosing W to be a function of x+ξ/2 and x−ξ/2 will be apparent shortly. This construction can be inverted as Note that if W(x 1 , x 2 ) is replaced by the quantum density matrix ρ(x 1 , x 2 ) and by , P (x, p) becomes the corresponding Wigner function [5] (see also Sec. 2). We can formally regard Eq. (17) as the inverse Wigner transform of the classical probability distribution with an arbitrarily chosen Planck's constant.
Taking the Fourier transform of Eq. (16) and using partial integration to evaluate the second term on the right hand side results in the following equation of motion for W: This equation can be made explicitly symmetric by switching to new variables (x 1 , where we also multiplied Eq. (19) by i . Eq. (20) is exact and holds for any value of . Now we can make a crucial simplification and take to be sufficiently small: in this case, we can see from Eq. (17) that in order for W to be nonzero, we also need to consider ξ = (x 1 − x 2 ) sufficiently small compared to . Namely, if P (x, p, t) does not significantly change when varying p over a scale on the order of /ξ, the integral over the exponential will average out to zero, and the only non-zero contributions to W will be for ξ similarly small. We can make an approximation and set ξV ( such that the equation of motion reduces to Similar to regular quantum mechanics, Eq. (21) can be simplified through an eigenvalue decomposition of W(x 1 , x 2 , t), interpreting W as an operator with (x 1 , x 2 ) as indices. Note that this operator is Hermitian: taking its transpose, i.e. exchanging x 1 and x 2 , corresponds to taking ξ → −ξ in Eq. (17), which is in turn equivalent to taking the complex conjugate. As such, it is guaranteed to be diagonalizable, with real eigenvalues. Then where the functions Ψ α (x, t) are the analogues of the time-dependent wave functions in the Schrödinger representation and w α are weights/eigenvalues. While the latter can in principle be time-dependent, it follows from Eq. (21) that they are time-independent, while the orthonormal wave functions satisfy the Schrödinger equation

Representation of observables through Hermitian operators.
The representation of the probability density as a (quasi-)density matrix immediately leads to the representation of observables as Hermitian operators acting on states. This mapping does not depend on any approximations, as we will now explore. For convenience, the time dependence of all distributions, states, and operators is also made implicit in this section.
First of all, we can choose the eigenstates to be normalized as From the normalization of the classical probability distribution, we then find that More generally, expressing the classical probability distribution P (x, p) through the function W, it is straightforward to find that From the first expression for x, one can see that W(x, x) = α w α |Ψ α (x)| 2 plays the role of the coordinate probability distribution. In fact, this correspondence holds for general observables that only depend on position, where for any function O(x) one has From Eq. (27) we see that the momentum is represented by the partial derivative such that the operatorsx (represented by multiplication by x) andp satisfy the usual commutation relation It is easy to check that this representation survives if we consider an expectation value of an arbitrary function O(p). Note that the definition of W in Eq. (17) implied choosing the coordinate representation of wave functions. Alternatively, we could have obtained the momentum representation where the momentum representation of operators would give p → p and x → i ∂ p .
One can similarly analyze more complicated observables involving products of x and p. For example, where we obtained the last equation by integrating by parts. We see that, as perhaps expected, This equation immediately extends to functions of the form pO(x), where It is straightforward to check that more general functions of the form p n O(x) correspond to so-called symmetrically-ordered operators (see e.g. Refs. [6,[17][18][19]), which can be defined by the recursion relation All symmetrically-ordered operators are explicitly Hermitian, which immediately follows from the recursion relation above, combined with the fact thatx andp are Hermitian. As an important example, for a time-independent Hamiltonian the energy of a system is conserved and given by which is independent of the limit → 0 necessary to obtain the Schrödinger equation.

Negative probabilities for classical systems. The uncertainty principle.
Much of the previous section exactly reproduced the language of quantum mechanics entirely within a classical formalism. The only approximation we made was in the equation of motion determing the dynamics, setting when acting on W(x 1 , x 2 ), as required to obtain Eq. (21) from the exact equation (20). This approximation is justified if the distribution P (x, p) is sufficiently smooth on the scale set by and is sufficiently small such that W(x 1 , x 2 , t) is only non-zero when |x 1 − x 2 | = |ξ| is small enough. Note that for harmonic potentials there are no approximations involved and Eq. (21) is exact for any choice of since then ( There seems to be an apparent contradiction with Heisenberg's uncertainty relation δxδp ≥ /2. In quantum mechanics this uncertainty relation is a direct consequence of the commutation relation (30). However, this commutation relation also holds classically, and even more, it holds for any choice of . Yet, the initial distribution P (x, p) can be chosen to be arbitrarily narrow and violate the uncertainty relation. To resolve this apparent paradox, it's necessary to conclude that W cannot be an exact density matrix: for a general distribution the weights w α entering Eq. (22), playing the role of probabilities, will not be positive. This allows us to interpret these weights as quasiprobabilities in the same way the non-positive quantum Wigner function W (x, p, t) is interpreted as a quasiprobability in phase space, since the weights are real and sum to unity, α w α = 1. We thus arrive at the interesting conclusion that, in the operator-state representation, it is now classical mechanics that leads to apparent negative probabilities. If we choose larger than , such that real quantum effects are neglected, but still small enough such that Eq. (20) holds, we can effectively realize density matrices with negative (quasi)-probabilities, which could lead to phenomena not possible in ordinary quantum mechanics.
As an immediate corollary, only classical distributions which satisfy the uncertainty relation δxδp ≥ /2 can have all non-negative weights w α . This can be exemplified from a Gaussian phase space distribution where the qualitative character of the eigenvalues will depend crucially on σ x σ p . First consider the distribution saturating this bound, where W(x 1 , x 2 ) can be easily obtained and shown to factorize as The quasi-density matrix has a single non-zero eigenvalue and the corresponding eigenstate is the ground state of a harmonic oscillator with frequency ω = /(2mσ 2 Now consider the case where σ x σ p > /2. The Wigner function for the Gibbs ensemble of a harmonic oscillator is exactly given by a Gaussian, which we can invert to show that in this case the quasi-density matrix following from Eq. (38) is given by (as shown in Appendix A.2) corresponding to the Gibbs distribution of a harmonic oscillator with frequency ω and ψ n (x) the n-th eigenstate, at inverse temperature β q , where the frequency and the inverse temperature 3 follow from the uncertainties as Since σ x σ p > /2, the second solution has a real and positive solution for β q , such that all weights w n = exp[−nβ q ω]/Z are strictly positive and have a clear interpretation as probabilities.
The final case we can consider is for σ x σ p < /2, i.e. when the distribution P (x, p) does not satisfy the uncertainty relation. This can be obtained by shifting β q in Eq. (42) to where we have used coth(x + iπ/2) = tanh(x). Given a Gaussian with σ x σ p < /2, these equations now have a real and positive solution forβ q . Introducing the same shift of β q in Eq. (42), the eigenstates of the quasi-density matrix violating the uncertainty relation remain those of a harmonic oscillator with frequency ω, but the weights will now become oscillatory and (see also Appendix A.2) Negative probabilities arise the moment the uncertainty relation is violated. Clearly, a narrow distribution P (x, p) with σ x σ p /2 will corresponds to a smallβ q and result in a very broad oscillatory distribution of the weights w n ∝ (−1) exp[−nβ q ω].
As an interesting observation, we note that a partition function of the form Z = 1/(1 − exp(β q ω)) naturally arises in the description of free bosons, whereas that of a free fermion is given by (1+exp(β q ω)). For the harmonic oscillator the eigenstate index n can be interpreted as an occupation number, leading to an average energy of ω( n + 1/2), and it can easily be checked that returning the Bose-Einstein and (minus the) Fermi-Dirac distributions respectively. While the bosonic Bose-Einstein distribution for oscillators with σ x σ p > /2 is not unexpected, the analogy between free fermions and the quantum oscillator at a complex temperature, in turn equivalent to a Gaussian classical probability distribution with σ x σ p < /2, is rather intriguing. At the moment we are not sure if this is a simple coincidence or if there is a deeper underlying reason.

Dynamics
As shown in Section 3.1, the equation of motion for W is exactly given by Introducing a discrete basis of eigenoperators of L, the coupled differential equations of classical mechanics will here lead to a solution of W described by dephasing eigenoperators of the superoperator L, familiar from quantum mechanics.
Denoting the complete set of orthonormal eigenoperators of L as O α (x 1 , x 2 ) with eigenvalues λ α , the classical dynamics is given by where the expansion coefficients of W are given by By making use of the fact that L is antisymmetric under exchange of This can easily be checked in a known limit: assuming the phase space distribution is smooth enough at all times such that we can approximate L[W] = [W,Ĥ], or that the potential is close to harmonic, the eigenoperators of L are simply products of eigenstates of the Hamiltonian. Labeling the eigenstates ofĤ as ψ n (x) with eigenvalue E n , the corresponding eigenoperators of L are given by O nm (x 1 , x 2 ) = ψ m (x 1 )ψ n (x 2 ) with eigenvalues λ nm = E n − E m . These clearly arise in pairs, since λ mn = E m − E n = −λ nm for the eigenoperator ψ m (x 2 )ψ n (x 1 ). Stationary states, which can also be obtained from the long-time average of W(x 1 , x 2 , t), here correspond to the zero-eigenvalue eigenoperators of L and reduce to the diagonal states ψ n (x 1 )ψ n (x 2 ) in this limit.

Stationary States
Both in quantum and classical mechanics a special role is played by stationary states/ stationary probability distributions. In quantum mechanics these states are defined as eigenstates of the Hamiltonian. While all possible stationary probability distributions are not classified in classical systems, an important class of such distributions is those corresponding to equilibrium statistical ensembles. Since we restrict ourselves to single-particle systems in this work we will focus on the two most common ensembles: canonical (fixed temperature) and microcanonical (fixed energy). As we will see the classical-quantum correspondence is most pronounced for canonical ensembles, such that we will focus on these in the following.
Following the previous Section, stationary distributions are necessarily eigenoperators of L with eigenvalue zero. It is possible to expand W in an arbitrary basis as where {ψ α (x)} are some complete set of orthonormal wave functions, which could be e.g. eigenstates of a quantum HamiltonianĤ. Plugging this expansion into Eq. (20), multiplying both parts of this equation by ψ γ (x 1 )ψ * δ (x 2 ) and integrating over x 1 and x 2 , we find the exact equation for the matrix elements of the stationary W αβ L αβ γδ W αβ = 0, where L is the superoperator with entries In the limit of sufficiently small this equation clearly reduces to the matrix form of the stationary von Neumann's equation In general, solutions of Eq. (50) are highly degenerate and the number of solutions generally corresponds to the number of eigenstates of the Hamiltonian, reflecting how each state results in a stationary distribution. However, different stationary distributions W(x 1 , x 2 ) are not expected to commute and different stationary contributions will generally have different eigenstates, such that the set of stationary eigenstates is not uniquely defined.

Microcanonical ensemble.
In one-dimensional systems, assuming that there are no spatially disconnected regions in phase space, any stationary distribution can be represented as a statistical mixture of microcanonical distributions [20]: where ρ(E) is the energy distribution function. The microcanonical distributions are characterized by an equal probability of occupying phase space points on the constant energy surface as Since P mc (x, p; E) is a function of the Hamiltonian of the system H(x, p) it automatically satisfies {H, P } = 0, such that it is necessarily stationary. Such microcanonical distributions naturally arise when considering long-time averages of classical trajectories [20].
In the following, we will analytically present numerical results for systems with a linear potential and a quadratic potential (harmonic oscillator). In these cases the operator L exactly reduces to the commutator with the Hamiltonian for an arbitrary . Therefore all zero-eigenvalue eigenoperators of L necessarily commute with the Hamiltonian, sharing the same eigenstates, which does not hold for more general potentials. Still, it is instructive to consider these simple examples to understand the behavior of the eigenvalues.

Linear potential.
For a linear potential V (x) = αx, the eigenstates of the quantum Hamiltonian can be expressed as Airy functions. As shown in Appendix A, the resulting quasi-density matrix can be obtained by combining two identities for the Airy functions, such that the transform of Eq. (54) for V (x) = αx can be explicitly written in its diagonal form as where we have introduced the customary length scale λ 3 = 2 /(2mα) and the Airy functions are the eigenstates of the Hamiltonian with linear potential 4 . For an eigenstate of the Hamiltonian with energyẼ, the corresponding eigenvalue of the quasi-density matrix of the microcanonical ensemble with energy E is given by Whereas it might be expected that a classical distribution with fixed energy E only contains contributions from quantum states with a similar energyẼ, quite the opposite happens: for a microcanonical state with classical energy E its quasi-density matrix contains contributions from almost all eigenstates of the Hamiltonian with quantum energiesẼ, where the eigenvalue is determined by the Airy function of E −Ẽ. For states withẼ E, the eigenvalue will be exponentially suppressed, whereas for states withẼ E the eigenvalues are highly oscillatory and only decaying as (Ẽ − E) −1/4 . Clearly, a large fraction of the latter eigenstates also have a negative eigenvalue, and the same oscillatory behaviour as for the Gaussian distribution (44) can be observed.
This distinction vanishes if we allow for sufficient uncertainty on the classical energy. Assuming a Gaussian uncertainty on the energy centered on E = 0 in Eq. (53) as the eigenstates of the corresponding W will remain unchanged (they do not explicitly depend on E), whereas the eigenvalues are now given by i.e. the Airy transform of the probability distribution of the microcanonical energy. The resulting eigenvalues are presented in the left panel of Fig. 2 for a Gaussian distribution with different widths. For small σ the resulting distribution still resembles the Airy function, albeit with a quicker decay, but for larger σ αλ all oscillations cancel out and we numerically obtain wẼ = ρ(Ẽ): all negative eigenvalues have effectively been averaged out to zero and the quantum and classical states agree.
Using known properties of the Hermite and Laguerre polynomials, we show in Appendix A.2 that for the harmonic oscillator the transform of the microcanonical ensemble can be expanded as in which L n are the Laguerre polynomials. The eigenstates are again the eigenstates of the quantum Hamiltonian, now labeled with a discrete index n, where the corresponding eigenvalue for the microcanonical ensemble with energy E is given by These exhibit the same qualitative behaviour as for the linear potential: given a microcanonical distribution with energy E, the eigenvalue distribution is peaked at the eigenstate of the Hamiltonian with the same energy. All other eigenvalues have strong positive and negative contributions, either exponentially decaying away from the peak at small n and correspondingly E n = ω(n + 1/2) < E or oscillating at E n > E. Adding some uncertainty on the energy of the microcanonical state, all oscillations cancel out, and we end up with an eigenvalue distribution centered around the corresponding eigenstate, as illustrated in the right panel of Fig. 2. The distribution of eigenvalues is now the Laguerre transform of the energy distribution, and we numerically observe that for sufficiently large uncertainty Note that the eigenvalues w n (E) are highly similar to the Wigner function in the positionmomentum space corresponding to a single n-th level of a quantum harmonic oscillator, as also argued in Appendix A.2.

General potentials.
Before continuing to canonical potentials, we briefly discuss how these previous results extend to more general potentials. Starting from a microcanonical ensemble with fixed energy E and potential V (x), it is straightforward to find the function W mc (x 1 , x 2 ) corresponding to the microcanonical distribution as where p E (x) = 2m(E − V (x)) is the classical momentum of the particle, v E (x) = p E (x)/m is the classical velocity, and θ(E − V (x)) is a step function that guarantees that W mc (x 1 , x 2 ) is nonzero only if the center-of-mass coordinate x = (x 1 + x 2 )/2 belongs to the allowed region. While this is a stationary state for general potentials, this operator is no longer expected to commute with the Hamiltonian, such that its eigenstates do not exactly correspond to the quantum eigenstates. Rather than exactly diagonalizing this operator, we can obtain a connection with the WKB approximation by assuming we are far away from the edges of the classically-allowed region, where E V (x) and consider the behaviour for small ξ = x 1 − x 2 (the region of W that is probed by local operators), and approximate in which we have introduced the classical action S E (x), where the lower limit for the integration can be chosen arbitrarily, and take p E (x) ≈ p E (x 1 )p E (x 2 ) to write, in the classicallyallowed region and for small ξ, where we have introduced orthonormalized states As an interesting observation, we note that the Bohr-Sommerfeld quantization of the action, is here equivalent to demanding that the approximation in (66) does not diverge when x 1 and x 2 are at opposite edges of the classically-allowed regions where both p E (x 1 ) and p E (x 2 ) go to zero. The two-dimensional case is introduced in Appendix B.

Canonical ensemble.
For the Hamiltonian given by Eq. (15) the canonical, or Gibbs or Maxwell-Boltzmann, probability distribution is given by in which the normalization constant or partition function is given by It straightforward to compute the function W(x + ξ/2, x − ξ/2) for this distribution by explicitly taking the Fourier transform according to Eq. (9) as Since this is a stationary state by construction and since L will reduce to the commutator with the Hamiltonian for sufficiently small 2 β, the eigenstates of W(x 1 , x 2 ) should approximate the eigenstates of the Hamiltonian. In the next Section we will show explicit examples of stationary eigenstates for particular potentials obtained in this way. We can advance analytically by deriving the equation for the stationary eigenstates ψ n (x) and quasi-probabilities w n . Writing out the eigenvalue equation and changing the integration variable x 2 to ξ = x 1 − x 2 , we find the following exact integral equation: As in the WKB method, it is convenient to define a complex action S n (x) as ψ n (x) = exp[iS n (x)/ ]. The above integral equation now reads This equation can be simplified in the limit of small , leading to essentially the same analysis as in Sec. 3.1. It is perhaps more interesting to note that the inverse temperature β can serve as an alternative different saddle point parameter, since in the limit β → 0 the first prefactor in the exponential diverges. Performing the Taylor expansion of the integrand in ξ up to the second order and integrating over ξ, which is equivalent to the saddle point approximation, we find Since the left-hand side of this equation is x-independent we must have where we have labelled the constant E n . As is well known from the WKB analysis [21] this equation is exactly identical to the stationary Schrödinger equation for the wave function with eigenvalue E n , We also immediately see that the quasi-probabilities are, up to a prefactor, simply the Boltzmann weights of the discrete energies E n There is a clear connection between the classical canonical distribution and quantum stationary states. If we take = and analyze the eigenstates of W, here the Fourier transform of the Gibbs distribution, the correct "quantum" stationary states are recovered. As we will show in the next Section, these include stationary states in a non-linear potential, tunneling states and random states in chaotic two-dimensional systems. Because the saddle point approximation is justified by the smallness of β and not = , the accuracy of classical eigenstates is independent from the accuracy of the WKB approximation, and as we will demonstrate one can very accurately recover both the ground and excited states. Moreover, as we numerically observe for smooth potentials, the difference between quantum and classical eigenstates remains surprisingly small even if β is on the order of one and only few states are effectively populated. We also note that within the saddle point approximation all probabilities w n are strictly positive -a small enough β such that the saddle point approximation holds implies a smooth enough phase space distribution such that all necessary uncertainty relations hold.

Examples of canonical stationary states
In this section we will analyze stationary states in several characteristic single particle systems increasing their complexity and compare them with corresponding quantum mechanical states. In particular, we will analyze a harmonic oscillator, a quartic potential, a double-well and periodic potential, and finally a two-dimensional non-linear potential.

Harmonic oscillator.
We will start from the harmonic oscillator, where all the eigenstates can be found analytically. The potential energy is then V (x) = 1 2 mω 2 x 2 such that the Boltzmann's distribution reads This is precisely the Gaussian distribution we analyzed earlier (see Eq. (38)) with σ x = 1/ βmω 2 and σ p = m/β such that σ p /σ x = mω and σ x σ p = 1/(βω). From Eq. (42) we see that this distribution maps to the equilibrium canonical ensemble, where ψ n (x) are the eigenstates of the quantum harmonic oscillator with the same parameters and → and for T = 1/β ≥ ω/2 whereas for T < ω/2 As expected, for T ω/2 we have β q ≈ β and in the opposite limit T ω/2 we havẽ β q ≈ 4/(β 2 ω 2 ).
We see that for the harmonic potential the stationary eigenstates ψ n (x) coincide at any temperature β with the eigenstates of a quantum harmonic oscillator if we set → . This also follows from the fact that the canonical distribution is stationary by construction and L is exactly the commutator with the Hamiltonian with = for a harmonic potential, such that W necessarily commutes with the Hamiltonian and they share a common eigenbasis. In this sense there is a precise correspondence between the classical Gibbs ensemble and the quantum Gibbs distribution for T ≥ ω/2 if we identify the quantum inverse temperature with β q according to the relations above. Interestingly, for T < ω/2 we still have the exact correspondence between the classical and quantum probability distributions but the quantum weights w n are no longer positive, with an additional oscillatory dependence on top of the exponential decay.

Quartic potential.
Now we move to a slightly more complicated case of a particle of mass m in a nonlinear quartic potential V (x) = 1 4 νx 4 . The exact quantum ground state energy for this model can be found numerically as This energy also provides a characteristic quantum energy scale for the system. We will now compare some eigenstates obtained from the classical Boltzmann's distribution with the quantum eigenstates at different values of β. For concreteness we will fix all the parameters , ν, m to be unity. In Fig. 3 we show examples of several wave functions describing the ground state, the first excited state, and the fifth excited state of the particle in the quartic potential. All plots illustrate a comparison between the exact quantum wave functions obtained by numerically solving the Schrödinger equation (full black lines) and the classical eigenstates obtained by diagonalizing W(x 1 , x 2 ) corresponding to the Gibbs distribution at three different values of β = 0.1, 0.5, 1. While at β = 1 there are clear visible differences between the classical and quantum states, the agreement between them is still strikingly good given that there is no single small parameter in the problem. The "high-temperature" classical eigenstates at β = 0.1 are visually indistinguishable from the quantum eigenstates. In Table 6.2 we list the energies of the classical eigenstates computed as usually as the expectation values of the Hamiltonian, The second column gives the (numerically) exact quantum-mechanical spectrum and the three following columns describe the energy spectrum following from the eigenstates of the classical Gibbs eigenstates computed at three different temperatures. As with the wave functions, the last column corresponding to β = 0.1 gives nearly exact results with about 0.01% accuracy. The lower temperature spectrum has a larger discrepancy with the quantum spectrum but is still pretty accurate. It is remarkable that for β = 1 even the tenth eigenstate, with energy about 15 times larger than the classical temperature and with a tiny occupation, is still reproduced reasonably well.
In the left panel of Fig. 4 we show the distribution of quasi-probabilities w n for the first twenty classical eigenstates corresponding to the four different temperatures of the Gibbs ensemble. We emphasize that these probabilities, together with the corresponding eigenstates ψ n (x), exactly represent the classical Gibbs ensemble. At inverse temperature β = 0.1 the State # Quantum Gibbs β = 1 Gibbs β = 0.5 Gibbs β = 0.  Table 1: Comparisons of energies of the first 10 eigenstates corresponding to a particle in a quartic potential. The first column corresponds to the exact quantum energies. The next three columns are the energies of the classical eigenstates obtained from the Gibbs distribution at three different temperatures. The parameters are the same as in Fig. 3.
distribution of quasi-probabilities almost exactly matches the quantum Gibbs distribution at the same temperature, as illustrated in the right panel of Fig. 4. While qualitatively the similarities with the quantum ensemble extend all the way to β = 1, one can clearly observe the emergence of negative weights with increasing β. From this plot it is clear that it is possible to generate classical probability distributions dominated by the ground state and yet still have very good agreement with the corresponding quantum ground state. In other words, despite the absence of any intrinsic quantization, classical Gibbs ensembles have excellent discrete representations where a small number of eigenstates can be occupied. As the temperature is further lowered and the classical distribution becomes too narrow, such that it is no longer possible to satisfy the uncertainty principle, then the discrete representation becomes broad again, with many states occupied and weights w n becoming strongly oscillatory (see the data for β = 5 in Fig. 4). This behavior of the quasi-probabilities is qualitatively very similar to those in the microcanonical ensembles and the Gaussian phase space distributions discussed in previous sections. In order to further quantify the behavior of the weights/eigenvalues/quasi-probabilities, we consider For positive eigenvalues, this corresponds to the Rényi entropy, which is bounded from below by zero and S α = 0 only for a factorizable distribution with a single nonzero w n = 1. When allowing for negative w n , this lower bound and the interpretation as entropy vanishes, but it is still instructive to consider how S α changes as β is varied. In order to avoid negative arguments for the logarithm, we consider α integer and even. Note that S 2 can be analytically obtained as The partition function integrals can be explicitly evaluated for the harmonic oscillator to return S 2 = − log(β ω/2) and for the quartic potential to return such that S 2 is positive for β −1 0.485 m 2/3 4/3 ν 1/3 , close to the ground-state energy of Eq. (82). For sufficiently small β the saddle-point approximation holds, such that all weights are positive and S 2 > S 4 > S 6 . Further increasing β, these entropies equal zero around (but not exactly at, see inset) the same value of β ≈ 2.06 4/3 ν 1/3 /m 2/3 . Near this point, the distribution is close to factorizable, with a dominant eigenvalue w 0 ≈ 0.990 and second and third largest eigenvalues 0.108 and −0.081. Further increasing β, all S α become negative with inverted ordering S 2 < S 4 < S 6 , strongly indicating negative eigenvalues. The same behavior would be observed for the harmonic oscillator, where S 2 becomes negative at β −1 = ω/2, where the distribution is exactly factorizable and the uncertainty relation is satisfied. More generally, given a fixed classical distribution it should also be possible to choose in such a way that the discrete representation is optimized, minimizing the number of non-negligible weights, and it would be interesting to identify such a choice that returns the actual Planck's constant as = . 6.3 Double-well potential and periodic potential. Next we consider a particle of mass m in a somewhat more complicated double-well potential, As before, we will choose m = 1 for the numerical analysis, although we now consider a smaller value of = = 0.1 as for larger values of , e.g. = 1, the potential is too weak to support tunneling states. In Fig. 6 we show the classical and quantum wave functions corresponding to the lowest symmetric and antisymmetric tunneling states (top) and to the second pair of symmetric and anti-symmetric tunneling states (bottom) at various values of β. Already for β = 1, the classical and the quantum states are visually indistinguishable.
To quantify the accuracy of the agreement between the quantum and classical tunneling states in Fig. 7, we plot the relative error in the tunneling gap between both pairs of tunneling states as a function of β. This relative error is defined as where the eigenstate index n = 1 corresponds to the two lowest tunneling states and n = 3 describes the second pair of states. The index β implies that the corresponding energies are obtained from the classical Gibbs distribution at the inverse temperature β and the index QM corresponds to the numerically calculated quantum eigenstates of the Hamiltonian. The mistake clearly increases with β and scales as β 2 , but even at β = 1 the relative error is still less than 0.5%, which is surprisingly accurate given that the tunneling splitting itself is a very small fraction of the actual eigenenergies.
If we would interpret this in the language of quantum mechanics, such an accurate representation implies that the lowest symmetric and antisymmetric tunneling states are maximally entangled. In the Fock basis of localized left-and right orbitals, these states read where |0 L,R and |1 L,R are the vacuum (no-particle) and the one-particle states corresponding to the left and right orbitals respectively. Both states are obviously maximally entangled. It is remarkably that such states are built into the classical Gibbs ensemble, and it would be interesting to see how this example generalizes to many identical particles. The number of minima in the potential can be systematically increased to consider, e.g., a periodic lattice. Taking a potential of the form with V 0,c constants, and where the last term represents an overall confining term helping to avoid dealing with boundary conditions and consider x ∈ [−40π, 40π]. In order to enhance dispersion, where tunneling again plays a crucial role, we choose a slightly smaller mass m = 0.5 while keeping = = 1 and β = 0.1. In Fig. 8   6.4 Two-dimensional coupled oscillator.
We now move to two-dimensional systems with a four-dimensional phase space, which we will take as (x, p x , y, p y ). For Hamiltonian evolution with the previously-obtained canonical distribution (71) at inverse temperature β readily extends to with Z = dx dy exp[−βV (x, y)] and where we can now define (x 1 , x 2 ) = (x+ξ x /2, x−ξ x /2) and (y 1 , y 2 ) = (y + ξ y /2, y − ξ y /2).
As an example, we can consider an asymmetric coupled two-dimensional oscillator where the asymmetry is tuned by δ and the coupling by ν. In Fig. 9, we again compare the eigenstates of W(x 1 , y 1 ; x 2 , y 2 ) with those of the quantum mechanical Hamiltonian, where we again choose β = 1 for m = 1 and = = 0.1. The same qualitative behaviour as for one-dimensional systems can be observed, not just at low-lying states, but also at higher excited chaotic states. We compare the ground state (n = 0), the fifth excited state (n = 5) and the 140th excited state (n = 140), where the latter is representative for general higherenergy states. The correspondence between the two is again excellent. Note that, while the effect of the coupling term on the ground and low-lying states might be small, it is crucial in the higher-energy states such as n = 140. The region where the wave function is nonzero corresponds to the classically-allowed region for a particle with a given energy, and the interaction is evidenced by the deformation of the edges of this region, where an ellipse would be obtained for two non-interacting oscillators with ν = 0.

Conclusions and Outlook
We highlighted how key concepts in classical mechanics can be reformulated in the language of Hermitian operators and states more familiar from quantum mechanics. This language naturally follows from applying the inverse Wigner-Weyl transform to the classical probability distribution P (x, p), mapping it to a function W(x 1 , x 2 ) = W * (x 2 , x 1 ) playing the role of the density matrix in the language of quantum mechanics. This function can in turn be diagonalized, with its eigenfunctions ψ n (x) playing a role similar to quantum wave functions. We showed that the correspondence with quantum mechanics is particularly striking if P (x, p) is described by the classical Gibbs ensemble. Then the corresponding classical eigenstates exactly match quantum eigenstates in the limit of high temperature. Surprisingly, this correspondence remains highly accurate even if all the parameters remain of the order of unity. In particular, we were able to accurately describe both ground and excited wave functions in a nonlinear quartic potential, a double-well potential containing tunneling states, a band structure in a periodic one-dimensional potential, and both low-energy states and highly excited states in a two-dimensional nonlinear (chaotic) potential.
Not only do these classical eigenstates correspond to quantum states at high temperatures, the eigenvalues of the quasi-density matrix, commonly interpreted as the probabilities to occupy the eigenstates, return the expected quantum Gibbs/Boltzmann factors. As the temperature is lowered, or more generally as the classical distribution becomes narrower, the classical weights of eigenstates start to acquire negative values and can now be interpreted only as quasi-probabilities. This is exactly dual to the Wigner function, where the phase space distribution following from a quantum density matrix can take non-positive values. Interestingly, there is always a minimum temperature set by the quantum uncertainty relation where the representation of W(x 1 , x 2 ) and hence P (x, p) through the eigenstates becomes 'maximally discrete', i.e. contains the fewest number of non-negligible components. This can be quantified through various measures, where we here consider an extension of Rényi entropies to negative weights. For example, given an oscillator at a temperature T = ω/2 ( is a free parameter in classical mechanics, playing the role of ) only a single (ground) state is occupied. At both higher and lower temperatures an increasing number of states are occupied. Interestingly, the classical partition function of oscillators with T > ω/2 maps to the quantum partition function of bosons with energy scale ω, while the classical partition function of oscillators with T < ω/2 maps exactly to the partition function of free fermions with the same energy scale ω. Whether this is a simple coincidence or if there is a deeper underlying reason remains to be understood.
In this work we focused only on unbounded potentials and hence did not emphasize the role of boundary conditions: obviously both classical and quantum probabilities have to vanish deep in the energetically-forbidden region. It is clear that understanding other, e.g. periodic, boundary conditions should prove to be very interesting. For example, if we consider a particle in a central potential the classical Gibbs probability distribution has to be a periodic function of the angular variable. This implies that the classical states can be both periodic and anti-periodic functions of the angle, with a striking similarity of these two possibilities to integer and half-integer spin. We left analyzing such possibilities to a future work. Similarly, the extension to multiple classical particles and the resulting particle statistics should also be highly interesting. In the presented mapping the number of particles play no role, so at least at high temperatures the many-particle stationary states should satisfy the correct Schrödinger equation. Many more questions such as adiabatic continuation, the manifestation of the discreteness of stationary states, the relaxation of interacting systems to equilibrium, linear response theory,... were left out of the present work, offering various directions for future research. One of the cornerstones of our current understanding of quantum thermalization is given by the eigenstate thermalization hypothesis [22], relating matrix elements of quantum states to thermal distributions, and the connection between quantum states and classical distributions could also be revisited in this context. Classically, stationary distributions distinct from the thermal and microcanonical ones are guaranteed to exist as time-averaged KAM trajectories [23][24][25], and it would similarly be interesting to check the correspondence between the eigenstates following from this classical distribution and the quantum Hamiltonian. From our discussion it should be clear that there is a close connection between the classical Liouville equation and the quantum Schrödinger equation, so it is inevitable that various quantum dynamical phenomena are encoded in the operator-state representation of classical systems.
Finally let us point out that the ideas presented here can go beyond classical mechanics. Given a classical probability distribution in position space, it is always possible to introduce momentum as an auxiliary degree of freedom, as is often done in e.g. annealing problems. The construction shown here can be seen as a way to map such a continuous to a discrete probability distribution, in the same way that stationary quantum mechanics presents a discrete representation of the Gibbs distribution. The parameter (or equivalently ) can then be chosen for convenience, e.g. minimizing the number of discrete components in such a representation.
The eigenvalues and eigenstates of the microcanonical ensemble for the harmonic oscillator can also be analytically obtained, where the eigenstates necessarily correspond to those of the quantum Hamiltonian. These are given by ψ n (x) = 1 and the Wigner function of these states is given by (see e.g. Ref. [28]) P n (x, p) = dξ 2π exp i pξ ψ n (x + ξ/2)ψ n (x − ξ/2) with L n the Laguerre polynomials. These satisfy the orthonormality relation δ (x − y) = e −(x+y)/2 ∞ n=0 L n (x)L n (y), which can be used to express Using the known transform of the oscillator states (100) then leads to and we have that for the harmonic oscillator the transform of the microcanonical ensemble can be expanded as (up to a normalization factor Z = 2π/ω) The eigenvalues are highly similar to the Wigner function, which can be understood by noting that and P n (x, p) is a function of p 2 2m + 1 2 mω 2 x 2 , immediately leading to the correct expression for the eigenvalues.
The relation (42) can similarly be obtained by making use of the known transform of the harmonic oscillator states (100), now using the generating function of the Laguerre polynomials [29] as given by its transform is given by where the summation can be explicitly evaluated from the generating function (106) by taking t = −e −β ω as P (x, p) = 1 πZ where we have switched to polar coordinates p = (p cos φ, p sin φ) in the second line and defined p( r) = 2m(E − V ( r)). Continuing, with J 0 a Bessel function of the first kind. Note that this expression was also obtained in M. Berry's original paper [30], introducing what is now known as Berry's conjecture.