Hello, bonjour, cama-i, ainngai,
Welcome to the the sixth week of Statistical Mechanics: Algorithms and Computations,
from the Physics Department of Ecole Normale supérieure.
We are right now in the middle of our three-weeks module on quantum physics and quantum statistical mechanics
And we have just started our second half of this course.
So let me start with a brief exposition
of where we have come from and where we stand right now.
This exposition will be in the light of what we already learnt
about the connection between sampling and integration.
So, five weeks ago we started by considering pebbles
on the Monte Carlo beach and on the heliport
the partition function was the integral in x from -1 to 1
and the integral in y from -1 to 1 of pi(x,y)
which was a constant equal to one.
We then moved to a discussion of hard disks
the partition function now was this integral, with the statistical weight pi(x)
which was equal to 1 if the configuration was legal, and 0 if it was illegal
last week we then moved to a discussion of quantum physics
The partition function now was the path integral
over dx_0, .., dx_(N-1)
of the statistical weight which was the density matrix between x_0 and x_1, beta/N
density matrix of x_1 and x_2, beta/N
density matrix of x_(N-1) and x_0, beta/N.
This is the path integral
We sampled it using a naive path integral Monte Carlo algorithm
that resembled what we already did for hard-disks.
But rather than picking one disk at random, we picked one bead at random,
then we moved it a little bit to the left or to the right,
and accepted it with the Metropolis acceptance probability.
This week, we will remain within this picture of quantum physics.
We will present a providential sampling algorithm that goes back to Paul Levy in 1940.
It is a rejection-free direct-sampling algorithm, something we never achieved for hard-disks.
Let me briefly explain what the Lévy algorithm achieves.
You pick a starting position x at tau=0
and a final position x' at tau=beta
and then you fill in the path between x and x' in one set.
You can do this with a number of slices that is as large as you want: 100, 1000, 1000000..
Here we have a free particle, but an analogous Levy construction
exists also for a particle in a harmonic trap.
Lévy's algorithm is really famous, and not just in quantum statistical mechanics.
It exists under various names in the physics of interfaces
in numerical analysis, in statistics, in financial mathematics,
and in electrical engineering.
The path you see here illustrates the Heisenberg uncertainty principle:
the founding principle of quantum mechanics, and one of the pillars of quantum statistical physics.
At low temperature (large beta) the path fluctuates wildly
and this means that the particle position x
is uncertain.
Next week, the Levy construction algorithm will be a crucial element
of many-body simulations, and then we will fathom
the second pillar of quantum statistical mechanics: the indiscernability of particles,
and their bosonic or fermionic nature
and the curious Bose-Einstein condensation
To discuss Bose-Einstein condensation, we will again mostly use the path integral language
but to understand why this is a condensation
we need to consider wavefunctions and energy levels.
And this is what we will do in this week's tutorial
where we will get down to the art of bosonic statistics
and we'll solve a model of finite (small) number of particles in a three-dimensional harmonic trap
In this week's homework session we consolidate what we learn in the lecture through practical calculations.
For simplicity, we remain within the framework of the harmonic potential
but the firework of sampling strategies and re-weighting schemes
applies to a much wider range of problems
than just to four little quantum particles in a harmonic trap.
Again, you must understand the great range of applications
of our concepts inside and outside physics.
Quantum paths are closely related to interfaces
as you see here for an interface between two liquids.
Quantum paths are also closely related to what is studied by financial analysts.
So ,let's get started
with week 6 of Statistical Mechanics: Algorithms and Computations.
We recall once more the final subject of last week's lecture
We arrived at a representation of the partition function
of a quantum particle in a harmonic potential
as a sum
or rather as a multiple integral over paths with a statistical weight giving the weight
of each path.
You see the partition function really is a sum over all paths.
We sampled this partition function
using a naive algorithm, naive_harmonic_path.py
Its key element was that we sampled a random bead k
and proposed a random move
from x_k to x_k + delta_x
and accepted or rejected the move with the Metropolis acceptance probability.
naive_harmonic_path satisfies detailed balance, is aperiodic and irreducible
but it painfully slow.
If we start at a configuration as the one shown here
the next configuration will be very little different from the one that we had already,
because we only move one bead and we move it only by a little bit.
If we made larger moves,
we would get into troubles with the 1/2 rule.
We would reject all the moves, and would not move either.
So you see that a simulation using naive_harmonic_path
would be very very slow in exploring all configurations in our system.
To make progress let us move back by two steps
first step: we take away the harmonic potential and consider a free particle,
second step: instead of sampling the value of k at each step, we keep k fixed.
We move the particle x_k
for fixed values of x_(k-1) and x_(k+1)
as programmed in naive_path_slice.py.
So we move x_k with the Metropolis algorithm,
and here is the histogram of the positions of x_k, it looks like a Gaussian,
and in fact IT IS a Gaussian,
because the probability pi(x_k)
is put together from the density matrix rho_free(x_(k-1), x_k)
and rho_free(x_k, x_(k+1))
The two are Gaussians, and their product is also a Gaussian for x_k.
So the distribution pi(x_k) is a Gaussian
with mean value (x_(k-1) + x_(k+1) ) / 2
and a variance sigma^2 = delta_tau / 2
Now, there is no reason in the world to use the Metropolis algorithm
as we know that the distribution of x_k is a Gaussian.
We should sample x_k from this Gaussian distribution with the correct mean value and variance.
You might be interested to implement this idea
in your version of naive_harmonic_path, as we did here for the free particle
but rather than doing it, let us slightly generalize this problem
and consider a position x_k
sandwiched between positions x' and x''
with intervals in tau: delta_tau' and delta_tau''
The probability distribution of x_k given x' and x''
is the product of two free density matrices:
the density matrix rho_free(x', x_k, delta_tau')
and the density matrix rho_free(x_k, x'', delta_tau'')
Both these density matrices are Gaussians and their product is also a Gaussian.
Putting all this together and noticing that in this distribution x_k is the variable
and x' and x'' are fixed
we find that the probability distribution pi(x_k) given x' and x''
is a Gaussian with mean value and variance as shown here.
Before continuing, please take a moment to download, run and modify the programs that we discussed in this section.
On the Coursera website, you'll find the program naive_path_slice.py
that treats the case of one variable x_k
sandwiched in between x' and x''
This is an innocent little program, but it should convince you
that the above formulas that we used are correct.
In any case, running this program, you'll understand that using the Metropolis algorithm in this situation is useless:
you might just as well sample x_k from the appropriate Gaussian.
Then - again - you'll find the naive_harmonic_path.py
In this algorithm, you could replace the Metropolis algorithm
for the x_k by the Gaussian direct sampling
but better wait for greater things to come.
In naive_path_sampling
we propose x_k with a flat distribution and accept it with a Gaussian,
The mismatch between proposed and accepted moves generates of course the rejection rate in the Metropolis algorithm.
We could modify the naive path sampling algorithm
by proposing x_k with the appropriate Gaussian distribution and there would be no rejections any more.
But we can put the conditional probability pi(x_k)
to much better use.
In fact, pi(x_k | x', x'')
gives the statistical weight of all paths that start at x'
that pass through x_k and that end up at x''.
Let's apply this idea for x_k = x_1
x' = x_0 and x'' = x_N.
We can sample the position x_1, given x_0 and x_N
Between the freshly sampled value of x_1 and x_N
we can then sample the value of x_2
and then of x_3 and then of x_4
and eventually the entire path
In Python, this gives the following little algorithm: levy_free_path.py
and look here at a configuration for a path with N=50000
produced by levy_free_path.py
It has no correlation with previous paths
and its construction has caused no rejection.
In the limit N going to infinity, x(tau)
is a continuous function of tau
As a curiosity, let us discuss for a moment
the return condition, that makes the Levy construction non-trivial.
The path has to return to the position x' at tau=beta
without this return condition, you would simply have a random walk,
and you would sample the position x_1 from a given position x_0,
then we would sample the position x_2 from the given position x_1 as a Gaussian.
Then you would sample x_3 from x_2, x_4 from x_3,
and so on, until x_N from x_(N-1).
The partition function of this problem would be written as Z_random_walk = integral dx_1 .. dx_N
of the density matrices rho_free(x_0, x_1, beta/N)
until rho_free(x_(N-1),x_N,beta/N)
And notice that we would sample over the position x_N, that means we integrate over the variable x_N.
This is implemented in the algorithm continuous_random_walk.py
but of course the probability to hit the position x' at tau=beta
the probability that x_N = x' is equal to zero.
Now, let me tell you of a simple algorithm (trivial_path.py)
that will rectify the situation.
What you can do is that you sample the Gaussian random walk that goes to x_N rather than x',
then you pull back the position x_N by the value of (x_N - x')
and you pull back all the intermediate positions
by a value proportional to (x_N-x') * tau/beta.
This algorithm produces output
which is indistinguishable from the output of the Levy construction
but I leave it to you as an exercise to prove this.
Before moving on to next section, please take a moment to download, run and modify the algorithms we discussed right now.
There is the algorithm levy_free_path.py
In this algorithm you really can modify many aspects
there is no reason to first sample x_1 than x_2, x_3 and so on
you can rewrite your algorithm so that it first samples x_(N/2)
and fills in the path in the lower and upper parts separately.
Then there is the algorithm continuous_random_walk.py
and the mysterious pull-back algorithm: trivial_free_path.py.
try to understand - at least empirically - that this algorithm gives exactly the same output as the Levy construction.
We have sampled so far free quantum paths
continuous non-differentiable objects that correspond to the free Hamiltonian.
This is possible because the free density matrix is a Gaussian and product and integrals of Gaussians are again Gaussians.
For the harmonic density matrix, in Trotter approximation, the same is true.
It is given by a product over three Gaussians: two for the Trotter formula and the potential, and one for the free density matrix.
So at high temperature (small beta)
the harmonic density matrix is given by an exponential of two terms: one in (x-x')^2 and one in (x+x')^2
Under matrix squaring, the harmonic density matrix keeps this form
and an exact calculation shows that the explicit formula
for the harmonic density matrix at all temperatures is given by this expression.
The diagonal term rho_harm(x, x, beta)
is given by exp(-x^2 * tanh(beta/2))
The fact that it is a Gaussian
allows us to derive an algorithm for the Levy construction of the particle in a harmonic trap
This is levy_harmonic_path.py
and it is as simple as the free Levy construction.
Again, we can now choose a configuration x at tau=0 and x' at tau=beta
and the algorithm fills in
a configuration, a path, between 0 and beta,
without rejections and without using Markov chain methods.
Of course, the path is now constrained through the harmonic potential
and its extension is fixed by the oscillator strength omega
and the temperature.
Before continuing, please take a moment to download, run and modify
the program we just discussed: levy_harmonic_path.py
Besides this program and the graphics version, you'll also find a short fact-sheet
that gives the relevant formulas and derives the harmonic density matrix through matrix squaring.
You also find on the website a multidimensional version of this program
We have no problem understanding that a three-dimensional quantum path can be sampled
by independent Levy constructions in x, y and z.
In this lecture, we studied the path integral picture of quantum statistical mechanics
and the providential path sampling algorithm that goes back to Levy in 1940.
In free space and in any dimension, this algorithm allows to sample a path from a position x at time tau=0
to x' at tau=beta.
Notice that in this picture we are discussing a quantum system in one dimension: x
but it looks like a classical line-shaped object in two dimensions
It is for this reason that one says that a quantum system in one dimension
corresponds to a classical system in two dimensions, or more generally:
a quantum system in dimension d corresponds to a classical system in (d+1).
As the temperature gets lower, beta gets higher
and the extension of this space in all dimensions becomes infinite.
As we discussed, the Levy construction also exists for particles in a harmonic potential
and you will rely heavily on this to create your own Bose-Einstein condensate
that resembles very much the experimental condensates that were created in this glass cell.
So now it is time for you to download and run the program that we discussed
and to create a quantum path like the one shown here.
Create your own quantum path, with hundreds and thousands of points,
and print it out on glossy paper, make yourself a frame and hang it up in your room.
As mentioned, there are many exciting subjects
that we will discuss in this week's tutorial and we will discover in this week's homework.
So thanks for your attention and see you again next week in Statistical Mechanics: Algorithms and Computations
for a week devoted to Bose-Einstein condensation.