Hello, bonjour, sdravo, ia orana,
welcome to the fifth tutorial of Statistical Mechanics: Algorithms and Computations,
from the Physics Department of Ecole Normale supérieure.
In this week's lecture, we made our first move into quantum mechanics
and quantum statistical mechanics.
We used two languages. The first language is the one of wavefunctions
and of quantized energy levels, both solutions of the Schroedinger equation.
For the harmonic potential, this is what we studied in harmonic_wavefunctions.py.
Besides this traditional approach to quantum mechanics, we also employed the language of density matrices.
This lead us directly to matrix squaring
to the Feynman path integral and to path integral Monte Carlo simulations.
To make the connection between wavefunctions on the one hand
and path integrals on the other, we must provide the missing link
the density matrix at small beta, that holds the path together.
This is like a spring constant
of our quantum mechanical elastic object
We will arrive at this connection between Schroedinger and Feynman in two steps
First, in what follows
Michael will actually compute the density matrix for a free particle
and derive and illustrate a result that was already presented in the lecture.
He will then extend this simplest of all calculations
using the path integral representation.
Then, Alberto will derive and explain
the Trotter decomposition.
Remember: among the three properties of the density matrix
this was the one which allowed to incorporate
arbitrary interaction at high temperature.
This simulation knows
about the free density matrix and the Trotter decomposition
but not about much else.
And certainly not about the solutions of the Schroedinger equation.
So, now in these simulations
let us spend a moment to speak about the axis.
In the wavefunction picture, the axis are the position x and the energy level E.
So things are clear.
But in the path integral picture, the axes here is again x, the position
but the y axis refers an "inverse temperature like" variable, tau.
In quantum statistical mechanics, we call this variable the imaginary time,
and this mysterious name adds to the fascination of path integrals.
Vivien, in the last part of this tutorial,
will tell us about the relationship between temperature and time.
He will convince us that the Trotter decomposition at high temperature
can also be used at small times.
This will allow you to simulate time evolution on the computer,
and to see some of the amazing effect of quantum physics, like the quantum tunneling effect.
But now, let us listen to Michael
explain to us the whence and the whither of the free density matrix
and how it combines with the language of path integrals.
As Werner mentioned in this week lecture
we are able to explicitly write down the density matrix of a particle
in an arbitrary potential V. We do this by use of the Trotter decomposition formula.
Unfortunately, though, this approach is limited to the high-temperature regime.
On the contrary, for a free particle (that is the case V=0)
we are actually able to calculate an exact result which is valid at any temperature.
So let us see how this is done.
As in this week's lecture, we consider the Schroedinger's equation,
and in the case of the free particle it reduces to the expression shown here.
This is simply a wave equation
and you can easily convince yourself that it is satisfied by any sine or cosine wave
as sine and cosine simply reproduce themselves
in the second derivative, up to a sign and a scaling factor.
When we put the particle into a box of size L
say with periodic boundary conditions
then we can only use sine or cosine waves
whose period is equal to the box length L.
Other waves simply don't fit.
The complete set of waves that do fit is shown here
The superscripts denote symmetric and antisymmetric wavefunctions
with respect to the center of the interval from 0 to L.
You see one ground-state
and the symmetric and antisymmetric solutions
for each excited state, with the energies E_n
= 2 pi^2 n^2 / L^2.
In Python, this gives the following program
it is analogous to the program harmonic_wavefunction.py
that was introduced by Werner in this week's lecture
and it computes the density matrix by summing over all states.
So let us have a look at the density matrix
by simply plotting the results of our little program.
The stripe on the diagonal of the matrix
indicated that the density matrix rho_free(x, x', beta)
is in fact a function of x-x' only.
This means that rho_free is translation invariant
as the individual values of x and x' do not matter.
only the distance x-x' counts.
rho_free is largest on the diagonal
In fact at high temperature (that is small beta)
it is all concentrated on the diagonal, where x=x'.
As we lower the temperature, we see that the off-diagonal part
become larger and larger with increasing beta.
This is the famous coherence, that is governed by the thermal De Broglie wavelength.
Now let us do a small calculation
but for simplicity, we will not use the sines and cosines
but rather the equivalent complex exponential wavefunctions.
You can actually check this equivalence by writing a small Python program.
When you compare the output of this program to the output of the previous program, you will see
that it arrives exactly at the same results as before.
However, we promised an exact expression
for the free density matrix, so let us now calculate it.
Using the complex exponential wavefunctions
the density matrix is given by this expression
where you notice that each term is in fact a complex number.
At first glance this seems a bit strange as we know that the finite result will be real.
This is because for each complex number in the sum
there's a complex conjugate counterpart and together they become real.
Having the particle in a box of size L
is no true limitation, as we can let L go to infinity
and replace the sum over n by an integral.
This may sound difficult, but it can be easily understood
when we introduce a dummy parameter delta_n = 1,
difference between two successive n values.
We then change the variables from n to lambda
with lambda = 2 pi n / L
Thus we have delta_lambda = 2 pi delta_n / L
This leads to this term by term equivalent sum
where we use this result
which follows when we take ctilde = c / sigma^2, from the Gaussian integral.
This Gaussian integral was already evaluated in Lecture 4, last week's lecture.
We finally arrive at the free density matrix
where we have reinserted Planck constant hbar
and the particle mass m.
The free density matrix is indeed a function of x-x'
and see how at large temperatures (that is small beta)
rho_free becomes infinitely peaked
at the diagonal x=x', as we noticed earlier.
Note in particular that the density matrix is positive, not only on the diagonal but everywhere.
For us Monte Carlo specialists
this is important as you can interpret the density matrix as a probability.
Now in this calculation
we made a transformation from a sum to an integral
when taking the limit L to infinity.
But what is the density matrix of a finite box with periodic boundary conditions, that is for finite L?
The above sum can be evaluated using the famous Poisson's sum rule
and the result is really simple.
The parameter w is called the winding number and you will see in a moment where this name comes from.
The density matrix rho(x,x',beta) is made up
of all the paths
(in Feynman's path integral picture) that go from x to x' directly
plus all the paths that wind around the box
once, that is they have the winding number w=1
plus all the paths that wind around the box twice and so on and so on for higher winding numbers.
Of course, we can also wind a path to the left, corresponding to negative winding numbers.
This concept of winding number will come in handy two weeks from now.
It is key to understand Bose-Einstein condensation.
So we have now obtained the free density matrix at all temperatures and it is a Gaussian of x-x'
with a normalization 1 / sqrt(2 pi beta)
and a variance of beta.
Alberto will determine a general expression
for all potentials at high temperatures, using the Trotter formula.
But before you go on, please take a moment or three to download, run and modify the programs discussed in this section.
On the Coursera website, you will find the program free_periodic_sine_cosine.py and its cousin
free_periodic_complex_exp.py
Take these programs apart to understand the inner workings of complex arithmetics
that lead to a final real-valued and positive result.
As explained by Michael, for a free particle we can compute exactly the density matrix
this is of great importance
as we can compute the general density matrix - so to speak - which grows from this seed, as we will discuss
For what follows
you must understand that rho(x,x',beta) is a matrix of two indices
x and x', but it is also an operator: exp(- beta H)
The matrix elements of the Hamiltonian H can be written in an arbitrary basis
and in this basis the matrix elements of rho write like this
This expression can be familiar to those of you who have studied a little bit of quantum mechanics
but for all of you, please retain that the density matrix
is an operator, the exponential of minus beta H
it is the direct generalization of the Boltzmann distribution
Let us now consider the density matrix of a system with the following Hamiltonian
H = Hfree + V
where Hfree is the Hamiltonian discussed by Michael
and V is the interaction potential, that for the harmonic oscillator is (1/2) x^2
Now, the Trotter decomposition corresponds to the following expression.
To see that this is true, we will expand both terms up to the order beta^2
but attention, Hfree and V are operators
the product V Hfree can be different from Hfree V.
Operators, matrices in general do not commute.
So, let's start from the second term of this equation
We can first expand exp( (-beta/2) V(x))
'a' corresponds to the linear term in beta, and 'b' to the second order term
Then we expand rho_free
'c' is the linear term and 'd' is the second order term,
finally we can also expand exp(-(beta/2) V(x'))
So let's reorganize these terms order by order in beta.
At the zero order we have one.
At the first order in beta, we have a+c+e
and then we need a little bit of work to write the second order term
This expression fully agrees with the direct expansion of the full density matrix
These two equations prove
that the Trotter expansion is correct up to beta^2.
Please note that this manipulations
are valid on operators and wavefunctions
in a finite box.
In Python we can implement the Trotter's decomposition using Michael's approach,
and plot the entire density matrix at high temperature.
Remember: the free density matrix
rho_free(x,x',beta)
is a function of x-x'
and this gives the nice stripe in this diagram.
Now let us add the harmonic potential V to the Trotter decomposition.
You see that for large x and x' the density matrix is suppressed
and this is what keeps close to the origin the quantum harmonic oscillator.
In this week's homework, you will use
both the Trotter decomposition
and the convolution properties
to compute the density matrix up to very low temperatures.
As you must square the matrix many times, the approximation errors will grow
However, the quality of this approximation (correct up to beta^2) is so high that it will survive
up to low temperatures.
The Trotter break-up is completely general:
it applies to an arbitrary potential
and (as will be shown by Vivien in a moment)
it can be even turned around to describe the real time evolution of the system.
But before understanding why the inverse temperature beta is an imaginary time,
please take a moment to download, run and modify
the program harmonic_trotter_movie.py
which implements the Trotter decomposition and provides a nice graphical output.
You can check this approximation
and actually in this week's homework
you will test this approximation against the output of harmonic_wavefunctions.py
In this section we examine a complementary problem
which is the time evolution of a quantum system
Werner, Michael and Alberto have explained to you
in the previous sequences
how to describe a quantum system at finite temperature
in equilibrium, by the use of the density matrix
However, in general, systems are described by
a wavefunction psi(t) which evolves in time.
So, if we fix this wavefunction psi(t=0)
to be equal to psi0 at initial time
it evolves in time according to the Schroedinger equation displayed here
which is i times dpsi/dt
= H psi
and here H is the Hamiltonian of the system
and we have also fixed the constant hbar to be equal to one.
To solve this equation, one may be tempted to replace
the time derivative by a first order approximation
which is displayed there
but this is not an option in our case because this would break the conservation of probability.
A better choice is to consider the operator solution
of this equation, which writes as follows
It is psi(t) = exp(-i H t)
times the initial vector psi(0)
You may indeed check, by differentiating with respect to t
that this expression is the solution of the Schroedinger equation
with initial condition psi(t=0)=psi0.
Now, let's have a closer look to this formula.
You note that by formally replacing
(i t) by beta
that is to say, by working in imaginary time
you obtain the very same expression as that of the density matrix
that we know for a quantum system at equilibrium at inverse temperature beta.
So this relation provides us a dictionary
a mathematical dictionary between time evolution and thermal equilibrium.
It is not only beautiful from the abstract point of view
but also very useful in practice
Indeed, you can apply what Alberto has told you
about the Trotter expansion
to devise an algorithm which allows you
to study very efficiently the evolution of quantum systems
Let us now focus our attention to a particle at position x
in a potential V(x) in dimension 1
If we set the mass of the particles to be m=1, its evolution
is governed by the Schroedinger equation
of Hamiltonian H = Hfree + V(x)
The idea is to decompose the full time interval t
into N time steps of duration delta_t = t / N
The full operator of evolution is now expressed exactly as a product of N
infinitesimal operator of evolution, as displayed here.
Each of these infinitesimal operators of evolution writes exp(-i delta_t H)
Now to understand the contribution
of the components Hfree and V(x) of this Hamiltonian
we can use Alberto's Trotter decomposition, which writes as follows
You recognize the same formula, where beta has been replaced by (i delta_t)
This approximation is valid up to order delta_t ^2
that is to say, the first error term is of order delta_t^3.
Let us examine each term of this product
the first and the last term are easy to understand
they consist in the multiplication by the same factor exp(-(i/2) delta_t V(x))
The free term is a bit enigmatic
but the mystery is lifted by transforming the wavefunction to Fourier space
for the spatial coordinate.
This consists in going from the real space psi(x,t)
to the momentum space psi_hat(p,t)
You will easily convince yourself - for instance by reading the short fact-sheet associated to this tutorial 5 -
that the action of the free operator
consists in a multiplication by a factor exp(- (i/2) p^2)
So, the action of this operator is just the multiplication by a Gaussian factor.
This provides a possible scheme to devise an algorithm which is the following:
first, define the initial wavefunction, then
discretize the time interval t into N time slices of duration delta_t,
define the Fourier procedure,
and apply the Trotter formula many times.
At each time step
first in real space, multiply by a factor exp(-(i/2) delta_t V(x)),
make a Fourier transformation,
multiply by exp(- (i/2) p^2 delta_t),
and go back to the real space multiplying by the same factor as previously.
This is the only step, you just have to iterate and it gives the following beautiful algorithm
quantum_time_evolution.py.
This method is really powerful.
To illustrate it, let's consider the fate of a quantum particle confined between two barriers:
an infinite barrier on its left and a finite barrier on its right.
The particle starts from some position between the two walls.
If the particle were classical, it would bounce forever between the two walls.
You observe that each time the wavefunction bounces on the wall on the right
a finite but tangible proportion of the wave actually goes through the wall.
This quantum tunneling effect is a genuine quantum effect and has no classical counterpart.
We have illustrated how our mathematical dictionary between
real time and imaginary time
has allowed us to observe the quantum tunneling effect in a simple one-dimensional situation.
However, have in mind that the very same effect
has deep physical consequences
from for instance the radioactive decay to the construction of the transistor in the CPUs of your computer,
or in more complex situation to the devise of the scanning tunneling microscope,
or even to the evaporation of black holes.
It is now time for you to download, run and modify
the beautiful program we have seen in this section: quantum_time_evolution.py.
Check it out.
See by yourself how a quantum particle can tunnel through a wall
and program maybe some other exciting potentials
Take also a look to the fact-sheet which explains to you how the Fourier transform from x
to p is implemented in practice.
In this tutorial
we zoomed in to the density matrix and its three main properties
the convolution, the free density matrix and the Trotter decomposition
we already perceive the plot of these three stage characters
that prepare for a great story and that will lead us to quite complicated situations
In what Alberto and Michael told us
the density matrix remains positive and real
allowing us to interpret it as a probability, and to do quantum Monte Carlo.
This positivity is preserved for many particles,
arbitrary interactions, and even for bosons
It is broken in the case of fermionic particles
In what Vivien told us
in the time-evolution all of this was broken also
the density matrix was complex-valued and it could be negative
The great algorithm produced only works in one-dimension
because quantum time evolution is a really hard nut to crack
So now it is time for you
to look at how these three characters that we described
intervene in practice
and we wish you good luck and a lot of fun with the homework session 5.