FORCETOC

In last week's lecture, we started our discussion of the density matrix with its three properties
  1. Convolution <math>\rho(x,x', 2\beta) = \int dx'' \rho(x,x'', \beta) \rho(x'', x', \beta)</math>
  2. High-temperature limit (Trotter formula) <math>\rho(x,x',\beta) \to \exp(-\beta/2 V(x)) \rho^{\text{free}}(x,x',\beta)\exp(-\beta/2 V(x'))</math>
  3. Free density matrix (at all temperatures)

These three elements were integrated in the matrix-squaring algorithm. In the present lecture, we take a second look at
matrix squaring (this time from an analytical viewpoint), then discuss the unity between sampling in classical Monte Carlo and
in Quantum Monte Carlo, and finally discuss a direct sampling approach to the free path integral, one of the most beautiful solutions
of a complicated sampling problem.

Density matrix, analytic computations and sampling

For a particle in a harmonic potential, the high-temperature Trotter approximation was already used in
last week's lecture, inside a naive sampling algorithm (algorithm naive-harmonic-path.py). Here we show that the
matrix squaring algorithm can also be used to derive an analytical result, the one for the density matrix of the
harmonic oscillator. We will need the results in this week's DM. See SMAC for a full discussion.

Analytic matrix squaring

Here, we start from the Trotter break-up:

<math>\rho(x,x',\beta) = \exp(-\beta x^2/4) \rho^{\text{free}}(x,x',\beta)\exp(-\beta x'^2/4)</math>.

It is simply rewritten as

<math>\rho(x,x', \beta)= c(\beta) exp\left[ -g(\beta) \frac{(x - x') ^2 }{2} - f(\beta) \frac{(x+x')^2}{2} \right]</math> with <math>c(\beta) = \frac{1}{\sqrt{2 \pi \beta}}</math> and <math>f(\beta)= \beta/4</math> and <math>g(\beta) = \frac{1}{\beta} + \frac{\beta}{4}</math>.

Products and convolutions of Gaussians are again Gaussians, and it follows that the density matrix is at all temperatures of the form

<math>\rho(x,x', 2\beta)= c(2\beta) exp\left[ -g(2\beta) \frac{(x - x') ^2 }{2} - f(2\beta) \frac{(x+x')^2}{2} \right]</math>

The final outcome, valid (and exact) at all temperatures, is

<math> g(\beta) = \frac{1}{2} \coth \frac{\beta}{2}</math>

<math> f(\beta) = \frac{1}{2} \tanh \frac{\beta}{2}</math>

<math> c(\beta) = \sqrt{\frac{1}{2 \pi \sinh \beta} }</math>

see <ref name="SMAC">W. Krauth ''Statistical mechanics: Algorithms and Computations'' (Oxford University Press, 2006) See also [http://www.smac.lps.ens.fr/ the book's website]</ref> SMAC sect. 3.2.2. This fundamental result was obtained by Felix Bloch (in a different way).
The above formula is derived differently form in the books by Feynman
<ref name="Feynman">R. P. Feynman ''Statistical mechanics'' (xxx)</ref>
, and by Landau & Lifshitz <ref name="LandauIII">L. D. Landau and Lifshitz ''Statistical mechanics'' (xxx)</ref>.

Observables, relationship between density matrices and path integrals

In our lectures on quantum mechanics, we discuss very little the computation of thermodynamic obervables. In fact, we mostly look at
snapshots, which correspond to the position operator. The calculation of the energies is subject of this week's homework (see DM 06).

We note that the path-integral representation is mathematically equivalent to the matrix-squaring procedure, but it corresponds to the
conceptual leap, discussed in a previous lecture () from integrals that one computes to integrals that one samples, very much as we have done for hard spheres.

  • A lady on the heliport (see ...).
  • Classical Hard spheres.
  • A quantum particle in a potential.

From naive sampling to the Lévy construction (free particles)


Inefficiency of the naive sampling algorithm


We analyze the naive sampling algorithm for a single particle (see <ref name="SMAC"/>, sect. 3.3.2, and especially fig. 3.12). What we learn in this
simple problem can be generalized to extremely complex problem of interacting quantum particles, the so-called N-body problem.

Lévy construction

The Lévy construction <ref name="Levy">P. Lévy ''Sur certains processus stochastiques homogènes'' [http://www.google.fr/url?sa=t&rct=j&q=p.%20l%C3%A9vy%20sur%20certains%20processus%20stochastiques%20homog%C3%A8nes%20composition%20mathematica%207%2C%20283%20(1940)&source=web&cd=1&ved=0CCQQFjAA&url=http%3A%2F%2Farchive.numdam.org%2Farticle%2FCM_19407283_0.pdf&ei=Qo6ZUOGbK4aZhQe80IGgBA&usg=AFQjCNHbWNnBFxKPVLotGKyjOwzcT_CLbw&cad=rja Composition mathematica 7, 283 (1940)]</ref> yields a rejection-free direct-sampling algorithm for non-interacting quantum particles. In a sense, it represents the exact computational solution of the free-particle propagator. See <ref name="SMAC"/>, sect. 3.3.2 for a detailed discussion. Note that the Lévy construction is also referred to as ''stochastic interpolation'' or as a ''Gaussian bridge''.

##========+=========+=========+=========+=========+=========+=========+=
    1. PROGRAM : levy.py
    2. PURPOSE : implement SMAC algorithm 3.5 (levy-free-path)
    3. LANGUAGE: Python 2.6
##========+=========+=========+=========+=========+=========+=========+
import numpy, pylab,time
from math import pi, exp, sqrt
from random import gauss
N = 10000
beta = 1.
Del_tau = beta / N
x = [0 for k in range(N + 1)]
y = range(N + 1)
for k in range(1,N):
Del_p = (N - k) * Del_tau
x_mean = (Del_p*x[k - 1] + Del_tau*x[N]) / (Del_tau + Del_p)
sigma = sqrt(1. / (1. / Del_tau + 1. / Del_p))
x[k] = gauss(x_mean,sigma)
pylab.figure(1,figsize=(8.,12.))
pylab.plot(x,y,'b-')
pylab.show()

This algorithm can be generalized for the harmonic oscillator (a particle in a harmonic trap) because, in that case also, the starting
density matrix (coming out of the Trotter formula) is made up of Gaussians only.
The approach can be generalized to free particles with hard walls see TD06 and, in combination with the Metropolis algorithm, the Lévy construction is an extremely efficient tool for the simulation of general quantum systems (and also, to simply understand Quantum Mechanics).

References