06_Quantum_Levy

In last week's lecture, we started our discussion of the density matrix with its three properties
 * 1) Convolution $$\rho(x,x', 2\beta) = \int dx \rho(x,x, \beta) \rho(x'', x', \beta)$$
 * 2) High-temperature limit (Trotter formula) $$\rho(x,x',\beta) \to \exp(-\beta/2 V(x)) \rho^{\text{free}}(x,x',\beta)\exp(-\beta/2 V(x'))$$
 * 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:

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

It is simply rewritten as

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

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

$$\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]$$

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

$$ g(\beta) = \frac{1}{2} \coth \frac{\beta}{2}$$

$$ f(\beta) = \frac{1}{2} \tanh \frac{\beta}{2}$$

$$ c(\beta) = \sqrt{\frac{1}{2 \pi \sinh \beta} }$$

see 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 , and by Landau & Lifshitz.

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, 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 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, 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.

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
 * 1) PROGRAM : levy.py
 * 2) PURPOSE : implement SMAC algorithm 3.5 (levy-free-path)
 * 3) LANGUAGE: Python 2.6
 * 1) LANGUAGE: Python 2.6

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=