05_density_matrix_path_integral

=Wave functions= Quantum systems are described by wave functions and eigenvalues, which are solutionsof the Schrödinger equation. In this lecture, we study for the case of the one-dimensional harmonic oscillator how exactly the probabilities of quantum physics combine with the statistical weights in the Boltzmann distribution, before moving on to more general problems.

Harmonic wave functions
As an example of wave functions of a Schrödinger equation, math H \psi_n = \left( -\frac{\hbar^2}{2 m} \frac{\partial^2}{\partial x^2} + V \right) \psi_n = E_n \psi_n math we consider the harmonic oscillator, a particle of mass //m// in a potential math V(x)= \frac{1}{2} m \omega^2 x^2 math

This program uses a dictionary (hash table) for the wave functions code format="python" import numpy, pylab from math import pi, exp, sqrt

grid = numpy.linspace(-5.0, 5.0, 500) psi = {} pylab.figure(1, figsize=(6.0, 10.0)) for x in grid: psi[x] = [pi ** ( - 0.25) * exp( - x **2 / 2.0)] psi[x].append(sqrt(2.0) * x * psi[x][0]) for n in range(2,50): psi[x].append(sqrt(2.0 / n) * x * psi[x][n - 1] - sqrt((n - 1.0) / n) * psi[x][n - 2]) for k in range(40): y = [psi[x][k] + k for x in grid] pylab.plot(grid,y,'r-') pylab.xlabel('$x$') pylab.ylabel('$\psi_n(x)$ (offset)') pylab.xlim(-5.0, 5.0) pylab.show code

=Density Matrix=

Definition and properties
For simplicity, we use in the following the position representation of operators and wave functions. Given a quantum system with an orthonormal set of wave functions and with energies, its density matrix at inverse temperature beta is defined as math \rho(x,x',\beta)=\sum_k e^{-\beta E_k}\psi_k(x)\psi^*_k(x'). math The density matrix is where statistical mechanics (the Boltzmann factor) meets quantum mechanics (the quantum-mechanical probability of the squared wave function).

The expectation value of an observable //O// is given by math \textrm{Tr}(\rho(\beta)O). math For example, the probability of finding the particle at position //x// is math \pi(x)=\textrm{Tr}(\rho(\beta)|x\rangle \langle x|)=\sum_k e^{-\beta E_k}|\psi_k(x)|^2 math

The density matrix has three important properties:

Convolution
Given two inverse temperatures $$\beta$$ and $$\beta'$$, math \rho(\beta)\rho(\beta')=\rho(\beta+\beta'). math This follows from the fact that the eigenstates constitute an orthonormal basis of the Hilbert space. In the position representation, and in a simplified case, we have math \rho(x, x', 2\beta) = \int dx \rho(x, x, \beta)\rho(x'', x', \beta) math

This exact property of the density matrix will allow us to compute the density matrix at low temperature (high \beta), when we know it at high temperature (low \beta ) only.

High-temperature limit
In the limit of high temperature, the density matrix of a quantum system described by a Hamiltonian math H= H^{\text{free}} + V math is given by a general expression known as the Trotter formula: math \rho(x,x', \beta) \xrightarrow[\beta \rightarrow 0]{} \exp (- \beta V(x)/2) \rho^{\text{free}} (x,x',\beta)\exp (- \beta V(x')/2) math This expression is very easy to prove by simply expanding all the exponentials (taking into account non-commutativity of quantum operators).

Density matrix for a free particle
The density matrix for a free particle is given by math \rho^{\text{free}}(x,x',\beta) = \sqrt{\frac{m}{2 \pi \hbar^2 \beta}} \exp \left[ - \frac{m (x-x')^2}{2 \hbar^2 \beta} \right] math

As we usually work in units with math m = \hbar = 1 math this means math \rho^{\text{free}}(x,x\beta) = \sqrt{\frac{1}{2 \pi \beta}} \exp \left[ - \frac{ (x-x')^2}{2 \beta} \right] math

Matrix square
In the following Python program, we implement the matrix squaring algorithm math \rho(x,x',2\beta) = \int dx \rho(x,x,\beta) \rho(x'',x',\beta) math for the harmonic oscillator, with potential //V//(//x//)=//x//^2/2. At high temperature, we initialize the density matrix following the Trotter formula, with half of the potential at //x// and the other half at //x'//.

code format="python" import numpy, pylab from math import pi, exp, sqrt, sinh, tanh

beta = 1.0 / 64 N = 100 grid = numpy.linspace(-5.0, 5.0, N) rho = numpy.zeros(shape=(N,N)) del_x = grid[1] - grid[0] for k in range(N): x = grid[k] for j in range(N): xp = grid[j] rho[k][j] = 1.0 / sqrt(2.0 * pi * beta) * exp(- beta / 4.0 * (x ** 2 + xp ** 2) -                   (x - xp) ** 2 / 2.0 / beta) for k in range(6): beta *= 2 rho = del_x*numpy.dot(rho, rho) y_arrow = [rho[k][k] for k in range(N)] pylab.plot(grid, y_arrow, 'ro') y2_arrow = [1.0 / sqrt(2.0 * pi * sinh(beta)) * exp( -x ** 2 * tanh(beta / 2.0)) for x in grid] pylab.title('harm. diag. density matrix at beta=1 matrix squaring, exact sol.') pylab.plot(grid, y2_arrow, 'b-') pylab.show code

700px|left|Comparing the exact solution for the harmonic oscillator density matrix to the solution obtained through matrix squaring.

=Path Integral=

code format="python" import pylab from math import exp from random import random, uniform, randint

def rho_free(x, xp, beta): return exp(-(x - xp) ** 2 / (2.0 * beta))

N = 20; beta = 4.0; delta = 0.1; del_tau = beta / N x = [0.0 for k in range(N)] data = [] y = [k/float(N)*beta for k in range(N)] for iter in range(1000000): k = randint(0, N-1) x_old = x[k] x_new = x_old + uniform(-delta, delta) xp = x[(k + 1) % N]   xm = x[k - 1] pi_old = rho_free(xm, x_old, del_tau) * rho_free(x_old, xp, del_tau) * exp(-del_tau * x_old ** 2 / 2.0) pi_new = rho_free(xm, x_new, del_tau) * rho_free(x_new, xp, del_tau) * exp(-del_tau * x_new ** 2 / 2.0) Upsilon = pi_new / pi_old if (random < Upsilon): x[k] = x_new if iter % 10 == 0: data.append(x[0]) pylab.hist(data, bins=20, range=[-2.0, 2.0], normed=True) pylab.show code

=References=