Semiclassical Approximation of the Wigner Function for the Canonical Ensemble

Abstract

The Weyl–Wigner representation of quantum mechanics allows one to map the density operator in a function in phase space—the Wigner function—which acts like a probability distribution. In the context of statistical mechanics, this mapping makes the transition from the classical to the quantum regimes very clear, because the thermal Wigner function tends to the Boltzmann distribution in the high temperature limit. We approximate this quantum phase space representation of the canonical density operator for general temperatures in terms of classical trajectories, which are obtained through a Wick rotation of the semiclassical approximation for the Weyl propagator. A numerical scheme which allows us to apply the approximation for a broad class of systems is also developed. The approximation is assessed by testing it against systems with one and two degrees of freedom, which shows that, for a considerable range of parameters, the thermodynamic averages are well reproduced.

Appendices

Appendix A: Wigner symbol of normal forms

In order to calculate the Wigner symbol of an operator of the form (47), it is sufficient to do so for the monomials ${\hat{o}}^n$, where ${\hat{o}} = {\hat{p}}^2+{\hat{q}}^2$. Our strategy will consist in finding a recurrence relation that allows us to calculate $o^{n+1}(p,q)$ in terms of $o^n(p,q)$. As the initial term $o^1(p,q) = o(p,q) = p^2+q^2$ is readily obtained, the problem is solved.

For that, we first observe that, as ${\hat{o}}^n$ is hermitian, $o^n(p,q)$ must be real. Using this fact, writing ${\hat{o}}^{n+1} = {\hat{o}} {\hat{o}}^{n}$, and applying Groenewold’s rule (11), we arrive at the recurrence relation

$$\begin{aligned} o^{n+1}(p,q) = \left( p^2+q^2 - \frac{\hbar ^2}{4}\nabla ^2 \right) o^n(p,q) \end{aligned}$$

(A1)

where $\nabla ^2 = \partial _p^2 + \partial _q^2$. This relation is further simplified if we introduce the coordinates $s,\phi$, defined by $p = \sqrt{s} \cos \phi , \ q = \sqrt{s} \sin \phi$, in terms of which the laplacian takes the form

$$\begin{aligned} \nabla ^2 = 4 \left( s \partial _s^2 + \partial _s \right) + \frac{1}{s}\partial _\phi ^2, \end{aligned}$$

(A2)

which allows us to rewrite (A1) as

$$\begin{aligned} o^{n+1}(s,\phi ) = \left[ s - \hbar ^2\left( s \partial _s^2 + \partial _s +\frac{1}{4s}\partial _\phi ^2\right) \right] o^n(s,\phi ) \end{aligned}$$

(A3)

Since $\partial _\phi o(s,\phi ) = 0$, and, as deduced from the recurrence relation, $\partial _\phi o^n(s,\phi ) = 0 \Rightarrow \partial _\phi o^{n+1}(s,\phi ) = 0$, we prove by induction that $\partial _\phi o^n(s,\phi ) = 0 \ \forall \ n$, which eliminates the derivative with respect to $\phi$ from (A3). This allows us to easily obtain the first terms in the recurrence relation, which, already expressed in terms of p, q, are given by

$$\begin{aligned} \begin{aligned} o^2(p,q)&= \left( p^2 + q^2 \right) ^2 - \hbar ^2 \\ o^3(p,q)&= \left( p^2 + q^2 \right) ^3 - 5\hbar ^2\left( p^2 + q^2 \right) \\ o^4(p,q)&= \left( p^2 + q^2 \right) ^4 - 14\hbar ^2\left( p^2 + q^2 \right) ^2 + 5 \hbar ^4 \\ \end{aligned} \end{aligned}$$

(A4)

We see that, in general, ${\hat{o}}^n(p,q)$ is a polynomial of order n in $(p^2+q^2)$, whose dominant term is $(p^2+q^2)^n$, while corrections proportional to even powers of $\hbar$ are also present.

Appendix B: Numerical Details

The calculations in this article were performed using the Julia language [50]. The package DifferentialEquations.jl [51] was used to solve the necessary differential equations in parallel. The calculations were performed on a 12th Gen Intel Core i5-12600K processor, which has 16 threads.

1.1 Morse System

The integrals related to the Morse system were performed using Gaussian quadrature. The integration region, in units of $\omega = \hbar = 1$, is given by

$$\begin{aligned} R = \left\{ (p,q) \in {\mathbb {R}}^2 \ \left| \ \chi p^2 + \frac{1}{4\chi } \left( 1-e^{-q} \right) ^2 < \frac{1}{4\chi } \right. \right\} \end{aligned}$$

(B5)

Introducing the variables ${\tilde{P}} = 2\chi p$ e $Q = 1-e^{-q}$, we obtain

$$\begin{aligned} \begin{aligned} R&= \left\{ \left( {\tilde{P}},Q\right) \in {\mathbb {R}}^2 \ \left| \ {\tilde{P}}^2 + Q^2 < 1 \right. \right\} \\ {}&= \left\{ \left( {\tilde{P}},Q\right) \in {\mathbb {R}}^2 \ \left| \ Q \in \left( -1,1 \right) ; \ {\tilde{P}} \in \left( -\sqrt{1-Q^2},\sqrt{1-Q^2} \right) \right. \right\} , \end{aligned} \end{aligned}$$

(B6)

which can be simplified by defining $P = {\tilde{P}}/\sqrt{1-Q^2}$. Then, we have

$$\begin{aligned} R = \left\{ \left( P,Q\right) \in {\mathbb {R}}^2 \ \left| \ Q \in \left( -1,1 \right) ; \ P \in \left( -1,1 \right) \right. \right\} . \end{aligned}$$

(B7)

The inverse transformation is then

$$\begin{aligned} {\left\{ \begin{array}{ll} p = \sqrt{1-Q^2}\dfrac{P}{2\chi } \\ q = -\ln \left( 1-Q \right) \end{array}\right. }, \end{aligned}$$

(B8)

which has jacobian determinant

$$\begin{aligned} \det \frac{\partial (p,q)}{\partial (P,Q)} = \frac{1}{2\chi } \sqrt{\frac{1+Q}{1-Q}}, \end{aligned}$$

(B9)

which is proportional to the weight function of a Gauss-Chebyshev quadrature of the 3$^\circ$ kind. We therefore use this quadrature rule to perform the integration over the Q coordinate, while a Gauss-Legendre quadrature is used to integrate over P. The advantage of Gaussian quadrature is that the integration points will be independent of $\theta$, and then a single set of points can be used to compute the thermodynamic quantities over a range of temperatures. In this work, we used a grid of $300 \times 300$ points to perform the integration, which corresponds to $9 \times 10^4$ trajectories.

1.2 Nelson System

In the semiclassical calculations for the Nelson system, different techniques were used for different set of parameters.

In the case of the energy, as well as the heat with $\mu = 1.5, 2$, we first performed the change of variables $(p_x,p_y,x,y) \mapsto (P_X,P_y,X,Y)$ with

$$\begin{aligned} {\left\{ \begin{array}{ll} P_x = \sqrt{\dfrac{\theta }{2}} p_x \\ P_y = \sqrt{\dfrac{\theta }{2}} p_y \\ X = \sqrt{\theta \mu } x \\ Y = \sqrt{\theta } (y - x^2/2) \end{array}\right. }. \end{aligned}$$

(B10)

This transformation has constant jacobian determinant and, in terms of the new variables, we have that the classical Boltzmann’s weight is simply

$$\begin{aligned} e^{-\beta H} = \exp \left[ - \left( P_x^2 + P_y^2 + X^2 + Y^2 \right) \right] . \end{aligned}$$

(B11)

The integration is then performed by an h-adaptive technique as described in [52, 53]. The Julia implementation can be found in [54]. We bounded the integration algorithm to use roughly $10^5$ integration points. We used the BS3 [55, 56] and Vern6 [56, 57] algorithms to solve the differential equations, and the tolerances varied between $10^{-2}$ and $10^{-6}$. For each $\mu$, the corresponding plot took around 40 seconds to 3 minutes to complete.

We found that the heat capacities with $\mu = 0.5, 1$ were much harder to integrate. In this case, we didn’t perform a change of variables and resorted to a Monte Carlo integration method, where the $10^7$ integration points were sampled from the classical Boltzmann’s distribution $e^{-\beta H({\varvec{x}})}/Z$ using the Metropolis-Hastings algorithm [58, 59]. In this case, for each $\mu$, the corresponding plot took around 3 hours to complete.

For the energy spectrum, which is used to calculate the quantum versions of the thermodynamic quantities, we used a grid of $160 \times 160$ points, where x spanned from $-4.5$ to 4.5, and y spanned from $-4$ to 5. We then approximated the laplacian of the time independent Schrödinger equation through a finite differences matrix over this grid. The discretizatation of this equation gives rise to a eigenvalue equation, which can be solve through standard linear algebra libraries.