An Introduction to Percolation and Many-body Physics by Computer

In the following we will see how computers are used to solve complex problems in science. Along the way we will need to learn a little statistics and about the numerical technique called the "Monte Carlo algorithm". Our goal is to solve the "wildfire problem" (also known as percolation): imagine a plain filled with patches of dry grass. A fire has started at one end. It has a certain probability of jumping from one patch of grass to another. The question we seek to answer is: What is the probability that the fire crosses the entire plain? To answer this we first need to take a little detour into random numbers.

If you just want to see the percolation applet, click here.

Random Numbers on a Computer

Our first task is to understand what randomness means on a computer. To aid in this click this link now. A window should pop up which contains an applet. When you click "go", the applet generates (=computes. There is a subtlety here) the numbers 1-6 randomly (representing a die). The histogram to the right keeps track of the total number of times each number is generated. This happens in real time so the graph will be modified as you watch it. Click on go now. When the applet is done you'll see that 13 of the hundred dice throws were 6's. We expect this number to be 100/6 = 16 2/3, but of course our process was random, and there is no reason to get exactly what we expect. If you click go again you'll get a different number of 6's.

Can we get a better estimate of 16 2/3? Yes, the more times we throw the dice (generate random numbers) the more likely it is that we'll up with an estimate of 0.1666. Try entering 1000 into the "total number of throws" box (you may have to hit enter while in the box to get it to register). Then hit "go". The number of 6's/1000 should be much closer to 0.1666 than it was when you only did 100 throws. ( I get 0.152 instead of 0.130 when I do it.)

This raises an interesting question. What if I had run the applet ten times with 100 throws each? How does this compare with running it once with 1000 throws? In the second case we know we get a better estimate of the true probability of throwing a 6. What about the first case? We can get a better estimate there too. Here's how: write down your answer for each of the ten runs and make a histogram of the results. It should look something like this.

If I repeat the experiment 100 times (and I did) I get the following histogram.

Exercise Look up the formula for a normal distribution. See if you can estimate the average and the standard deviation by comparing the formula to the graph above. You may want to plot them on the same graph and adjust the strength, mean, and standard deviation of the normal distribution so that it looks like the data.

Monte Carlo Algorithms

Now you know how random numbers are calculated on a computer and a little about how to analyse the results of statistical experiments. Let's move on to computing something. This example involves calculating the area of a circle. In this case we know that the answer is pi R² (I have chosen R=1 in this example), but in general it may be an area that is too complicated to evaluate conveniently and we want to develop a computational method. Imagine throwing darts at a square of size 2 by 2 with a circle in it of diameter 2. If all of our darts land in the square what is the probability that one lands in the circle? The answer is given by the ratio of the area of the circle to the square. In this case pi (D/2)² divided by D². Thus if we keep track of the number of times a dart has landed in the circle, we can get an estimate of the area of a circle! To see this in action, close the dice window and open this window.

The box on the right is the square which we are throwing darts at. When you click go, each dart throw is shown as a red dot if it lands in the circle and blue if it lands outside the circle. The area of the circle is then calculated as the number of darts landing (the third box) divided by the total number of darts times the area of the box (in these units it is 4). This should be pi=3.141592654... My first run with 100 darts gives 3.08, not too bad. With 1000 darts I get 3.16. You can make a histogram of results just as with the dice example, and analysing the histogram will allow you to make a better estimate of pi and of your error in the estimate.

Throwing darts is the basis of the Monte Carlo method. In practice it is wasteful to throw darts into the square so we always try to find a shape which is close to the one we want to calculate the area of (things get more complicated in many dimensions.) But remember, we must know the area of the reference shape! Otherwise we don't know what to multiply by. Choosing a good reference shape is called importance sampling and it is a vital part of the Monte Carlo method.

Exercise Evaluate the area of the circle 20 times with 100 throws each. Make a histogram of your results. Make a rough estimate of the error in your best estimate of the area.

Exercise Write a computer code which calculates the volume of a sphere.

Exercise Write a Monte Carlo code which evaluates the following integral.

Exercise Evaluate the probability for a dart landing in a D-dimensional super-sphere of radius 1 if it is embedded in a D-dimensional super-cube with sides of length 2. You should get (take D to be even). Plot your answer vs. D. You should see that it is exceedingly unlikely to hit the super-sphere once D gets large.

Percolation

We are now ready to tackle the wildfire problem. Our approach will be to make a simple stochastic (random) model of the system. We will then analyse it using an applet. Along the way we will see many surprises.

A Simple Model of Percolation

Here's the model: divide a square into smaller squares and color each square red with a certain probability, p. For example, we may have a 4x4 grid of squares and color each square red with probability 0.6. We will let the computer do this. In practice this means that the computer generates a random number between zero and one for each square. If the random number is smaller than 'p', the square is marked, otherwise it is left unmarked.

Here is an example of doing this on a 6x6 grid with p=1/2.

Eighteen squares are colored, this is right since we expect about 1/2 of the 36 squares to be marked. It is luck that it is exactly 1/2 in this case! Each square is marked separately so anywhere between 0 and 36 marked squares are possible, with 18 being the most likely. The applet has an algorithm which indicates whether marked squares are next to each other -- in other words it figures which red squares are in connected clumps. For convenience it colors the largest such clump blue. In the figure this is the 9 square clump in the upper right. An example on a grid of 120x120 squares is shown below.

Notice that things get a lot more complicated when there are a lot of squares!

If a connected area goes from top to bottom we say that the system has "percolated". This will be indicated by the applet as a green clump. Here is an example from a 120x120 grid:

The Big Question

Here's what we want to find out: for a given probability of marking each square what is the probability of percolating? Also, how does this change with the grid size? The answer will amaze you!

Before we start, let's sharpen our understanding of the model by looking at some special cases. First, if p=0, then no squares will ever be marked and the probability of percolation will be zero. If p=1 then all the squares will be marked and the probability of percolation will be unity. In between, we expect some intermediate answer.

The 2x2 Case

Now consider a 2x2 grid. It is a fairly simple matter to prove that the probability of percolation is p²(2-p²).

Exercise Prove this

Let's use the computer to check this result. Close the pi window and open the percolation window now. Select a 2x2 grid using the pull down menu on the top right (it should already show 2x2). Make sure that the next pull down menu reads 'percolation' and that the "prob." window reads 0.3. Now click go.

What you'll see is a series of 2x2 grids flashing by (two per second). Each one of these represents the result of one experiment where the squares are marked with probability,p, in this case 0.3. Each one of these is called a configuration. The wide "total number of configurations" box shows how many configurations that the computer will generate. The smaller boxes to the left indicate what is happening. From top to bottom: the number of configurations which have been generated at that time, the number of configurations which have successfully percolated, and the ratio (labelled "stat"). When I ran this I found 3 of 20 configurations had percolated, giving a ratio of 0.15.

If you now hit plot you'll get a little graph. The x-axis is the probability you ran at (0.3 in this case) and the y-axis is the probability of percolation (0.15 in this case). You should see a blue dot at the point (0.3,0.15). This is the experiment you just ran. You should also see a blue curve. This is the analytic result given above. Notice that your point doesn't lie exactly on the curve -- its stochastic! As discussed above, the more configurations you generate the closer you will get to the theoretical curve.

Now let's fill out the whole curve. if you click on auto the applet will automatically cycle through various probabilities and plot the results. If you find the pace too slow, put a smaller number into the "pause" box (this is the pause between configurations measured in milliseconds).

Exercise Increase the number of configurations and generate a new curve. It should be more accurate.

Exercise Think of a method to associate an error with each point you calculate.

Larger Grids and Phase Transitions

What happens as the grid size increases? The limit of large size systems is of interest to physicists because it is a good approximation to real systems which contain huge numbers of degrees of freedom (atoms, molecules, etc). In fact we are often interested in the bulk limit where the system size goes to infinity.

Let's get an idea of what happens as the grid size (we'll call it L for an LxL grid) increases by running the auto experiment for a 4x4 grid. To do this select 4x4 in the top pull-down menu. Compare the 2x2 and 4x4 plots. The 4x4 plot should look steeper in the middle section than the 2x2. Repeat this for larger systems.

Exercise What is happening as L increases? Try to put all of the plots onto a single graph with each curve labelled by L. Speculate on the infinite L shape of the curve.

What you will find is that the curve gets sharper and sharper as L is increased. In fact the infinite L curve is a step function: it is zero up to a special value of p called p_c -- the critical probability. Above p_c the percolation probability is 1. This is an example of a phenomenon called a phase transition. For infinite systems there is a magical probability below which you never percolate and above which you always percolate.

p_c and L-extrapolation

Our task in this section is to find an accurate estimate of p_c. If you run auto for the 120x120 system, you can get a rough estimate of the critical probability. When I do it, it looks like p_c is between 0.5 and 0.6. To do better than this we need to run lots of cases around p=0.55 checking to see if the percolation probability is (near) 0 or 1.

We'd like to look at the percolation probability on an infinite system, but of course this is impossible. The best we can do is watch what happens to the probability as we increase L. Imagine fixing p at 0.5 and calculation the percolation probability for L=2,4,6,8,etc. Then plot the probabilities vs. L and try to see where is is going, extrapolate to L=infinity. Better yet, plot the percolation probability vs. 1/L. Then the L=infinity corresponds to the intercept at the origin.

The applet does this for you: choose a probability (say 0.3) and select "extrapolate" from the second pull down menu. Then click go. You will get a graph of the perc prob vs 4/L for L=4,10,20,40,60. For p=0.3 it is pretty clear that perc prob(L=infinity) = 0. Try it for p=0.7. It is clear that perc prob(L=infinity) = 1.

Exercise See how accurately you can determine p_c. Compare to my answer.

You will find that it is pretty hard to do the extrapolation near p=0.6. This is because you are near the phase transition and the capabilities of stochastic methods break down there. To make progress you may need to use larger lattices (120x120) and run for a lot of configurations.

Exercise Make your own 1/L graph including larger lattices with many configurations.

Systematic Error

Before we move on, there is another subtlety which we shouldn't ignore. Why did I plot the perc probs vs. 4/L? Why not 1/L² or exp(-L)? Sometimes physicists have theoretical reasons to choose one form over another. But often we don't know and have to guess. The error induced in our conclusions this way is called systematic error. There is little we can do about it except try to estimate it reliably. The errors which arise because we are doing numerical experiments with random numbers are called statistical errors.

Exercise Make some extrapolation graphs using 1/L² as a plotting variable. Does it do better or worse than 1/L?

Percolation and Fractals

We now wish to analyse the structure of percolating clusters more carefully. To start, let's take a look at one on a 120x120 grid.

Notice that the percolating cluster (in green) is not very smooth. It is filled with holes and its border wanders all over the place. If you look in the red "hole" at the center bottom, you will see that it is full of holes too. In fact if you blow up any region of this plot (and we had run it for a very large grid size) you'll see a pattern that looks very similar to what you were looking at before. (You may have seen demonstrations of this with something called the Mandelbrot set.) This is an example of the phenomenon of scale invariance which occurs at phase transitions. Scale invariance means that the thing you are looking at looks the same at all magnifications.

The jaggedness of the border of the percolating region may remind you of something else you've probably heard of -- fractals. A fractal is a surface which is so bumpy that it is not really a line and not really an area -- it is something in between. A typical example is a coast line which looks jagged and bumpy no matter what scale you draw the map at.

More to come on the definition of the fractal dimension and an applet demonstration of it

The fractal dimension is expected to be 91/48 = 1.90.

Scaling Laws, Critical Exponents, and Universality

We continue our analysis of the character of percolating clusters by examining the fraction of squares which are in the percolating cluster. We call this fraction F. It turns out that F obeys a simple formula near the critical probability:

Here is a number known as a critical exponent. It tells us how the fraction of squares in the percolating cluster goes to zero as we approach p_c from above. Its value is 5/36 = 0.139.

The applet provides a way to see this scaling law. Select "F calculator" from the second pull down menu. Choose a lattice size and then click auto. The applet will cycle through various probabilities and calculate the average number of spins in each percolating cluster, F. In the bulk limit we expect F to be zero below p_c. If you run on a 40x40 or larger lattice you should confirm this. Above p_c it is less clear what will happen, except that F=1 when p=1. What you will find is that F approaches zero rapidly as p approaches p_c from above. This means that even though a cluster has percolated, the fraction of spins in it is essentially zero at the critical point! Percolating clusters must be very wispy things near p_c. Of course, this is a reflection of their fractal nature.

Exercise How would you measure ?

It turns out that critical exponents have a very general validity. Imagine that we solve a variety of two-dimensional percolation problems, say with a square lattice (as we have been doing), a triangular lattice, and a honeycomb lattice. Now calculate for each theory. You'll find that it doesn't change! This is an example of a phenomenon called universality. What it means in this case is that the structure of the percolating cluster doesn't depend on the details of the model that you use to study percolation. This is a very deep property of physical systems near phase transitions. It underlies a concept called the renormalization group in particle physics and the epsilon-expansion in condensed matter physics. In 1982 Ken Wilson won the Nobel Prize by exploiting these concepts to invent a way to make calculations with phase transitions.

Exercise Think about the percolation problem in two dimensions with a triangular lattice. Do you expect p_c to be larger or smaller than it is for the square lattice? Write a code to simulate it. What fractal dimension does its percolating clusters have?

Exercise Examine percolation for a three dimensional square lattice (a grid of cubes). What happens to p_c now?

Exercise Examine universality by calculating different quantities in several different two-dimensional models. Why isn't p_c universal? Can you think of any other quantities which are universal?

This site is now a part of the Merlot Mathematics Digital Library

Further Reading

• Introduction to Computational Physics, Giardino.
• Introduction to Percolation Theory, D. Stauffer.
• Numerical Recipes in C, W.H. Press et al.