# Wave functions

Quantum systems are described by wave functions and eigenvalues, which are solutions of 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,
$H \psi_n = \left( -\frac{\hbar^2}{2 m} \frac{\partial^2}{\partial x^2} + V \right) \psi_n = E_n \psi_n$
we consider the harmonic oscillator, a particle of mass m in a potential
$V(x)= \frac{1}{2} m \omega^2 x^2$

This program uses a dictionary (hash table) for the wave functions
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()
 Output of the above Python program: Wave functions of the Schrödinger equation for the harmonic oscillator.

# 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
$\rho(x,x',\beta)=\sum_k e^{-\beta E_k}\psi_k(x)\psi^*_k(x') .$
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
$\textrm{Tr}(\rho(\beta)O).$
For example, the probability of finding the particle at position x is
$\pi(x)=\textrm{Tr}(\rho(\beta)|x\rangle \langle x|)=\sum_k e^{-\beta E_k}|\psi_k(x)|^2$

The density matrix has three important properties:

### Convolution

Given two inverse temperatures $\beta$ and $\beta'$,
$\rho(\beta)\rho(\beta')=\rho(\beta+\beta').$
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
$\rho(x, x', 2\beta) = \int dx'' \rho(x, x'', \beta)\rho(x'', x', \beta)$

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
$H= H^{\text{free}} + V$
is given by
a general expression known as the Trotter formula:
$\rho(x,x', \beta) \xrightarrow[\beta \rightarrow 0]{} \exp (- \beta V(x)/2) \rho^{\text{free}} (x,x',\beta)\exp (- \beta V(x')/2)$
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
$\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]$

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

## Matrix square

In the following Python program, we implement the matrix squaring algorithm
$\rho(x,x',2\beta) = \int dx'' \rho(x,x'',\beta) \rho(x'',x',\beta)$
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'.

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()
 Output of the above Python program: Comparing the exact solution for the harmonic oscillator density matrix to the solution obtained through matrix squaring.

700px|left|Comparing the exact solution for the harmonic oscillator density matrix to the solution obtained through matrix squaring.
<br clear="all" />

# Path Integral

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()
 Probability $pi(x)$ in the harmonic potential at beta=1, obtained through path-integral simulation.