Hello, bonjour, ni hao, hola, welcome to the second week of Statistical Mechanics: Algorithms and Computations from the Physics department of Ecole Normale Superieure. This week's subject is the physics of hard disks. Remember, hard disks: two-dimensional idealization of billiard balls, not the hard disk in your computer. These hard disks - or at least their idealization - were studied with two approaches. Molecular dynamics, the solution of Newtonian deterministic mechanics, and Monte Carlo, the solution of Boltzmann statistical mechanics. As discussed in the lecture, for hard disks, Boltzmann statistical mechanics means equal weights for all legal configurations. We also discussed that Boltzmann and Newton are equivalent for all thermodynamic properties, that means for all properties that do not explicitly depend on time. In this tutorial, we will study direct sampling for hard disks. We will find out that although the program direct_disks.py is not a very powerful program it is a VIP, a very important program, as it contains important connections between physics and computation. In what follows, Michael will scrutinize the equal probability principle and the tabula rasa rule. Then Alberto will explain the profound link between computation and physics, or rather, between the acceptance probability of the direct sampling algorithm and the partition function of the physical system. Finally, Vivien will explore an analytical method for studying the hard disks system, the virial expansion. The virial expansion has existed a long time before numerical calculations became available. Ludwig Boltzmann himself, the founder of statistical mechanics wrote a famous paper in 1874 on the virial expansion of the hard sphere system. So you see, we also keep to the center of the road of statistical mechanics. In this week's lecture, Lecture 2, Werner explained that the equiprobability principle is one of the pillars of statistical physics. This principle states that the probability pi(a) of a configuration is a function of the energy of the configuration, that is pi(a) equals pi of the energy of a. For hard disks, all allowed configurations have the same energy E=0. This gives the same statistical weight to each configuration pi(a) = pi(b) = pi(c). The Monte Carlo algorithm must sample these configurations with equal probability, this was achieved by the direct sampling Monte Carlo algorithm. Remember that the algorithm direct_disks_box.py first places one disk at a randomly chosen position inside the box, then a second one, then a third one, but in the case of an overlap applies the tabula rasa rule and wipes out the entire configuration, and starts afresh. But why don't we simply take the last disk (the one that caused the overlap) and try to put it a second time or a third time until it finally finds its place? This algorithm is called random sequential deposition, and it is important for adhesion and catalysis, but not for equilibrium statistical physics. Its configurations are not equally probable as we will now see, using two arguments. To do so, let us first consider a discrete hard rods model, a one dimensional discrete version of hard disks. On the line shown here, suppose that the rod centers can be placed on the sites 0, 1, 2, 3 and 4. Now assume that there also is an overlap condition, just like for the two-dimensional hard disks. We say that two adjacent rod centers must differ by at least three, this means that there must be at least two free sites in between two adjacent rods. For example we cannot place rods on the sites 0 and 2, but there's no problem to place them on the sites 1 and 4. Let us do a simple simulation of random sequential deposition, where we impose that we end up with two rods: the red rod and the blue rod. We place the red rod, then we place the blue rod, but we have to pull it off again because there's an overlap and we put it again, we pull it off again, we put it again until we succeed. Here you can see a second simulation. We now analyze the algorithm: we first place the red rod and it'll fall with the probability of 1/5 on any of the sites 0, 1, 2, 3, 4. Rather, with a probability of 1/4 because there is no two-disks configuration with one of the disks on site 2. Now, if the red rod is on site 0, the blue rod can be placed with a probability of 1/2 on the sites 3 and 4. This gives the allowed configurations a and b. The same reasoning is valid for the configurations e and f. These four allowed configurations have the probability 1/4 times 1/2, equal 1/8. In contrast, if the red rod is placed either on site 1 or on site 3, there's only one site left for the blue rod to go. This leads to the configurations c and d, which have the probabilities 1/4 times 1, equals to 1/4. Clearly, we have violated the equiprobability principle. What would happen if we apply - instead of random sequential deposition - the tabula rasa strategy of the Monte Carlo algorithm? In Python, this would give the program shown here as you will discover in the following little quiz. Question 1: how many configurations (legal and illegal ones) of red and blue rods can be generated by the program direct_discrete.py? The answer is 25, as there are five positions for the red rod and five positions for the blue rod to be placed. Question 2: how many legal configurations are among these 25 configurations? Remember that there must be at least two free sites in between the centers of the red and the blue rods. The answer is 6, as only the 6 configurations a, b, c, d, e and f (that we discussed earlier) survive the tabula rasa wipe out. Question 3: what is the probability of the legal configurations a,...,f? The legal configurations a,...,f have the probabilities pi(a) = pi(b) = .. = 1/6. The direct sampling Monte Carlo algorithm generates the 25 legal and illegal configurations with equal probabilities 1/25. The legal configurations a,..,f have equal probabilities before the tabula rasa wipe out and consequently have equal probability after. The argument we just used for the discrete model applies also for the four hard disks in the two dimensional box, if we reprogram direct_disks_box.py so that it places the four disks at randomly chosen positions without checking for an overlap. Now, there's no doubt that these four disks configurations are uniformly distributed, in an eight-dimensional space as we have two coordinates (namely x and y) for each of the four hard disks. In fact, as we do not check for overlaps, each point in this eight-dimensional space is sampled with equal probability. It is only after we have generated all these configurations that we cut out pieces. In this program, the vector L is a uniform random vector in an eight-dimensional hypercube from sigma to (1-sigma) in all dimensions. We then cut out pieces from this hypercube with the overlap condition. This is like in the discrete model, where we threw away 19 of the 25 configurations. In the next section, Alberto will deepen our understanding of the eight-dimensional space picture for hard disks, but no longer with walls, but with periodic boundary condition in the x and y directions. Before following these exciting developments you should download, run and modify the relevant programs. On the Coursera website, you will find programs for the discrete rods model. First, the program for sequential random deposition which also outputs a nice histogram, and as usual also a random sequential deposition program with nice graphics. You can also find the program for direct sampling of the four hard disks in the box, if you have not already downloaded it in Lecture 2. And even the slower program which first samples all the positions, and checks for overlaps only in the end. Make sure that you understand that direct_disks_box and direct_disks_box_slow are statistically speaking strictly identical in their output. Let's move to direct sampling of hard disks with periodic boundary conditions. As you can see here, when a disk is close to the right hand side of the box, it will appear also on the left hand side of the box. If you have two disks they overlap not only if they're close together but also if for example one disk is on the lower right corner and the other disk is in the upper right corner In a system with periodic boundary conditions, walls play no role Here we show the multirun Python version with periodic boundary conditions. Many legal configurations are generated Check out for yourself that the function dist() correctly implements the periodic boundary conditions. Now we consider more than four disks. Let's say N=16 disks, and we fix also the density eta, which is the fraction of space occupied by disks. Here, we show the output for a density eta=0.3, only 30% of the space is occupied by disks. For our box of side 1, the radius of the hard disks can be obtained by the relation 16 x pi x sigma^2 equal to 0.3. Something chilling is happening: we typically have to wipe off our table 100000 times before having a single legal configuration. This means that our algorithm is exact but is very bad, because the acceptance ratio is very small. For a system of 16 disks with a density of 0.3 the acceptance ratio is of the order of 5 x 10^(-6). Michael spoke a few minutes ago about the eight-dimensional space of configurations of four disks These configurations are described by the eight coordinates of the disks (x_1, y_1), (x_2, y_2), (x_3, y_3), (x_4, y_4). Here we have 16 disks, which means that the dimension of the space is 32 Furthermore, if our box has a volume (or an area) equal to V the 32-dimensional space has a volume V^16 or V^N, where N is the number of disks. If the system has a density equal to zero, which means that the radius of the hard disks is equal to zero, all the points of the configurational space are legal configurations. So the number of legal configurations, which is the partition function of the system, is V^N. At density 0, the acceptance ratio is 1. We found that at density 0.3 the acceptance ratio is 5 x 10^(-6). This means that the number of legal configurations is 200000 smaller than the volume of the original 32-dimensional space (V^N). The number of legal configurations at density eta - which is the partition function of the hard disks - is equal to the product of the partition function at density 0 times the acceptance ratio at density eta. The partition function is a central object in equilibrium statistical mechanics and is a central object in Statistical Mechanics: Algorithms and Computations. In a few moments, Vivien will discuss this in more details. and will provide an approach (an analytical approach) to compute the acceptance ratio and the partition function for hard disks. Before doing this, please take a moment to download, run and modify the programs we discussed in this section On the Coursera website, you will find the algorithm for direct sampling with periodic boundary conditions. This algorithm allows you to pick out legal configurations from the hypercube of dimension V^N. These configurations are as rare as diamonds, and if you wish they are almost as beautiful. Note that all coordinates are now random numbers between 0 and 1 because no wall has to be considered. You can inspect these legal configurations, the diamonds of statistical physics using, as usual, the graphic version of the program. Let us go one step further and determine the acceptance ratio of the direct sampling algorithm at any density eta. You may be tempted to run the program direct_disks_multirun for several runs and maybe on many computers for different values of eta. But I'd rather look at the following configuration of N=16 points in a box with periodic boundary conditions. At zero density, that is to say for zero radius it clearly leads to a valid configuration of hard disks and this is also the case for a density equal to 0.0001 In contrast, at higher density the radius is larger, and the points lead to a configuration of hard disks which is illegal. In fact, from these 16 points in the box you can obtain a hard disks configuration for all values of the radius which are smaller than the maximal radius, equal to half of the minimal distance between two of these points. In this present configuration, the value sigma_max of this maximal radius is equal to 0.025 and it corresponds to a maximal density eta_max = 0.03143 You may consider the probability distribution p(eta_max) associated to this maximal density and from it - as displayed on this formula - you can reconstitute the probability distribution pi_accept(eta) of the acceptance ratio of the direct sampling algorithm. You can find this formula on the fact-sheet associated to this class session. For 16 disks you obtain the following graph which has been obtained from a few millions realizations of the program. It describes the acceptance rate as a function of the density eta and you observe that it decays faster than exponentially with eta. We will now obtain an analytical approximation of the acceptance rate and of the partition function of hard disks. It will be excellent in the limit of small densities, and valid for any value of N To do so, let's consider the following N=16 points in the box with periodic boundary conditions This configuration is described by the set of vectors x_0, x_1, .., x_(N-1) each of these vectors contains the two spatial coordinates of the point, x and y. For relatively small radius, this configuration will be allowed In general this configuration will be valid provided that the distances between disk 0 and disk 1, disk 0 and disk 2 up to N-1 and N-2 are all larger than 2 sigma. If one of these conditions is not valid, then the configuration is illegal. Following Alberto a few minutes ago, we will rewrite this condition in integral form, as follows. This integral runs over the 2 N components of the N vectors of the configuration and Upsilon is equal to 1 if there is an overlap, and 0 otherwise. It is important to note that if one of these brackets is equal to zero then the contribution of the configuration vanishes in the integral This expression of the partition function is exact, but no one has been able to compute it analytically. Let's consider all the brackets, each of them is equal to 1 if there is no overlap, and 0 otherwise. We can expand and multiply out all these brackets For instance if one picks the ones in each of these brackets and then one computes the integral, one obtains V^N. This is in fact the partition function - as Alberto showed - of the ideal gas of hard disks of radius equal to zero. To go to the next order, let's consider all the brackets there are (N-1)N/2 of them, so what you can do is to pick one of the Upsilon in a single term, and to pick one in all the other brackets. Doing so, you obtain a four dimensional integral, which is displayed here. Thanks to the invariance by translation and the periodic boundary condition, it reduces to V times a two-dimensional integral. We observe - as displayed in this figure - that the excluded region corresponds to a disk of radius 2sigma. This allows to compute the integral as follows. Now, we can come back to the full expression of the partition function, and if we stop the evaluation at this term we obtain the following expression. For large N, it behaves as V^N times the exponential of - (N-1) times eta. This means that the acceptance probability decays exponentially with both the density and the number N of particles. Remember that we have a computer program called direct_disks_any.py, which allows to evaluate numerically the acceptance rate. We can thus compare it to our prediction, and this is displayed on the following picture As you observe the analytical prediction matches very well the numerical results in the low density regime. Let's now remark that the partition function is proportional to V^N To obtain an extensive quantity, in physics we take the logarithm of this quantity and you obtain the following result. Let us now explain why it corresponds to a virial expansion. To do so, we can first differentiate log(Z(eta)) with respect to V and then multiply by V/N. You obtain the following expression the first term is a constant equal to 1 and it corresponds to the partition function V^N of the ideal hard disks gas of radius zero, that Alberto explained. The next term is proportional to 1/V In fact, we obtain this term in 1/V because we consider only the term with one Upsilon in the expansion. We could go on in this expansion of the brackets and consider pairs of Upsilon, triplets and so on This would yield an expansion - as previously - in powers of 1/V However, this is not the state of the art of statistical mechanics for the 21st century. Indeed, Ludwig Boltzmann in 1874 computed the 4th term of the virial expansion, that is to say the term in 1/V^3. It was for three-dimensional hard spheres. This computation was really complicated, it was subject to years of debate, but in the end he was right and others were not. The virial expansion was once believed to provide a systematic access to the thermodynamics of liquids and gases, in our case hard disks. In the same as the expansion of the exponential function exp(x) = 1 + x + x^2/2 + x^3/6 and so on. This expansion is valid for all values of x, real or complex. However, for the virial expansion this is not the case any more. What is proven today is that this expansion is valid only up to some finite values of the density eta. However it cannot predict the existence of a phase transition. Those phase transitions will be the subject of our next week. Now, please, download run and modify the programs we have seen during this session. You may also have a look to the fact-sheet that corresponds to the computation of the virial expansion we have just seen. On the Coursera website you may also have a look on the following program. It is sweet and short, it allows to determine the acceptance ratio for any density eta. Have fun going through this those are 15 lines, but not easy to understand. In conclusion, we have moved in this tutorial from the tabula rasa rule to the discussion of the high-dimensional abstract spaces of our hard disk configurations. We moved towards the center of statistical physics: the partition function and the virial expansion. Our discussion revealed a beautiful connection between algorithms and the physical description of our systems. But our analytical approach remained perturbative. This is a true limitation We will see next week that our algorithms have the potential to carry us beyond this horizon and this limitation of the perturbative approach This will allow us to study phase transitions, when the liquid of hard disks becomes solid, and this is a major revolution in physics. In the meantime, have fun with the homework session of this week and see you again next week at Statistical Mechanics: Algorithms and Computations.