Path-Integral Monte Carlo

Path-integral Monte Carlo (PIMC) samples the imaginary-time paths that represent a thermal density matrix. For positive path weights, equilibrium quantum averages become classical expectations over an enlarged configuration space. The extra beads resolve quantum fluctuations; they also introduce strongly correlated sampling modes.

Required background. Path-integral discretization defines the action and slice error; MCMC defines acceptance and equilibration.

Sampling the closed-path weight

For one particle in $d$ dimensions, use $M$ beads, $\epsilon=\beta/M$, periodic bead index $x_M=x_0$, and

$$S_M=\sum_j\left[ \frac{m|x_{j+1}-x_j|^2}{2\hbar^2\epsilon} +\epsilon V(x_j)\right].$$

The target is proportional to $e^{-S_M}$. A symmetric single-bead displacement needs only the two adjacent springs and the local potential to compute $\Delta S_M$; accept with $\min(1,e^{-\Delta S_M})$. This locality makes the update cheap, but as $\epsilon$ decreases the springs stiffen and long-wavelength path motion becomes slow.

A simultaneous translation of all beads leaves spring terms unchanged and moves the centroid. Staging, bisection, or normal-mode moves additionally sample correlated internal displacements. If a proposal draws from the exact free conditional bridge, its kinetic density cancels from the Metropolis–Hastings ratio, leaving the potential-action change. That cancellation requires the actual forward and reverse proposal probabilities.

Diagonal observables and thermal energy

For a position-diagonal observable, an unbiased finite-$M$ estimator is the bead average,

$$\widehat O=\frac1M\sum_j O(x_j).$$

Different beads on the same path are correlated; they are not $M$ independent Monte Carlo samples. Block whole-path measurements over Monte Carlo time.

Differentiating the complete finite-$M$ partition function at fixed $M$ gives the thermodynamic energy estimator,

$$E_{\rm th}= \frac{dM}{2\beta} -\frac{mM}{2\hbar^2\beta^2}\sum_j|x_{j+1}-x_j|^2 +\frac1M\sum_jV(x_j).$$

The first two terms nearly cancel and their fluctuations can become large at fine discretization. Virial or centroid-virial estimators reorganize this cancellation and can reduce variance under their differentiability assumptions. Compare estimators before investing only in longer runs.

Oscillator verification at two levels

For $m=\hbar=\omega=1$, the exact thermal energy is $\frac12\coth(\beta/2)$, and $\langle x^2\rangle$ has the same numerical value. At a finite number of beads, compare the sampled bead variance first with the finite-$M$ Gaussian covariance from discretization. Agreement tests the sampler without requiring the continuum approximation to be good.

Then increase $M$ at fixed $\beta$ and compare the extrapolated result with the exact thermal expression. Finally vary $\beta$ if ground-state behavior is the intended target. These tests isolate sampling, Trotter, and temperature errors.

A centroid-only move samples the overall position but leaves all relative bead coordinates fixed. A single-bead-only algorithm can explore all modes but slowly. Mixing both local and collective moves is useful only if their stationary weights are each correct.

Exchange and other limitations

For indistinguishable bosons, path closure includes particle permutations. A sampler that never changes permutation sector misses exchange even if every bead moves. Worldline methods address sector-changing updates and winding observables.

For fermions, permutation parity introduces signed weights; the sign problem is a statistical obstruction, not a small change to a positive sampler. Real-time path weights are generally oscillatory as well. Imaginary-time data do not directly provide a stable real-frequency spectrum without an additional inverse problem.

Exercises

Derive why translating every bead by the same vector leaves the spring action unchanged.

Solution

Each difference becomes $(x_{j+1}+a)-(x_j+a)=x_{j+1}-x_j$. Only potential terms change. Therefore a symmetric centroid-shift proposal has acceptance determined by $\epsilon\sum_j[V(x_j+a)-V(x_j)]$. This move changes the zero Fourier mode and leaves all nonzero modes unchanged.

References

• Ceperley, David M., “Path integrals in the theory of condensed helium,” Reviews of Modern Physics 67, 279–355 (1995), doi:10.1103/RevModPhys.67.279 — density matrices, path sampling, and estimators.
• Hastings, W. K., “Monte Carlo sampling methods using Markov chains and their applications,” Biometrika 57, 97–109 (1970), doi:10.1093/biomet/57.1.97 — stationary sampling and proposal correction.