Quantum Monte Carlo in Classical Phase Space with the Wigner-Kirkwood Commutation Function. Results for the Saturation Liquid Density of $^4$He

This paper presents a Metropolis Monte Carlo algorithm capable of handling complex phase space weights in quantum statistical mechanics and demonstrates its accuracy by successfully calculating the saturation liquid density of $^4$He near the $\lambda$-transition using a third-order Wigner-Kirkwood expansion.

The Gist

Imagine you are trying to simulate a crowded dance floor where the dancers are tiny, invisible particles called atoms. In the "classical" world (like normal people dancing), you can predict exactly where everyone will be and how fast they are moving. But in the quantum world (where these atoms actually live), things get weird: the dancers are fuzzy, they can be in two places at once, and they don't like to be too close to each other because of a fundamental rule of the universe called the Heisenberg Uncertainty Principle.

This paper is about a new way to simulate these quantum dancers using a computer, specifically for Helium-4 (a type of helium gas that becomes a super-fluid liquid at very cold temperatures).

Here is the breakdown of what the author, Phil Attard, did and found:

1. The Problem: The "Fuzzy" Dance Floor

For a long time, simulating quantum particles was like trying to film a dance floor in slow motion by taking thousands of pictures of every single step. It was incredibly expensive and slow.

• The Old Way: One famous method (by Ceperley) treated the particles as if they were walking through time, taking many tiny steps. It was accurate but required a supercomputer to simulate just 64 atoms.
• The New Approach: Attard developed a way to simulate these particles on a "classical" dance floor (where positions and speeds are clear) but adds a special "ghost" rule to account for the quantum fuzziness. This allowed him to simulate 5,000 atoms on a regular personal computer.

2. The Secret Sauce: The "Commutation Function"

The main trick in this paper is a mathematical tool called the Wigner-Kirkwood commutation function.

• The Analogy: Imagine the classical dance floor has a rule that says, "If you get too close to your neighbor, you must pay a fine." In the quantum world, this "fine" isn't just a number; it's a complex, wavy rule that makes the particles act "fuzzier" and keeps them further apart than they would be in a normal crowd.
• The Innovation: Attard didn't just use a simple rule; he expanded this rule into a series of steps (like a recipe with ingredients). He tested the recipe using the first, second, and third ingredients (orders of the expansion).
  • Order 0 (No quantum rules): The atoms clump together too tightly. The liquid is way too dense (about 3 times denser than real life).
  • Order 2 (Adding some quantum rules): The atoms spread out a bit. The density drops by half, getting closer to reality.
  • Order 3 (The full recipe): The atoms spread out just right. The simulated density matches the measured density of real liquid helium almost perfectly.

3. The Results: A Perfect Match

The paper reports that by using this "third-order" recipe, the computer simulation of 5,000 helium atoms created a liquid droplet that was the exact same density as real liquid helium found in nature.

• Why this matters: Before this, if you tried to simulate a large, uniform block of liquid helium on a computer, it would fall apart (cavitate) because the atoms were too crowded. By adding these quantum "fuzziness" rules, the simulation stays stable at the real density, which is a huge achievement.

4. What Happened to the "Symmetrization"?

In quantum mechanics, identical particles (like helium atoms) are so alike that swapping them doesn't change anything. This is called "symmetrization."

• The Paper's Stance: The author admits he did not include this specific rule in this particular simulation. He focused entirely on the "fuzziness" (the commutation function) because it was the main cause of the density error. He says, "I'll tackle the swapping rule in my next paper." He argues that for the temperatures he studied (near the transition point), the fuzziness was the most important factor to get right first.

5. A Few Glitches and Limits

• The "Hard Core": Sometimes, the math got so wild that the computer thought two atoms were on top of each other (which is impossible). To fix this, the author put in a "hard core" rule: "If atoms get closer than X distance, the computer rejects the move." This kept the simulation from crashing.
• The "Solid-Like" Droplet: At the coldest temperatures tested, the liquid droplet in the simulation started to look a bit like a solid crystal (the atoms lined up in rows). The author notes this might be an artifact of the simulation setup (like the walls of the container or the size of the droplet) rather than real helium, which stays liquid even at absolute zero unless squeezed hard.

Summary

Phil Attard created a new, faster way to simulate quantum liquids on a regular computer. By adding a specific mathematical "fuzziness" rule (the third-order Wigner-Kirkwood expansion), he managed to make a virtual bottle of liquid helium that is just as dense as real liquid helium. This proves that you don't always need a supercomputer to simulate quantum matter; you just need the right mathematical recipe.

Technical Summary

Technical Summary: Quantum Monte Carlo in Classical Phase Space with the Wigner-Kirkwood Commutation Function

Problem Statement
The simulation of quantum condensed matter remains a significant computational challenge. While Feynman path integral methods (e.g., Ceperley, 1995) effectively handle the non-commutativity of position and momentum operators (the Wigner-Kirkwood commutation function), they often require substantial computational resources, such as thousands of temperature slices, limiting system sizes. Conversely, classical phase space algorithms handle wave function symmetrization efficiently but traditionally struggle with non-commutativity. The specific challenge addressed here is improving the classical phase space algorithm to accurately account for the Wigner-Kirkwood commutation function without incurring the computational burden of path integrals, specifically for liquid $^4$He near the $\lambda$-transition.

Methodology
The paper presents a Metropolis Monte Carlo algorithm designed for complex phase space weights, applicable to general quantum statistical mechanics. The core of the method involves a subsystem of $N$ identical bosons in a volume $V$, where the phase space probability density $\wp(\Gamma)$ includes a classical Hamiltonian term, a symmetrization function, and the Wigner-Kirkwood commutation function $e^{W(\Gamma)}$.

1. Formalism: The commutation function $W(\Gamma)$ is defined via the relation $e^{-\beta H(\Gamma)}e^{W(\Gamma)} = e^{p \cdot q / i\hbar} e^{-\beta \hat{H}(q)} e^{-p \cdot q / i\hbar}$. Since $W(\Gamma)$ is complex ($W = W_r + iW_i$), the algorithm focuses on the real part for averaging, noting that the imaginary part averages to zero in equilibrium. To prevent rapid oscillations that cancel out, a restriction $|W_i(\Gamma)| < \pi/2$ is imposed, effectively preventing particles from approaching too closely, consistent with the Heisenberg uncertainty relation.
2. Algorithm: A trial move involves changing both position and momentum for one particle. The move is accepted if the ratio of new to old weights satisfies the Metropolis criterion:
   $$\frac{e^{-\beta \Delta H(\Gamma)} e^{\Delta W_r(\Gamma)} \cos W_i(\Gamma_{new})}{\cos W_i(\Gamma_{old})} \ge r$$
   where $r$ is a random number in $[0,1]$. Moves violating the $|W_i| < \pi/2$ constraint are rejected.
3. Expansion: The commutation function is expanded in powers of the inverse temperature $\beta$. The paper derives terms up to the fourth order but reports numerical results using expansions terminated at the third order ($n_{max}^W = 3$). The expansion coefficients involve gradients of the pair potential.
4. Simulation Details: Simulations were performed on Lennard-Jones $^4$He using 1,000 and 5,000 atoms. A hard-core cutoff was implemented to prevent numerical overflow in the repulsive core region. The study focused on the saturation curve, allowing a liquid droplet to condense within a vapor phase under periodic boundary conditions.

Key Contributions

• Generalized Algorithm: The derivation of a Metropolis algorithm capable of handling complex phase space weights, specifically addressing the Wigner-Kirkwood commutation function within a classical phase space framework.
• Higher-Order Expansion: The provision of a fourth-order temperature expansion for the commutation function, with numerical implementation up to the third order.
• Handling of Constraints: The implementation of a specific rejection criterion based on the imaginary part of the commutation function to maintain numerical stability and physical consistency (preventing particles from violating the uncertainty principle).

Results
Simulations were conducted for Lennard-Jones $^4$He near the $\lambda$-transition ($k_B T / \epsilon \approx 0.45 - 0.8$).

• Density: Classical simulations ($n_{max}^W = 0$) yielded a saturated liquid density approximately three times the measured experimental value. Including the second-order commutation term reduced the density by a factor of two. Including the third-order term ($n_{max}^W = 3$) brought the simulated saturation liquid density into agreement with measured values.
• Kinetic Energy: In the classical and second-order cases, the kinetic energy per particle remained at the classical value of $3k_B T/2$. In the third-order case, the kinetic energy was reduced below this value, with the reduction increasing as temperature decreased.
• Structure: The radial distribution function $g(r)$ showed that as the order of the commutation expansion increased, the peak of the distribution shifted to larger separations. This shift is attributed to the delocalization of bound particles due to the Heisenberg uncertainty relation.
• Phase Behavior: At $k_B T / \epsilon = 0.5$, the system exhibited liquid-like behavior. At $k_B T / \epsilon = 0.45$, the density profile suggested a solid-like state, though the authors note this may be an artifact of the Lennard-Jones potential, the truncation of the expansion, or the Laplace pressure in the nano-droplet, as real $^4$He does not solidify at saturation pressure down to absolute zero.
• Validation: The consistency between the energy derivative and the directly simulated heat capacity served as a sensitive test of the mathematical expressions and their implementation, confirming the correctness of the algorithm after extensive debugging.

Significance
The paper claims that the third-order approximation for the Wigner-Kirkwood commutation function enables the classical phase space simulation of Lennard-Jones $^4$He at approximately the measured saturation liquid density within a homogeneous system. Without this commutation function, the system would cavitate into a liquid drop at the much higher bare Lennard-Jones density ($\rho_{sat}^{LJ}\sigma^3 \approx 0.9$). The ability to simulate the system at the experimentally observed density represents a significant advance in the efficiency and accuracy of quantum Monte Carlo methods for condensed matter, preserving the advantages of classical phase space algorithms (such as handling symmetrization) while correcting for non-commutativity. The authors note that methods for including the symmetrization function are available but were not implemented in this specific study, which focused on delineating the role of the commutation function.