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.