Geometric phase transition


The goal of this program is to demonstrate some of the important properties of percolation, especially near the geometrical phase transition. The key idea is the non-existence of a spanning path for p < pc and the existence of a spanning path for p ≥ pc, where pc is the percolation threshold.

    Percolation is a geometrical model that exhibits a continuous phase transition. The simplest example of a percolation model is formulated on a lattice. Assume that every lattice site can be in one of two states, "occupied" or "empty." Each site is occupied independently of its neighbors with probability p. This model is called site percolation. The nature of percolation is related to the properties of the clusters of occupied sites. Two occupied sites belong to the same cluster if they are linked by a path of nearest-neighbor bonds joining occupied sites.

    We generate a random number for each lattice site and occupy a site if its random number is less than p. Because each site is independent, the order that the sites are visited is irrelevant. If p is small, there are many small clusters. As p is increased, the size of the clusters increases. If p ≅ 1, most of the occupied sites form one large cluster that extends from one end of the lattice to the other. Such a cluster is said to "span" the lattice and is called a spanning cluster. What happens for intermediate values of p, for example between p = 0.5 and p = 0.7? It has been shown that in the limit of an infinite lattice there exists a well defined threshold probability pc such that

    The essential characteristic of percolation is connectedness. Because the connectedness exhibits a qualitative change at a well defined value of p, the transition from a state with no spanning cluster to a state with one spanning cluster is an example of a geometrical phase transition. It is convenient to define an order parameter, a quantity that vanishes for p < pc and is nonzero for p ≥ pc. A convenient choice of the order parameter for percolation is P, the probability that an occupied site is part of the spanning cluster. We can estimate P for a given configuration from its definition:

P = (number of sites in the spanning cluster)/(total number of occupied sites).

    For p < pc there is no spanning cluster and P = 0. (There are configurations that span a lattice for p < pc, such as a column of occupied sites, but these configurations have such a low probability of occurring in the limit of an infinite lattice that they can be ignored.) At p = 1, P has its maximum value of unity because only the spanning cluster exists. These properties suggest that P is a reasonable choice for the order parameter. To calculate P we need to average over all possible configurations for a given value of p.

     The simulation shows the lattice for each trial as well as plots of S(p), the mean cluster size, Pspan, the probability of a spanning cluster, and P as a function of p. The program also plots the mean cluster size distribution ns as a function of cluster size s.


  1. The sites of a lattice are occupied with probability p. (The occupancy of a given site is independent of the occupancy of all other sites.)
  2. The Newman-Ziff algorithm is used to identify the clusters.
  3. The probability of a site being in the spanning cluster, the probability of a spanning cluster, and the mean cluster size are computed as a function of p
  4. The cluster size distribution is computed at the chosen value of p.
  5. Repeat steps 1–4 and average over the trials and plot results.


  1. Run the program and look at the configurations of occupied sites. The program only shows configurations near the value of p inputted by the user. Choose p = pc ≅ 0.5927 and note that clusters appear at many different sizes. Because the system is close to the percolation threshold, there is sometimes a spanning cluster and sometimes not. Increase L and qualitatively describe the distribution of cluster sizes.
  2. After about 100 trials stop the program and look at the cluster size distribution ns(p =pc) as a function of s (a log-log plot). Do you see a linear plot for at least part of the data? What functional form does this linear dependence imply? Use Data Tool menu item under Tools for a linear fit to the data: First generate a log-log plot of the data in Data Tool by using Data Builder and defining functions for log(x) and log(y). To see the defined functions on the plot, drag the columns just created by the Data Builder to the columns labeled horiz and vert (and then "uncheck" the boxes that display the original data as needed). Then, you can select the linear part of the data for a linear curve fit. The data should be in the form ns ∼ Asτ. The exact result for this slope for an infinite system is τ = 187/91 ≅2.05. How does your result for τ compare?
  3. A spanning cluster is defined in the program to be one that touches both the left and right edge of the lattice. In the infinite lattice limit clusters that span horizontally will always span vertically. For L = 128 estimate the value of pc as the value of p for which half of the trials span. How does your estimate compare with the value 0.5927?
  4. Choose L = 128 and do at least 100 trials (1000 is better). Copy the data for the mean cluster size S(p) and the probability P(p) that an occupied site is in the spanning cluster into Data Tool. Use Fit Builder to fit P to the form P ∼ (p - pc)β and mean cluster size S(p) to S(p) ∼ |p-pc| with pc \approx 0.5927. (Alternatively, copy the data into a spreadsheet or plotting program and find the best linear fit to log-log plots of each quantity.) The exact results for γ and β for an infinite lattice is γ = 43/18 and β = 5/36. The quantities τ, β, and γ are examples of critical exponents.
  5. A better way to determine the critical exponents β and γ is to use finite size scaling. The underlying assumption of finite size scaling is that there is only one important length in the system near p = pc, the correlation length ξ ∼ |p - pc|. We write P ∼ (p - pc)β ∼ ξ-β/ν. For a finite system we replace ξ by L and write P ∼ L-β/ν. Similarly, S ∼ Lγ/ν. Hence, a log-log plot of P versus L at p = pc should be linear with a slope of -β/ν. Use the program to generate configurations at p = pc for L = 10, 20, 40, 80, and 160, and determine the ratios β/ν and γ/ν. Use the exact result ν = 4/3, and compare your results with the exact results for β = 5/36 and γ = 43/18.


Java Classes

Updated 28 December 2009.