Percolation Theory with Python Guide
Percolation Theory with Python Guide
Anders Malthe-Sørenssen
Percolation
Theory Using
Python
Lecture Notes in Physics
Founding Editors
Wolf Beiglböck
Jürgen Ehlers
Klaus Hepp
Hans-Arwed Weidenmüller
Volume 1029
Series Editors
Roberta Citro, Salerno, Italy
Peter Hänggi, Augsburg, Germany
Betti Hartmann, London, UK
Morten Hjorth-Jensen, Oslo, Norway
Maciej Lewenstein, Barcelona, Spain
Satya N. Majumdar, Orsay, France
Luciano Rezzolla, Frankfurt am Main, Germany
Angel Rubio, Hamburg, Germany
Wolfgang Schleich, Ulm, Germany
Stefan Theisen, Potsdam, Germany
James Wells, Ann Arbor, MI, USA
Gary P. Zank, Huntsville, USA
The series Lecture Notes in Physics (LNP), founded in 1969, reports new develop-
ments in physics research and teaching - quickly and informally, but with a high
quality and the explicit aim to summarize and communicate current knowledge in
an accessible way. Books published in this series are conceived as bridging material
between advanced graduate textbooks and the forefront of research and to serve
three purposes:
Dr Lisa Scalone
[Link]@[Link]
Anders Malthe-Sørenssen
© The Editor(s) (if applicable) and The Author(s) 2024. This book is an open access publication.
Open Access This book is licensed under the terms of the Creative Commons Attribution 4.0 Inter-
national License ([Link] which permits use, sharing, adaptation,
distribution and reproduction in any medium or format, as long as you give appropriate credit to the
original author(s) and the source, provide a link to the Creative Commons license and indicate if changes
were made.
The images or other third party material in this book are included in the book’s Creative Commons
license, unless indicated otherwise in a credit line to the material. If material is not included in the book’s
Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the
permitted use, you will need to obtain permission directly from the copyright holder.
The use of general descriptive names, registered names, trademarks, service marks, etc. in this publication
does not imply, even in the absence of a specific statement, that such names are exempt from the relevant
protective laws and regulations and therefore free for general use.
The publisher, the authors and the editors are safe to assume that the advice and information in this book
are believed to be true and accurate at the date of publication. Neither the publisher nor the authors or
the editors give a warranty, expressed or implied, with respect to the material contained herein or for any
errors or omissions that may have been made. The publisher remains neutral with regard to jurisdictional
claims in published maps and institutional affiliations.
This Springer imprint is published by the registered company Springer Nature Switzerland AG
The registered company address is: Gewerbestrasse 11, 6330 Cham, Switzerland
This textbook was developed for the course Fys4460 Disordered media and
percolation theory. The course was developed in 2004 and taught every year since at
the University of Oslo. The original idea of the course was to provide an introduction
to basic aspects of scaling theory to a cross-disciplinary group of students. Both
geoscience and physics students have successfully taken the course.
This book follows the underlying philosophy that learning a subject is a hands-
on activity that requires student activities. The course that used the book was
project driven. The students solved a set of extensive computational and theoretical
exercises and were supported by lectures that provided the theoretical background
and group sessions with a learning assistant. The exercises used in the course have
been woven into the text, but are also given as a long project description in an
appendix. This textbook provides much of the same information that was provided
in the lectures.
I believe that in order to learn a subject such as scaling, the student needs to gain
hands-on experience with real data. The student should learn to generate, analyze
and understand data and models. The focus is not to generate perfect data. Instead,
we aim to teach the student how to make sense of imperfect data. The data presented
in the book and the data that students may generate using the supplied programs are
therefore not from very long simulations, but instead from simulations that take a
few minutes on an ordinary computer. The experience from this course has been
that students learn most effectively by being guided through the process of building
models and generating data. Some details of the computer programs have therefore
been provided in the text, and we strive to use a similar notation in the computer
code and in the mathematics in order to make the transfer from mathematics to
computational modeling as simple as possible.
Another aspect of the book is that it tries to be complete in exposition and worked
examples. Not only are the theoretical arguments carried out in detail. The computer
codes needed to generate data are provided in such a form that they can be run and
can generate the data in the examples. This provides students with a complete set of
worked examples that contain theory, modeling (the transfer from theory to model),
implementation, analysis and the resulting connection between theory and analysis.
In the full course, this textbook was only one half of the curriculum. For
the first 10 years the first part of the course focused on random walks and the
last part focused on random growth processes. For the second 10 years of the
v
vi Preface
This textbook was written using doconce—a document translation and format-
ting tool that allows simple integration of text, mathematics and computer code
developed by Hans Petter Langtangen.
1 Introduction to Percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1 Basic Concepts in Percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
1.2 Percolation Probability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
1.3 Spanning Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
1.4 Percolation in Small Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
1.5 Further Reading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2 One-Dimensional Percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.1 Percolation Probability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
2.2 Cluster Number Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.3 Spanning Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2.4 Correlation Length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
3 Infinite-Dimensional Percolation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.1 Percolation Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
3.2 Spanning Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
3.3 Average Cluster Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.4 Cluster Number Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
4 Finite-Dimensional Percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.1 Cluster Number Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
4.2 Consequences of the Scaling Ansatz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.3 Percolation Thresholds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
5 Geometry of Clusters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.1 Geometry of Finite Clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2 Characteristic Cluster Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
5.3 Geometry of the Spanning Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
5.4 Spanning Cluster Above pc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
6 Finite Size Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
6.1 General Aspects of Finite Size Scaling. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
6.2 Finite Size Scaling of P (p, L) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3 Average Cluster Size . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
6.4 Percolation Threshold . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
ix
x Contents
7 Renormalization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
7.1 The Renormalization Mapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102
7.2 Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108
8 Subset Geometry. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
8.1 Singly Connected Bonds . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
8.2 Self-Avoiding Paths on the Cluster . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
8.3 Renormalization Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
8.4 Deterministic Fractal Models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
8.5 Lacunarity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
9 Flow in Disordered Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
9.1 Introduction to Disorder . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135
9.2 Conductivity and Permeability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
9.3 Conductance of a Percolation Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138
9.4 Scaling Arguments for Conductance and Conductivity . . . . . . . . . . . . 147
9.5 Renormalization Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151
9.6 Finite Size Scaling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 152
9.7 Internal Distribution of Currents . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
9.8 Real Conductivity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 158
10 Elastic Properties of Disordered Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
10.1 Rigidity Percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 163
11 Diffusion in Disordered Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
11.1 Diffusion and Random Walks in Homogeneous Media . . . . . . . . . . . . 171
11.2 Random Walks on Clusters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
12 Dynamic Processes in Disordered Media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
12.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 195
12.2 Diffusion Fronts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 196
12.3 Invasion Percolation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 198
References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 209
Index . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 211
Introduction to Percolation
1
In this chapter we motivate the study of disordered media through the example of
a porous system. The basic terms in percolation theory are introduced, and you
learn how to generate, visualize and measure on percolation systems in Python.
We demonstrate how to find exact results for small systems in two dimensions by
addressing all configurations of the system, and show that this approach becomes
unfeasible for large system sizes.
Percolation is the study of connectivity of random media and of other properties
of connected subsets of random media [8, 30, 37]. Figure 1.1 illustrates a porous
material—a material with holes, pores, of various sizes. This is an example of a
random material with built-in disorder. In this book, we will address the physical
properties of such media, develop the underlying mathematical theory and the
computational and statistical methods needed to discuss the physical properties of
random media. In order to do that, we will develop a simplified model system,
a model porous medium, for which we can develop a well-founded mathematical
theory, and then afterwards we can apply this model to realistic random systems.
Fig. 1.1 Illustration of a porous material from a simulation of nanoporous silicate (SiO.2 ). The
colors inside the pores illustrates the distance to the nearest part of the solid
be .Vp = (1 − p)V . We usually call the relative volume of the pores, the porosity,
.φ = Vp /V , of the material. The solid is called the matrix and the relative volume
Fig. 1.2 Illustration of an array of .4 × 4 random numbers, and the set sites for different values of
p
The resulting matrices are shown in Fig. 1.2 for various values of p. The left
figure illustrates the random values in the matrix, z, and the right figures the set sites
for various values of p. You can think of this process as similar a to changing the
water level in a landscape (.zi ) and observing what parts of a landscape is below
water (.zi ≤ p).
from one side to another (from the top to the bottom, from the left to the right, or
both) the percolation threshold. For a given realization of the matrix, there is well-
defined value for .pc , but another realization would give another realization of .pc .
We therefore need to either use statistical averages to characterize the properties of
the percolation system, or we need to refer to a theoretical—thermodynamic—limit,
such as the value for .pc in an infinitely large system. When we use .pc here, we will
usually refer to the thermodynamic value.
In this book, we will develop theories describing various physical properties
of the percolation system as a function of p. We will characterize the sizes of
connected regions, the size of the region connecting one side to another, the size
of the region that contributes to transport (fluid, thermal or electrical transport), and
other geometrical properties of the system. Most of the features we study will be
universal, that is, independent of many of the details of the system. From Fig. 1.2
we see that .pc depends on the details. For example, it depends on the definition of
connectivity. It would also depend on the type of lattice used: square, triangular,
hexagonal, etc. The value of .pc is specific. However, many other properties of
the system are general. For example, how the conductivity of the porous medium
depends on p near .pc does not depend on the type of lattice or the choice of
connectivity rule. It is universal. This means that we can choose a system which
is simple to study in order to gain intuition about the general features, and then
apply that intuition to the special cases afterwards.
While the connectivity or type of lattice does not matter, some things do matter.
For example, the dimensionality matters: The behavior of a percolation system is
different in one, two and three dimensions. However, the most important changes in
behavior occur between one and two dimensions, where the difference is dramatic,
whereas the difference between two and three dimensions is more of a degree that
we can easily handle. Actually, the percolation problem becomes simpler again in
higher dimensions. In two dimensions, it is possible for a path to around a hole
and still have connectivity. But it is not possible to have connectivity of both the
pores and the solid in the same direction at the same time. This is possible in
three dimensions: A two-dimensional creature would have problems with having
a digestive tract, as it would divide the creature in two, but in three dimensions
this is fully possible. Here, we will therefore focus on two and three-dimensional
systems.
We will first address percolation in one and infinite dimensions, since we can
solve the problems exactly in these cases. We will then address percolation in two
dimensions, where there are no exact solutions. However, we will see that if we
assume that the distribution of cluster sizes has a particular scaling form, we can
still address the problem in two dimensions and make powerful predictions. We will
also see that close to the percolation threshold the porous medium has a self-affine
scaling structure—it is a fractal. This property has important consequences for the
1.1 Basic Concepts in Percolation 5
physical properties of random systems. We will also see how this is reflected in a
systematic change of scales, a renormalization procedure, which is a general tool
that can applied to rescaling in many areas.
The resulting matrix is illustrated in Fig. 1.3. However, this visualization does not
provide us with any insight into the connectivity of the sites in this system. Let us
instead analyze the connected regions in the system.
Definitions
• two sites are connected if they are nearest neighbors (there are 4 nearest
neighbors on a square lattice)
• a cluster is a set of connected sites
• a cluster is spanning if it spans from one side to the opposite side
• a cluster that is spanning is called the spanning cluster
• a system is percolating if there is a spanning cluster in the system
Fortunately, there are built-in functions in python that finds connected regions
in an image.1 The function [Link] finds clusters based on a
given connectivity. For example, with a coordination number .Z = 4, that is nearest
neighbor connectivity, we find
from [Link] import measurements
lw, num = [Link](m)
1 Notice
that these program lines will need the libraries numpy, [Link] and
measurements to be loaded. We will assume you have done this in the following.
6 1 Introduction to Percolation
Fig. 1.3 Illustration of the index array for a .10 × 10 system for .p = 0.45
This function returns the matrix lw, which for each site in the original array tells
what cluster it belongs to. Clusters are numbered sequentially, and each cluster is
given an index number. All the sites with the same index number belong to the
same cluster. The resulting array is shown in Fig. 1.3, where the index number for
each site is shown and a color is used to indicate the clusters. Notice that there is a
distribution of cluster sizes, but no cluster is large enough to reach from one side to
another, and as a result the system does not percolate.
In order to visualize the clusters effectively, we give the various clusters different
colors.
[Link](lw, origin=’lower’)
Unfortunately, this colors the clusters gradually from the bottom up. This is a
property of the underlying algorithm: Clusters are indexed starting from the top-
left of the matrix (which is the bottom-left of the image). Hence, clusters that are
close to each other will get similar colors and can therefore be difficult to discern
unless we shuffle the colormap. We can fix this by shuffling the labeling:
b = [Link]([Link]() + 1)
[Link](b)
shuffledLw = b[lw]
[Link](shuffledLw, origin=’lower’)
The resulting image is shown to the right in Fig. 1.3. (Notice that in these figures
we have reversed the ordering of the y-axis. Usually, the first row is in the top-left
corner in your plots, but when we use the keyword lower the first row is in the
bottom-left).
It may also be useful to color the clusters based on the size of the clusters, where
size refers to the number of sites in a cluster. We can do this using
area = [Link](m, lw, index=[Link]([Link]() + 1))
areaImg = area[lw]
1.1 Basic Concepts in Percolation 7
[Link](areaImg, origin=’lower’)
[Link]()
Let us now study the effect of p on the set of connected clusters. We vary the
value of p for the same underlying random matrix, and plot the resulting images:
import numpy as np
import [Link] as plt
from [Link] import measurements
[Link](figsize=(10,8))
L = 100
pv = [0.2,0.3,0.4,0.5,0.6,0.7]
z = [Link](L,L)
for i in range(len(pv)):
p = pv[i]
m = z<p
lw, num = [Link](m)
area = [Link](m, lw, index=[Link]([Link]() + 1))
areaImg = area[lw]
[Link](2,3,i+1)
tit = ’p=’+str(p)
[Link](areaImg, origin=’lower’)
[Link](tit)
Figure 1.4 shows the clusters for a .100×100 system for p ranging from .0.2 to .0.7
in steps of .0.1. We see that the clusters increase in size as p increases. At .p = 0.6
there is one large cluster spanning the entire region. We have a percolating cluster,
and we call the cluster that spans the system the spanning cluster. The transition is
very rapid from .p = 0.5 to .p = 0.6. We therefore look at this region in more detail
Fig. 1.4 Plot of the clusters in a .100 × 100 system for various values of p
8 1 Introduction to Percolation
Fig. 1.5 Plot of the clusters in a .100 × 100 system for various values of p
in Fig. 1.5. We see that the size of the largest cluster increases rapidly as p reaches a
value around .0.6, which corresponds to .pc for this system. At this point, the largest
cluster spans the entire system. For the two-dimensional system illustrated here we
know that in an infinite lattice the percolation threshold is .pc ≃ 0.5927.
The aim of this book is to develop a theory to describe how this random porous
medium behaves close to .pc . We will characterize properties such as the density of
the spanning cluster, the geometry of the spanning cluster, and the conductivity and
elastic properties of the spanning cluster. We will address the distribution of cluster
sizes and how various parts of the clusters are important for particular physical
processes. We start by characterizing the behavior of the spanning cluster near .pc .
When does the system percolate? When there exists a path connecting one side to
another. This occurs at some value .p = pc . However, in a finite system, like the
system we simulated above, the value we find for .pc will vary with each realization.
It may be slightly above or slightly below the .pc we find in an infinite sample. Later,
we will develop a theory to understand how the effective .pc in a finite system varies
from the thermodynamic .pc . But already now we realize that as we perform different
numerical experiments, we will measure various values of .pc . We can characterize
this behavior by introducing a probability .Π (p, L):
are present both on the left and on the right side of the lattice. In that case, they are
spanning (We could also have included a test for clusters spanning from the top to
the bottom, but this does not change the statistics significantly). In this case, there is
a spanning cluster—the system percolates. We count how many times, .Ni , a system
percolates for a given .pi and then divide by the total number of experiment, N,
to estimate the probability for percolation for a given .pi , .Π (pi , L) ≃ Ni /N . We
implement this as follows. First, we generate a sequence of 100 .pi values from .0.35
to .1.0:
p = [Link](0.35,1.0,100)
Then we prepare an array for .Ni with the same number of elements as .pi :
nx = len(p)
Pi = [Link](nx)
We will then loop over all samples, and for each sample we generate a new random
matrix. The for each value of .pi we perform the cluster analysis as we did above.
We use the function [Link] to label the clusters. Then we find
the intersection between the labels on the left and the right side of the system and
store in perc_x. If the length of the set of intersections is larger than zero, there
is at least one percolating cluster, and we find the label of the spanning cluster(s) in
perc:
lw,num = [Link](z)
perc_x = np.intersect1d(lw[0,:],lw[-1,:})
perc = perc_x[[Link](perc_x>0)]
Now, we are ready to implement this into a complete program. For a given value
of p, we count in how many simulations .Np (p) there is a path spanning from one
side to another and estimate .Π̄ (p) ≃ Np (p)/N , where N is the total number of
simulations/samples. This is implemented in the following program:
import numpy as np
import [Link] as plt
from [Link] import measurements
p = [Link](0.4,1.0,100)
nx = len(p)
Ni = [Link](nx)
P = [Link](nx)
N = 1000
L = 100
for i in range(N):
z = [Link](L,L)
for ip in range(nx):
m = z<p[ip]
lw, num = [Link](m)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
10 1 Introduction to Percolation
perc = perc_x[[Link](perc_x>0)]
if (len(perc)>0):
Ni[ip] = Ni[ip] + 1
Pi = Ni/N
[Link](p,Pi)
[Link](’$p$’)
[Link](’$\Pi$’)
The resulting plot of .Π (p, L) is seen in Fig. 1.6. The figure shows the resulting
plots as a function of system size L. We see that as the system size increases,
.Π (p, L) approaches a step function at .p = pc .
The probability .Π (p, L) describes the probability for there to be a spanning cluster,
but what about the spanning cluster itself, how can we characterize it? We see from
Fig. 1.4 that the spanning cluster grows quickly around .p = pc . Let us therefore
characterize the cluster by its size, .MS , or by its density, .P (p, L) = MS /L2 , which
corresponds to the probability for a site to belong the spanning cluster.
1.0
L = 50
0.8 L = 100
L = 200
0.6
Π(p, L )
0.4
0.2
0.0
0.4 0.5 0.6 0.7 0.8 0.9 1.0
p
1.0
L = 50
0.8 L = 100
L = 200
0.6
P (p, L )
0.4
0.2
0.0
0.4 0.5 0.6 0.7 0.8 0.9 1.0
p
Fig. 1.6 Plot of .Π (p, L), the probability for there to be a connected path from one side to anther,
as a function of p for various system sizes L, and .P (p, L), the density of the spanning cluster
1.3 Spanning Cluster 11
Density of the Spanning Cluster The probability .P (p, L) for a site to belong
to a spanning cluster is called the density of the spanning cluster.
We can measure .P (p, L) by counting the mass .Mi of the spanning cluster as
a function of .pi for various values of .pi . We can find the mass of the spanning
cluster by finding a cluster that spans the system (there may be more than one) as
we did above, and then measure the number of sites in the cluster using area =
[Link](m, lw, perc).
We do this in the same program as we developed above. For each .pi , we see if
a cluster is spanning from one side to another, and if it is, we add the mass of this
cluster to .MS (pi ). We implement these features in the following program, which
measures both .Π (p, L) and .P (p, L) for a given value of L:
import numpy as np
import [Link] as plt
from [Link] import measurements
p = [Link](0.4,1.0,100)
nx = len(p)
Ni = [Link](nx)
P = [Link](nx)
N = 1000
L = 100
for i in range(N):
z = [Link](L,L)
for ip in range(nx):
m = z<p[ip]
lw, num = [Link](m)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x>0)]
if (len(perc)>0):
Ni[ip] = Ni[ip] + 1
area = [Link](m, lw, perc[0])
P[ip] = P[ip] + area
Pi = Ni/N
P = P/(N*L*L)
[Link](2,1,1)
[Link](p,Pi)
[Link]("$\Pi(p)$")
[Link](2,1,2)
[Link](p,P)
[Link]("P(p)")
[Link]("p")
The resulting plot of .P (p, L) is shown in the bottom of Fig. 1.6. We see that
P (p, L) changes rapidly around .p = pc , and that it grows slowly—approximately
.
sites that are present, which again is proportional to p. We will now develop a theory
for the observations of .Π (p, L), .P (p, L) and other features of the percolation
system. First, we see what insights we can gain from small, finite systems.
develop a more advanced strategy than for .L = 1. Our strategy will be to list all
possible outcomes, find the probability for each outcome, and then use this to find
the probability for the various physical properties we are interested in. The possible
configurations are shown in Fig. 1.7.
Our plan is to use a basic result from probability theory: If we want to calculate
the probability of an event A, we can do this by summing the probability of A given
B multiplied by the probability for B over all possible outcomes B (as long as the
set of outcomes B span the space of all outcomes and are mutually exclusive, that
is, that they have no intersection). In this case:
P (A) =
. P (A|B)P (B) , (1.1)
B
where we have used the notation .P (A|B) to denote the conditional probability of
A given that B occurs. We can use this to calculate properties such as .Π (p, L) and
.P (p, L) by summing over all possible configurations c of the system:
Π (p, L) =
. Π (p, L|c)P (c) , (1.2)
c
where .Π (p, L|c) is the value of .Π for the particular configuration c, and .P (c) is
the probability of this configuration.
The configurations for .L = 2 have been numbered from .c = 1 to .c = 16
in Fig. 1.7. However, configurations that are either mirror images or rotations of
each other will have the same probability and the same physical properties since
percolation can take place both in the x and the y directions. It is therefore only
necessary to group the configurations into 6 different classes, k, as illustrated in the
bottom of Fig. 1.7, but we then need to include the multiplicity, .gk , for each class
1.4 Percolation in Small Systems 13
Fig. 1.7 The possible configurations for a .L = 2 site percolation lattice in two-dimensions. The
configurations are indexed using the cluster configuration number c
Table 1.1 lists the classes, the number of configurations in each class, the probability
of one such configuration in a class, and the value of .Π (p, L|k) for this class.
We should check that we have actually listed all possible configurations. In
2
general, the number of configurations for an .L × L system is .2L . The total number
of configurations is .1 + 4 + 2 + 4 + 4 + 1 = 16, which is equal to .24 as it should.
We have therefore included all the configurations.
1
L=1
0.8 L=2
L=3
L=4
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
1
L=1
0.8 L=2
L=3
0.6 L=4
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
We can then find the probability for .Π by direct calculation of the sum:
which we can simplify further if we want. The shape of .Π (p, L) for .L = 1, and
L = 2 is shown in Fig. 1.8.
.
These two results gives us an indication about how the percolation probability
Π (p, L) is approaching the step function when .L → ∞.
.
There are good general introduction texts to percolation theory such as the popular
books by Stauffer and Aharony [37], by Sahimi [30], by Christensen and Moloney
[8], and the classical book by Grimmet [14]. Mathematical aspects are addressed
by Kesten [21] and phase transitions in general are introduced by e.g. Stanley [35].
Applications of percolation theory are found in many fields such as in geoscience
[22], porous media [18] or social networks [33] and many more. We encourage you
to explore these books for a more theoretical introduction to percolation theory.
Exercises
Exercise 1.3 (Percolation in Small Systems in 3d) In this exercise we will study
the three-dimensional site percolation system for small system sizes.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
One-Dimensional Percolation
2
The percolation problem can be solved exactly in two limits: in the one-dimensional
and the infinite dimensional cases. Here, we will first address the one-dimensional
system. While the one-dimensional system does not allow us to study the full com-
plexity of the percolation problem, many of the concepts and measures introduced
to study the one-dimensional problem can generalized to higher dimensions.
Π (p, L) = pL .
. (2.1)
This shows that the percolation threshold is .pc = 1 in one dimension. However,
the one-dimensional system is anomalous, and in higher dimensions we will
always have .pc < 1, so that we can study the system both above and below .pc .
Unfortunately, for the one-dimensional system we can only study the system below
.pc .
In the simulations in Fig. 1.4 we saw that the percolation system was characterized
by a wide distribution of clusters—regions of connected sites. The clusters have
varying shape and size. If we increase p to approach .pc we saw that the clusters
increased in size until they reached the system size. We can use the one-dimensional
system to learn more about the behavior of clusters as p approaches .pc .
Figure 2.1 illustrates a realization of an .L = 16 percolation system in one
dimension below .pc = 1. In this case there are 5 clusters of sizes: 1,1,4,2,1
measured as the number of sites in each cluster. The clusters are numbered, indexed,
from 1 to 5 as we did for the numerical simulations in two dimensions. How can we
characterize the clusters in a system? In percolation theory we characterize cluster
sizes by asking a particular question: If you point at a (random) site in the lattice,
what is the probability for this site to belong to a cluster of size s?
It is common to use the notation .sn(s, p) for this probability for a given site to
belong to a cluster of size s. Why is it divided into two parts, s and .n(s, p)? Because
we must divide the question into two parts: (1) What is the probability for a given
site to be a specific site in a cluster of size s, and (2) how many such specific sites
are there in a cluster? What do we mean by a specific site? For the cluster with
index 3 in Fig. 2.1 there are 4 sites. We could therefore ask the question: What is
the probability for a site to be the left-most site in a cluster of size s? This is what
1 2 3 3 3 3 4 4 5
Empty sites
Fig. 2.1 Realization of a .L = 16 percolation system in one dimension. Separate clusters are
illustrated by the indexes, shown as numbers inside the sites. Occupied sites are marked with gray
squares
2.2 Cluster Number Density 19
we mean with a specific site. We could ask the same question about the second left-
most, the third left-most and so on. We call the probability for a site to belong to a
specific site in a cluster of size s (such as the left-most site in the cluster) the cluster
number density, and we use the notation .n(s, p) for this. To find the probability
.sn(s, p) for a site to belong to any of the s sites in a cluster of size s we must sum
the probabilities for each of the specific sites. This is illustrated for the case of a
cluster of size 4:
because each of these probabilities are the same. What is the probability for a site to
be the left-most site in a cluster of size s in one dimension? In order for it to be in
a cluster of size s, the site must be present, which has probability p, and then .s − 1
sites must also be present to the right of it, which has probability .ps−1 . In addition,
the site to the left must be empty (illustrated by an X in Fig. 2.1 bottom part), which
has probability .(1 − p) and the site to the right of the fourth site (illustrated by an
X in Fig. 2.1 bottom part), which also has probability .(1 − p). Since the occupation
probabilities for each site are independent, the probability for the site to be the left-
most site in a cluster of size s is the product of these probabilities:
Cluster Number Density The cluster number density .n(s, p) is the probabil-
ity for a site to be a particular site in a cluster of size s. For example, in one
dimension, .n(s, p) can be interpreted as the probability for a site to be the
left-most site in a cluster of size s.
Let us check that this is indeed the case for the one-dimensional expression for
n(s, p) by calculating the sum:
.
∞
∞
∞
. sn(s, p) = sps (1 − p)2 = (1 − p)2 p sps−1 , (2.6)
s=1 s=1 s=1
which gives
∞
∞
. sn(s, p) = (1 − p)2 p sps−1 = (1 − p)2 p (1 − p)−2 = p . (2.8)
s=1 s=1
Since .P = 0 when .p < 0 we see that the probability is normalized. We can use
similar tricks to calculate moments of any order.
In order to gain further insight into the distribution of cluster sizes, let us study
Fig. 2.1 in more detail. There are 3 clusters of size .s = 1, one cluster of size .s = 2,
and one cluster of size .s = 4. We could therefore introduce a histogram of cluster
sizes, which is what we would do if we studied the cluster distribution numerically.
Let us write .Ns as the number of clusters of size s so that .N1 = 3, .N2 = 1, .N3 = 0
and .N4 = 1.
How can we estimate .sn(s, p), the probability for a given site to be part of a
cluster of size s, from .Ns ? The probability for a site to belong to cluster of size
s can be estimated by the number of sites belonging to a cluster of size s divided
by the total number of sites. The number of sites belonging to a cluster of size s
2.2 Cluster Number Density 21
is .sNs , and the total number of sites is .Ld , where L is the system size and d is
the dimensionality. (Here, .d = 1). This means that we can estimate the probability
.sn(s, p) from
sNs
sn(s, p) =
. , (2.9)
Ld
where we use a bar to show that this is an estimated quantity and not the actual
probability. We divide by s on both sides, and find
Ns
n(s, p) =
. . (2.10)
Ld
This argument and the result are valid in any dimension, not only for .d = 1. We
can also see why this quantity is called the cluster number density: it is the number
of clusters divided by the volume measured in number of sites. We have therefore
found a method to estimate the cluster number density:
Ns
n(s, p) =
. . (2.11)
Ld
N1 3
n(1, p) =
. = , (2.12)
L1 16
N2 1
n(2, p) =
.
1
= , (2.13)
L 16
N3 0
n(3, p) =
. = , (2.14)
L1 16
N4 1
n(4, p) =
.
1
= , (2.15)
L 16
22 2 One-Dimensional Percolation
which is our estimate of .n(s, p) based on this single realization. We check the
consistency of the result by ensuring that the estimated probabilities also are
normalized:
3 1 1 9
. sn(s, p) = 1 · +2· +3·0+4· = =p, (2.16)
s
16 16 16 16
where .p is estimated from the number of present sites divided by the total number
of sites.
In order to produce good statistical estimates for .n(s, p), we must sample from
many random realization of the system. If we sample from M realizations, and
then measure the total number of clusters of size s, .Ns (M), summed over all the
realizations, we estimate the cluster number density from
Ns (M)
n(s, p) =
. . (2.17)
MLd
Notice that all simulations are for finite L, and we would therefore expect deviations
due to a finite L as well as due to the finite number of samples. However, we expect
the estimated .n(s, p; L) to approach the underlying .n(s, p) as M and L approaches
infinity.
In Fig. 2.2 we have plotted .n(s, p) for various values of p. In order to compare the
s-dependence of the plot directly for various p-values we plot
50
0
-50
-100
-150
-200 p=0.900
p=0.990
-250 p=0.999
-300
0 1 2 3 4 5 6
50
0
-50
-100
-150
-200 p=0.900
p=0.990
-250 p=0.999
-300
-3 -2 -1 0 1 2 3
Fig. 2.2 (Top) A plot of .n(s, p)(1 − p)2 as a function of s for various values of p for a one-
dimensional percolation system shows that the cut-off increases as a function of s. (Bottom) When
the s axis is rescaled by .sξ to .s/sξ , all the curves fall onto a common scaling function, that is,
.n(s, p) = (1 − p) F (s/sξ )
2
What we are seeing in Fig. 2.2 is therefore the exponential cut-off curve, where the
cut-off .sξ (p) increases as .p → 1. We call it a cut-off because the value of .n(s, p)
decays very rapidly (exponentially fast) when s is larger than .sξ .
How does .sξ depend on p?. We see from (2.21) that as p approaches .pc = 1,
the characteristic cluster size .sξ will diverge. The form of the divergence can be
determined in more detail through a Taylor expansion:
1
sξ = −
. (2.22)
ln p
where we have used that .ln(1 − x) = −x + O(x 2 ), which is simply the Taylor
expansion of the logarithm, where .O(x 2 ) is term that is on the order of .x 2 . As a
result
1 1
sξ ≃
. = = |p − pc |−1 . (2.24)
1−p pc − p
This shows that the divergence of .sξ as p approaches .pc is a power-law with
exponent .−1. This power-law behavior is general in percolation theory:
sξ ∝ |p − pc |−1/σ ,
. (2.25)
The value of the exponent .σ depends on the lattice dimensionality, but it does
not depend on the details of the lattice. It would, for example, be the same also for
next-nearest neighbor connectivity.
The functional form we have found is also an example of a data collapse. We
see that if we plot .(1 − p)−2 n(s, p) as a function of .s/sξ , all data-points for various
values of p should fall onto a single curve:
as illustrated in Fig. 2.2. We call this a data-collapse. We have one behavior for small
s and then a rapid cut-off when s reaches .sξ . We can rewrite .n(s, p) so that all the
.sξ dependence is in the cut-off function by realizing that since .sξ ≃ (1 − p)
−1 we
−2
have that .(1 − p)2 = sξ . This gives
2
s s
. n(s, p) = sξ−2 e−s/sξ = s −2 e−s/sξ = s −2 F . (2.27)
sξ sξ
where .F (u) = u2 e−u . We will see later that this form for .n(s, p) is general—it is
valid for percolation in any dimension, although with other values for the exponent
.−2 and other shapes of the cut-off function .F (u). In percolation theory, we call this
exponent .τ :
n(s, p) = s −τ F (s/sξ ) ,
. (2.28)
2.2 Cluster Number Density 25
Let us now test the measurement method and the theory through a numerical study
of the cluster number density. According to the theory developed above we can
estimate the cluster number density .n(s, p) from
Ns (M)
n(s, p) =
. , (2.29)
L2 M
where .Ns (M) is the number of clusters of size s measured in M realizations of the
percolation system. We generate a one-dimensional percolation system and index
the clusters using
import numpy as np
from [Link] import measurements
L = 20
p = 0.90
z = [Link](L)
m = z<p
lw, num = [Link](m)
Now, lw contains the indices for all the clusters. We can extract the size of the
clusters by summing the number of elements for each label:
labelList = [Link]([Link]() + 1)
area = [Link](m, lw, labelList)
We need to collect all the areas of all the clusters for many realizations, and then
calculate the number of clusters of each size s based on this long list of areas. This
is all brought together by continuously appending the area-array to the end of an
array allarea that contains the areas of all the clusters.
import numpy as np
import [Link] as plt
from [Link] import measurements
nsamp = 1000
L = 1000
p = 0.90
26 2 One-Dimensional Percolation
allarea = [Link]([])
for i in range(nsamp):
z = [Link](L)
m = z<p
lw, num = [Link](m)
labelList = [Link]([Link]() + 1)
area = [Link](m, lw, labelList)
allarea = [Link](allarea,area)
n,sbins = [Link](allarea,bins=int(max(allarea)))
s = 0.5*(sbins[1:]+sbins[:-1])
nsp = n/(L*nsamp)
sxi = -1.0/[Link](p)
nsptheory = (1-p)**2*[Link](-s/sxi)
[Link](s,nsp,’o’,s,nsptheory,’-’)
[Link](’$s$’)
[Link](’$n(s,p)$’)
This script also calculates .Ns using the histogram function with L bins to
ensure that there is at least one bin for each value of s:
n,sbins = [Link](allarea,bins=int(max(allarea)))
s = 0.5*(sbins[1:]+sbins[:-1])
For comparison with theory, we calculate values from the theoretically predicted
expression .n(s, p), which is .n(s, p) = (1 − p)2 exp(−s/sξ ), where .sξ = −1/ ln p.
This is calculated for the same values of s as used for the histogram using:
sxi = -1.0/[Link](p)
nsptheory = (1-p)**2*[Link](-s/sxi)
When we use the histogram-function with many bins, we risk that many of the
bins contain zero elements. To remove these elements from the plot, we can use the
nonzero function from numpy to find the indices of the elements of n that are
non-zero:
i = [Link](n)
And then we only plot the values of .n(s, p) at these indices. The values for the
theoretical .n(s, p) are calculated for all values of s, and the two are plotted in the
same plot:
[Link](s[i],nsp[i],’o’,s,nsptheory,’-’)
The resulting plot is shown in Fig. 2.3. We see that the measured results and the
theoretical values fit nicely, even though the theory is for an infinite system size, and
the simulations where performed at .L = 1000. We also see that for larger values of
s there are fewer observed values. It may therefore be a good idea to make the bins
2.2 Cluster Number Density 27
0.01
0.008
0.006
0.004
0.002
0
0 20 40 60 80 100 120
10-2
10-4
10-6
10-8
0 20 40 60 80 100 120
Fig. 2.3 Plot of the predicted .n(s, p), based on .M = 1000 samples of a .L = 1000 system with
.p = 0.9, and the theoretical .n(s, p) curve on a linear scale (top) and a semilogarithmic scale
(bottom). The semilogarithmic plot shows that .n(s, p) follows an exponential curve
used for the histogram larger for larger values of s. We will return to this when we
measure the cluster number density in two-dimensional systems in Chap. 4.
Since we have an exact expression for the cluster number density, .n(s, p) we can use
it to calculate the average cluster size. However, what do we mean by the average
cluster size in this case? In percolation theory it is common to define the average
cluster size as the average size of a cluster connected to a given (random) site in
our system. That is, we will use the cluster number density, .n(s, p), as the basic
distribution for calculating the moments.
1 2 (1 − p)2 2 s
.S= s n(s, p) = s p. (2.32)
p s p s
(1 − p)2 d d (1 − p)2 d d s
= p p ps = p p p. (2.33)
p s
dp dp p dp dp s
(1 − p)2 d d 1 d p
= p p = (1 − p)2 . (2.34)
p dp dp 1 − p dp (1 − p)2
1 2p 1+p
= (1 − p)2 ( + )= , (2.35)
(1 − p)2 (1 − p)3 1−p
where we have used the trick introduced in (2.7) to move the derivation out through
the sum.
This shows that we can write
1+p Γ
S=
. = , (2.36)
1−p |p − pc |γ
with .γ = 1 and .Γ (p) = 1 + p. That is, the average cluster size, S, also diverges as
a power-law when p approaches .pc . The exponent .γ = 1 of the power-law is again
universal. That is, it depends on features such as dimensionality, but not on details
such as the lattice structure.
The density of the spanning cluster, .P (p; L), can be found using similar approaches.
The spanning cluster only exists for .p ≥ pc . The discussion for .P (p; L) is therefore
not that interesting for the one-dimensional case. However, we can still introduce
some of the general notions.
The behavior of .P (p; ∞) for .L → ∞ in one dimension is
0 when p < 1
P (p; ∞) =
. . (2.37)
1 when p = 1
What is the relation between .P (p; L) and the distribution of cluster sizes? The
distribution of the sizes of finite clusters is described by .sn(s, p), which is the
2.4 Correlation Length 29
which is valid for percolation in any dimension, since we have not assumed anything
about the dimensionality in this argument.
We can use this relation to find the density of the spanning cluster from the cluster
number density .n(s, p) through
P (p) = p −
. sn(s, p) . (2.39)
s
This illustrates that the cluster number density .n(s, p) is a fundamental property,
which can be used to deduce many aspects of the percolation system.
From the simulations in Fig. 1.4 we see that the size of the clusters increases as .p →
pc . We expect a similar behavior for the one-dimensional system. We have already
seen that the mass (or area) of the clusters diverges as .p → pc . The characteristic
cluster size .sξ characterizes the mass (or area) of a cluster. How can we characterize
the extent of a cluster?
To characterize the linear extent of a cluster, we find the probability for two
sites at a distance r to belong to the same cluster. This probability is called the
correlation function, .g(r):
For one-dimensional percolation, two sites a and b are part of the same cluster if
and only if all the points in between a and b are occupied. If r denotes the number
of points between a and b (not counting the start and end positions) as illustrated in
30 2 One-Dimensional Percolation
r
a b
L
Fig. 2.4 An illustration of the distance r between two sites a and b. The two sites a and b are
connected if and only if all the sites between a and b are occupied
g(r) = pr = e−r/ξ ,
. (2.40)
where .ξ = − ln1p is called the correlation length. The correlation length diverges
as .p → pc = 1. We can again find the way in which it diverges from a Taylor
expansion in .(1 − p) when .p → 1
ξ = ξ0 (pc − p)−ν ,
. (2.42)
Exercises
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Infinite-Dimensional Percolation
3
If we start from the center in Fig. 3.1 and move along a branch, we find .(Z − 1) new
neighbors from each of the branches. To get a spanning cluster, we need to ensure
that at least one of the .Z − 1 sites are occupied on average. That is, the occupation
(a) (b)
Branch 1 t
oin
al p
entr
C
Branch 2 Branch 3
Fig. 3.1 Four generations of the Bethe lattice with coordination number .Z = 3. (a) Illustration of
a two-dimensional embedding of the lattice. (b) Illustration of how a central node is connected to
three branches
.p(Z − 1) ≥ 1 , (3.1)
1
pc =
. . (3.2)
Z−1
We will use a standard trick to find the density .P (p) of the spanning cluster when
p > pc . The technique is based on starting from a “central” site, and then address
.
clusters with probability . s sn(s, p). This means that
p=P +
. sn(s, p) , (3.3)
s
The sum . s sn(s, p) is the probability that the site is part of a finite cluster,
which means that it is not connected to infinity. We introduce .Q(p) to denote the
probability that a branch does not lead to infinity. The concept of a central point and
a branch is illustrated in Fig. 3.1.
Deriving an Equation for .Q(p) If the probability that a site is not connected
to infinity in a particular direction is Q, then the probability that the site is
not connected to infinity in any direction is .QZ . The probability that the site is
connected to infinity is therefore .1 − QZ . In addition, we need to include the
probability p that the site is occupied. The probability that a given site is connected
to infinity, that is, that it is part of the spanning cluster, is therefore
P = p(1 − QZ ) .
. (3.4)
.Q
Z−1 . The probability Q for the branch not to be connected to infinity is therefore
Q = (1 − p) + pQZ−1 .
. (3.5)
Solving to Find .Q(p) We could solve this equation for general Z. However, for
simplicity we will restrict ourselves to .Z = 3, which is the smallest Z that gives a
behavior different from the one-dimensional system. We recall from (3.2) that for
.Z = 3, .pc = 1/2. We insert .Z = 3 in (3.5), which gives
Q = 1 − p + pQ2 ,
. (3.6)
pQ2 − Q + 1 − p = 0 .
. (3.7)
36 3 Infinite-Dimensional Percolation
First, we notice that for .(2p − 1) = 0, that is for .p = 1/2 (which is .pc for .Z = 3),
we find that .Q = 1/2p = 1. Therefore, no branch propagates to infinity for .p =
pc = 1/2. Second, for .(2p − 1) < 0, which corresponds to .p < 1/2 = pc , we get
the solution
1 ± (1 − 2p)
.Q= , (3.9)
2p
which has two solutions: .Q = 1 or .Q = (1 − p)/p. The second solution .(1 − p)/p
is greater than 1 when .p < 1/2. We therefore conclude that for .p < 1/2 and for
.p = 1/2, no branch propagates to infinity. This means that there is no infinite cluster
(no spanning cluster) for .p ≤ 1/2 = pc . Third, for .(2p−1) > 0, which corresponds
to .p > 1/2 = pc , we get the solution
1 ± (2p − 1)
.Q= , (3.10)
2p
Finding .P (p) We insert .Q = (1 − p)/p back into the equation for .P (p) to find
the behavior of .P (p) for .p > pc = 1/2:
1−p 3
P = p(1 − Q3 ) = p(1 − (
. ) ). (3.11)
p
This result is illustrated in Fig. 3.2. We notice that when .p → 1, we have that
.(1 − p)/p → 0. Therefore .P ∝ p in this limit. To address the behavior when
.p → pc = 1/2, we use the identity .(1 − a ) = (1 − a)(1 + a + a ) and that
3 2
1−p 3 1−p 1−p 1−p 2
.P = p(1 − ( ) )=p 1− 1+ + . (3.12)
p p p p
3.2 Spanning Cluster 37
0.8
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
50
40
30
20
10
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig. 3.2 (Top) A plot of .P (p) as a function of p for the Bethe lattice with .Z = 3. The tangent
at .p = pc is illustrated by a straight line. (Bottom) A plot of the average cluster size, .S(p), as a
function of p for the Bethe lattice with .Z = 3. The average cluster size diverges when .p → pc =
1/2 both from below and above
P ≃ 6(p − pc ) + O (p − pc )2
. . (3.15)
38 3 Infinite-Dimensional Percolation
Since we are only interested in the leading term, we have found that for .p > pc =
1/2 we can approximate .P (p) as:
P (p) ≃ B(p − pc )β ,
. (3.16)
We will use a similar method to find the average cluster size, .S(p). We start by
defining .T (p) as the average number of sites connected to a given site on a specific
branch, such as in branch 1 in Fig. 3.1. The average cluster size S is then given as
S = 1 + ZT ,
. (3.17)
where the 1 represents the central point, T is the average number of sites on each
branch and Z is the number of branches. We will again attempt to find a self-
consistent solution for T , starting from a center site. The average cluster size T
is found from summing the probability that the next site k is empty, .(1 − p),
multiplied with the contribution to the average, in this case (0), plus the probability
that the next site is occupied, p, multiplied with the contribution in this case, which
is the contribution from the site (1) and the contribution of the remaining .Z − 1
subbranches. In total:
T = (1 − p) 0 + p (1 + (Z − 1)T ) ,
. (3.18)
This expression diverges when .1 − p(Z − 1) = 0, that is, for .p = 1/(Z − 1), which
we recognize as .pc = 1/(Z − 1). We insert this in (3.17) as find that S is
1+p 1+p pc (1 + p)
S = 1 + ZT =
. = = , (3.20)
1 − (Z − 1)p 1 − ppc pc − p
3.4 Cluster Number Density 39
which is illustrated in Fig. 3.2. This is a special case for .S(p), which in general can
be written on the general form
Γ
S=
. . (3.21)
(pc − p)γ
For the Bethe lattice with coordination number Z, we have found that .pc = 1/(Z −
1) and that the exponent is .γ = 1. where our argument determines .pc = 1/(Z − 1),
and the exponent .γ = 1. The average cluster size S therefore diverges as a power-
law when p approaches .pc . The exponent .γ characterizes the behavior. This is a
general results. The value of .γ depends on the dimensionality, but not on the details
of the lattice. Here, we notice that .γ does indeed not depend on Z.
In order to find the cluster number density for the Bethe lattice, we start by
developing a more general way to find the cluster number density. To find the cluster
number density for a given s, we need to find all possible configurations, .c(s), of
clusters of size s, and sum up their probability:
n(s, p) =
. ps (1 − p)t (c) . (3.22)
c(s)
This was simple in one dimension, because there was only one possible configu-
ration for a given s. Just like in one dimension we include the term .ps , because
we know that we must have all the s sites of the cluster present. In addition, we
need to include a term that takes into account that all the neighboring sites must be
empty. In general, we introduce the notation .t (c) for the number of neighbors for
configuration c. In one dimension, t is always 2, but in higher dimensions different
clusters may have different numbers of neighbors, as illustrated for the case of a
two-dimensional system in Fig. 3.3. We therefore need to include the term .(1 − p)t
to account for the probability for the t neighbors to be unoccupied. Based on this,
we realize that we may sum over all values of t. However, we then need to include
the effect that there are several clusters that can have the same t. We will then have
to introduce a degeneracy factor .gs,t which gives the number of different clusters
that have size s and a number of neighbors equal to t. The cluster number density
can then be written as
.n(s, p) = p gs,t (1 − p)t .
s
(3.23)
t
Degeneracy for the Bethe Lattice For the Bethe lattice, there is a particularly
simple relation between the number of sites, s, and the number of neighbors, t. We
can see this by looking at the first few generations of a Bethe lattice grown from a
central seed. For .s = 1, the number of neighbors is .t1 = Z. To add one more site,
we have to remove one neighbor from what we had previously, and then we add
.Z − 1 new neighbors, that is, for .s = 2 we have .t2 = t1 + (Z − 2). Consequently
we get an iterative relation for the number of neighbors .tk when we have k sites:
tk = tk−1 + (Z − 2) ,
. (3.24)
and therefore:
ts = s(Z − 2) + 2 .
. (3.25)
Cluster Number Density The cluster number density, given by the sum over all t,
is therefore reduced to only a single term for the Bethe lattice
For simplicity, we will write .gs = gs,ts . In general, we do not know .gs , but we will
show that we still can learn quite a lot about the behavior of .n(s, p). The cluster
density can therefore be written as
n(s, p) = gs ps (1 − p)2+(Z−2)s .
. (3.27)
1
Z=3
0.8 Z=4
Z=5
0.6
0.4
0.2
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
Fig. 3.4 A plot .f (p) = p(1 − p)Z−2 , which is a term in the cluster number density .n(s, p) =
gs [p(1 − p)Z−2 ]s (1 − p)2 for the Bethe lattice. We notice that .f (p) has a maximum at .p = pc ,
and that the second derivative, .f '' (p), is zero in this point. A Taylor expansion of .f (p) around
.p = pc will therefore have a second order term in .(p − pc ) as the lowest-order term—to lowest
order it is a parabola at .p = pc . It is this second order term which determines the exponent .σ ,
which consequently is independent of Z
Taylor Expansion Around .p = pc Let us now address .n(s, p) for p close to .pc
for a general value of Z. In this range, we will do a Taylor expansion of the term
.f (p) = p(1 − p)
Z−2 , which is raised to the power s in the equation for .n(s, p) in
(3.28). The shape of .f (p) as a function of p is shown in Fig. 3.4. The maximum
of .f (p) occurs for .p = pc = 1/(Z − 1). This is also easily seen from the first
derivative of .f (p).
which shows that .f ' (pc ) = 0. (We leave it to the reader to show that .f '' (pc ) < 0.)
The Taylor expansion of .f (p) around .p = pc is then:
1
f (p) = f (pc ) + f ' (pc )(p − pc ) + f '' (pc )(p − pc )2 + O((p − pc )3 ) ,
. (3.33)
2
where we already have found the first order term, .f ' (pc ) = 0. We can therefore
write
1
f (p) ≃ f (pc ) − f '' (pc )(p − pc )2 = A(1 − B(p − pc )2 ) .
. (3.34)
2
42 3 Infinite-Dimensional Percolation
Cluster Number Density We will now insert this Taylor expansion back into the
expression for the cluster number density:
When p is close to .pc , .(p − pc )2 is small number, and we can use the first order of
the Taylor expansion of .ln(1 − x) ≃ −x + O(x 2 ), to get
n(s, pc ) = gs As (1 − pc )2 .
. (3.38)
sξ = B −1 (p − pc )−2 .
. (3.41)
This implies that the characteristic cluster size .sξ diverges as a power-law with
exponent .1/σ = 2 as p approaches .pc . The general scaling form for the
characteristic cluster size .sξ is
sξ ∝ |p − pc |−1/σ ,
. (3.42)
where the exponent .σ is universal, which means that is does not depend on lattice
details such a Z, but it does depend on lattice dimensionality. It will therefore have
a different value for two-dimensional percolation.
Scaling Ansatz for .n(s, pc ) We can use our knowledge of the behavior of P and S
when p approaches .pc to constrain the scaling behavior of .n(s, pc ). We know that
3.4 Cluster Number Density 43
we can find S and P from the cluster number density. The average cluster size S as
p approaches .pc is
Γ
S=
. , (3.43)
pc − p
which should diverge when p approaches .pc . We rewrite S in terms of .n(s, p):
S=
. s 2 n(s, p) , (3.44)
which should diverge as p approaches .pc . We rewrite the sum as an integral in the
limit of .p = pc
∞
. S= s 2 n(s, pc )ds . (3.45)
0
This integral must diverge. We can therefore conclude that .n(s, pc ) is not an
exponential, otherwise the integral then would converge. We therefore assume that
the cluster number density has a power-law shape, that is, we introduce the scaling
ansatz:
n(s, pc ) ≃ Cs −τ .
. (3.46)
This expression is only valid in the limit when .s ⪢ 1. We introduce this scaling
ansatz into the expression for P :
. sn(s, p) = p − P , (3.47)
s
1. The integral
∞ ∞
. s 2 n(s, pc )ds = s 2 Cs −τ ds (3.49)
0 0
1. The integral
∞ ∞
. sn(s, pc )ds = sCs −τ ds (3.50)
0 0
Scaling Behavior of .S(p) When p Is Close to .pc We can sum up our arguments
so far in the relation
Let us use this expression to calculate S, for which we know the exact scaling
behavior, and then again use this to find the value for .τ
∞
S=C
. s 2−τ e−s/sξ → C s 2−τ e−s/sξ ds . (3.53)
s 1
We now make a rough estimate. This is useful, since it is in the spirit of this book,
and it also provides the correct behavior. We assume that
∞ sξ
.S=C s 2−τ e−s/sξ ds ∼ C s 2−τ ds ∼ sξ3−τ , (3.54)
1 1
where we have used the sign .∼ to denote that the expressions have the same scaling
behavior. We can do it slightly more elaborately:
∞
S≃C
. s 2−τ e−s/sξ ds . (3.55)
1
sξ ∼ (p − pc )−2 ,
. (3.58)
and that
S ∼ (p − pc )−1 .
. (3.59)
5
τ=
. . (3.60)
2
This relation indeed satisfies the exponent relations we found above, since .2 <
5/2 ≤ 3. A plot of the scaling form is shown in Fig. 3.5.
Preliminary Scaling Theory for Cluster Number Density We have now devel-
oped a preliminary scaling theory for the cluster number density. In the coming
chapters, we will demonstrate that similar scaling relations also are valid for
percolation in other dimension. We have found that in the vicinity of .pc , we do
not expect deviations until we reach large s, that is, until we reach a characteristic
cluster size .sξ that increases as .p → pc . The general scaling form for the cluster
number density is
s
n(s, p) = n(s, pc )F (
. ), (3.61)
sξ
0
(p-pc )=0.100
-2
(p-pc )=0.010
-4 (p-pc )=0.001
-6
-8
-10
-12
-14
-16
-18
-20
0 1 2 3 4 5 6 7 8
where
n(s, pc ) = Cs −τ ,
. (3.62)
and
sξ = s0 |p − pc |−1/σ .
. (3.63)
P (p) ∼ (p − pc )β ,
. (3.64)
ξ ∼ |p − pc |−ν ,
. (3.65)
and
S ∼ |p − pc |−γ ,
. (3.66)
with a possible non-trivial behavior for higher moments of the cluster density.
Exercises
Exercise 3.1 (P (p) for Z = 4) Find P (p) for Z = 4 and determine β for this
value of Z.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Finite-Dimensional Percolation
4
intersect except in a single point, because there are no loops in the Bethe lattice. This
is not the case in two and three dimensions, where loops always will be possible.
When there are loops present, we cannot use the arguments we used for the Bethe
lattice, because a branch cut off at one point may be connected again further out.
For the Bethe lattice, we could also estimate the multiplicity .g(s, t) of the clusters,
the number of possible clusters of size s and surface t, since t was a function of
s. In a two- or three-dimensional system this is not similarly simple, because the
multiplicity .g(s, t), that is the number of different cluster configurations with size s
and surface t, is not simple even in two dimensions, as illustrated in Fig. 4.1.
This means that the solution methods used for the one dimensional and the
infinite dimensional systems cannot be extended to address two-dimensional or
three-dimensional systems. However, several of the techniques and observations we
have made for the one-dimensional and the Bethe lattice systems, can be used as the
basis for a generalized theory that can be applied in any dimension. Here, we will
therefore pursue the more general features of the percolation system, starting with
the cluster number density, .n(s, p).
We have found that the cluster number density plays a fundamental role in our
understanding of the percolation problem, and we will use it here as our basis for
the scaling theory for percolation.
When we discussed the Bethe lattice, we found that we could write the cluster
number density as a sum over all possible configurations of cluster size, s:
n(s, p) =
. ps (1 − p)tj , (4.1)
j
where j runs over all different configurations, and .tj denotes the number of
neighbors for this particular configuration. We can simplify this by rewriting the
sum to be over all possible number of neighbors, t, and include the degeneracy .gs,t ,
the number of configurations with t neighbors:
.n(s, p) = gs,t ps (1 − p)t . (4.2)
t
The values of .gs,t can be found for smaller values of s. However, while this may
give us interesting information about the smaller cluster, and therefore for smaller
values of p, it does not help us to develop a theory for the behavior for p close to
.pc .
In order to address the cluster number density, we will need to study the
characteristics of .n(s, p), for example by generating numerical estimates for its
scaling behavior, and then propose a general scaling form which will be tested in
various settings.
4.1 Cluster Number Density 49
s=4, t=10 s=4, t=10 s=4, t=9 s=4, t=9 s=4, t=9 s=4, t=9
s=4, t=8 s=4, t=8 s=4, t=9 s=4, t=9 s=4, t=9 s=4, t=9
s=4, t=8 s=4, t=8 s=4, t=8 s=4, t=8 s=4, t=8 s=4, t=8
s=4, t=8
s=3, t=8 s=3, t=8 s=3, t=7 s=3, t=7 s=3, t=7 s=3, t=7
s=1, t=4
Fig. 4.1 Illustration of the possible configurations for two-dimensional clusters of size .s =
1, 2, 3, 4
50 4 Finite-Dimensional Percolation
Ns
n(s, p; L) =
. , (4.3)
M · Ld
where .Ns is the total number of clusters of size s measured for M simulations in a
system of size .Ld and for a given value of p. We perform these simulations just as
we did in one dimension, using the following program:
import numpy as np
import [Link] as plt
from [Link] import measurements
M = 2000
L = 200
p = 0.58
allarea = [Link]([])
for i in range(M):
z = [Link](L,L)
m = z<p
lw, num = [Link](m)
labelList = [Link]([Link]() + 1)
area = [Link](m, lw, labelList)
allarea = [Link](allarea,area)
n,sbins = [Link](allarea,bins=int(max(allarea)))
s = 0.5*(sbins[1:]+sbins[:-1])
nsp = n/(L*nsamp)
i = [Link](n)
[Link](figsize=(12,4))
[Link](1,2,1)
[Link](s[i],nsp[i],’o’)
[Link](’$s$’)
[Link](’$n(s,p)$’)
[Link](1,2,2)
[Link](s[i],nsp[i],’o’)
[Link](’$s$’)
[Link](’$n(s,p)$’)
The resulting plot of .n(s, p; L) for .L = 200 is shown in Fig. 4.2. Unfortunately,
this plot is not very useful. The problem is that there are too many values of s for
which we have little or no data at all. For small values of s we have many clusters
for each value of s and the statistics is good. But for large values of s, such as for
clusters of size .s = 104 and above, we have less than one data point for each value
of s. Our measured distribution .n(s, p; L) is therefore a poor representation of the
real .n(s, p; L) in this range.
4.1 Cluster Number Density 51
0.00010
10− 1
linear bins
logarithmic bins
0.00008 −3
10
0.00006
n(s, p )
n(s, p )
10− 5
0.00004
10− 7
0.00002
10− 9
0.00000
200 400 600 800 1000 1200 1400 100 101 102 103 104
s s
Fig. 4.2 Plot of .n(s, p; L) estimated from .M = 1000 samples for .p = 0.58 and .L = 200. (Left)
Direct plot. (Right) Log-log plot with linear and logaritmic binning
The problem with the measured results in Fig. 4.2 occur because we have chosen a
very small bin size for the histogram. For small values of s we want to have a small
bin size, since the statistics here is good, but for large values of s we want to have
larger bin sizes. This is often solved by using logarithmic binning: We make the bin
edges .a i , where a is the basis for the bins and i is bin number. If we chose .a = 2 as
the basis for the bins, the bin edges will be .20 , 21 , 22 , 23 , . . ., that is .1, 2, 4, 8, . . ..
We then count how many events occur in each such bin. If we number the bins
using the index i, then the edges of the bins are .si = a i , and the width of bin i is
.Δsi = si+1 − si . We then count how many events, .Ni , occur in the range from .si
to .si + Δsi . The average number of clusters, .N̄i in each bin in the interval .Δsi is
.N̄i = Ni /Δsi for a single realization and .N̄i = Ni /(Δsi M) for M realizations. The
estimate for the cluster number density in the middle point of the bin, that is for
.s̄i = (si+1 + si )/2, is
N̄i Ni
n(s̄i , p; L) =
. = . (4.4)
Ld MΔsi Ld
A common mistake is to forget to divide by the bin size .Δsi when the bin sizes are
not all the same! We implement method by generating an array of all the bin edges.
First, we find an upper limit to the bins, that is, we find an .im so that
We can for example round the right hand side up to the nearest integer
a = 1.2
logamax = [Link]([Link](max(allarea))/[Link](a));
52 4 Finite-Dimensional Percolation
where allarea corresponds to all the s-values. We can then generate an array of
indices from 1 to this maximum value
logbins = a**[Link](0,logamax)
And we can further generate the histogram with this set of bin edges
nl,nlbins = [Link](allarea,bins=logbins)
The complete code for this analysis is found in the following script
a = 1.2
logamax = [Link]([Link](max(s))/[Link](a))
logbins = a**[Link](0,logamax)
nl,nlbins = [Link](allarea,bins=logbins)
ds = [Link](logbins)
sl = 0.5*(logbins[1:]+logbins[:-1])
nsl = nl/(M*L**2*ds)
[Link](sl,nsl,’.b’)
The resulting plot for .a = 1.2 is shown in Fig. 4.2. Notice that the logarithmically
binned plot is much easier to interpret than the linearly binned plot. You should,
however, always reflect on whether your binning method may influence the resulting
plot in some way, since there may be cases where your choice of binning method
may affect the results you get. Although this is not expected to play any role in
your measurements in this book. We will in the following adapt logarithmic binning
strategies whenever we measure a dataset which is sparse.
law behavior over some range of s-values, but drops rapidly for larger s-values. This
is an effect of the characteristic cluster size, which also can be visually observed in
Figs. 1.4 and 1.5, where we see that the characteristic cluster size increases as p
approaches .pc . How can we characterize the characteristic cluster size based on this
measurement of .n(s, p)? We could measure .sξ directly from the plot, by drawing
4.1 Cluster Number Density 53
a 0
-2 p=0.45
p=0.50
-4 p=0.54
p=0.57
-6 p=0.58
-8
-10
-12
-14
0 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5
b 104
4
0
0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8
Fig. 4.3 (a) Plot of .n(s, p; L) as a function of s for various values of p for a .512 × 512 lattice.
(b) Plot of .sξ (p) measured from the plot of .n(s, p) corresponding to the points shown in circles in
(a)
a straight line parallel to but below .n(s, pc ), as illustrated in Fig. 4.3. When the
measured, .n(s, p) intersects this drawn line, .n(s, p) has fallen by a constant factor
below .n(s, pc ). We define this as .sξ , and we measure it by reading the values from
the s-axis. The resulting set of .sξ values are plotted as a function of p in Fig. 4.3.
We see that .sξ increases and possibly diverges as p approaches .pc . This is an effect
we also found in the one-dimensional and the infinite-dimensional case, where we
found that
sξ ∝ |p − pc |−1/σ
. (4.7)
where .σ was 1 is one dimension. We will now use this to develop a theory for
both .n(s, p; L) and .sξ based on our experience from one and infinite dimensional
percolation.
54 4 Finite-Dimensional Percolation
When we develop a theory, we realize that we are only interested in the limit .p →
pc , that is .|p − pc | ⪡ 1, and .s ⪢ 1. In this limit, we expect .sξ to mark the cross-
over between two different behaviors. There is a common behavior for .n(s, p) for
all p-values for small s, up to a cut-off, .sξ , as we also observe in Fig. 4.3: The curves
for different p-values are approximately equal for small s.
Based on what we observed in one-dimension and infinite-dimensions, we expect
and propose the following scaling form for .n(s, p):
s
n(s, p) = n(s, pc )F (
. ), (4.8)
sξ
n(s, pc ) = Cs −τ ,
. (4.9)
sξ = s0 |p − pc |−1/σ .
. (4.10)
Based on the methods presented in this book, we have estimated the exponents for
various systems and listed them in Table 4.1. You can find an up-to-date list of all the
exponents in the wikipedia article on percolation thresholds at [Link]
org/wiki/Percolation_critical_exponents.
We will often simplify the scaling form by writing it on the form:
n(s, p) = s −τ F̂ ((p − pc )s σ ) ,
. (4.12)
where we have introduced the function .F̂ (u) = F (uσ ). These forms are equivalent,
but in some cases this form produces simpler calculations.
This scaling form should in particular be valid for both the 1d and the Bethe
lattice cases—let us check this in detail.
Table 4.1 Values for scaling exponents for percolation in 1, 2, 3, 4 and infinite dimensions [8,37]
d .β .τ .σ .γ .ν D .μ .Dmin .Dmax .DB
1 2 1 1 1
2 .0.14 .2.05 .0.4 .2.4 .1.33 .1.89 .1.3 .1.1 .1.4 .1.6
3 .0.4 .2.2 .0.45 .1.8 .0.9 .2.5 .2.0 .1.3 .1.6 .1.7
-1
-2
-3 p=0.45
p=0.50
p=0.54
-4 p=0.57
p=0.58
-5
-5 -4 -3 -2 -1 0 1
Fig. 4.4 A plot of .n(s, p)s τ as a function of .|p − pc |1/σ s shows that the cluster number density
satisfies the scaling ansatz of (4.11)
In the case of one-dimensional percolation, we know that we can write the cluster
density exactly as
where .F (u) = u2 e−u . This is indeed on the general scaling form with .τ = 2.
For the Bethe lattice we found that the cluster density was approximately on the
form
n(s, p) = s −τ F (s/sξ ) .
. (4.16)
56 4 Finite-Dimensional Percolation
While the scaling ansatz has a simple form, it has powerful consequences. Here,
we address the theoretical consequences of the scaling ansatz, and demonstrate how
we can use the scaling ansatz in theoretical arguments. The methods we introduce
here are important methods in scaling theories, and we will use them in theoretical
arguments throughout this text.
Let us demonstrate how we can use the scaling ansatz to calculate the scaling of the
average cluster size, S, and how this can be used to provide limits for the exponent
.τ .
We assume that we study systems where p is close to .pc so that the cluster number
density .n(s, p) is wide and that its drop-off (crossover) .sξ is rapid. The sum over
s will then be a sum with many non-negligible terms and we can approximate this
sum by an integral over s instead:
∞
∞
S(p) =
. s n(s, p) ≃
2
s 2 n(s, p) ds . (4.18)
s=1 1
We know that the function .F (s/sξ ) goes very rapidly to zero when s is larger than .sξ ,
and that it is approximately a constant when s is smaller than .sξ . We will therefore
approximate .F (u) by a step function which is a constant up to .u = 1 and then 0
for .u > 1. Consequently, we only integrate up to .sξ , over a region where .F (s/sξ ) is
4.2 Consequences of the Scaling Ansatz 57
approximately a constant:
∞ sξ
S(p) =
. s 2−τ F (s/sξ ) ds ≃ Cs 2−τ ds . (4.20)
1 1
3−τ
γ =
. . (4.23)
σ
Consequences for .τ We have demonstrated that the average cluster size diverges
when .p → pc , which implies that the exponent .γ must be positive. In turn, this
implies that
3−τ
γ >0 ⇒
. >0 ⇒ 3>τ . (4.24)
σ
We have therefore found a first bound for .τ : .τ < 3. As an exercise, you can check
that this relation holds for the one-dimensional system and the Bethe lattice.
We may use a similar argument to find the behavior of .P (p) from the cluster number
density, which will give us further scaling relations between exponents and another
bound on the exponent .τ .
This equation expresses that a site picked at random is occupied with probability p
(right hand side), and that this site must either be in a finite cluster, with a probability
58 4 Finite-Dimensional Percolation
corresponding to the sum . s sn(s, p), or in the infinite cluster with probability
.P (p). We can therefore find .P (p) from
P (p) = p −
. sn(s, p) . (4.26)
s
Finding the Sum Using the Scaling Ansatz We can find the sum over .sn(s, p)
when p is close to .pc by transforming the sum to an integral and inserting the
scaling ansatz .n(s, p) = s −τ F (s/sξ ):
∞
∞
. sn(s, p) ≃ ss −τ F (s/sξ ) . (4.27)
s=1 s=1
Here, .F (s/sξ ) is approximately a constant when .s < sξ and goes very rapidly to
zero when .s > sξ , so we integrate up to .sξ assuming that .F (s/sξ ) is a constant C
up to .sξ , giving
∞ sξ
. s 1−τ
F (s/sξ )ds ≃ Cs 1−τ ds = c1 + c2 sξ2−τ . (4.29)
1 1
We insert this back into the expression for .P (p) in (4.26) getting:
P (p) = p −
. sn(s, p) ≃ p − c1 − c2 sξ2−τ . (4.30)
s
2−τ ≤0 ⇒ 2≤τ .
. (4.31)
Relating the Exponents .β and .τ We can rewrite the expression in (4.30) for .P (p)
and insert .sξ = s0 |p − pc |−1/σ , getting:
2−τ
P (p) ≃ p − c1 − c2 sξ2−τ ≃ (p − pc ) + c2 |p − pc |−1/σ
. (4.32)
We realize that when .p → pc the linear term .(p − pc ) will be smaller than the term
|p − pc |(τ −2)/σ . And we remember that .P (p) ∝ (p − pc )β . This gives us a new
.
τ −2
.β= . (4.33)
σ
We have therefore again demonstrated the power of the scaling ansatz by both
calculating bounds for .τ and by finding relations between the scaling exponents.
While the exponents are universal and independent of the details of the lattice but
dependent on the dimensionality, the percolation threshold, .pc , depends on details
of the system such as the lattice type and the type of percolation. We typically
discern between site percolation, where neighboring sites on a lattice are connected
if present, and bond percolation, where the presence of bonds between the sites
determines the connectivity. Table 4.2 provides basic values for the percolation
thresholds. These results have been measured with code in this book and therefore
have limited precision. You can find an updated set of percolation threshold for
various models on the Wikipedia page for percolation at [Link]
wiki/Percolation_threshold.
Exercises
We have then produced the array lw that contains labels for each of the connected
clusters.
(a) Familiarize yourself with labeling by looking at lw, and by studying the second
example in the python help system on the image analysis toolbox.
We can examine the array directly by mapping the labels onto a color-map,
using imshow.
[Link](lw)
You can also extract information about the clusters using the skimage.
measure module. This provides a powerful set of tools that can be used to
characterize the clusters in the system. For example, you can determine if a
system is percolating by looking at the extent of a cluster. If the extent in any
direction is equal to L, then the cluster is spanning the system. We can use this
to find the area of the spanning cluster or to mark if there is a spanning cluster:
import skimage
props = [Link](lw)
spanning = False
for prop in props:
if ([Link][2]-[Link][0]==L or
[Link][3]-[Link][1]==L):
# This cluster is percolating
area = [Link]
spanning = True
break
Exercises 61
(b) Using these features, write a program to calculate P (p, L) for various p for the
two-dimensional system.
(c) How robust is your algorithm to changes in boundary conditions? Could you do
a rectangular grid where Lx ⪢ Ly ? Could you do a more complicated set of
boundaries? Can you think of a simple method to ensure that you can calculate
P for any boundary geometry?
(a) Write a program to find P (p, L) and Π (p, L) for L = 2, 4, 8, 16, 32, 64, 128.
Comment on the number of samples you need to make to get a good estimate
for P and Π .
(b) Test the program for small L by comparing with the exact results from above.
Comment on the results?
P (p, L) ∼ (p − pc )β .
. (4.34)
Use the data from above to find an expression for β. For this you may need that
pc = 0.59275.
Your task is to determine the distribution function fZ (z) for this distribution.
(a) Find the cumulative distribution, that is, P (Z > z). You can then find the actual
distribution from
dP (Z > z)
fZ (z) =
. . (4.35)
dz
(b) Generate a method to do logarithmic binning in python. That is, you estimate
the density by doing a histogram with bin-sizes that increase exponentially in
size.
Exercise 4.5 (Cluster Number Density n(s, p)) We will generate the cluster
number density n(s, p) from the two-dimensional data-set.
(a) Find the average (finite) cluster size S(p) for p close to pc , for p above and
below pc .
(b) Determine the scaling exponent S(p) ∼ |p − pc |−γ .
(c) In what ways can you generate S (k) (p)? What do you think is the best way?
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Geometry of Clusters
5
We have seen how we can characterize clusters by their mass, s. As p approaches .pc ,
the typical cluster size s increases as well as the characteristic cluster diameter. In
this chapter we will discuss the geometry of clusters, and by geometry we will mean
how the number of sites in a cluster is related to the linear size of the cluster. We
will introduce several measures to characterize the spatial extent, the characteristic
radius .Rs , of clusters of size s. We will measure .Rs to motivate that it is proportional
to .s 1 /D, where D is a new exponent characterizing the dimension of clusters.
We will demonstrate that the percolation system is characterized by two lengths,
the system size L and a characteristic cluster size .ξ , and that the system shows
fractal, self-similar scaling when the characteristic length diverges. We develop
scaling theories for .P (s, L) for .p > pc and lay the foundations for a geometrical
understanding and description of the spanning cluster.
We have so far studied the clusters in our model porous material, the percolation
system, through the distribution of cluster sizes, .n(s, p), and properties that can be
found from .n(s, p), such as the average cluster size, S and the characteristic cluster
size, .sξ . However, clusters with the same mass, s, can have very different shapes.
Figure 5.1 illustrates three clusters all with .s = 20 sites. Notice that the linear and
the compact clusters are unlikely, but possible realizations. How can we characterize
the spatial extent of these clusters?
There are many ways to define the extent of a cluster. We could, for example,
define the maximum distance between any two points in a cluster i (.Rmax,i ) to be
the extent of the cluster, or we could use the average distance between two points in
the cluster. However, it is common to use the standard deviation of the position of
the sites in a cluster, which we recognize as the radius of gyration of a cluster. The
radius of gyration .Ri for a cluster i of size .si with sites at .rj for .j = 1, . . . , si , is
Fig. 5.1 Illustrations of three clusters all with .s = 24. The red circle illustrates the radius of
gyration for the clusters
defined as
1
si
q
Ri2 =
. rj − rcm,i , (5.1)
si
j =1
where .rcm,j is the center of mass of cluster i. The values for .Ri of the clusters in
Fig. 5.1 are illustrated by circles.
As we see from Fig. 5.1, clusters of the same size s can have different radii. How
can we then find a characteristic size for a given cluster size s? We find that by
averaging .Ri2 over all clusters of the same size s:
Rs2 = 〈Ri2 〉i .
. (5.2)
Let us see how this can be done analytically in one dimension and numerically in
two dimensions.
We can use the one-dimensional percolation system to gain insight into how we
expect .Rs to depend on s. In one dimension, a cluster of size s can only be realized
in one way, as a line of length s. If the cluster runs from 1 to s, the center of mass is
at .s/2, and the sum over all sites runs from 1 to s:
1
s
Rs2 =
. (i − s/2)2 , (5.3)
s
i=1
5.1 Geometry of Finite Clusters 65
where we assume that s is so large that we only need to address the leading term in
s, and that we do not have to treat even and odd s separately. This can be expanded
to
1 2
s
s2 1 s(s + 1)(2s + 1) s(s + 1) s2
Rs2 =
. [ i − is + ] = −s +s ∝ s2 ,
s 4 s 6 2 4
i=1
(5.4)
where we have used that . si=1 s 2 = s(s + 1)(2s + 1)/6 and . si=1 s = s(s + 1)/2
and where we only have kept the leading term in s. This shows that .Rs2 ∝ s 2 , which
means that .s ∝ Rs in one dimension. This is indeed what we expected. The extent
of the cluster is proportional to s because the cluster is s sites long.
For the one-dimensional system we found that .s ∝ Rs . How does this generalize to
higher dimensions? We start by measuring the behavior for a given value of p for a
finite system of size L. Our strategy is: (i) to generate clusters on a .L × L lattice;
(ii) for each cluster, i, of size .si , we will find the center of mass and the radius of
gyration, .Ri2 ; and (iii) for each value of s we will find the average radius, .Rs2 , by
a linear average. For larger values of s we will collect the data in bins, using the
logaritmic binning approach we developed to measure .n(s, p).
1
si
rcm,i =
. ri,j , (5.5)
si
j =1
We generate a lattice, ensure that each cluster are marked with the index of the
cluster, and find the center of mass cm using a built-in command:
L = 400
p = 0.58
z = [Link](L,L)
m = z<p
lw, num = [Link](m)
cm = measurements.center_of_mass(m, lw, labelList)
Second, we calculate the radius gyration by running through all the sites ix,iy in
the lattice. For each site, we find the index i of the cluster that it belongs to from i
= lw[ix,iy]. If the site belongs to a cluster, that is if i>0, we add the sum of
66 5 Geometry of Clusters
the square of the distance from the site to the center of mass to the radius of gyration
for cluster i:
dr = [Link]([ix,iy])-cm[i]
rad2[i] = rad2[2] + [Link](dr,dr)
After running through all the site, we divide by the mass (size) .si of each cluster to
find the radius of gyration according to the formula
1
si
2
2
.Ri = ri,j − rcm,i , (5.6)
si
j =1
def radiusofgyration(m,lw,L):
labelList = [Link]([Link]() + 1)
area = [Link](m, lw, labelList)
cm = measurements.center_of_mass(m, lw, labelList)
rad2 = [Link](int([Link]()+1))
for ix in range(L):
for iy in range(L):
ilw = lw[ix,iy];
if (ilw>0):
dr = [Link]([ix,iy])-cm[ilw]
dr2 = [Link](dr,dr)
rad2[ilw] = rad2[ilw] + dr2
rad = [Link](rad2/area)
return area,cm,rad2
M = 20 # Nr of samples
L = 400 # System size
p = 0.58 # p-value
allr2 = [Link]([])
allarea = [Link]([])
for i in range(M):
z = [Link](L,L)
m = z<p
lw, num = [Link](m)
area,rcm,rad2 = radiusofgyration(m,lw,L)
allr2 = [Link](allr2,rad2)
allarea = [Link](allarea,area)
[Link](allarea,allr2,’k.’)
We use this function to calculate the average radius of gyration for each cluster size
s for M different lattice realizations, and plot the results using the following script:
5.1 Geometry of Finite Clusters 67
M = 20 # Nr of samples
L = 400 # System size
p = 0.58 # p-value
allr2 = [Link]([])
allarea = [Link]([])
for i in range(M):
z = [Link](L,L)
m = z<p
lw, num = [Link](m)
area,rcm,rad2 = radiusofgyration(m,lw,L)
allr2 = [Link](allr2,rad2)
allarea = [Link](allarea,area)
[Link](allarea,allr2,’k.’)
[Link]("$s$")
[Link]("$R_s^2$")
Scaling Behavior of the Radius of Gyration The resulting plots for several values
of p are shown in Fig. 5.2. We see that there is an approximately linear relation
between .Rs2 and s in this double-logarithmic plot, which indicates that there is a
power-law relationship between the two:
.Rs2 ∝ s x . (5.7)
How can we interpret this relation? Equation (5.7) relates the radius .Rs and the area
(or mass) of the cluster. We are more used to the inverse relation:
s ∝ RsD ,
. (5.8)
105
p=0.40
p=0.45
104
p=0.50
p=0.55
103 p=0.59
102
R s2
101
100
10-1
100 101 102 103 104 105
s
Fig. 5.2 Plot of .Rs2 as function of s for simulations on two-dimensional systems with .L = 400.
The largest cluster for each value of p is illustrated by a circle. The dotted line shows the curve
.Rs ∝ s
2 2/D for .D = 1.89
68 5 Geometry of Clusters
where .D = 2/x is the exponent relating the radius to the mass of a cluster. This
corresponds to our intuition from geometry. We know that for a cube of size L, the
mass (or volume) of the cube is .M = L3 . For a square of length L, the mass (or area)
is .M = L2 , and similarly for a circle .M = π R 2 , where R is the radius of the circle.
For a line of length L, the mass is .M = L1 . We see a general trend, .M ∝ R d , where
R is a characteristic length for the object, and d describes the dimensionality of the
object. If we extend this intuition to the relation in (5.8), which is an observation
based on Fig. 5.2, we see that we may interpret D as the dimension of the cluster.
However, the value of D is not an integer. We have indicated the value of .D = 1.89
with a dotted line in Fig. 5.2. This non-integer value of D may seem strange, but it
is fully possible, mathematically, to have non-integer dimensions. This is a feature
frequently found in fractal structures, and percolation clusters as p approaches .pc
are indeed good examples of self-similar fractals. We will return to this aspect of
the geometry of the percolation system in Sect. 5.3.
Characteristic Cluster Radius The largest cluster and its corresponding radius
of gyration is indicated by a circle for each p value in Fig. 5.2. We see that as p
approaches .pc , both the mass and the radius of the largest cluster increases. Indeed,
this corresponds to the observation we have previously made for the characteristic
cluster size, .sξ . We may define a corresponding characteristic cluster radius, .Rsξ
through
1/D
sξ ∝ RsDξ
. ⇐⇒ Rsξ ∝ sξ . (5.9)
This length is a characteristic length for the system for a given value of p,
corresponding to the largest cluster size or the typical cluster size in the system. In
Sect. 5.2 we see how we can relate this length directly to the cluster size distribution.
We have already found that the characteristic cluster size .sξ diverges as a power law
as p approaches .pc :
sξ ≃ s0 (p − pc )−1/σ ,
. (5.10)
when .p < pc . The behavior is similar when .p > pc , but the prefactor .s0 may
be different. How does .Rsξ behave when p approaches .pc ? We can find this by
combining the scaling relations for .sξ and .Rsξ from (5.9):
1/D
∝ (p − pc )−1/σ ∝ (p − pc )−1/σ D ,
1/D
Rsξ ∝ sξ
. (5.11)
where we introduce the symbol .ν = 1/(σ D). For two-dimensional percolation, the
exponent .ν is a universal number, just like .σ and D. By universal, we mean that it
5.2 Characteristic Cluster Size 69
does not depend on details such as the lattice type or the connectivity of the lattice,
although it does depend on the dimensionality of the system. We know the exact
value of .ν in two dimensions, .ν = 4/3.
The argument we have provided here is an example of a scaling argument. In
these arguments we are only interested in the exponent in the scaling relation, which
gives us the functional form, and not in the values of the prefactors.
We have now defined the characteristic size of a cluster of size s through .Rs . In
addition, we introduced a characteristic cluster length .Rsξ , which characterizes the
whole system and not only clusters of a particular size s. However, there are several
ways we can introduce a length scale to describe the typical cluster size in a system.
Here, we will introduce two such measures, the average radius of gyration R and
the characteristic length .ξ .
We have now defined the characteristic length .Rsξ through the definition of the
characteristic cluster size, .sξ , and the scaling relation .s ∝ RsD . However, it may be
more natural to define the characteristic length of the system as the average radius
and not the cut-off radius. We introduced the radius of gyration for clusters of size
s by averaging the radius of gyration .Ri over all clusters i of size s:
Rs2 = 〈Ri2 〉i ,
. (5.12)
This gives us the radius of gyration .Rs , which we found to scale with cluster mass s
as .s ∝ RsD .
1
S=
. s sn(s, p) , ZS = sn(s, p) . (5.13)
ZS s s
120
L= 64
L=128
100 L=256
L=512
80
R (p)
60
40
20
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
p
Fig. 5.3 A plot of the average radius of gyration R as a function of p and L. We observe that R
increases as .p → pc , but is limited in magnitude by the finite system size L
1 2 2
R=
. Rs s n(s, p) , ZR = s 2 n(s, p) , (5.15)
ZR s s
where we notice that the normalization sum .ZR is the average cluster size, S, .ZR =
S. Figure 5.3 shows a plot of the average R as a function of p for various systems
sizes L. We see that R diverges as p approaches .pc . How can we develop a theory
for this behavior?
A Scaling Form for R We know that the cluster number density .n(s, p) has the
approximate scaling form
n(s, p) = s −τ F s/sξ , sξ ∝ |p − pc |−1/σ .
. (5.16)
We can use this to calculate the average radius of gyration, R, when p is close to
pc . We find the scaling behavior of the average radius of gyration by replacing the
.
5.2 Characteristic Cluster Size 71
where we have inserted .Rs2 ∝ s 2/D . While this expression only is valid when .s <
sξ , we can apply it here since .F (s/sξ ) goes rapidly to zero when .s > sξ , and
therefore only the .s < sξ values will contribute significantly to the integral. We
change variables to .u = s/sξ , getting:
2/D+3−τ ∞
sξ 1/sξ u2/D+2−τ F (u) du
R2 ∝
.
∞ . (5.19)
sξ3−τ 1/sξ u2−τ F (u) du
∞ 2/D+2−τ
2/D 0 u F (u) du 2/D
∝ sξ ∞ 2−τ ∝ sξ , (5.20)
0 u F (u) du
where the lower limit .1/sξ goes to zero for large .sξ and the two integrals over
F (u) simply are numbers and therefore have been included in the constant of
.
proportionality.
2/D
The Characteristic Lengths Are Proportional This shows that .R 2 ∝ sξ . We
2/D
found above that .Rsξ ∝ Therefore, .R ∝ Rsξ ! These two characteristic
sξ .
lengths therefore have the same behavior. They are only different by a constant
of proportionality, .R = c Rsξ . We can therefore use either length to characterize the
system—they are effectively the same up to a prefactor. This is not only true for
these two lengths, but all lengths have the same asymptotic scaling behavior close
to .pc . For example, Fig. 5.4 illustrates the radius of gyration of the largest cluster
with a circle and the average radius of gyration, R, by the length of the side of the
square. As p increases, both the maximum cluster size and the average cluster size
increases in concert.
Correlation Length
We can also define the typical size of a cluster from the correlation function. We
recall that the correlation function .g(r, p) is the probability that an occupied site i
at .ri is connected to a site j at .rj , where .r = rj − ri and .r = |r|. The correlation
function is only a function of the relative position of the two sites, .r, which we
usually only write as r, because we assume that the correlation function is isotropic.
We define the correlation length, .ξ , as the average squared distance between two
72 5 Geometry of Clusters
Fig. 5.4 Illustration of the largest cluster in .512 × 512 systems for .p = 0.55, .p = 0.57, and .p =
0.59. The circles illustrate the radius of gyration of the largest cluster, and the boxes show the size
of the average radius of gyration, .R = 〈Rs 〉. We observe that both lengths increase approximately
proportionally as p approaches .pc
connected sites:
2
r r g(r; p)
.ξ =
2
. (5.21)
r g(r; p)
where the sum is over all relative positions .r, that is, over all space. The denominator
is a normalization sum, which corresponds to the average cluster size, S. This length
is called the correlation length. However, to gain insight into this length, we will first
address the correlation function, its scaling behavior and its relation to the average
cluster size S.
g(r) = pr = e−r/ξ ,
. (5.22)
where .f (u) is a function that decays rapidly when u is larger than 1. We will assume
that this behavior can be generalized to higher dimension. That is, we expect the
correlation function to decay rapidly beyond a length, .ξ , that corresponds to the
typical extent of clusters in the system.
dimension. However, we can still measure it from our simulations, although such
measurements typically are computationally intensive. How can we measure it?
We can loop through all sites i and j and find their distance .rij . We estimate the
probability for two sites at a distance .rij to be connected by counting how many
of the sites that are a distance .rij apart that are connected and compare it to how
many sites in total that are a distance .rij apart. This is done through the following
program, which returns the correlation function g(r) estimated for a lattice lw
which contains the cluster indexes for each site, similar to what is returned by
the lw, num = [Link](m) command. We write a subroutine
perccorrfunc to find the correlation function for a given lattice lw, and then we
use this function to find the correlation function for several values of p and L:
import numpy as np
import [Link] as plt
from [Link] import measurements
from numba import jit
@jit
def perccorrfunc(m,lw,L):
r = [Link](2*L) # Positions
pr = [Link](2*L) # Correlation function
npr = [Link](2*L) # Nr of elements
for ix1 in range(L):
for iy1 in range(L):
lw1 = lw[ix1,iy1]
if (lw1>0):
for ix2 in range(L):
for iy2 in range(L):
lw2 = lw[ix2,iy2]
if (lw2>0):
dx = (ix2-ix1)
dy = (iy2-iy1)
rr = [Link](dx,dy)
# Find corresponding box
nr = int([Link](rr)+1)
pr[nr] = pr[nr] + (lw1==lw2)
npr[nr] = npr[nr] + 1
pr = pr/npr
return r,pr
m = z<p
lw, num = [Link](m)
r,g = perccorrfunc(m,lw,L)
pr[:,ip] = pr[:,ip] + g
rr[:,ip] = rr[:,ip] + r
pr = pr/M
r = r/M
Figure 5.5 shows the resulting plots of the correlation function .g(r; p) for various
values of p for an .L = 400 system. First, we notice that the scaling is rather poor.
We will understand this as we develop a theory for .g(r; p) below. The plot shows
that there also in two dimensions appear to be a cross-over length .ξ , which we call
the correlation length, beyond which the correlation function falls rapidly to zero.
For .r < ξ the correlation function appears to approximately be a power-law. Based
on our experience with percolation systems, we suggests the following functional
form
g(r; p) = r −x f (r/ξ ) ,
. (5.24)
where the cross-over function .f (u) falls rapidly to zero when .u > 1 and
is approximately constant when .u < 1. We expect that as p approaches .pc ,
100
10−1
10−2
g(r; p)
10−3
10−4
p = 0.57
10−5 p = 0.58
p = 0.59
r−η
10−6
101 102 103
r
Fig. 5.5 A plot of .g(r; p) as a function of r for various values of p. When p approaches .pc the
function approaches a power-law behavior .g(r) ∝ r −η , which is illustrated by a dashed line with
.η = 0.208. Notice the rapid cross-over for large r, which occurs at a characteristic length .ξ
5.2 Characteristic Cluster Size 75
the correlation length .ξ grows to infinity, and the correlation function .g(r; pc )
approaches a power-law .r −x for all values of r up to a length limited by the system
size L.
How can we relate .g(r; p) to .S(p)? We notice that .S(p) is the average number
of sites in a cluster, that is, the average number of sites connected to a given
occupied site. We can therefore find .S(p) by summing the probability for a site
to be connected for all possible relative positions .r:
S(p) =
. g(r; p) . (5.25)
r
where the integral is over a volume in space corresponding to all relative positions
r. We change integration variables to the radial distance r and the solid angle .Ω
.
S(p) ≃
. g(r)dr d = g(r)r d−1 dr dΩ , (5.27)
where the integral is from .r = 0 to .r → ∞ and over all solid angles .Ω.
Average Cluster Size at .pc We know that when .p → pc , then the average cluster
size, S, diverges. At .p = pc , the scaling form for the correlation function in (5.24)
is .g(r; p) ∝ r −x . The condition that the integral in (5.27) must diverge therefore
provides bounds for the exponent x:
∞
S(pc ) ≃
. r −x r d−1 dr dΩ = c r −x+d−1 dr . (5.28)
1
This integral diverges when .−x + d > 0, that is, for .x < d. It is common to
introduce an exponent .η so that .x = (d − 2 + η). The condition for .S(pc ) to diverge
is then .2 − η > 0. The correlation function at .pc is then:
g(r; pc ) ∝ r −(d−2)−η .
. (5.29)
76 5 Geometry of Clusters
This is consistent with what we found for the one-dimensional system where .x = 0
and .η = 1.
Finding .S(p) for p below .pc We use the scaling form of .g(r; p) in (5.24) to
calculate the integral in (5.27) for .p < pc :
r ∞ r
S(p) ≃
. r −(d−2)−η f ( )r d−1 dr dΩ = c r 1−η f ( ) dr . (5.30)
ξ 1 ξ
∞ 1
S(p) ≃ c
. ξ 2−η u1−η f (u) du = cξ 2−η u1−η du = cξ 2−η (1 − ξ −(2−η) ) .
1/ξ 1/ξ
(5.31)
Because .2 − η > 0, we see that .ξ −(2−η) approaches zero as p approaches .pc and .ξ
grows. The right-hand term is therefore approximately 1, and we get:
S(p) ∝ ξ 2−η .
. (5.32)
where we have introduced a new critical exponent, .ν, and related it to other
exponents through .ν = γ /(2 − η). For percolation in two dimensions, .ν = 4/3,
whereas in three dimensions it is .ν = 0.9.
Finding the Correlation Length from .n(s, p) The correlation length .ξ can be
found from the correlation function by
2
r r g(r)
.ξ =
2
, (5.34)
r g(r)
where . r g(r) = S(p). You can check this by inserting .g(r) = r −(d−2)−η f (r/ξ )
and calculating the integral.
We recall that the average radius of gyration is
2 2
s Rs s n(s, p)
.R =
2 2
, (5.35)
s s n(s, p)
5.2 Characteristic Cluster Size 77
averaged over all clusters of size i. Here, . s s 2 n(s, p) = S(p). Notice the
similarity between these two definitions. This indicates that the correlation length
can be defined in terms on .n(s, p). The common definition is
2 2
s 2Rs s n(s, p)
.ξ =
2 2
= 2R 2 . (5.36)
s s n(s, p)
The important aspect is that the two lengths .ξ and R are proportional to each other.
And we already round that .R ∝ Rsξ . This means that they all have the same scaling
behavior as p approaches .pc . This means that
Correlation Length All the lengths R, .Rsξ and .ξ has the same scaling
behavior.
The correlation length .ξ scales as
ξ ∝ |p − pc |−ν
. when p → pc , (5.38)
L. This is not surprising, since we cannot observe clusters that are larger than the
system size L. Figure 5.6 illustrates .ξ ∝ |p − pc |ν as we would expect it in an
infinite system. However, for a finite system, the curve for .ξ(p) is cut off at a length
proportional to L. This means that for p in some region around .pc , as illustrated in
Fig. 5.6, we have that .ξ > L. In this range we cannot determine if the system is at
.pc or not, because the system is not large enough for us to make this distinction.
ξ− (pc − p)−ν = ξ ≃ L ,
. (5.39)
78 5 Geometry of Clusters
pc p
Fig. 5.6 Illustration of the behavior of .ξ when p approaches .pc . In a finite system of size L, the
system will be percolating for all p in the range indicated by the arrows
Fig. 5.7 Illustration of the largest cluster in .512 × 512 systems with .p > pc , for .p = 0.593,
.p= 0.596, and .p = 0.610. The circles illustrate the radius of gyration of the largest cluster. We
observe that the radius of gyration increases as p approaches .pc
when .p < pc , and a similar expression for .p > pc . This means that when we
observe spanning we can only be sure that p is within a certain range of .pc :
|p − pc | = cL−1/ν .
. (5.41)
The correlation length .ξ is therefore the natural length to characterize the system.
At distances smaller than .ξ , the system behaves as if it is at .p = pc . However, at
distances much larger than .ξ , the system is essentially homogeneous. As we can
observe in Fig. 5.7 the system becomes more and more homogeneous when p goes
away from .pc . We will now address this feature in more detail when .p > pc .
5.3 Geometry of the Spanning Cluster 79
How can we develop a scaling theory for the spanning cluster, that is, a theory
for how the mass of the spanning cluster depends on system size L and the
characteristic cluster size, .ξ ? We know that as p is increased from below towards
.pc , the characteristic cluster size .ξ diverges. The mass of a characteristic cluster of
size .ξ is expected to follow the scaling relation .sξ ∝ ξ D . For a given value of p, we
can therefore choose the system size L to be equal to .ξ , .L = ξ(p). In this case, a
cluster of size .ξ would correspond to a cluster of size L, and it would be a spanning
cluster in this system. For this system of size .L = ξ , we therefore expect the mass
of the spanning cluster to be .M(p, L) ∝ ξ D ∝ LD . This suggests that the mass
of the spanning cluster in a system close to or at .pc depends on the system size L
according to .M(pc , L) ∝ LD .
The Density of the Spanning Cluster The density, .P (p, L), of the spanning
cluster at .p = pc therefore has the following behavior:
M(pc , L)
P (pc , L) =
. ∝ LD /Ld ∝ LD−d . (5.42)
Ld
Because we know that .P (pc , L) does not diverge when .L → ∞, we deduce that
D < d. The value of D in two-dimensional percolation is .D = 91/48 ≃ 1.90.
.
This implies that the density of the spanning cluster depends on the system
size, L. Indeed, since .D < d, we see that the density decreases with system size.
This may initially seem surprising, since we may be used to thinking of density
as a material property. However, we recognize this behavior from bodies that are
embedded in dimensions higher than themselves, such as for a thin sheet or a thin
rod embedded in three dimensions.
Embedded, Regular Bodies For example, consider a thin, flat sheet of thickness
h, and dimensions .L × L , placed in a three-dimensional space. If we cut out a
volume of size .L × L × L, so that .L ⪡ L , the mass of the sheet inside that volume
is
.M = hL2 , (5.43)
hL2
ρ=
. = hL−1 . (5.44)
L3
It is only in the case when we use a two-dimensional volume .L × L with a third
dimension of constant thickness H larger than h, that we recover a constant density
.ρ independent of system size.
80 5 Geometry of Clusters
Fig. 5.8 Illustration of three generations of the Sierpinski gasket starting from an equilateral
triangle
Dimension For an object with mass M and linear size L, we define the
dimension of the object as D, if .M(L) = cLD , that is, if the mass is
proportional to L to the power of D.
ln 3
D ln 2 = ln 3
. ⇒ D= ≃ 1.58 . (5.46)
ln 2
Thus, the dimension of this object is also not an integer, but lies between 1 and 2.
This indicates that the Sierpinski gasket is an object that is somewhere between a
plane and a line. It fills space less efficiently that a plane, but more efficiently than
a line.
We have presented a robust method for finding the dimension of iteratively
defined objects. However, the method also gives us the correct dimension for e.g. a
cube: For a cube of size L, if we double the size of the cube, that is .L' = 2L, the
mass is increased by a factor of 8, .M ' = M(2L) = 8M(L), which gives a dimension
of .d = ln 8/ ln 2 = 3.
L
N(δ) = ( )d ,
. (5.47)
δ
and for a fractal structure we expect
L
N(δ) = ( )D ,
. (5.48)
δ
We leave it as an exercise for the reader to address what happens when .δ → 1, and
when .δ → L.
Let us now return to the discussion of the mass .M(p, L) of the spanning cluster for
p > pc in a finite system of size L. The behavior of the percolation system for .p >
.
pc is illustrated in Fig. 5.7. We notice that the correlation length .ξ diverges when
p approaches .pc . At lengths larger than .ξ , the system is effectively homogeneous
because there are no holes significantly larger than .ξ . For such systems, there are
two effectively types of behavior, depending on whether L is larger than or smaller
than the correlation length .ξ .
Systems Effectively at .pc When .L ⪡ ξ and .p > pc , we are in the situation where
we cannot discern a system at p from a system at .pc because the size of the holes
(empty regions described by .ξ when .p > pc ) in the spanning cluster is much larger
than the system size. We say that the system is effectively at .pc . When we look at
82 5 Geometry of Clusters
Fig. 5.9 Illustration of the spanning cluster in a .512 × 512 system at .p = 0.595 > pc . In this
case, the correlation length is .ξ = 102. The system is divided into regions of size .ξ . Each such
region has a mass .M(p, ξ ) ∝ ξ D , and there are .(L/ξ )d ≃ 25 such regions in the system
this system through a window of size L, that is, in a finite system size, we do not
know if the spanning cluster in that system is an infinite cluster or just a cluster
that is larger than the system size L. However, because the mass of finite clusters
all scale as .s ∝ RsD up to a length .ξ that is much larger than L, the mass of the
spanning cluster of size L will also follow this scaling relation:
M(p, L) ∝ LD when L ⪡ ξ .
. (5.49)
Systems Away from .pc In the other case, when .L ⪢ ξ and .p > pc , the typical
size of a hole in the percolation cluster is .ξ , as illustrated in Fig. 5.9. This means that
on lengths much larger than .ξ , the percolation cluster is effectively homogeneous.
That is, it has holes up to a size .ξ , but not much larger than that. It looks like a Swiss
cheese with a finite size of the holes. We can therefore divide the .L × L system into
.(L/ξ ) regions of size .ξ . In each such region of size .𝓁 = ξ , the system is effectively
d
at the percolation threshold .pc . The mass of the spanning cluster in this region is
therefore given by (5.49), so that, .m ∝ 𝓁D , where .𝓁 = ξ , and therefore .m ∝ ξ D .
Consequently, the total mass of the spanning cluster is the mass of one such region
multiplied with the number of regions:
Scaling Behavior of the Mass Above .pc We have therefore derived the behavior
of the mass, .M(p, L), of the spanning cluster for .p > pc for system sizes L that are
Exercises 83
LD when L ⪡ ξ
M(p, L) ∝
. . (5.51)
ξ D−d Ld when L ⪢ ξ
This behavior can be rewritten in what we call the standard scaling form, with a
scaling behavior and a cut-off function:
L
M(p, L) = LD Y ( ) ,
. (5.52)
ξ
constant u ⪡ 1
Y (u) =
. . (5.53)
ud−D u ⪢ 1
We will use this function form many times when we discuss the behavior of finite
system sizes in the next chapter.
Exercises
(a) Find the mass M(L) of the percolating cluster at p = pc as a function of L, for
L = 2k , k = 4, . . . , 11.
(b) Plot log(M) as a function of log(L).
(c) Determine the exponent D.
1 1 1
Rs2 =
. 〈 (ri − rcm )2 = (ri − rj )2 ,
s 2 s2
i ij
where the sum over ij are over all pairs of sites in a cluster of size s.
Hint Show that both expressions are equal to i ri · ri − rc m · rcm .
(a) Write a program to find the correlation function, g(r, p, L) for L = 256.
(b) Plot g(r, p, L) for p = 0.55 to p = 0.65 for L = 256.
(c) Find the correlation length ξ(p, L) for L = 256 for the p-values used above.
(d) Plot ξ as a function of p − pc , and determine ν.
84 5 Geometry of Clusters
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Finite Size Scaling
6
In this chapter we will introduce the theory of finite size scaling and demonstrate
how we can apply the theory to improve our measurements of the properties
of percolation clusters. Usually, we attempt to measure properties of percolation
system in the largest possible system we can simulate. Here, we demonstrate that if
the system behaves according to simple scaling relations, it is instead much better
to systematically vary the system size and the interpolate to infinite system sizes.
This approach is generally called finite size scaling and we provide a thorough
introduction to the theory and its applications to understand the scaling of the
density of the spanning cluster, .P (p, L), the average cluster size, .S(p, L), and the
percolation probability .Π (p, L).
How can we utilize a disadvantage, such as a finite system size, to an advantage?
Usually, we have found a finite system size to be a hassle in simulations. We would
like to find the general behavior, but we are limited by the largest finite system
size we can afford to simulate. It may be tempting to put all our resources into one
attempt—to make one simulation in a really large system. However, this is usually
not a good strategy. Because we will then know that our results are limited by the
system size, but we do not know to what degree the finite system size affects our
result.
Instead, we will follow a different strategy: the strategy of finite size scaling.
We will systematically increase the system size, measure the quantities we are
interested in, and then try to extrapolate to an infinite system size. This has several
advantages: It allows us to understand and estimate the errors in our predictions,
and it allows us to use simulations of smaller systems. Indeed, it turns out that it is
more important to do simulations in smaller systems, than only to try to simulate
that largest system possible. However, for this to be effective, we need to have a
theoretical understanding of finite size scaling [7].
The methods we develop here are powerful and can be generalized to many
other experimental and computational situations. In many experiments it is also
tempting to try to perform the perfect experiment by reducing noise or measurement
errors. For example, we may perform an experiment where we need to make the
experimental system as horizontal as possible, because deviations from a horizontal
system would introduce errors. Instead of trying to make the system as horizontal as
possible, we may instead systematically vary the orientation, and then extrapolate
to the case when the system is perfectly horizontal. This allows us to control the
uncertainty. Of course, we cannot vary all possible uncertainties in an experiment or
a simulation, but this alternative mindset provides us with a new tool in our toolbox,
and a new way to deal with uncertainties.
In practical situations, we will always be limited by finite system sizes. If you
measure the size of earthquakes in the Earth’s crust, your results are limited by the
thickness of the crust or by the extent of a homogeneous region. If you simulate
a molecular system, you are definitely limited by the number of atoms you can
include in your simulation. Thus, better insight into how we can systematically vary
the system size and use this to gain insight are general tools of great utility.
Here, you will learn how to systematically vary system size L in order to
find good estimates for exponents and percolation thresholds. Indeed, my hope is
that you will see that finite size scaling is a powerful tool that can be used both
theoretically and computationally. To introduce this tool, we need to address specific
examples that can help build our intuition and shape our mindset. We will therefor
start from a few examples, such as the finite size scaling for the density of the
spanning cluster, .P (p, L), and then apply the method to a new case, the percolation
probability .Π (p, L).
We have found that a percolation system is described by three length-scales: the size
of a site, the system size L, and the correlation length .ξ . Finite size scaling addresses
the change in behavior of a system as we change the system size L. Typically, we
divide the behavior into two categories:
• When the system size L is much smaller than the correlation length .ξ , .L ⪡ ξ ,
the system appears to be on the percolation threshold.
• When L is much larger than .ξ , .L ⪢ ξ , the geometry is essentially homogeneous
at lengths longer than .ξ .
We will then address the behavior close to .pc . In the case of percolation, we
usually assume that the behavior is a power-law in .p − pc . For example, the mass
.M(p; L) of the spanning cluster:
M(p) ∝ (p − pc )−x ,
. (6.1)
term and a cut-off function, as you have seen several times in this book:
x L
M(p, L) = L ν f
. , (6.2)
ξ
Measuring .P (p, L) for finite L Let us now apply this methodology to study
the behavior of the density of the spanning cluster, .P (p, L), for finite system sizes.
First, we generate a plot of .P (p, L) for various values of L using the following
program:
import numpy as np
import [Link] as plt
from [Link] import measurements
LL = [25,50,100,200]
p = [Link](0.4,0.75,50)
nL = len(LL)
nx = len(p)
Ni = [Link](nx)
P = [Link]((nx,nL),float)
for iL in range(nL):
L = LL[iL]
N = int(2000*25/L)
for i in range(N):
z = [Link](L,L)
for ip in range(nx):
m = z<p[ip]
lw, num = [Link](m)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x>0)]
if (len(perc)>0):
Ni[ip] = Ni[ip] + 1
area = [Link](m, lw, perc[0])
P[ip,iL] = P[ip,iL] + area
P[:,iL] = P[:,iL]/((L*L)*N)
for iL in range(nL):
L = LL[iL]
[Link](p,P[:,iL],label="L = "+str(L))
88 6 Finite Size Scaling
0.21
L = 25
L = 50 0.20
0.6
L = 100
L = 200 0.19
P (pc ; L )
P (p, L )
0.4
0.18
0.2 0.17
0.16
0.0
0.4 0.5 0.6 0.7 50 100 150 200
p L
(a) (b)
Fig. 6.1 (a) Plot of .P (p, L). (b) Plot of .P (pc ; L) as a function of L
[Link](’$P(p,L)$’)
[Link](’$p$’)
[Link]()
The resulting plot of .P (p, L) is shown in Fig. 6.1. We see that as L increases,
P (p, L) approaches the shape expected in the limit when .L → ∞. We can see
.
P (p) ∝ (p − pc )β ∝ ξ −β/ν .
. (6.3)
This is valid in the limit when L → ∞, that is, when L ⪢ ξ . In the limit when
L ⪡ ξ , which eventually will occur as p approaches pc and ξ diverges, we see
from Fig. 6.1 that P (pc , L) depends on L. In this case, we have previously found
that
M(pc , L) LD
P (p, L) ≃ P (pc , L) =
. ∝ ∝ LD−d ∝ L−β/ν . (6.4)
L2 Ld
Combined, we therefore have the behavior
ξ −β/ν when L ⪢ ξ
P (p, L) ∝
. . (6.5)
L−β/ν when L ⪡ ξ
6.2 Finite Size Scaling of P (p, L) 89
Finite Size Scaling Ansatz The fundamental idea of finite size scaling is then to
assume a particular form of a function that encompasses this behavior both when
ξ ⪡ L and ξ ⪢ L, by rewriting the expression for P (p, L) as
Where we have assumed that the only relevant length scales are L and ξ , and that
the function therefore only can depend on a ratio between these two length scales.
How must the function f (u) behave for this general form to reduce to Eqs. (6.3) and
(6.4)?
First, we see that when ξ ⪢ L the function f (L/ξ ) should be a constant, that
is, f (u) is a constant when u ⪡ 1. Second, we see that when ξ ⪡ L, we need
the function f (L/ξ ) to cancel all the L-dependency in order to find the relation in
Eq. (6.3):
This will occur if and only if f (u) is a power-law, that is, f (u) = ua . In order to
cancel the L-dependency, the power-law exponent for the L-term must be zero:
Indeed, we could have used this in order to find the exponent in the relation ξ −β/ν .
It would simply have been enough to assume that P (p, L) ∝ ξ x for some exponent
x in the limit of ξ ⪡ L.
In order to satisfy these conditions, the scaling form of P (p, L) must therefore
be
where
const. for u ⪡ 1
f (u) =
. (6.11)
uβ/ν for u ⪢ 1
Testing the Scaling Ansatz We can now test the scaling ansatz by plotting
P (p, L) according to the ansatz, following a strategy similar to what we developed
for n(s, p). We rewrite the scaling function P (p, L) = L−β/ν f (L/ξ ) in terms of
|p − pc | by inserting ξ = ξ0 |p − pc |−ν :
L β/ν P ( p, L )
0.8
0.6
0.4
0.2
0.0
− 10 −5 0 5
L 1/ν ( p − pc)
Consequently, if we plot L1/ν (p − pc ) along the x-axis and Lβ/ν P (p, L) along
the y-axis, we expect the data from simulations for various L-values to fall onto a
common curve, the curve f (u). This is illustrated in Fig. 6.2, which shows that the
measured data is consistent with the scaling ansatz. We call such as plot a scaling
data collapse plot.
We can characterize the distribution of cluster sizes using moments of the cluster
number distribution. The k-th moment .Mk (p, L) is defined as:
∞
Mk (p, L) =
. s k n(s, p; L) . (6.17)
s=1
We have already introduced the second moment, .M2 (p, L), which we called the
average cluster size, .S(p, L).
∞
S(p, L) = M2 (p, L) =
. s 2 n(s, p; L) . (6.18)
s=1
Now, let us see if we can apply the finite-size scaling approach to develop a scaling
theory for .S(p, L). First, we will measure .S(p, L), and then develop and test a
scaling theory.
How would we measure .S(p, L)? We recall that we measure the cluster number
density from
Ns
n(s, p; L) =
. , (6.19)
Ld
where .Ns is the number of clusters of size s. Thus we can estimate .S(p, L) from:
∞
∞
Ns
. S(p, L) = s 2 n(s, p; L) = s2 . (6.20)
Ld
s=1 s=1
We realize that we can perform this sum by summing over all possible s and then
including how many clusters we have for a given s, or we can alternatively sum over
all the observed clusters .si . (Try to convince yourself that this is the same by looking
at a sequence of clusters of sizes .1, 2, 1, 5, 1, 2.). Thus, we can estimate the second
moment from the sum:
.S(p, L) = si2 /L2 . (6.21)
i
92 6 Finite Size Scaling
And similarly by summing over .sik for the k-th moment. We implement this in the
following program:
import numpy as np
import [Link] as plt
from [Link] import measurements
LL = [25,50,100,200]
p = [Link](0.4,0.75,50)
nL = len(LL)
nx = len(p)
S = [Link]((nx,nL),float)
for iL in range(nL):
L = LL[iL]
M = int(2000*25/L)
for i in range(M):
z = [Link](L,L)
for ip in range(nx):
m = z<p[ip]
lw, num = [Link](m)
labelList = [Link]([Link]() + 1)
area = [Link](m, lw, labelList)
# Remove spanning cluster by setting area to zero
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x>0)]
if (len(perc)>0):
area[perc[0]] = 0
S[ip,iL] = S[ip,iL] + [Link](area*area)
S[:,iL] = S[:,iL]/(L**2*M)
# Plotting the results
[Link](figsize=(6,4))
for iL in range(nL):
L = LL[iL]
lab = "$L="+str(L)+"$"
[Link](p,S[:,iL],label=lab)
[Link](’$S(p,L)$’)
[Link](’$p$’)
[Link]()
How can we understand these plots and how can we develop a theory for .S(p, L)?
We previously found that S diverges as p approaches .pc :
S(p) = S0 |p − pc |−γ ,
. (6.22)
6.3 Average Cluster Size 93
1250 L = 25 1200
L = 50
1000 L = 100 1000
L = 200
750 800
S ( pc ; L )
S ( p, L )
600
500
400
250
200
0 0
0.4 0.5 0.6 0.7 50 100 150 200
p L
(a) (b)
Fig. 6.3 (a) Plot of .S(p, L). (b) Plot of .S(pc ; L) as a function of L
where the exponent is .γ = 43/18 for .d = 2. Following the approach for finite-size
scaling introduced above, we introduce the finite size L through a scaling function
.f (L/ξ ), giving us a finite-size scaling ansatz (our hypothesis):
L
S(p, L) = S0 |p − pc |−γ f
. . (6.23)
ξ
Now, we see from Fig. 6.3 that when .p = pc , .S(pc , L) does not diverge, but is
limited by L, as we would expect for a finite system. Thus we know that in the
limit when .p → pc , .S(p, L) can only depend on L. This implies that the function
.f (L/ξ ) in this limit must be so that the .ξ in .f (L/ξ ) cancels the .ξ
γ /ν in front of it.
γ /ν
L
. S(p, L) ∝ ξ γ /ν
∝ Lγ /ν . (6.25)
ξ
L − γ/ν S ( p, L )
0.06
0.04
0.02
0.00
− 10.0 − 7.5 − 5.0 − 2.5 0.0 2.5 5.0 7.5
L 1/ν ( p − pc)
The resulting plot is shown in Fig. 6.4, which indeed demonstrates that the measured
data is consistent with the scaling theory. Success!
We can measure the percolation probability for a set of finite system sizes using the
methods we developed previously. Here, we have implemented the measurement in
the following program which is very similar to the program developed to measure
.P (p, L)
import numpy as np
import [Link] as plt
from [Link] import measurements
LL = [25,50,100,200]
p = [Link](0.4,0.75,50)
nL = len(LL)
nx = len(p)
6.4 Percolation Threshold 95
Ni = [Link]((nx,nL),float)
Pi = [Link]((nx,nL),float)
for iL in range(nL):
L = LL[iL]
N = int(2000*25/L)
for i in range(N):
z = [Link](L,L)
for ip in range(nx):
m = z<p[ip]
lw, num = [Link](m)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x>0)]
if (len(perc)>0):
Ni[ip,iL] = Ni[ip,iL] + 1
Pi[:,iL] = Ni[:,iL]/N
for iL in range(nL):
L = LL[iL]
lab = "$L="+str(L)+"$"
[Link](p,Pi[:,iL],label=lab)
[Link](’$\Pi(p,L)$’)
[Link](’$p$’)
[Link]()
The resulting plot of .Π (p, L) for various values of L is shown in Fig. 6.5.
Let us now assume that we do not a priori know .pc or any of the scaling exponents.
How can we use this data-set to estimate the value for .pc ?
The simplest approach may be to estimate .pc as the value for p that makes
.Π (p, L) = 1/2. This corresponds to intersection between the horizontal line
1.0 L = 25
L = 50
0.8 L = 100
L = 200
0.6
Π( p, L )
0.4
0.2
0.0
0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75
p
0.600
L = 25
0.575 0.59300
L = 50
L = 100
0.550 0.59275
L = 200
0.525 0.59250
Π( p, L )
p1/ 2
0.500 0.59225
0.475 0.59200
0.450 0.59175
0.425 0.59150
0.400
0.590 0.591 0.592 0.593 0.594 0.595 25 50 75 100 125 150 175 200
p L
(a) (b)
Fig. 6.6 (a) Plot of .Π (p, L). (b) Plot of .p1/2 as a function of L
.Π = 1/2 and the curves in Fig. 6.5. This is illustrated in Fig. 6.6. Here, we
have also plotted .p1/2 as a function of L, where .p1/2 is the value for p so that
.Π (p1/2 , L) = 1/2. These values for .p1/2 are calculated by a simple interpolation
as illustrated in the following program. (Notice that as usual in this book, we do not
aim for high precision in this program. The simulations are for small system sizes
and few samples, but are meant to illustrate the principle and be reproduceable for
you.)
for iL in range(nL):
ipc = [Link](Pi[:,iL]>0.5) # Find first i where Pi>0.5
# Interpolate from ipc-1 to ipc to find intersection
ppc = p[ipc-1] + (0.5-Pi[ipc-1,iL])*\
(p[ipc]-p[ipc-1])/(Pi[ipc,iL]-Pi[ipc-1,iL])
Pic = 0.5
[Link](LL[iL],ppc,’o’)
[Link](’$L$’)
[Link](’$p_{1/2}$’)
From Fig. 6.6 we see that as L increases the value for .p1/2 gradually approaches
pc . Well, we cannot really see that it is approaching .pc , but we guess that it will.
.
We apply the same method as before to develop a theory for .Π (p, L). First. we
notice that at .pc .Π (pc , L) does not either diverge or go to zero. This means that
6.4 Percolation Threshold 97
Π (p, L) cannot be a function of .ξ alone, but instead must have the scaling form:
.
L
Π (p, L) = ξ 0 f
. . (6.29)
ξ
We introduce a new function .Φ(u) = f ξ0 u1/ν :
Π (p, L) = Φ L1/ν (p − pc ) .
. (6.31)
How can we use this theory to estimate .pc ? We follow a technique similar to what
we used above: We find the value .px that makes .Π (px , L) = x. Above, we did this
for .x = 1/2, but we can do this more generally. Actually, as .L → ∞, we expect
any such .px to converge to .pc . We notice from above that .px is a function of L:
.px = px (L).
where it is important to realize that the right hand side, .Cx , is a number which only
depends on x and not on L. We can therefore rewrite this as
px − pc = Cx L−1/ν .
. (6.34)
If we know .ν, we see that this gives a method to estimate the value of .pc .
Figure 6.7 shows a plot of .p1/2 − pc as a function of .L−1/ν for .ν = 4/3. We can
use this plot to extrapolate to find .pc in the limit when .L → ∞ as indicated in the
plot. The resulting value for .pc extrapolated from .L = 50, 100, 200 is .pc = 0.5935,
which is surprisingly good given the small system sizes and small sample sizes used
for this estimate. (The best known value is .pc = 0.5927). This demonstrates the
power of finite size scaling.
98 6 Finite Size Scaling
p1/ 2
.L → ∞ is .pc = 0.5935
0.5920
0.5915
However, this approach depends on us knowing the value for .ν. What if we did not
know neither .ν nor .pc ? How can we estimate both from the scaling ansatz? One
alternative is to generate plots of .px as a function of .L−1/ν for several values of x.
Then we adjust the values of .ν until we get a straight line, in that case we can read
of the intersect with the .px axis as the value for .pc .
However, we can do even better by noticing a trick: For two x values .x1 and .x2 ,
we get
and we can therefore plot .log(dp) as a function of .log(L) to get .ν, and then use
this to estimate .pc . As an exercise, the reader is encouraged to demonstrate that this
scaling ansatz is valid for .d = 1, and in this case find .Cx explicitly.
Exercises
(a) Show that the scaling ansatz for Π (p, L) is valid for d = 1.
(b) Find an explicit expression for Cx for d = 1.
(a) Find px for x = 0.3 and x = 0.8 for L = 25, 50, 100, 200, 400, 800. Plot px as
a function of L.
Exercises 99
(b) Plot log (p0.8 − p0.3 ) as a function of log(L) to estimate the exponent ν. How
does it compare to the exact result?
In the following, please use the exact value ν = 4/3. The scaling theory also
predicted that
px = pc + Cx L−1/ν .
. (6.37)
(a) Develop a finite size scaling ansatz/theory for n(s, pc , L). You should provide
arguments for the behavior in the various limits.
(b) Plot n(s, pc , L) as a function of s for L = 100, 200, 400, 800.
(c) Demonstrate the validity of the scaling theory by producing a data-collapse plot
for n(s, pc , L).
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Renormalization
7
We can conclude that when .p < pc , the system behaves as a uniform, unconnected
system and when .p > pc , the system is uniform and connected.
The argument we have sketched above is the essence of the renormalization
group argument. It is only exactly at .p = pc that an iterative rescaling is a non-
trivial fix point: the system iterates onto itself because it is a self-similar fractal.
When p is away from .pc , rescaling iterations will make the system progressively
more homogeneous, and effectively bring the rescaled p towards either 0 or 1.
In this chapter we will provide an introduction to the theoretical framework for
renormalization. This is a powerful set of techniques, introduced for equilibrium
critical phenomena by Kadanoff [19] in 1966 and by Wilson [39] in 1971. Wilson
later received the Nobel prize for his work on critical phenomena.
We see that the system appears similar throughout the iteration, and the final single
site is occupied. This is because the system is a self-similar fractal and does not
change significantly through the iterations. What happens when .p > pc ? In this
case, we see that the system becomes more homogeneous through each iteration
and eventually the whole system is filled. This means that the effective percolation
probability becomes higher through the iterations. Similarly, when .p > pc , the
system becomes more homogeneous, but also more empty, as the iterations proceed.
c=1 c=2 c=3 c=4 c=5 c=6 c=7 c=8 c=9 c=10 c=11 c=12 c=13 c=14 c=15 c=16
k=1
Fig. 7.1 Illustration of a renormalization rule for a site percolation problem with a rescaling
.b= 2. The top row shows the 16 configurations c. The middle row show the 6 classes k with
multiplicities .g(k)
7.1 The Renormalization Mapping 103
This means that the effective percolation probability becomes lower through the
iterations.
p' = R(p) ,
. (7.1)
where the renormalization function .R(p), which provides the mapping, depends on
the details of the rule used for renormalization.
Selecting a Renormalization Rule There are many choices for the mapping
between the original and the renormalized lattice. We have illustrated a particular
mapping with a rescaling .b = 2 in Fig. 7.1. Such a mapping describes how each of
the .42 = 16 possible configurations c of the .2 × 2 system is mapped onto a .1 × 1
single site through a function .f (c), where .f (c) is 1 if the new site is occupied and
0 if it is empty. The renormalization mapping is then
R(p) =
. P (c)f (c) , (7.2)
c
k is called the multiplicity .g(k) of the class. Expressed in terms of the classes k, the
renormalization mapping is
.R(p) = g(k)P (k)f (k) . (7.3)
k
For the particular mapping provided in Fig. 7.1, the renormalization mapping
becomes
= p4 + 4p3 (1 − p) + 4p2 (1 − p2 ) .
This illustrates a particular rule, but there are many possible rules. Usually, we
want to ensure that important aspects of the percolation system is preserved by the
mapping. For example, we would want the mapping to conserve connectivity. That
is, we would like to ensure that
L
Π (p, L) = Π (p' ,
. ). (7.5)
b
However, even though we may ensure this on the level of the mapping, this does
not ensure that the mapping actually conserves connectivity when applied to a large
cluster. It may, for example, connect clusters that were unconnected in the original
lattice, or disconnect clusters that were connected, as illustrated in Fig. 7.3.
b
Fig. 7.3 Illustration of a single step of renormalization on an .8 × 8 lattice of sites. We see that the
renormalization procedure introduces new connections: the blue cluster is now much larger than
in the original. However, the procedure also removes previously existing connections: the original
yellow cluster is split into two separate clusters
7.1 The Renormalization Mapping 105
a map exists and study its qualitative features. Then we will address detailed
properties of the renormalization mapping through two worked examples.
For any choice of mapping, the rescaling will result in a change in the correlation
length .ξ :
1
ξ ' = ξ(p' ) =
. ξ(p) . (7.6)
b
We will use this relation to address the behavior of fixpoints of the mapping.
Fixpoint A fixpoint of a mapping .R(p) is a point .p∗ that does not change
when the mapping is applied. That is
p∗ = R(p∗ ) .
. (7.7)
Trivial Fixpoints At a fixpoint, the iteration relation for the correlation length
becomes:
ξ(p∗ )
ξ(p∗ ) =
. . (7.8)
b
The only possible solutions for this equation are that .ξ = 0 or .ξ = ∞. We call
the case when .ξ = 0 a trivial fixed point. There are two trivial fixed points for any
renormalization mapping at .p = 0 and at .p = 1.
Stable and Unstable Fixpoints Let us assume that there exists a non-trivial
fixpoint .p∗ , and let us address the behavior for p close to .p∗ . We notice that for any
finite .ξ , iterations by the renormalization relation will reduce .ξ . That is, both for
∗ ∗
.p < p and for .p > p iterations will make .ξ smaller. This implies that iterations
will take the system further away from the non-trivial fixpoint, where the correlation
length is infinite. The non-trivial fixpoint is therefore an unstable fixpoint. Similarly,
for p close to a trivial fixpoint, where .ξ = 0, iterations will decrease .ξ , and the
renormalized system will move closer to the fixpoint in each iteration. The trivial
fixpoint is therefore stable.
R(p)
1
p’>p
p’<p
p* 1 p
Fig. 7.4 Illustration the renormalization mapping .p ' = R(p) as a function of p. The non-trivial
fixpoint .p ∗ = R(p ∗ ) is illustrated. Two iterations sequences are illustrated by the lines with arrows.
Let us look at the path starting from .p > p ∗ . Through the first application of the mapping, we read
off the resulting value of .p ' . This value will then be the input value for the next application of
the renormalization mapping. A fast way to find the corresponding position along the p axis is to
reflect the .p ' value from the line .p ' = p shown as a dotted line. This gives the new p value, and
the mapping is applied again producing yet another .p ' which is even further from .p ∗ . With the
drawn shape of .R(p) there is only one non-trivial fixpoint, which is unstable
We are now ready for a more quantitative argument for the effect of iterations
through the renormalization mapping .R(p). First, we notice that the non-trivial
fixpoint corresponds to the percolation threshold of the renormalization model, since
the correlation length is diverging for this value of p. (This does not imply that .p∗
is equal to .pc . As we shall see, .p∗ depends on the choice of .R(p)).
We will now assume that .R(p) is differentiable, which it should be since .R(p)
is based on sums of polynomials of p and .1 − p. Let us study the behavior close to
∗ '
.p through a Taylor expansion of the mapping .p = R(p). First, we notice that
because .p' = R(p) and .p∗ = R(p∗ ). The Taylor expansion of .R(p) for a p close
to .p∗ is:
p' − p∗ ≃ Λ(p − p∗ ) ,
. (7.11)
We see that the value of .Λ characterizes the fixpoint. For .Λ > 1 the new point .p'
will be further away from .p∗ than the initial point p. Consequently, the fixpoint is
unstable. By a similar argument, we see that for .Λ < 1 the fixpoint is stable. For
.Λ = 1 we call the fixpoint a marginal fixpoint.
Let us now assume that the fixpoint is indeed the percolation threshold. In this
case, when p is close to .pc , we know that the correlation length is
for the renormalized point. We will now use (7.11) for .p∗ = pc , giving
p' − pc = Λ(p − pc ) .
. (7.14)
1
ξ(p' ) =
. ξ(p) . (7.17)
b
Consequently, we have found that
b = Λν .
. (7.18)
This implies that the exponent .ν is a property of the fixpoint of the mapping .R(p).
We can find .ν from
ln b
ν=
. , (7.19)
ln Λ
7.2 Examples
Let us first address the one-dimensional percolation problem using the renormal-
ization procedure. We have illustrated the one-dimensional percolation problem in
Fig. 7.5. We generate the renormalization mapping by ensuring that it conserves
connectivity. The probability for two sites to be connected over a distance b is .pb
when the occupation probability for a single site is p. A renormalization mapping
that conserves connectivity is therefore:
p' = Π (p, b) = pb .
. (7.20)
p∗ = (p∗ )b ,
. (7.21)
0.8
0.6
p’
0.4
0.2
0
0 0.2 0.4 0.6 0.8 1
p
Fig. 7.6 Illustration of a renormalization rule for a one-dimensional site percolation system with
.b=3
We can also apply the theory directly to find the exponent .ν. The renormalization
relation is .p' = R(p) = pb . We can therefore find .Λ from:
∂R
Λ= = b(p∗ )b−1 = b , (7.22)
∂p p∗
.
where we are now studying the unstable fixpoint .p∗ = 1. We can therefore
determine .ν from (7.19):
ln b
ν=
. =1. (7.23)
ln Λ
We notice that b was eliminated in this procedure, which is essential since we do
not want the exponent to depend on details such as the size of renormalization cell.
The result for the scaling of the correlation length is therefore
1
ξ∝
. , (7.24)
1−p
when .1 − p ⪡ 1.
Let us now use this method to address a renormalization scheme for two-
dimensional site percolation. We will use a scheme with .b = 2. The possible
configurations for a .2 × 2 lattice are shown in Fig. 7.7.
110 7 Renormalization
c=1 c=2 c=3 c=4 c=5 c=6 c=7 c=8 c=9 c=10 c=11 c=12 c=13 c=14 c=15 c=16
k=1
Fig. 7.7 Possible configurations for a .2 × 2 site percolation system. The top row indicates various
configurations and the middle row how the configurations are categorized into 7 classes k, where
.g(k) is the number of configurations in class k
where .f (k) = 1 if class k is mapped onto an occupied site and .f (k) = 0 if class k
is mapped onto an empty site. The renormalization relation is illustrated in Fig. 7.8.
We will now follow steps 3 and 4. First, in step 3, we determine the fixpoints of
the renormalization relation. That is, we find the solutions to the equation
The trivial solution .p∗ = 0 is not of interest. Therefore we divide by .p∗ to produce
(p∗ )2 + p∗ − 1 = 0 .
. (7.28)
1.0 R ( p)
p
0.8
0.6
R ( p), p
0.4
0.2
0.0
Fig. 7.8 Plot of the renormalization relation .p ' = R(p) = p 4 + 4p 3 (1 − p) + 2p 2 (1 − p)2 for a
two-dimensional site percolation problem
We have therefore found an estimate of .pc by setting .pc = p∗ . This does not
produce the correct value for .pc in a two-dimensional site percolation system, but
the result is still reasonably correct. We can similarly estimate the exponent .ν by
calculating .R ' (p∗ ).
We will now use the same method to address percolation on site percolation on a
triangular lattice. A triangular lattice is a lattice where each point has six neighbors.
In solid state physics, the lattice is known as the hexagonal lattice because of its
hexagonal rotation symmetry. Site percolation on the triangular lattice is particularly
well suited for renormalization treatment, because a coarse grained version of the
lattice
√ is also a triangular lattice, as illustrated in Fig. 7.9, with a lattice spacing
.b = 3 times the original lattice size.
We will use the majority rule for the renormalization mapping. That is, we will
map a set of three sites onto an occupied site if a majority of the sites are occupied,
meaning that two or more sites are occupied. Otherwise, the renormalized site is
empty. This mapping is illustrated in Fig. 7.9. This mapping does, as the reader may
easily check, on the average conserve connectivity. The renormalization mapping
becomes
p∗ = 3(p∗ )2 − 2(p∗ )3 .
. (7.31)
112 7 Renormalization
We observe that the trivial fixpoints .p∗ = 0 and .p∗ = 1 indeed satisfy (7.31). The
non-trivial fixpoint is .p∗ = 1/2. We are pleased to observe that this is the exact
solution for .pc for site percolation on the triangular lattice.
We can use this relation to determine the scaling exponent .ν. First, we calcu-
late .Λ:
3
Λ = R ' (p∗ ) = 6p(1 − p)|p= 1 =
. . (7.32)
2 2
As a result we find the exponent .ν from
1 ln Λ ln 3/2
. = = √ ≃ 1.355 , (7.33)
ν ln b ln 3
which is very close to the exact result .ν = 4/3 for two-dimensional percolation.
7.2 Examples 113
(a) (d)
k=1 k=2 k=3 k=4 k=5
Fig. 7.10 (a) Illustration of a renormalization scheme for bond percolation on a square lattice
in two dimensions. The rescaling factor is .b = 2. (b) In general, the renormalization involves
a mapping from 8 to two bonds. However, we will consider percolation only in the horizontal
direction. This simplifies the mapping, to the figure shown in (c). For this mapping, the classes are
shown and enumerated in (d)
13
p' = R(p) = Π =
. n(c)P (c)Π |c , (7.34)
c=1
where we have used k to denote the various classes, .P (k) is the probability for
one instance of class k, .n(k) is the number of different configurations due to
symmetry consideration in class k, and .Π |k is the spanning probability given that
the configuration is in class k. The resulting relation is
p' = R(p).
. (7.35)
= p + p (1 − p) + 4p (1 − p) + 2p (1 − p)
5 4 4 3 2
. (7.36)
+2p3 (1 − p)2 + 4p3 (1 − p)2 + 2p2 (1 − p)3. (7.37)
= 2p − 5p + 2p + 2p .
5 4 3 2
(7.38)
114 7 Renormalization
The fixpoints for this mapping are .p∗ = 0, .p∗ = 1, and .p∗ = 1/2. The fixpoint
.p
∗ = 1/2 provides the exact solution for the percolation threshold on the bond
13
Λ = R ' (p∗ ) =
. . (7.39)
8
The corresponding estimate for the exponent .ν is
ln b
. ν= ≃ 1.428 , (7.40)
ln Λ
which should be compared with the exact result of .ν = 4/3 for two-dimensional
percolation.
Exercises 115
Exercises
(a) Find all 28 possible configurations for the 2 × 2 × 2 renormalization cell for
three-dimensional site percolation.
(b) Determine a renormalization scheme - what configurations map onto an occu-
pied site?
(c) Find the renormalized occupation probability p ' = R(p).
(d) Plot R(p) and f (p) = p.
(e) Find the fixpoints p∗ so that R(p∗ ) = p ∗ .
(f) Find the rescaling factor Λ = R ' (p∗ ).
(g) Determine the exponent ν = ln Λ/ ln b.
z
y
x
116 7 Renormalization
# Coarsening procedure
import numpy as np
import [Link] as plt
from [Link] import measurements
def coarse(z,f):
#function zz = coarse(z,f)
# The original array is z
# The transfer function is f given as a vector
# with 16 possible places
# f applied to a two-by-two matrix should return
# the renormalized values
#
# The various values of f correspond to the following
# configurations of the two-by-two region that is renormalized,
# where I have used X to mark a present site, and 0 to mark an
# empty sites
#
# 0 00 4 00 8 00 12 00
# 00 X0 0X XX
#
# 1 X0 5 X0 9 X0 13 X0
# 00 X0 0X XX
#
# 2 0X 6 0X 10 0X 14 0X
# 00 X0 0X XX
#
# 3 XX 7 XX 11 XX 15 XX
# 00 X0 0X XX
#
nx = [Link](z)[0]
ny = [Link](z)[1]
if (nx%2==1): # Must be even number
raise ValueError(’nx must be even’)
if (ny%2==1): # Must be even number
raise ValueError(’ny must be even’)
nx2 = int(nx/2)
ny2 = int(ny/2)
zz = [Link]((nx2,ny2),float) # Generate return matrix
Exercises 117
x = [Link]((2,2),float)
for iy in range(0,ny,2):
for ix in range(0,nx,2):
x = z[ix,iy]*1 + z[ix,iy+1]*2 + \
z[ix+1,iy]*4 + z[ix+1,iy+1]*8
xx = f[int(x)]
zz[int((ix+1)/2),int((iy+1)/2)] = xx
return zz
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Subset Geometry
8
So far, we have studied the geometry of the percolation system. Now, we will
gradually address the physics of processes that occur in a percolation system.
We have addressed one physics-like property of the system, the density of the
spanning cluster, and we found that we could build a theory for the density P as a
function of the porosity (occupation probability) p of the system. In order to address
other physical properties, we need to have a clear description of the geometry of
the percolation system close to the percolation threshold. In this chapter, we will
develop a simplified geometric description that will be useful, indeed essential,
when we discuss physical process in disordered media. We will introduce various
subsets of the spanning cluster—sets that play roles in specific physical processes.
We will start by introducing singly connected bonds, the backbone and dangling
ends and provide a simplified image of the spanning cluster in terms of the blob
model for the percolation system [2, 9, 16, 34].
We will start with an example of a subset of the spanning cluster, the set of singly
connected sites (or bonds). This will demonstrate what we mean by a subset and
how the subset is connected to a physical problem.
Singly Connected Site A singly connected site is a site with the property that
if it is removed, the spanning cluster will no longer be spanning.
We can relate this to a physical property: If we study fluid flow in the spanning
cluster, all the fluid has to go through the singly connected sites. These sites are
also often referred to as red sites, because if we were studying a set of random
Fig. 8.1 Illustration of the spanning cluster, the singly connected bonds (red), the backbone (blue),
and the dangling ends (green) for a .256 × 256 bond percolation system at .p = pc . (Figure from
Martin Søreng)
resistors, the highest current would have to go through the singly connected bonds,
and they would therefore heat up and become “red”. Several examples of subsets of
the spanning cluster, including the singly connected bonds, are shown in Fig. 8.1.
Scaling Hypothesis We have learned that the spanning cluster may be described
by the mass scaling relation .M ∝ LD , where D is termed the fractal dimension of
the spanning cluster. Here, we will make a daring hypothesis, which we will also
substantiate: We propose that subsets of the spanning cluster obey similar scaling
relations.
For example, we propose that the mass of the singly connected sites (.MSC ) has
the scaling form
MSC ∝ LDSC ,
. (8.1)
where we call the dimension .DSC the fractal dimension of the singly connected
sites. Because the set of singly connected sites is a subset of the spanning cluster,
we know that .MSC ≤ M. It therefore follows that
DSC ≤ D .
. (8.2)
8.2 Self-Avoiding Paths on the Cluster 121
Based on this simple example, we will generalize the approach to other subsets of
the spanning cluster. However, first we will introduce a new concept, a self-avoiding
path on the spanning cluster.
The study of percolation is the study of connectivity, and many of the physical
properties that we are interested in depends on various forms of connecting paths
on the spanning cluster between two opposite edges. We can address the structure
of connected paths between the edges by studying self-avoiding paths (SAPs). A
Self-Avoiding Path (SAP) is a set of connected sites that correspond to the sites on
the path of a walk on the spanning cluster that does not intersect itself going from
one side to the opposite side.
Minimal Path
The shortest path between the two edges is called the shortest SAP between the two
edges. (Notice, that there may be more than one path the satisfy this criterion. We
chose one of these paths randomly). We call this the minimal path and denote its
length .Lmin . The length here refers to the number of sites in the path, which we also
call the mass of the path, .Mmin = Lmin . We will use mass instead of length in the
following to describe the paths.
We assume that mass of the minimal path also scales with the system size
according to the scaling form:
Where we have introduced the scaling exponent of the minimal path to be .Dmin .
Similarly, we call the longest SAP between the two edges the longest path with a
mass .Mmax . Again, we assume that the mass has a scaling form, .Mmax ∝ LDmax . We
notice that .Mmin ≤ Mmax . Consequently, a similar relation holds for the exponents
.Dmin ≤ Dmax .
We also introduce the term the average path, meaning the average mass (length)
of all possible SAPs going between opposite sides of the system, .〈MSAP 〉 ∝ LDSAP .
The dimension .DSAP will lie between the dimensions of the minimal and the
maximal path.
122 8 Subset Geometry
Backbone
Intersection of All Self-Avoiding Paths The notion of SAPs can also be used
to address the physical properties of the cluster, such as we saw for the singly
connected bonds. The set of singly connected bonds is the set of intersections
between all SAPs connecting the two sides. That is, the singly connected bonds
is the set of points that any path must go through in order to connect the two sides.
From this definition, we notice that the dimension .DSC < Dmin , and as we will see
further on, .DSC = 1/ν which is smaller than 1 for two-dimensional systems.
Union of All Self-Avoiding Paths Another useful set is the union of all SAPs that
connect the two edges of the cluster. This set is called the backbone with dimension
DB .
This set has a simple physical interpretation for a random porous material, since it
corresponds to the sites that are accessible to fluid flow if a pressure is applied across
the material. The remaining sites are called dangling ends. The backbone are all the
sites that have at least two different paths leading into them, one path from each side
of the cluster. The remaining sites only have one (self-avoiding) path leading into
them, and we call this set of sites the dangling ends. The spanning cluster consists
of the backbone plus the dangling ends, as illustrated in Fig. 8.2. The dangling ends
are therefore pieces of the cluster that can be cut away by the removal of a single
bond.
Fig. 8.2 Illustration of the spanning cluster consisting of the backbone (red) and the dangling ends
(blue) for a 512 × 512 site percolation system for (a) p = 0.58, (b) p = 0.59, and (c) p = 0.61
8.2 Self-Avoiding Paths on the Cluster 123
Generally, we will find that the dimension of the backbone, .DB , is smaller than the
dimension of the spanning cluster. For example, in two dimensions, we find that
.DB ≃ 1.6, whereas .D ≃ 1.89. This has implications for the relative size of the
where .M0 and .M0,B are constant prefactors. To see what happens when .L → ∞,
we divide by M:
Since .DB ≤ D, we see that the fraction .MDE /M goes to a constant (one) as L
approaches infinity. Consequently, we have found that .MDE ∝ M ∝ LD . This also
implies that as the system size goes to infinity most of the mass is in the dangling
ends. This means that the backbone occupies a smaller and smaller portion of the
total mass of the system as the system size increases.
We can provide a better argument for why the various subsets should scale with the
system size L to various exponents. We notice that the following relation between
the masses must be true:
where the first inequality .L1 ≤ Mmin follows because even the minimum path must
be at least of length L to go from one side to the opposite side.
Now, if this is to be true for all values of L, it can be argued that because all
the masses are between two scaling relations, .L1 and .Ld , also the scaling of the
124 8 Subset Geometry
Let us now try to formulate our geometric description of the spanning cluster into a
model of the spanning cluster [36]. We have found that the spanning cluster can be
subdivided first into two parts: the backbone and the dangling ends. The backbone
may again be divided into two parts: a set of blobs where the are several parallel
paths and a set of sites, the singly connected sites, that connect the blobs to each
other and the blobs to the dangling ends. Thus, we have ended up with a model with
three components:
Each of the blobs and the dangling ends will again have a similar substructure of
dangling ends, blobs with parallel paths, and singly connected bonds as illustrated
in Fig. 8.3. This cartoon image of the clusters provides very useful intuition about
the geometrical structure of percolation clusters, which we will use when we address
the physics of disordered systems in the next chapters.
The exponents can be calculated either by numerical simulations, where the masses
of the various subsets are measured as a function of system size at .p = pc , or by the
renormalization group method. Numerical results based on computer simulations
using the code provided in this book are listed in Table 8.1. You can find up-to-date
Table 8.1 A list of known exponent for the various subset types in two dimensions
d .DSC .Dmin .Dmax .DB D .DDE
results for exponents in the percolation system at the Wikipedia page: [Link]
[Link]/wiki/Percolation_critical_exponents.
We will now use the renormalization group approach to address the scaling exponent
for various subsets of the spanning cluster at .p = pc . For this, we will here
use the renormalization procedure for bond percolation on a square lattice in
two dimensions following Hong and Stanley [17], where we have found that the
renormalization procedure produces the exact result for the percolation threshold,
∗
.pc = p = 1/2, which is a fixpoint of the mapping.
Our strategy will be to assume that all the bonds have a mass .M = 1 in the
original lattice, and then find the mass .M ' in the renormalized lattice, when the
length has been rescaled by b. For a property that displays a self-similar scaling, we
will expect that
M ' ∝ bDx M ,
. (8.8)
where .Dx denotes the dimension for the particular subset we are looking at. We can
use this to determine the fractal exponent .Dx from
ln M ' /M
Dx =
. . (8.9)
ln b
We will do this by calculating the average value of the mass of the H-cell, by
taking the mass of the subset we are interested in for each configuration, .Mx (c),
and multiplying it by the probability of that configuration, summing over all
configurations:
. 〈M〉 = Mx (c)P (c) . (8.10)
c
We have now calculated the average mass in the original 2 by 2 lattice, and this
should correspond to the average renormalized mass, .〈M ' 〉 = p' M ' , which is the
mass of the renormalized bond, .M ' multiplied with the probability for that bond to
be present .p' . That is, we will find .M ' from:
p' M ' =
. M(c)P (c) , (8.11)
c
126 8 Subset Geometry
Table 8.2 Numerical exponents for the exponent describing various subsets of the spanning
cluster defined using the set of Self-Avoiding Walks going from one side to the opposite side of the
cluster. The last line shows the exponents found from numerical simulations in a two-dimensional
system
c P(c) MSC Lmin LAV G Lmax MBB M
1 p5 (1 − p)0 0 2 2.5 3 5 5
2 p4 (1 − p)1 0 2 2 2 4 4
4 2p3 (1 − p)2 2 2 2 2 2 3
5 2p3 (1 − p)2 3 3 3 3 3 3
6 4p3 (1 − p)2 2 2 2 2 2 3
7 2p2 (1 − p)3 2 2 2 2 2 2
〈Mx 〉 26/25 34/25 36.5/25 39/25 47/25 53/25
ln 13 ln 17 ln 36.5 ln 39 47
ln 16 53
ln 16
Dx 8
ln 2
8
ln 2
16
ln 2
16
ln 2 ln 2 ln 2
Dx 0.7004 1.0875 1.1898 1.2854 1.5546 1.7279
Dx,n 3/4 1.13 1.4 1.6 91/48
We have found that we can calculate the behavior of infinite-dimensional and one-
dimensional systems exactly. However, for finite dimensions such as for .d = 2 or
.d = 3 we must rely on numerical simulations and renormalization group arguments
to determine the exponents and the behavior of the system. However, in order to
learn about physical properties in systems with scaling behavior, we may be able
to construct simpler models that contain many of the important features of the
percolation cluster. For example, we may be able to introduce deterministic, iterative
fractal structures that reproduce many of the important properties of the percolation
cluster at .p = pc , but that are deterministic systems. The idea is that we can use
such a system to study other properties of the physics on fractal structures.
8.4 Deterministic Fractal Models 127
Fig. 8.4 Illustration of first three generations of the Mandelbrot-Given curve. The length is scaled
by a factor .b = 3 for each iteration, and the mass of the whole structure is increased by a factor of
8. The fractal dimension is therefore .D = ln 8/ ln 3 ≃ 1.89
.m(l) = 8 , and the linear size of the cluster is .L(l) = 3 . If we assume a scaling on
l l
ln 8
D=
. ≃ 1.89 . (8.12)
ln 3
This is surprisingly similar to the fractal dimension of the percolation cluster. We
can also look at other dimensions, such as for the singly connected bonds, the
minimum path, the maximum path and the backbone.
Single Connected Bonds Let us first address the singly connected bonds. In the
zero’th generation, the system is simply a single bond, and the length of the singly
connected bonds, .LSC is 1. In the first generation, there are two bonds that are
singly connecting, and in the second generation there are four bonds that are singly
connecting. The general relation is that
.LSC = 2l , (8.13)
where l is the generation of the structure. The dimension, .DSC , of the singly
connected bonds is therefore .DSC = ln 2/ ln 3 ≃ 0.63, which should be compared
with the exact value .DSC = 3/4 for two-dimensional percolation.
Minimum Path The minimum path will for all generations be the path going
straight through the structure, and the length of the minimal path will therefore
be equal to the length of the structure. The scaling dimension .Dmin is therefore
.Dmin = 1.
128 8 Subset Geometry
Maximum Path The maximum path increases by a factor 5 for each iteration. The
dimension of the maximum path is therefore .Dmax = ln 5/ ln 3 ≃ 1.465.
Backbone We can similarly find that the mass of the backbone increases by a
factor 6 for each iteration, and the dimension of the backbone is therefore .DB =
ln 6/ ln 3 ≃ 1.631.
Model System This deterministic iterative fractal can be used to perform quick
calculations of various properties on a fractal system, and may also serve as a useful
hierarchical lattice on which to perform simulations when we are studying processes
occurring on a fractal structure.
8.5 Lacunarity
The fractal dimension describes the scaling properties of structures such as the
percolation cluster at .p = pc . However, structures that have the same fractal
dimension, may have a very different appearance. As an example, let us study
several variations of the Sierpinski gasket introduced in Sect. 5.3. As illustrated in
Fig. 8.5, we can construct several rules for the iterative generation of the fractal
that all result in the same fractal dimension, but have different visual appearance.
The fractal dimension .D = ln 3/ ln 2 for both of the examples in Fig. 8.5, but by
increasing the number of triangles that are used in each generation, the structures
become more homogeneous. How can we quantify this difference?
1:
2:
Fig. 8.5 Two versions of the Siepinski gasket. In version 1, the next generation is made from 3
of the structures from the last generation, and the spatial rescaling is by a factor .b = 3. In version
2, the next generation is made from 9 of the structures from the last generation, and the spatial
rescaling is by a factor .b = 6. The resulting fractal dimension is .D2 = ln 9/ ln 4 = ln 32 / ln 22 =
ln 3/ ln 2 = D1 . The two structures therefore have the same fractal dimension. However, version 1
have large fluctuations that version 2
8.5 Lacunarity 129
sizes. We can characterize and measure the fractal dimension of a fractal structure
using box-counting, as explained in Sect. 5.3. The structure, such as the percolation
cluster, is divided into boxes of size .𝓁. In each box, i, there will be a mass .mi (𝓁).
The fractal dimension is found by calculating the average mass per box of size .𝓁:
However, there will variations in the masses .m(𝓁) in the boxes, characterized by a
distribution .P (m, 𝓁), which gives the probability for mass m in a box of size .𝓁. We
can characterize this distribution by its moments:
where this particular scaling form implies that the structure is unifractal: the scaling
exponents for all the moments are linearly related.
We can determine the unknown scaling exponent x from the scaling of the zero’th
moment, that is, from the normalization of the probability density: .〈m0 〉 = 1 implies
that .D(0 + 1) + x = 0, and therefore, that .x = −D. The scaling ansatz for the
distribution of masses is therefore
m
P (m, 𝓁) = 𝓁−D f (
. ). (8.21)
𝓁D
And we found that the moments can be written as
.〈m 〉 = 𝓁 x k f (x)dx = Ak 𝓁kD ,
k D(k+1)−D
(8.22)
130 8 Subset Geometry
This number describes another part of the mass distribution relation than the scaling
relation, and can be used to characterize fractal set. For the percolation problem,
this number is assumed to be universal, independent of lattice type, but dependent
on the embedding dimensionality of the system.
Exercises
Exercise 8.1 (Singly Connected Bonds) Use the example programs from the text
to find the singly connected bonds.
(a) Run the programs to visualize the singly connected bonds. Can you understand
how this algorithms finds the singly connected bonds? Why are some of the
bonds of a different color?
(b) Find the mass, MSC , of the singly connected bonds as a function of system size
L for p = pc and use this to estimate the exponent DSC : MSC ∝ LDSC .
(c) Can you find the behavior of PSC = MSC /Ld as a function of p − pc ?
(a) Use the following programs to generate and illustrate of the singly connected
bonds for a 100 × 100 system. Check that the illustrated bonds correspond to
the singly connected bonds.
import numpy as np
import numba
@[Link](cache=True)
def walk(z):
# Left turning walker
# Returns left: nr of times walker passes a site
# First, ensure that array only has one contact point at
# left and right : topmost points chosen
nx = [Link][0]
ny = [Link][1]
i = [Link](z[0,:] > 0)
ix0 = 0 # starting row for walker is always 0
iy0 = i[0][0] # starting column (first element where
# there is a matching column which is zero)
# First do left-turning walker
directions = [Link]((4,2))
directions [0,0] = -1 # west
directions [0,1] = 0
directions [1,0] = 0 # south
directions [1,1] = -1
directions [2,0] = 1 # east
directions [2,1] = 0
directions [3,0] = 0 # north
directions [3,1] = 1
nwalk = 1
ix = ix0
iy = iy0
direction = 0 # 0 = west, 1 = south, 2 = east, 3 = north
left = [Link]((nx,ny))
right = [Link]((nx,ny))
while (nwalk >0):
132 8 Subset Geometry
left[ix,iy] = left[ix,iy] + 1
# Turn left until you find an occupied site
nfound = 0
while (nfound==0):
direction = direction - 1
if (direction < 0):
direction = direction + 4
# Check this direction
iix = ix + int(directions[direction,0])
iiy = iy + int(directions[direction,1])
if (iix >= nx):
nwalk = 0 # Walker escaped
nfound = 1
iix = nx
ix1 = ix
iy1 = iy
# Is there a site here?
elif(iix >= 0):
if(iiy >= 0):
if (iiy < ny):
if (z[iix,iiy]>0): # site present
ix = iix
iy = iiy
nfound = 1
direction = direction + 2
if (direction > 3):
direction = direction - 4
#left
nwalk = 1
ix = ix0
iy = iy0
direction = 1 # 1=left, 2 = down, 3 = right, 4 = up
while(nwalk >0):
right[ix,iy] = right[ix,iy] + 1
# ix,iy
# Turn right until you find an occupied site
nfound = 0
while (nfound==0):
direction = direction + 1
if (direction > 3):
direction = direction - 4
# Check this directionection
iix = ix + int(directions[direction,0])
iiy = iy + int(directions[direction,1])
if (iix >= nx):
if (iy >= iy1):
nwalk = 0 # Walker escaped
nfound = 1
iix = nx
# Is there a site here?
elif(iix >= 0):
if(iiy >= 0):
if (iiy < ny):
Exercises 133
(e) (Advanced) Develop a theory for the behavior of PH (p, L), the probability for
a site to belong to the Hull as a function of p and L for p > pc .
(f) (Advanced) Measure the behavior of PH (p, L) as a function of p for L =
512 × 512.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Flow in Disordered Media
9
We have now developed the tools to address the statistical properties of the geometry
of a disordered system such as a model porous medium: the percolation system. In
the following chapters, we will apply this knowledge to address physical properties
of disordered systems and to study physical processes in disordered materials.
We have learned that the geometry of a disordered system displays fractal scaling
close to the percolation threshold. Material properties such as the density of singly
connected sites, or the backbone of the percolation cluster, display self-similar
scaling. The backbone is the part of the spanning cluster that participates in fluid
flow. The mass, .MB , of the backbone scales with the system size, L, according to
the scaling relation .MB = LDB , where .DB is smaller than the Euclidean dimension.
The density of the backbone therefore decreases with system size. This implies that
material properties which we ordinarily would treat as material constants, depend on
the size of the sample. In this part we will develop an understanding of the origin of
this behavior, and show how we can use the tools from percolation theory to address
the behavior in such systems.
The behavior of a disordered system can in principle always be addressed by
direct numerical simulation. For example, for incompressible, single-phase fluid
flow through a porous material, the effective permeability of a sample can be found
to very good accuracy from a detailed numerical model of fluid flow through the
system. However, it is not practical to model fluid flow down to the smallest scale
in more applied problems. We would therefore need to extrapolate from the small
scale to the large scaling. This process, often referred to as up-scaling, requires that
we know the scaling properties of our system. We will address up-scaling in detail
in the following chapters.
We may argue that a system at the percolation threshold is anomalous and that
any realistic system, such as a geological system, would be far away from the
percolation threshold. In this case, the system will only display an anomalous, size-
dependent behavior up to the correlation length, and over larger lengths the behavior
will be that of a homogeneous material. We should, however, be aware that many
physical properties are described by broad distributions of material properties, and
this will lead to a behavior similar to the behavior close to the percolation threshold,
as we will discuss in detail in this part. In addition, several physical processes ensure
that the system is driven into or is exactly at the percolation threshold. One such
example is the invasion-percolation process, which gives a reasonable description of
oil-water emplacement processes such as secondary oil migration. For such systems,
the behavior is well described by the scaling theory we have developed.
In this and following chapters, we will first provide an introduction to the scaling
of material properties such as conductivity (Chap. 9), elasticity (Chap. 10) and
diffusion (Chap. 11). Then we will demonstrate how processes occurring in systems
with frozen disorder, such as a porous material, often lead to the formation of fractal
structures (Chap. 12).
We will start our studies of physics in disordered media by addressing flow, either
in the form of incompressible fluid flow in a random, porous system or in the form
of electric current in a random, porous materials. First, let us address the similarities
between these two flow phenomena.
We can also use fluid flow in porous medium as our basic physical system. If done
in the limit of slow, incompressible fluid flow these two systems are practically
identical. For fluid flow in a porous medium of length L and cross-sectional area
A, the system is described by Darcy’s law which provide a relation between, .Φ, the
amount of fluid volume flowing through a given cross-sectional area, A, per unit
time and the pressure drop .Δp across the sample:
kA Δp
.Φ= , (9.2)
η L
where k is the called permeability of the material and is a property of the material
geometry, and .η is the viscosity of the fluid. Again, we would like a description so
that k is material property, and all the information about the geometry of the material
goes into the permeability of the sample through the length L and the cross-sectional
area A. Generalized to a d-dimensional system, the relation is
kLd−1 k
Φ=
. Δp = Ld−2 Δp . (9.3)
ηL η
From this, we see that the electric conductivity problem in this limit is the same
as the Darcy-flow permeability problem, where .Δp/L corresponds to the voltage
138 9 Flow in Disordered Media
Let us first address the conductance of a .Ld percolation system. The system may be
either a site or a bond percolation system, but many of the concepts we introduce
are simpler to explain if we just consider a bond percolation system.
We will start with a simplified system: a network of bonds that are present with
probability p. We assume that all bonds have the same conductance, which we can
set to 1 without loss of generality. The bonds are removed with probability .1 − p,
and we model this by setting the conductance of a removed bond to be zero.
The conductance of the .Ld sample is found by solving the flow problem illustrated
in Fig. 9.1. A potential difference V is applied across the whole sample, and we
find (measure) the resulting current I . We find the conductance from Ohm’s law (or
(a) (b)
Ii,j
Vi Vj
Fig. 9.1 (a) Illustration of flow through a bond percolation system. The bonds shown in red are
the singly connected bonds: all the flux has to go through these bonds. The bonds shown in blue
are the rest of the backbone: The flow only takes place on the singly connected bonds and the
backbone, the remaining bonds are the dangling ends, which do not participate in fluid flow. (b)
Illustration of the potentials .Vi and .Vj in two adjacent sites and the current .Ii,j from site i into
site j
9.3 Conductance of a Percolation Lattice 139
I
I = GV
. ⇒ G= . (9.4)
V
Local Potentials and Currents Let us look at the network in more detail.
Figure 9.1b illustrates a small part of the whole system. The two adjacent sites i
and j are connected with a bond of conductance .Gi,j . If the bond is present (with
probability p in the percolation system), the conductance is .Gi,j is 1, otherwise it is
zero.
The current from site i to site j is related to the difference in potential between
the two sites:
Ii,j = Gi,j Vi − Vj ,
. (9.5)
where we notice that the current is positive if the potential is higher in site i than in
site j .
In electromagnetism this is called Kirchhoff’s rule for currents. We can rewrite this
in terms of the local potentials .Vi instead by inserting (9.5) in (9.6):
. Gi,k (Vi − Vk ) = 0 . (9.7)
k
This provides us with a set of equations for all the potentials .Vi , which we must solve
to find the potentials and hence the currents between all the sites in a percolation
system.
Finding Currents and Potentials We can use this to find all the potentials for
a percolation system. Let us address a two-dimensional system of size .L × L.
The potential in a position .(x, y) on the lattice is .V (x, y), where x and y are
140 9 Flow in Disordered Media
This is effectively a set of .Ld equations for .Vi . In addition we have the boundary
conditions that .V (0, j ) = V and .V (L − 1, j ) = 0 for .j = 0, 1, . . . , L − 1. This
defines the system as a tri-diagonal set of linear equations that can be solved easily
numerically.
Computational Methods
50 50
100 100
150 150
200 200
250 250
300 300
350 350
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350
Current BB and DE
0 0
50 50
100 100
150 150
200 200
250 250
300 300
350 350
0 50 100 150 200 250 300 350 0 50 100 150 200 250 300 350
Fig. 9.2 Plots of the spanning cluster, the potential, .V (x, y), the absolute value of the current
flowing into each site, and the singly connected bonds, the backbone and the dangling ends
function coltomat transforms the resulting array of potentials back into a matrix
form, .V (x, y). The following programs are used to calculate the potentials and
currents and visualize the results. The resulting plots are shown in Fig. 9.2.
import numpy as np
import [Link] as plt
from [Link] import spdiags, dia_matrix, coo_matrix
from [Link] import spsolve
from [Link] import measurements
V_out = 0.
# Calls MK_EQSYSTEM .
B,C = MK_EQSYSTEM (A , X , Y )
# Kirchhoff ’ s equations solve for V
V = spsolve(B, C)
# The pressure at the external sites is added
# ( Boundary conditions )
V = [Link]((V_in*[Link](X),V,V_out*[Link] (X)))
# Calculate Ceff
# second-last X elements of V multiplied with second-last
# elem. of A, these are the second last column of the
# system gives the conductivity of the system per row
Ceff = [Link]((V[-1-2*X:-1-X]-V_out).T,A[-1-2*X:-1-X,1]) \
/ ( V_in - V_out )
return V , Ceff
def MK_EQSYSTEM (A , X , Y ):
# Total no of internal lattice sites
sites = X *( Y - 2)
# Allocate space for the nonzero upper diagonals
main_diag = [Link](sites)
upper_diag1 = [Link](sites - 1)
upper_diag2 = [Link](sites - X)
# Calculates the nonzero upper diagonals
main_diag = A[X:X*(Y-1), 0] + A[X:X*(Y-1), 1] + \
A[0:X*(Y-2), 1] + A[X-1:X*(Y-1)-1, 0]
upper_diag1 = A [X:X*(Y-1)-1, 0]
upper_diag2 = A [X:X*(Y-2), 1]
main_diag[[Link](main_diag == 0)] = 1
# Constructing B which is symmetric , lower=upper diagonals
B = dia_matrix ((sites , sites)) # B *u = t
B = - spdiags ( upper_diag1 , -1 , sites , sites )
B = B + - spdiags ( upper_diag2 ,-X , sites , sites )
B = B + B.T + spdiags ( main_diag , 0 , sites , sites )
# Constructing C
C = [Link](sites)
# C = dia_matrix ( (sites , 1) )
C[0:X] = A[0:X, 1]
C[-1-X+1:-1] = 0*A [-1 -2*X + 1:-1-X, 1]
return B , C
def sitetobond ( z ):
# Function to convert the site network z(L,L) into a
# (L*L,2) bond network
# g [i,0] gives bond perpendicular to direction of flow
9.3 Conductance of a Percolation Lattice 143
while (len(perc)==0):
ncount = ncount + 1
if (ncount >100):
break
z=[Link](lx,ly)<p
lw,num = [Link](z)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x > 0)]
print("Percolation attempt", ncount)
zz = [Link]((lw == perc[0]))
# zz now contains the spanning cluster
zzz = zz.T # Transpose
g = sitetobond ( zzz ) # Generate bond lattice
V, c_eff = FIND_COND (g, lx, ly) # Find conductivity
x = coltomat ( V , lx , ly ) # Transform to nx x ny lattice
V = x * zzz
g1 = g[:,0]
g2 = g[: ,1]
z1 = coltomat( g1 , lx , ly )
144 9 Flow in Disordered Media
z2 = coltomat( g2 , lx , ly )
# Plot results
[Link](figsize=(16,16))
ax = [Link](2,2,1)
[Link](zzz, interpolation=’nearest’)
[Link]("Spanning cluster")
[Link](2,2,2, sharex=ax, sharey=ax)
[Link](V, interpolation=’nearest’)
[Link]("Potential")
We can now use this program to measure the conductance .G(p, L) of the system
and how it varies with both p and L. The idea is to calculate G from .G = I /V ,
where we select V and calculate I using the program. We find I as the sum of
9.3 Conductance of a Percolation Lattice 145
all the currents escaping (or entering) the system. In the program, we have set the
potential on the left side to be 1. We recall that we describe positions with the index
.j = x + yL, where the left side corresponds to .x = 0 and therefore .j = yL, which
we write at iL with .i = y is the position along the y-axis. The potentials along
the left side are therefore .ViL = 1 for .i = 0, 1, . . . , L − 1. The conductance from
site iL into a site on the right, that is, into .iL + 1, is .Gi L, iL + 1, which is 1. The
current into the system from the left side, that is from site iL into site .iL + 1, is
.IiL,iL+1 = GiL,iL+1 (ViL − ViL+1 ). The total current I into the system is therefore:
L−1
L−1
L−1
I=
. IiL,iL+1 = GiL,iL+1 (ViL − ViL+1 ) = (ViL − ViL+1 ) .
i=0 i=0 i=0
(9.14)
We use the following program to find the conductance, .G(p, L), for an .L×L system
for .L = 400, as well as the density of the spanning cluster .P (p, L).
import numpy as np
import [Link] as plt
from [Link] import measurements
Lvals = [400]
pVals = [Link](np.log10(0.58), np.log10(0.85), 20)
C = [Link]((len(pVals),len(Lvals)),float)
P = [Link]((len(pVals),len(Lvals)),float)
nSamples = 600
G = [Link](len(Lvals))
for iL in range(len(Lvals)):
L = Lvals[iL]
lx = L
ly = L
for pIndex in range(len(pVals)):
p = pVals[pIndex]
ncount = 0
for j in range(nSamples):
ncount = 0
perc = []
while (len(perc)==0):
ncount = ncount + 1
if (ncount > 1000):
print("Couldn’t make percolation cluster")
break
z=[Link](lx,ly)<p
lw,num = [Link](z)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x > 0)]
if len(perc) > 0: # Found spanning cluster
area = [Link](z, lw, perc[0])
P[pIndex,iL] = P[pIndex,iL] + area # Find P(p,L)
zz = [Link]((lw == perc[0])) # zz=spanning
zzz = zz.T
146 9 Flow in Disordered Media
The resulting behavior for .L = 400 and .M = 600 different realizations is shown
in Fig. 9.3. We observe two things from this plot: First we see that the behaviors of
.G(p, L) and .P (p, L) are qualitatively different around .p = pc : .P (p, L) increases
very rapidly as .(p − pc )β where .β is less than 1. However, it appears that .G(p, L)
increases more slowly. Indeed, from the plot it looks as if .G(p, L) increases as
.(p − pc ) with an exponent x that is larger than 1. How can this be? Why does the
x
density of the spanning cluster increase very rapidly, but the conductance increases
much slower? This may be surprising, but we will develop an explanation for this
observation in the following.
G
0.8
P
0.6
G, P
0.4
0.2
0.0
0.60 0.65 0.70 0.75 0.80 0.85
p
Fig. 9.3 Plots of the conductance .G(p, L) and the density of the spanning cluster .P (p, L) for
.L= 400
9.4 Scaling Arguments for Conductance and Conductivity 147
We will now use the same scaling techniques we introduced to find the behavior of
P (p, L) to develop a theory for the conductance .G(p, L). First, we realize that we
.
L
G(ξ, L) = ( )d−2 G(ξ, ξ ) ,
. (9.15)
ξ
148 9 Flow in Disordered Media
G(ξ, ξ )
g(ξ, L) = L−(d−2) G(ξ, L) =
. . (9.16)
ξ d−2
What is .G(ξ, ξ )? A system with correlation length equal to the system size is
indistinguishable from a system at .p = pc . The conductance .G(ξ, ξ ) is therefore
the conductance of the spanning cluster at .p = pc in a system of size .L = ξ .
Let us therefore find the conductance of a finite system of size L at the percolation
threshold.
What is the conductance, .G(∞, L), of the spanning cluster at .p = pc ? We know that
the spanning cluster consists of the backbone and the dangling ends, and that only
the backbone will contribute to the conductivity. The backbone can be described
by the blob model (see Sect. 8.2 for a discussion of the blob model): The backbone
consists of blobs of bonds in parallel, and links of singly connected bonds between
them.
Starting from a scaling hypothesis for the conductance, we will derive the
consequences of this assumption, and then test these consequences to see if the
data is consistent with our assumption. This will then corroborate the hypothesis.
Our scaling hypothesis will be that the conductance of a system of size L at .pc can
be described by a scaling exponent .ζ̃R :
G(∞, L) ∝ L−ζ̃R .
. (9.17)
Finding Bounds for the Scaling Behavior In many cases, we cannot find the
scaling exponents directly, but we may be able to find bounds for the scaling
exponents. We will pursue this approach here. We will find bounds for the scaling
of .G(∞, L), and use them to determine bounds for the exponent .ζ̃R .
Lower Bound for the Scaling Exponent First, we know that the spanning cluster
consists of blobs in series with the singly connected bonds. This implies that the
resistivity .R = 1/G of the spanning cluster is given as the resistivity of the singly
connected bonds .RSC plus the resistivity of the blobs, .Rblob since resistances are
added for a series of resistances:
This implies that .R > RSC . The singly connected bonds are connected in series,
one after another. Their total resistance is the sum of the resistances of each bond,
which is the resistance of a single bond, multiplied with the number of sites, .MSC .
Because .Rblob is positive, we see from (9.18) that
Because .MSC ∝ LDSC , and we have assumed that .R ∝ Lζ̃R , we find that
Upper Bound for the Scaling Exponent We can find an upper bound by exam-
ining the minimal path. The resistance of the spanning cluster will be smaller than
or equal to the resistance of the minimal path, since the spanning cluster will have
some regions, the blobs, where there are bonds in parallel. Adding parallel bonds
will always lower the resistance. Hence, the resistance is smaller than or equal to
the resistance of the minimal path. Since the minimal path is a series of resistances
in series, the total resistance of the minimal path is the mass of the minimal path
multiplied by the resistance of a single bond. Consequently, the resistance of the
spanning cluster is smaller than the mass of the minimal path, .Mmin , which we
know scales with system size, .Mmin ∝ LDmin . We have therefore found an upper
bound for the exponent
and therefore
ζ̃R ≤ Dmin .
. (9.22)
Upper and Lower Bound Demonstrate the Scaling Relation We have therefore
demonstrated (or proved) the scaling relation
Because this scaling relation shows that the scaling of R is bounded by two power-
laws in L, we have also proved that the resistance R is a power-law, and that the
exponents are within the given bounds. We notice that when the dimensionality of
the system is high, the probability of loops will be low, and blobs will be unlikely.
In this case
We use this to find an expression for .G(ξ, ξ ), which is the conductance of the
spanning cluster at .p = pc in a system of size .L = ξ , by inserting .L = ξ in
(9.25):
G(ξ, ξ ) ∝ ξ −ζ̃R .
. (9.26)
We insert this in (9.16) in order to establish the behavior of the conductivity, g, for
p > pc , finding that
.
G(ξ, ξ )
g=
. ∝ ξ −(d−2+ζ̃R ). (9.27)
ξ d−2
∝ (p − pc )ν(d−2+ζ̃R ) ∝ (p − pc )μ (9.28)
μ = ν(d − 2 + ζ̃R ) .
. (9.29)
We notice that for two-dimensional percolation, any value of .ζ̃R larger than .1/ν
will lead to .μ > 1, which was what was observed in Fig. 9.3. The exponent .μ is
therefore larger than 1, which is significantly different from the exponent .β, which
is less than 1, which describes the mass of the spanning cluster.
Another theoretical approach to address and understand the behavior of the system
is through the renormalization calculation. Here, we will use the renormalization
approach for a square bond lattice in order to estimate the exponent .ζ̃R .
In order to apply the renormalization approach, we calculate the average
resistance .〈R〉 of a .2×2 cell. We use the H-cell approach and only look at percolation
in the horizontal direction. The various configurations c and their degeneracy .g(c)
is illustrated in Table 9.1. (The degeneracy is the number of configurations in the
same class). We assume that the resistance of a single bond is .R0 . The average
resistance, .〈R〉, of the renormalized cell is then the probability of the renormalized
cell to be occupied, .p' , multiplied with the resistance of the renormalized cell, .R ' ,
so that .p' R ' = 〈R〉. Using the scaling relation for the resistance, .R ∝ Lζ̃R , we can
determine the exponent from
ln R '
.ζ̃R = . (9.30)
ln b
where all the values are calculated for .p∗ , which we recall is .p∗ = 1/2 for this
scheme. The renormalization scheme and the values used are shown in Table 9.1,
where we use .p' = p∗ = 1/2 to calculate .R ' . The resulting value for the
renormalized resistance is
1
R' =
. g(c)P (c)R(c). (9.31)
p' c
5
1 1 5
= ' 1+1+4· +2·2+2·3+4·2+2·2 . (9.32)
p 2 3
≃ 1.917 . (9.33)
5 p3 (1 − p)2 2 3
6 p3 (1 − p)2 4 2
7 p2 (1 − p)3 2 2
152 9 Flow in Disordered Media
ln 1.917
. ζ̃R ≃ ≃ 0.939 . (9.34)
ln 2
This value is consistent with the scaling bounds set by the scaling relation in (9.24).
g ∝ (p − pc )μ ∝ ξ −μ/ν ,
. (9.36)
with the exponent .μ given as .μ = ν d − 2 + ζ̃R .
How can we use this scaling behavior as a basis for a finite-size scaling ansatz?
We extend the behavior of the infinite system to the finite size system by the
introduction of a finite size scaling function .f (L/ξ ):
L
g(ξ, L) = ξ −μ/ν f ( ) .
. (9.37)
ξ
We find the behavior of the scaling function, .f (u), by addressing the limiting
cases. When .ξ → ∞, we know that .g(ξ, L) will only depend on L, which means
that .f (L/ξ ) must cancel the .ξ −μ/ν term, that is, .f (u) ∝ u−μ/ν for .u ⪡ 1.
Similarly, when .ξ ⪡ L, we know that .g(ξ, L) will only depend on .ξ , which means
that .f (L/ξ ) must be a constant. The scaling function .f (u) therefore has the form
How does the scaling ansatz correspond to the observations? We can use the
program we have developed to measure the conductivity as a function of both
9.6 Finite Size Scaling 153
p and system size L. The following program has been modified for this type of
measurement:
import numpy as np
import [Link] as plt
from [Link] import measurements
from [Link] import ListedColormap
Lvals = [25,50,100,200,400]
pVals = [Link](np.log10(0.58), np.log10(0.85), 20)
C = [Link]((len(pVals),len(Lvals)),float)
P = [Link]((len(pVals),len(Lvals)),float)
nSamples = 600
mu = [Link](len(Lvals))
for iL in range(len(Lvals)):
L = Lvals[iL]
for pIndex in range(len(pVals)):
p = pVals[pIndex]
ncount = 0
for j in range(nSamples):
ncount = 0
perc = []
while (len(perc)==0):
ncount = ncount + 1
if (ncount > 1000):
print("Couldn’t make percolation cluster")
break
z=[Link](L,L)<p
lw,num = [Link](z)
perc_x = np.intersect1d(lw[0,:],lw[-1,:])
perc = perc_x[[Link](perc_x > 0)]
if len(perc) > 0:
zz = [Link]((lw == perc[0]))
# zz now contains the spanning cluster
zzz = zz.T
# # Generate bond lattice from this
g = sitetobond ( zzz )
# # Generate conductivity matrix
Pvec, c_eff = FIND_COND(g, L, L)
C[pIndex,iL] = C[pIndex,iL] + c_eff
C[pIndex,iL] = C[pIndex,iL]/nSamples
for iL in range(len(Lvals)):
L = Lvals[iL]
[Link](pVals,C[:,iL],label="L="+str(L))
[Link](r"$p$")
[Link](r"$g(p,L)$")
[Link]()
The results for .L = 25, 50, 100, 200, 400 are shown in Fig. 9.4. Here, we plot
both the raw data, .g(p, L), and the behavior of .g(pc , L) as a function of L on a
.log − log-scale, showing that .g(pc , L) indeed scales as a power-law with L.
154 9 Flow in Disordered Media
0.6 L = 25 − 1.2
L = 50
0.5 L = 100 − 1.4
L = 200
0.4 L − 1.6
log10 (g( pc , L ))
= 400
g( p, L )
0.3 − 1.8
0.2 − 2.0
0.1 − 2.2
− 2.4
0.0
0.60 0.65 0.70 0.75 0.80 0.85 1.4 1.6 1.8 2.0 2.2 2.4 2.6
p log10 L
160
140
120
100
L μ/ ν g( p, L )
80
60
40
20
0
0 5 10 15 20
L 1/ ν ( p − pc )
Fig. 9.5 Finite-size data scaling collapse for .g(p, L) showing the validity of the scaling ansatz
Scaling Data Collapse We can also test the scaling ansatz by plotting a finite-size
scaling data collapse. We expect that the conductivity will behave as
In Fig. 9.5 we demonstrate that this scaling form is valid by getting a data collapse
when we plot .Lμ/ν g(p, L) as a function of .L1/ν (p − pc ).
Estimating the Exponent .μ from the Data We can also use the results from the
simulations to measure .μ directly by plotting .g(p, L) as a function of .(p − pc )
9.7 Internal Distribution of Currents 155
1.3
− 0.25 L = 25
L = 50
− 0.50 L = 100 1.2
L = 200
− 0.75
L = 400
log10 g ( p, L )
1.1
− 1.00
μ
1.0
− 1.25
− 1.50 0.9
− 1.75
0.8
− 2.00
− 1.8 − 1.6 − 1.4 − 1.2 − 1.0 − 0.8 − 0.6 100 200 300 400
log10 |p − pc | L
Fig. 9.6 (a) Plot of .g(p, L) for increasing values of L. (b) Plot of the exponent .μ calculated by a
linear fit for increasing system sizes L
and fitting a linear function on a log-log plot. We do this for increasing values of L
in Fig. 9.6 (Notice that the curves for small values of L clearly are not linear, and
we should, ideally, have fitted the linear curve to only the part of the curve that is
approximately linear. We will not address methods to do this here, but you should
develop such methods in your own research.).
When we solve the flow problem on a percolation cluster, we find a set of currents
Ib = Ii,j for each bond .b = (i, j ) on the backbone. For all other bonds, the currents
.
will be identically zero. How can we describe the distribution of currents on the
backbone?
For electrical flow, the conservation of energy is formulated in the expression:
RI 2 =
. rb Ib2 , (9.41)
b
where R is the total resistance of the system, I is the total current, .rb is the resistivity
of bond b and .Ib is the current in bond b. We can therefore rewrite the total resistance
R as
Ib 2 2
R=
. rb ( ) = r b ib , (9.42)
I
b b
156 9 Flow in Disordered Media
where we have introduced the fractional current .ib = Ib /I . We have written the
total resistance as a sum of the square of the fractional currents in each of the bonds.
1 2q 1
〈i 2q 〉 =
. ib = i 2q n(i)di . (9.44)
MB MB
b
There is, unfortunately, no general way to simplify this relation, since we do not
know whether the function .n(i) has a simple scaling form.
What can we say about the scaling exponents .y(q) for moment q?
• For .q → ∞, the only terms that will be important in the sum are the terms where
.ib = 1, because all other terms will be zero. The bonds with .ib = 1 are the singly
connected bonds: all the current passes through these bonds. Therefore, we have
. (ib2 )∞ ∝ Ly(∞) ∝ MSC ∝ LDSC , (9.47)
b
• For .q = 1, we find from (9.42) that the sum is given as the total resistance of the
cluster
. (ib2 )1 = R ∝ Lζ̃R , (9.48)
b
2
b ib
Mq = 〈i 2q 〉 =
. ∝ Ly(q)−DB . (9.49)
MB
Because each term in the sum . b (ib )2q is monotonically decreasing in q, the sum
is also monotonically decreasing. We can therefore illustrate the curve .y(q) as in
Fig. 9.7, where we see that .y(q) is a non-trivial function of q. This is in contrast
to a unifractal distribution. We have seen unifractals in e.g. cluster number density.
The moments of the cluster number density has the form .Mq ∝ ξ xq with .x =
γ (β + 1)/ν. This means that all the moments are effectively described by a single
exponent, x. We call such distributions unifractal, whereas distributions where the
relationship is non-linear, such as for .y(q), we call multifractal.
In real resistor-networks, the case is even more complex, because the resistivity
is due to impurities, and the impurities diffuse. Therefore, the fluctuations in the
resistivity will also have a time-dependent part. This is the origin of thermal noise
in the circuit. If we keep the total current I constant, fluctuations in the resistivity
will lead to fluctuations in the voltage.
y(q)
DB
DSC
1 2 q
Fig. 9.7 Illustration of the exponents .y(q) characterizing the scaling of the moments of the
distribution of fractional currents, as a function q, the order of the moment
158 9 Flow in Disordered Media
1 with probability p
Gb = Gi,j =
. . (9.50)
0 with probability 1 − p
G2 with probability p
Gb = Gi,j =
. . (9.51)
G1 with probability 1 − p
∞ with probability p
Gb = Gi,j =
. . (9.52)
1 with probability 1 − p
In this case, we expect the conductance to diverge when p approaches .pc from
below, and that the conductance is infinite when .p > pc . It can be shown that the
behavior for the random superconductor network is similar to that of the random
resistor network, but that the exponent describing the divergence of the conductance
(and consequently conductivity) when p approaches .pc is s:
G ∝ (pc − p)−s ,
. (9.53)
Combining the Two Approaches How can we address both these problems? For
any system with a finite smallest conductance, .G< , we can always use the smaller
conductance as the unit for conductance, and write the functional form for the
conductance of the whole system as
G( G1 G2
G1 , G1 , p) G2
G(G1 , G2 , p) = (
. ) = G( , p) , (9.54)
G1 G1
(G1
G2 )
G = G2 (p − pc ) f± (
.
μ
), (9.55)
(p − pc )y
infinite conductance is canceled from the prefactor. That is, we need .f− (u) ∝ u.
We insert this into (9.55), getting
G1
G2
G ∝ G2 (p − pc )μ
. ∝ G1 |p − pc |μ+y . (9.56)
(p − pc )y
Finite .G2 and .G1 When .p → pc the conductance G should approach a constant
number when both .G2 and .G1 are finite. However, .p → pc corresponds to the
argument .x → +∞ in the function .f± (x). The only way to ensure that the total
conductance is finite, is to require that the two dependencies on .(p − pc ) cancel
exactly. We achieve this by selecting
f± (x) ∝ x μ/(μ+s) .
. (9.57)
G1
G2
G = G2 |p − pc | (
.
μ
)μ/(μ+s) , (9.58)
|p − pc |μ+s
which results in
G1 μ+s
μ
G = G2 (
. ) . (9.59)
G2
Exercises
PB (p) = P0 (p − pc )x ,
. (9.62)
and find an expression for the exponent x. You can assume that the fractal
dimension of the backbone, DB , is known.
(b) Assume that the functional form of PB (p) when p → pc+ and ξ ⪡ L is
PB (p) = P0 (p − pc )x ,
. (9.63)
Determine the exponent x by numerical experiment. If needed, you may use that
ν = 4/3.
Exercise 9.2 (Flow on Fractals) Use the example programs from the text to study
fluid flow in a percolation system.
(a) Run the example programs provided in the text to visualize the currents on the
spanning cluster.
(b) Modify the program to find the backbone and the dangling ends of the spanning
cluster.
(c) Use the program to find the singly connected bonds in the spanning cluster.
Exercise 9.4 (Current Distribution) Use the example programs from the text to
find the currents Ib in each bonds b on a spanning cluster at p = pc , p = 0.585,
and p = 0.60.
Exercises 161
Exercise 9.5 (Bivariate Porous Media) Rewrite the programs in the text to study
a bivariate distribution of conductances. That is, for each site, the conductance is 1
with probability p and g0 < 1 with probability 1 − p.
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Elastic Properties of Disordered Media
10
There are various physical properties that we may be interested in for a disordered
material. In the previous chapter, we studied flow problems in disordered materials
using the percolation system as a model disordered material. In this chapter we will
address mechanical properties of the disordered material.
We will address the behavior of the disordered material in the limit of fractal
scaling. In this limit we expect material properties such as Young’s modulus to
display a non-trivial dependence on system size. That is, we will expect material
properties such as Young’s modulus to have an explicit system size dependence. We
will use the terminology and techniques already developed to study percolation to
address the mechanical behavior of disordered systems such as the coefficients of
elasticity [4, 11, 20, 28, 40]
What are the elastic properties of a percolation system? First, we need to decide on
how to convert a percolation system into an elastic system. We will start by modeling
an elastic material as a bond lattice, where each bond represents a local elastic
element. The element will in general have resistance to stretching and bending.
Systems with only stretching stiffness are termed central force lattices. Here, we
will address systems with both stretching and bending stiffness.
Models for Stretching and Bending Stiffness We can formulate the effect of
bending and stretching through the elastic energy of the system. The energy will
have terms that depend on the elongation of bonds—these will be the terms that
are related to stretching resistance. In addition, there will be terms related to the
bending of bonds. Here we will introduce the bending terms through the angles
between bonds. For any two bonds connected to the same site, there will be an
j
ui i
uj
energy associated with changes in the angle of the bond. This can be expressed as
1 1
U=
. kij (ui − uj )2 + κij k φij2 k , (10.1)
2 2
ij ij k
where U is the total energy, the sums are over all particle pairs ij or all particle
triplets ij k. The force constant is .kij = k for bonds in contact and zero otherwise,
and .κij k = κ for triplets with a common vertice, and zero otherwise. The vector
.ui gives the displacement of node i from its equilibrium position. The various
Elastic Modulus Let us address the effective elastic behavior of the percolation
system using a material property such as Young’s modulus, E, or the shear modulus,
G. Let us consider a three-dimensional sample with cross-sectional area .A = L2 and
length L. Young’s modulus, E, relates the tensile stress, .σzz , applied normal to the
surface with area A to the elongation .ΔL in the z-direction.
Fz ΔLz
σzz =
. =E , (10.2)
A L
We can therefore write the relation between the force .Fz and the elongation .ΔLz as
EA EL2
Fz =
. ΔL = ΔL = Ld−2 EΔL . (10.3)
L L
We recognize this as a result similar to the relation between the conductance and
the conductivity of the sample, and we will call .K = Ld−2 E the compliance of the
system. We recognize this as being similar to the spring constant of a spring.
10.1 Rigidity Percolation 165
Elastic Properties When .p < pc What happens to the compliance of the system
as a function of p? When .p < pc there are no connecting paths from one side
to another, and the compliance will therefore be zero. It requires zero force .Fz
to generate an elongation .ΔLz in the system. Notice that we are only interested
in the infinitesimal effect of deformation. If we compress the sample, we will of
course eventually generate a contacting path, but we are only interested in the initial
response of the system.
where .τ is an exponent describing the elastic system. We will now use our
knowledge of percolation to show that this behavior is indeed expected, and to
determine the value of the exponent .τ .
Let us address the Young’s modulus .E(p, L) of a percolation system with occu-
pation probability p and a system size L. We could also write E as a function
of the correlation length .ξ = ξ(p), so that .E = E(ξ, L). Young’s modulus is
in general related to the compliance through .E(ξ, L) = K(ξ, L)Ld−2 . We can
therefore address the compliance of the system and then calculate Young’s modulus.
Dividing the System into Boxes of Size .ξ We will follow an approach similar to
what we used to derive the behavior of .P (p, L). First, we address the case when the
correlation length .ξ ⪡ L. In this case, we can subdivide the .Ld system into boxes
166 10 Elastic Properties of Disordered Media
of linear size .ξ as illustrated in Fig. 10.2. There will be .(L/ξ )d such boxes. On this
scale the system is homogeneous. Each box will have a compliance .K(ξ, ξ ), and
the total compliance will be .K(ξ, L).
K(ξ, ξ )
E(ξ, L) = L−(d−2) K(ξ, L) =
. . (10.6)
ξ d−2
In order to progress further we need to find the compliance .K(ξ, ξ ). This is the
compliance of the percolation system at .p = pc when the system size L is equal
to the correlation length .ξ . We are therefore left with the problem of finding the
compliance of the spanning cluster at .p = pc as a function of system size L.
10.1 Rigidity Percolation 167
Again, we expect from experience that the compliance will scale with the system
size with a dimension .ζ̃K :
K ∝ Lζ̃K .
. (10.7)
We will follow our now standard approach: We assume a scaling behavior, establish
a set of bounds for K, which will also serve as a proof of the scaling behavior of K,
and then use this result to develop a general theory for .K(p, L).
Energy, Force and Elongation of the System We will use arguments based on the
total energy of the system. The total energy of a system subjected to a force .F = Fz
resulting in an elongation .ΔL is:
1
U=
. K(ΔL)2 , (10.8)
2
where the elongation .ΔL is related to the force F through, .ΔL = F /K.
Consequently,
1 F 1 F2
U=
. K( )2 = . (10.9)
2 K 2 K
We can therefore relate the elastic energy of a system subjected to the force F
directly to the compliance of that system.
Upper Bound for the Compliance Our arguments will be based on the geomet-
rical picture we have of the spanning cluster when .p = pc . The cluster consists
of singly connected bonds, blobs, and dangling ends. The dangling ends do not
influence the elastic behavior, and can be ignored in our discussion. It is only the
backbone that contribute to the elastic properties of the spanning cluster. We can
find an upper bound for the compliance by considering the singly connected bonds.
The system consists of blobs and singly connected bonds in series. The compliance
must include the effect of all the singly connected bonds in series. However, adding
the blobs in series as well will only contribute to lowering the compliance. We will
therefore get an upper bound on the compliance, by assuming all the blobs to be
infinitely stiff, and therefore only include the effects of the singly connected bonds.
Let us therefore study the elastic energy in the singly connected bonds when the
cluster is subjected to a force F . The energy, U , can be decomposed in a stretching
part, .Us , and a bending part, .Ub : .U = Us + Ub .
168 10 Elastic Properties of Disordered Media
For a singly connected bond from site i to site j , the change in length, .δ𝓁ij , due
to the applied force F is .δ𝓁ij = F /k, where k is the force constant for a single bond.
The energy due to stretching, .Us , is therefore
1 1 F 2 1 MSC 2
. Us = kδ𝓁2ij = k( ) = F , (10.10)
2 2 k 2 k
ij ij
1 1 ri F 2 1
Ub =
. κ(δφij )2 = κ( ) = 2
MSC RSC F2 , (10.11)
2 2 κ 2κ
ij ij
where
1 2
2
RSC
. = ri , (10.12)
MSC
ij
where the sum is taken over all the singly connected bonds.
The elastic energy of the singly connected bonds is therefore:
1 R2
.USC = ( + SC )MSC F 2 , (10.13)
2k 2κ
and the compliance of the singly connected bonds is
F2 1
KSC =
. = . (10.14)
2U (1/k + RSC
2 /κ)M
SC
Lower Bound for the Compliance We can make a similar argument for a lower
bound for the compliance K of the system. The minimal path on the spanning cluster
provides the minimal compliance. The addition of any bonds in parallel will only
make the system stiffer, and therefore increase the compliance. We can determine
the compliance of the minimal path by calculating the elastic energy of the minimal
path. We can make an identical argument to what we did above, but we need to
replace .MSC with the mass, .Mmin , of the minimal path, and the radius of gyration
2 2
.R
SC with the radius of gyration of the bonds on the minimal path .Rmin .
2
Kantor [20] has provided numerical evidence that both .Rmin and .RSC 2 are
proportional to .ξ . When we are studying the spanning cluster at .p = pc
2
10.1 Rigidity Percolation 169
this corresponds to .Rmin and .RSC being proportional to L. This shows that the
dominating term for the energy is the bending and not the stretching energy when p
is approaching .pc .
Kmin ≤ K ≤ KSC ,
. (10.15)
where we have found that when .L ⪢ 1, .Kmin ∝ L−(Dmin +2) and .KSC ∝
L−(DSC +2) . That is:
Because .K(L) is bounded by two power-laws in L (for all values of L), we have also
demonstrated that .K(L) also is a power-law in L with an exponent .ζ̃K satisfying the
relation
This scaling relation gives us .K(pc , L). We use this expression to find .K(ξ, ξ ), the
compliance of a system of size .ξ from (10.6):
K(ξ, ξ ) ξ ζ̃K
E(ξ, L) =
. ∝ ∝ ξ ζ̃K −(d−2) . (10.18)
ξ d−2 ξ d−2
Bounds on the Exponent .τ The resulting bounds on the scaling exponents are:
3.41 ≤ τ ≤ 3.77 ,
. (10.22)
170 10 Elastic Properties of Disordered Media
Similarity Between the Flow and the Elastic Problems We see that the bounds
are similar to the bounds we found for the exponent .ζ̃R . This similarity lead Sahimi
[29] and Roux [27] to conjecture that the elastic coefficient E and the conductivity
g is related through
E
. ∝ ξ −2 . (10.23)
g
τ = μ + 2ν = (d + ζ̃R )ν .
. (10.24)
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Diffusion in Disordered Media
11
A typical example of a random walk is the random motion of a small dust particle
due to random collisions with air molecules, a process called Browian motion.
Random walks are general processes that we often use as physical, theoretical or
conceptual models.
n
rn = r0 +
. ui , (11.1)
i=1
where .ui is step i. We will usually assume that the steps .ui are independent and
isotropically distributed.
Fig. 11.1 Plots of 10 random walks of size .n = 100 (left) and .n = 1000 (right)
import numpy as np
import [Link] as plt
n = 100
u = 2*[Link](n,2)-1
r = [Link](u,axis=0)
[Link](r[:,0],r[:,1])
We can develop a theory for the position .rn as a function of the number of steps n.
For simplicity, we start the walker at the origin, so that .r0 = 0. First, we see find the
average position after n steps:
n
n
〈rn 〉 = 〈
. un 〉 = 〈ui 〉 = 0 , (11.2)
i=1 i=1
where we have used that since .ui are isotropic, .〈ui 〉 = 0. This is not surprising, the
random walker has the same probability to walk in all directions and therefore does
not get anywhere on average.
However, from Fig. 11.1 we see that the extent of the path increases with the
number of steps n. We can characterize this using the same measures we used to
11.1 Diffusion and Random Walks in Homogeneous Media 173
describe the geometry of the percolation clusters, by measuring .rn2 . We find the
average value of .rn2 using
〈rn2 〉 = 〈rn · rn 〉
= 〈( ui ) · ( uj )〉
i j
=〈 ui · uj 〉
i j
. (11.3)
=〈 ui · uj 〉 + 〈 ui · ui 〉
i=j i/=j
= 〈ui · ui 〉 + 〈u · u 〉
i i
i i/=j =0
= nδ , 2
The Dimension of the Random Walk Here we have demonstrated that the size of
the random walk, measured as .r 2 , is proportional to the number of elements in the
random walk. This is similar to the way we measured the size of a cluster using the
radius of gyration of the cluster. Indeed, it can be shown that these two definitions
give the same relation .rn2 = b2 n, where b is a constant of unit length that describes
the distribution of a single step. We realize that n is the number of elements in the
random walk, corresponding to s, the number of sites in a cluster. We have therefore
found that .rn = bn1/2 , or similarly, that .n = (rn /b)2 ∝ r Dw . This implies that
the dimension, .Dw , of the random walk always is .Dw = 2, independent of the
embedding dimension d. This means that for .d = 1 the random walk will overfill
space. Indeed, we expect it to step on top of itself repeatedly. For .d = 2 the random
walk will just fill space since .Dw = d, whereas for .d = 3 and higher dimensions the
random walk will fill a diminishing portion of space. Just like the spanning cluster
174 11 Diffusion in Disordered Media
had a smaller scaling exponent than the spatial dimension, and hence the density of
the spanning cluster decreased for larger systems.
We can also describe the motion of the random walker through the probability
density .P (r, t), where .P (r, t) dr dt is the probability for the random walker to be in
the volume .r dr in the time period t to .t + dt.
For a random walker on a grid, the probability to be at a grid position i is given
as .Pi (t). The probability for the walker to be at a position i at the time .t = t + δt is
then
.Pi (t + δt) = Pi (t) + [σj,i Pj (t) − σi,j Pi (t)]δt , (11.5)
j
where the sum is over all neighbors j of the site i. The term .σi,j is the transition
probability from site i to site j . The first term in the sum represents the probability
that the walker during the time period .δt walks into site i from site j , and the second
term represents the probability that the walker during the time period .δt walks from
site i to one of the neighboring sites j .
When .δt → 0 this equation approaches a differential equation
∂Pi
. = [σj,i Pj (t) − σi,j Pi (t)] . (11.6)
∂t
j
If we assume that the transition probability is equal for all the neighbors, so that
σi,j = 1/Z, where Z is the number of neighbors, the differential equation simplifies
.
to
∂P
. = D∇ 2 P , (11.7)
∂t
which we recognize as the diffusion equation, where the diffusion constant D is
related to the transition probabilities .σi,j and Z.
The general solution to this equation is
1 1 − 1 ( r )2
e−r /2Dt =
2
P (r, t) =
. e 2 |R| , (11.8)
(2π Dt)d/2 (2π ) |R|
d/2 2
√
where we have introduced .|R| = Dt.
It can be shown that the moments of this distribution are
〈r k 〉 = Ak R(t)k ∝ t k/2 ,
. (11.9)
11.2 Random Walks on Clusters 175
which displayes the same relationship as found above between extent, .〈r 2 〉, and
time, t, where the time is .t = nδt and .δt is the time a single step takes.
ix = [Link](L)
iy = [Link](L)
where L is the system size. If this site is empty, the walk stops immediately and its
length is zero:
if not cluster[ix,iy]:
return
176 11 Diffusion in Disordered Media
Storing the Trajectory of the Walker We store the trace of the walker
in two arrays (we need both to handle periodic boundary conditions later):
walker_map which consists of the positions ix,iy of the walker for each
step, and displacement, which consists of the positions relative to the initial
position of the walker.
Random Selection of Next Step How do we select where the walker can move?
The walker is restricted to move to nearest neighbor sites that are present. There are
several approaches:
• We may select a direction at random and try to move in this direction. If the
walker cannot move in this direction it stays put for this step, and then tries again
in the next step. In this case, the walker may have many steps without any motion.
• We may find all the directions the walker can possibly move in, and then select
one of these directions at random. In this case the walker will move onto a new
site in each step.
Both these methods effectively produce the same behavior. We will select the second
method. We therefore need to create a list of the possible directions to move in. In
order to make this list, we have a list called directions of possible movement
directions:
For each step, we need to collect all the possible steps into a list called
neighbor_arr. This is done by the following loop:
neighbor = 0
for idir in range([Link][1]):
dr = directions[:,idir]
iix = ix + dr[0]
iiy = iy + dr[1]
if 0<=iix<L and 0<=iiy< L and cluster[iix,iiy]:
neighbor_arr[neighbor] = idir
neighbor += 1
If this list is empty, that is, if neighbor is zero, there are no possible places to
move. This means that the walker has landed on a cluster of size .s = 1. In this case,
we stop and return with .n = 1.
11.2 Random Walks on Clusters 177
Finally, we select one of the neighbor directions at random, move the walker
into this site, update walker_map and displacement and repeat the process.
Preparing the Function We put this into a function and use the numba library to
speed up simulation times.
import numba
import numpy as np
@[Link](cache=True)
def percwalk(cluster, max_steps):
"""Function performing a random walk on the spanning cluster
Parameters
----------
cluster : [Link]
Boolean array with 1’s signifying a present site
max_steps : int
Maximum number of walker steps to perform.
Returns
-------
walker_map : [Link]
A coordinate map of walk, x in [0] and y in [1]
displacement : [Link]
A coordinate map relative pos., x in [0] and y in [1]
num_steps : int
Number of steps performed.
"""
walker_map = [Link]((2, max_steps))
displacement = np.zeros_like(walker_map)
directions = [Link]((2, 4), dtype=np.int64)
neighbor_arr = [Link](4, dtype=np.int64)
# X-dir: east and west, Y-dir: north and south.
directions[0, 0] = 1
directions[1, 0] = 0
directions[0, 1] = -1
directions[1, 1] = 0
directions[0, 2] = 0
directions[1, 2] = 1
directions[0, 3] = 0
178 11 Diffusion in Disordered Media
directions[1, 3] = -1
# Initial random position
Lx, Ly = [Link]
ix = [Link](Lx)
iy = [Link](Ly)
walker_map[0, 0] = ix
walker_map[1, 0] = iy
step = 0
if not cluster[ix, iy]: # Landed outside the cluster
return walker_map, displacement, step
while step < max_steps-1:
# Make list of possible moves
neighbor = 0
for idir in range([Link][1]):
dr = directions[:,idir]
iix = ix + dr[0]
iiy = iy + dr[1]
if 0<=iix<Lx and 0<=iiy<Ly and cluster[iix,iiy]:
neighbor_arr[neighbor] = idir
neighbor += 1
if neighbor == 0: # No way out, return
return walker_map, displacement, step
# Select random direction from 0 to neighbor-1
randdir = [Link](neighbor)
dir = neighbor_arr[randdir]
ix += directions[0, dir]
iy += directions[1, dir]
step += 1
walker_map[0, step] = ix
walker_map[1, step] = iy
displacement[:,step]=displacement[:,step-1]+\
directions[:,dir]
return walker_map, displacement, step
Testing the Function Let us test the newly generated function on a few simplified
cases. First, we try it on a system with .p = 1, that is, on a homogeneous system.
import numpy as np
import [Link] as plt
L = 50
p = 1
z = [Link](L,L)<p
[Link](z,origin="lower")
walker_map, displacement, steps = percwalk(z,200)
# walker_map is oriented as row-column (ix, iy)
[Link](walker_map[1,:steps],walker_map[0,:steps])
Walks from 10 such simulations are shown in Fig. 11.2. This looks reasonable
and nice, but we do notice that quite a few of these walks reach the boundaries of
the system. We may wonder how this finite system size affects the behavior and
statistics of the system.
11.2 Random Walks on Clusters 179
30
20
10
0
0 10 20 30 40
Measuring .r 2 (t) for a Random Walker The function percwalk returns the
displacements, .rn , for the walking starting from .r0 = 0. We find .rn2 and visualize
the result for a single walk:
L = 50
p = 1
z = [Link](L,L)<p
walker_map, displacement, steps = percwalk(z,200)
r2 = [Link](displacement**2,axis=0)
t = [Link](len(r2))
[Link](t,r2)
The resulting plot is shown in Fig. 11.3. We do not really learn much from this plot.
We need to collect more statistics. We need to generate many different walks and
then average over all the walks to find a statistically better measure for .r 2 (t).
350
40 300
250
30
200
r2
20 150
100
10
50
0 0
0 10 20 30 40
0 100 200 300 400 500
t
Fig. 11.3 (a) Trajectory of a random walk for a (homogeneous) system with .L = 50 and .p = 1.
(b) Plot of the corresponding .r 2 (t)
Notice a few details: If the number of steps is 1 or smaller it means that the walker
landed either on an empty size (.n = 0) or on a single site (.n = 1). We do not
want to include these in our statistics since they provide little information about
the behavior of the random walker. We use this program to collect statistics from
.M = 500 random walks of length .n = 10,000 steps on a .L = 100 system:
p = 1.0
L = 100
max_steps = 10000
num_walkers = 500
num_systems = 100
displacements = find_displacements(p,L,num_systems,\
num_walkers,max_steps)
dr1 = displacements[1:]
t = [Link](len(dr1))
[Link](t,dr1)
The resulting plot in Fig. 11.4 shows that the system indeed behaves as we expect—
for small values of t. However, as t increases, we see that the effect of the finite
system size L starts to affect the results. This is because the random walker is limited
by the wall and eventually we will be limited the .L × L system. This problem will
also arise when we study the percolation system. How can we reduce this problem?
log10 (r 2 )
− 0.5
− 1.0
− 1.5
− 2.0
− 2.5
that this solves the problem with the boundaries. This aspect will be even more
important when we study percolation systems in non-uniform media.
@[Link](cache=True)
def percwalk(cluster, max_steps):
"""Function performing a random walk on the spanning cluster
Parameters
----------
cluster : [Link]
Boolean array with 1’s signifying a site present
max_steps : int
Maximum number of walker steps to perform.
Returns
-------
walker_map : [Link]
A coordinate map of walk, x in [0] and y in [1]
displacement : [Link]
A coordinate map of relative pos, x in [0] and y in [1]
num_steps : int
Number of steps performed.
"""
walker_map = [Link]((2, max_steps))
displacement = np.zeros_like(walker_map)
directions = [Link]((2, 4), dtype=np.int64)
neighbor_arr = [Link](4, dtype=np.int64)
# X-dir: east and west, Y-dir: north and south.
directions[0, 0] = 1
directions[1, 0] = 0
directions[0, 1] = -1
directions[1, 1] = 0
directions[0, 2] = 0
directions[1, 2] = 1
directions[0, 3] = 0
182 11 Diffusion in Disordered Media
directions[1, 3] = -1
# Initial random position
Lx, Ly = [Link]
ix = [Link](Lx)
iy = [Link](Ly)
walker_map[0, 0] = ix
walker_map[1, 0] = iy
step = 0
# Check if we landed outside the spanning cluster
if not cluster[ix, iy]:
# Return the map with starting pos and nr of steps
return walker_map, displacement, step
while step < max_steps-1:
# Make list of possible moves
neighbor = 0
for idir in range([Link][1]):
dr = directions[:,idir]
iix = ix + dr[0]
iiy = iy + dr[1]
# Periodic BC
if iix>=Lx:
iix = iix-Lx
if iix<0:
iix = iix+Lx
if iiy>=Ly:
iiy = iiy-Ly
if iiy<0:
iiy = iiy+Ly
if cluster[iix, iiy]:
neighbor_arr[neighbor] = idir
neighbor += 1
if neighbor == 0: # No way out, return
return walker_map, displacement, step
# Select random direction from 0 to neighbor-1
randdir = [Link](neighbor)
dir = neighbor_arr[randdir]
ix += directions[0, dir]
iy += directions[1, dir]
step += 1
walker_map[0, step] = ix
walker_map[1, step] = iy
displacement[:,step]=displacement[:,step-1]+\
directions[:,dir]
return walker_map, displacement, step
We now have all the tools to start studying the behavior of a random walker on top
of a percolation system. We select .p = pc and drop the random walker on a random
position on the lattice. The resulting set of walks from such a simulation can be seen
in Fig. 11.5.
11.2 Random Walks on Clusters 183
80
60
40
20
0
0 20 40 60 80
− 1.5
log(r 2 )
− 2.0
− 2.5
− 3.0
0 1 2 3 4
log(n )
L = 100
p = 0.5927
z = [Link](L,L)<p
[Link](z,origin="lower")
for i in range(30):
walker_map, displacement, steps = percwalk(z,10000)
[Link](walker_map[1,:steps],walker_map[0,:steps],’o’)
Understanding Behavior for .p > pc We then simulate a larger set of walks for
p = 0.45, 0.50, 0.55, pc . The resulting plots of .r 2 (t) are shown in Fig. 11.6. We see
.
that when .p < pc , .r 2 (t) ∝ t x for some time, but then after some time, .r 2 (t) crosses
over to a constant instead. How can we understand this behavior?
184 11 Diffusion in Disordered Media
Long-Term Behavior When .p < pc For a single walker that lands on a cluster
of size s, we expect that the walker will be limited to walk on this cluster and
therefore cannot reach positions that are much further away than .Rs . Thus, after
a long time, we expect .r 2 (t) ∝ Rs2 . If we repeat this experiment many times, each
time dropping the walker onto a random occupied point in the system, we need to
take the average over all clusters of size s and over all starting positions, to find
the average of .r 2 (t) for all these different walks. If we drop the walker at a random
position, the probability for that walker to land on a cluster of size s is .sn(s, p),
and the contribution from this cluster to .r 2 (t) after a long time is .Rs2 . Therefore, the
average .〈r 2 (t)〉 for the walker is:
. 〈r 2 〉 ∝ Rs2 = Rs2 sn(s, p) . (11.11)
s
We approximate this sum by an integral and replace .n(s, p) by the scaling ansatz
n(s, p) = s −τ F (s/sξ ), getting
.
∞
. Rs2 = Rs2 ss −τ F (s/sξ )ds . (11.12)
1
We realize that the function .F (s/sξ ) falls to zero very rapidly when .s > sξ and it is
effectively constant below that, we therefore replace the integral with an integral up
to .sξ :
sξ
2
. Rs = Rs2 ss −τ ds . (11.13)
1
We now insert that .Rs2 ∝ s 2/D and perform the integral, getting:
2/D+2−τ 2/D
. Rs2 ∝ sξ ∝ sξ sξ2−τ . (11.14)
where we recognize the first factor as .ξ 2 ∝ (p − pc )−2ν and the second factor from
(4.33) as .(p − pc )β so that
. Rs2 ∝ (p − pc )β−2ν . (11.15)
We notice that in this case the average is of .Rs2 over .sn(s, p), but when we calculated
the correlation length in (5.18) the average was of .Rs2 over .s 2 n(s, p), and this is the
reason for the appearance of the exponent .β − 2ν and not simply .−2ν as we got for
the correlation length.
to be that .r 2 (t) ∝ t 2k for some exponent 2k. We notice that as p approaches .pc ,
the crossover time .t0 increases. All the curves for various p-values appear to have
similar, or possibly the same behavior for .t < t0 .
In Fig. 11.6 we notice that the exponent 2k is not 1, as we found for the
homogeneous case. It is clearly lower than 1. If we measure it, we find that
.2k ≃ 0.66 and .k ≃ 0.33. We call this behavior anomalous diffusion because the
mean squared distance .r 2 (t) does not grow linearly with time, but with an exponent
different than 1. What can we say about the crossover time .t0 ? We will return to this
after examining the case when .p > pc .
Diffusion at p = pc
From Fig. 11.6 we also see that for .p = pc the random walk follows .r 2 (t) ∝ t 2k .
This behavior is as expected. For times shorter than .t0 , the walker behaves as if
it is on .pc , whereas after a long time, .t > t0 , we start noticing that the walker
is restricted when it diffuses on the finite clusters. Another way to think of this
is that the crossover time .t0 increases as .p → pc , and diverges at .p = pc . The
exponent k is a universal exponent for diffusion on percolation systems. It does not
depend on the lattice structure or the rules for connectivity, but it does depend on
the embedding dimension d.
We can use the same computational approach to study the behavior of the random
walker when .p > pc . The resulting plots for .p = 0.8, 0.75, 0.70, 0.65 and .pc are
shown in Fig. 11.7. The plots show that when .p > pc , for short times the .r 2 (t) curve
follows the behavior for .p = pc with .r 2 (t) ∝ t 2k as shown by the solid line. But
for a crossover time .t0 , the behavior changes and crosses over to a behavior where
.r (t) ∝ t , that is, it crosses over to the behavior of a homogeneous system. How
2 1
− 1.0
− 1.5
− 2.0
− 2.5
− 3.0
〈r 2 〉 = D(p)t , r ⪢ ξ .
. (11.16)
D(p) ∝ g(p) ∝ (p − pc )μ .
. (11.17)
We therefore expect that when .p > pc , and the time is larger than a crossover time
t0 (p), that the behavior is scaling with exponent .μ, identical to that of conductivity.
.
And for a time shorter than the crossover time, the behavior is identical to the
behavior at .p = pc . We can understand this in the same way as above: When .t < tc
the walker does still not experience that the characteristic clusters are limited by a
finite characteristic length .ξ .
Scaling Theory
Let us develop a scaling theory for the behavior of .〈r 2 〉. We will assume that when
the time is smaller than a cross-over time, the behavior is according to a power-law
with exponent 2k, and that when the time is larger than the cross-over time, the
behavior is either that of diffusion with diffusion constant .D(p) for .p > pc , or it
reaches a constant plateau for the case when .p < pc .
Let us introduce a scaling ansatz with these properties:
. 〈r 2 〉 = t 2k f [(p − pc )t x ] . (11.18)
Notice that we could have started from any of the end-points, such as from the
assumption that
t
〈r 2 〉 = (pc − p)β−2ν G1 ( ) ,
. (11.19)
t0
or
t
〈r 2 〉 = (p − pc )μ G2 ( ) .
. (11.20)
t0
11.2 Random Walks on Clusters 187
We have two unknown exponents k and x that must be determined from independent
knowledge. We will assume that the function .f (u) has the behavior
⎧
⎨ const. for |u| ⪡ 1
.f (u) = uμ for u ⪢ 1 (11.21)
⎩
(−u)β−2ν for u ⪡ −1
Let us now address the various limits in order to determine the scaling exponents k
and x in terms of known exponents.
Scaling Behavior in the Limit .p > pc First, we know that when .p > pc , that is
when .u ⪢ 1, we have that
〈r 2 〉 ∝ (p − pc )μ t ,
. (11.22)
(p − pc )μ t ∝ t 2k f ((p − pc )t x ) ∝ t 2k [(p − pc )t x ]μ .
. (11.23)
2k = 1 − μx ,
. (11.24)
or
1 − μx
k=
. . (11.25)
2
Scaling Behavior in the Limit .p < pc Similarly, we know that the behavior in the
limit of .u ⪡ −1 should be proportional to .(pc − p)β−2ν . Consequently, the scaling
ansatz gives
2k + x(β − 2ν) = 0 .
. (11.27)
Solving to Find the Exponents We solve the two equations for x and k, finding
1 μ
k=
. [1 − ], (11.28)
2 2ν + μ − β
188 11 Diffusion in Disordered Media
log(t − 2k r 2 (t))
− 4.00
− 4.25
− 4.50
− 4.75
and
1
x=
. . (11.29)
2ν + μ − β
Our argument therefore shows that the scaling ansatz is indeed consistent with the
limiting behaviors we have already determined, and it allows us to make a prediction
for k and x.
Testing the Scaling Ansatz We can test the scaling function by a direct plot of the
simulated result. The scaling relation states that .r 2 (t) = t 2k f [(p − pc )t x ], which
means that .r 2 (t)t −2k = f [(p − pc )t x ]. If we therefore plot .r 2 (t)t −2k on one axis
and .(p − pc )t x on the other axis, all the data for the various values of p should
fall onto a common curve corresponding to the function .f (u). This is illustrated in
Fig. 11.8, which shows that the scaling ansatz is in good correspondence with the
data. Indeed, the plot also shows that the assumptions about the shape of the scaling
function .f (u) are correct.
We can write this relation in the same way as we wrote the behavior of an ordinary
random walk,
.t ∝ r dw , (11.31)
where .dw is the dimension of the random walk. We have therefore found that
1 μ
dw =
. =2+ , (11.32)
k ν − β2
11.2 Random Walks on Clusters 189
which is a number larger than 2. This means that for a given time, the walk remains
more compact, which is consistent with our intuition.
Defining the Cross-Over Time We have introduced a cross-over time, .t0 , which is
defined so that
(p − pc )t0x ≃ 1 ,
. (11.33)
which gives
t0 ∝ |p − pc |−1/x ∝ |p − pc |−(2ν+μ−β) .
. (11.34)
Interpreting the Crossover Time How can we interpret this relation? We could
decompose the relation to be:
|p − pc |β−2ν
t0 ∝
. , (11.35)
|p − pc |μ
where we know that the average radius of gyration for clusters are
[Rs2 ] ∝ |p − pc |β−2ν ,
. (11.36)
[Rs2 ]
t0 (p) ∝
. , (11.37)
D
where D is the diffusion constant. Why is this time not proportional to .ξ 2 /D,
the time it takes to diffuse a distance proportional to the correlation length? The
difference comes from the particular way we devised the experiment: the walker
was dropped onto a randomly selected occupied site.
Interpreting the Behavior for .p > pc Let us now address what happens when
p > pc . In this case, the contributions to the variance of the position has two main
.
terms: one term from the spanning cluster and one term from the finite clusters.
P '
[〈r 2 〉] = Dt =
. D t + Rs2 , (11.38)
p
where the first term, .P /pD ' t is the contribution from the random walker on the
infinite cluster. This term consists of the diffusion constant .D ' for a walker on the
spanning cluster, and the prefactor .P /p which comes from the probability for the
walker to land on the spanning cluster: For a random walker placed randomly on an
occupied site in the system, the probability for the walker to land on the spanning
cluster is .P /p, and the probability to land on any of the finite clusters is .1 − P /p.
190 11 Diffusion in Disordered Media
The second term is due to the finite cluster. This term reaches a constant value for
large times t. The only time dependence is therefore in the first term, and we can
write:
P '
Dt =
. Dt, (11.39)
p
Dp μ−β
D' =
. ∝ (p − pc )μ−β ∝ ξ − ν ∝ ξ −θ . (11.40)
P
where we have introduce the exponent
μ−β
θ=
. . (11.41)
ν
Interpreting the Crossover Time for .p > pc We have therefore found an
interpretation of the cross-over time .t0 , and, in particular for the appearance of the
.β in the exponent. We see that the cross-over time is
|p − pc |β−2ν ξ2
t0 ∝
. ∝ ' . (11.42)
|p − pc | μ D
The interpretation of .t0 is therefore that .t0 is the time the walker needs to travel a
distance .ξ when it is diffusing with diffusion constant .D ' on the spanning cluster.
How does the random walker behave on the spanning cluster? We have found that
for .p > pc and for .t > t0 the mean square displacement increases according to
〈r 2 〉 = D ' t ∝ (p − pc )μ−β t ,
. (11.43)
Interpretation of .t0 for Walks on the Spanning Cluster We expect the relations
to be valid up to the point .(t0 , ξ 2 ), where both descriptions should provide the same
result. Therefore we expect
'
ξ ∝ t02k ∝ D ' t0 ,
. (11.45)
11.2 Random Walks on Clusters 191
ξ2 (p − pc )−2ν
t0 ∝
. ∝ ∝ (p − pc )−(2ν+μ−β) . (11.46)
D' (p − pc )μ−β
Consequently, the value of .t0 is the same for diffusion on the spanning cluster
as for diffusion on any cluster including the spanning cluster. In general, we can
interpret .t0 as the time it takes for the walker to diffuse to the end of the cluster
when .p < pc , and the time it takes to diffuse to a distance .ξ on the spanning cluster
when .p > pc .
Interpretation of k for Walks on the Spanning Cluster Let us check the other
exponent, .k ' . We find that
'
ξ 2 ∝ (p − pc )−2(2ν+μ−β)k ,
. (11.47)
which is not the same as we found in (11.28) for all clusters. We find that .k ' is
slightly larger than k.
this gives
'
r 1/k ∝ t ∝ r dw ,
. (11.50)
where .dw can be interpreted as the dimension of the random walk. For the case of
random walkers on the spanning cluster at .p = pc we have therefore found that4
μ−β
dw = 2 +
. . (11.51)
ν
The fractal dimension is larger than 2. This corresponds to the walker getting stuck
on the percolation cluster, and the structure of the walk is therefore more dense or
compact.
192 11 Diffusion in Disordered Media
We can use the theory we have developed so far to address the behavior of the
diffusion constant with time. Fick’s law can generally be formulated as
. 〈r 2 〉 = Dt , (11.52)
or, equivalently, we can find the diffusion constant for Fick’s law from:
∂ 2
.D= 〈r 〉 . (11.53)
∂t
Now, we have established that for diffusion on the spanning cluster for .p = pc , the
diffusion is anomalous. That is, the relation between the square distance and time is
not linear, but a more complicated power-law relationship
'
. 〈r 2 〉 ∝ t 2k . (11.54)
As a result, we find that the diffusion constant .D ' for diffusion on the spanning
cluster defined through Fick’s law is
∂ 2k ' '
D' ∝
. t ∝ t 2k −1 . (11.55)
∂t
We can therefore interpret the process as a diffusion process where .D decays with
time.
In the anomalous regime, we find that
'
r ∝ tk ,
. (11.56)
We could therefore also say that the diffusion constant is decreasing with distance.
The reverse is also generally true: Whenever .D depends on the distance, we will
end up with anomalous diffusion.
Exercises 193
We can rewrite the dimension, .dw , of the walk to make the relation between the
random walker and the dimensionality of the space on which it is moving more
obvious:
μ β
dw = 2 − d +
. +d − , (11.59)
ν ν
where we recognize the first term as
μ
. ζ̃R = 2 − d + , (11.60)
ν
and the second term as the fractal dimension, D, of the spanning cluster:
β
D=d−
. . (11.61)
ν
We have therefore established the relation
dw = ζ̃R + D .
. (11.62)
This relation is actually generalizable, so that for a random walker restricted to only
walk on the backbone, the dimension of the walker is
dw,B = ζ̃R + DB .
. (11.63)
Exercises
Exercise 11.1 (Random Walks on the Spanning Cluster) In this exercise we will
use and modify the program percwalk from the text to study random walks in
percolation systems, and on the spanning cluster in particular. We want to find the
dimension dw of a two-dimensional random walk on the spanning cluster.
(a) Find the distance 〈r 2 〉 as a function of the number of steps N for random walks
on the spanning cluster for p = pc .
(b) Find the dimension, dw of the walk, from the relation 〈r 2 〉 ∝ N 2/dw .
(c) Find the distribution P (r, N) for the position r as a function of the number of
steps N for a random walker on the percolation cluster.
(d) (Advanced) Can you produce a data-collapse for the distribution P (r, N).
(e) (Advanced) Can you determine the functional form of the distribution P (r, N).
Is it a Gaussian?
(a) Find the distance 〈r 2 〉 as a function of the number of steps N for random walks
on the spanning cluster for p < pc and for p > pc .
(b) Plot log〈R 2 〉 as a function of N for various values of p.
(c) Can you find the behavior of the correlation length ξ from this plot?
(d) Discuss the behavior of the characteristic cross-over time t0 based on the plot.
(a) Find the distance 〈R 2 〉 as a function of the number of steps N for random walks
on the spanning cluster for p = pc .
(b) Find the dimension, dw of the walk, from the relation 〈R 2 〉 ∝ N 2/dw .
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
Dynamic Processes in Disordered Media
12
12.1 Introduction
So far, we have studied the behavior and properties of systems with disorder, such as
a model porous material in the form of a percolation system. That is, we have studied
properties that depend on the existing disorder of the material. In this chapter, we
will start to address dynamical processes that generate percolation-like disordered
structures, but where the structures evolve, develop, and change in time.
The first dynamic problem we will address is the formation diffusion fronts, and
we will demonstrate that the front of a system of diffusing particles can be described
as a percolation system.
The second dynamic problem we will address is the slow displacement of one
fluid by another in a porous medium. We will in particular demonstrate that the
invasion percolation process generates a fractal structure similar to the percolation
cluster by itself - it is a process that drives itself to a critical state, similar to the
notion of Self-Organized Criticality [3]. We will then address how we can study
similar processes in the gravity field, and, in particular, the influence of stabilizing
and destabilizing mechanisms.
The first dynamical problem we will address is the structure of a diffusion front [31]
on a square lattice. One example of such a process is the two-dimensional diffusion
of particles from a source at .x = 0 into the .x > 0 plane, when particles are not
allowed to overlap. The system of diffusing particles is illustrated in Fig. 12.1.
Exact Solution For this problem we know the exact solution for the concentration,
c(x, t), of particles, corresponding to the occupation probability .P (x, t). The
.
This solution produces the expected variance .〈x 2 〉 = Dt, where D is the diffusion
constant for the particles. There is no y (or z) dependence for the solution.
(a) (b)
200
50 210
100 220
150 230
200 240
250 250
300
260
350
270
400
280
450 xc 290
500
0 50 100 150 200 250 300 350 400 450 500 70 80 90 100 110 120 130 140 150
Fig. 12.1 Illustration of the diffusion front. Particles are diffusing from a source at the left
side. We address the front separating the particles connected to the source from the particles not
connected to the source. (a) The average distance is given by .xc shown in the figure. The width
of the front, .ξ(xc ), as .xc is also illustrated. The different clusters are colored to distinguish them
from each other. The close-up in figure (b) illustrates the finer details of the diffusion fronts and
the local cluster geometries
12.2 Diffusion Fronts 197
.P (xc , t) = pc , (12.3)
defines the center of the diffusion front: the front separating the particles that are
connected to the source
√ from the particles that are not connected to the source. We
notice that .xc (t) = Dt.
Width of the Diffusion Front What is the width of the diffusion front? For a given
time t, the occupation probability decreases with distance from .xc . The correlation
length will therefore depend on the distance .δx = x − xc to .xc . We expect that
a cluster may be connected to the front if it is within a distance .ξ of .xc . Particles
that are further away than the local correlation length, .ξ , will not be connected over
such distances, and will therefore not be connected. Particles that are closer to .xc
than .ξ will typically by connected through some connecting path. We will therefore
introduce .ξ as the width of the front, corresponding to the distance at which the
local correlation length, due to the occupation probability .P (x, t), is equal to the
distance from .xc . The local correlation length at a position x, .ξ(x), is given as
The distance w at which .ξ(xc + w) = w gives the width of the front. We can write
this self-consistency equation for w as
w = ξ0 |P (x + w, t) − pc |−ν .
. (12.5)
√
where we recognize that .xc ∝ Dt gives
dP 1 1
.
∝√ ∝ . (12.7)
dx xc Dt xc
198 12 Dynamic Processes in Disordered Media
which gives
ν/(1+ν)
w ∝ xc
. . (12.9)
The width of the front therefore scales with the average position of the front, and
the scaling exponent is related to the scaling exponent of the correlation length for
the percolation problem.
Time Development What happens in this system with time? Since .xc is increasing
with time, we see that the relative width decreases:
ν/(1+ν) 1
w xc − 1+ν
. ∝ ∝ xc . (12.10)
xc xc
This effect will also become apparent under renormalization. Applying a renor-
malization scheme with length b, will result in a change in the front width by a
factor .bν/(1+ν) , but along the y-direction the rescaling will simply be by a factor b.
Successive applications will therefore make the front narrower and narrower.
We will now study the slow injection of a non-wetting fluid into a porous medium
saturated with a wetting fluid. In the limit of infinitely slow injection, this process is
termed invasion percolation for reasons that will soon become obvious [12, 38].
Γ
Pc =
. , (12.11)
ϵ
where .Γ is the interfacial surface tension, and .ϵ is the characteristic size of the
pore-throats in the porous medium. However, there will be some disorder present in
the porous medium corresponding to local variation in the characteristic pore sizes
.ϵ. This will lead to a distribution of the capillary pressure thresholds, .Pc , needed
to invade a particular pore. We will assume that the medium can be described as
12.3 Invasion Percolation 199
Fig. 12.2 Illustration of the invasion percolation process in which a non-wetting fluid is slowly
displacing a wetting fluid. The left figure shows the interface in a pore throat: the pressure in the
invading fluid must exceed the pressure in the displaced fluid by an amount corresponding to the
capillary pressure .Pc = Γ /ϵ, where .Γ is the interfacial surface tension, and .ϵ is a characteristic
length for the pore throat. The right figure illustrates the invasion front after injection has started.
The fluid may invade any of the sites along the front indicated by small circles. The site with the
smallest capillary pressure threshold will be invaded first, changing the front and exposing new
boundary sites
a set of pores connected with pore throats with a uniform distribution of capillary
pressure thresholds, and we will assume that the capillary pressure thresholds are
not correlated but statistically independent. We can therefore rescale the pressure
scale, by subtracting the minimum pressure threshold and dividing by the range
of pressure thresholds, and describe the system as a matrix of critical pressures .Pi
required to invade a particular site.
Modeling the Fluid Displacement Process The fluid displacement process can be
modeled by assuming that all the sites on the left side of the matrix are in contact
with the invading fluid. The pressure in the invading fluid is increased slowly,
until the fluid invades the connected site with the lowest pressure threshold. This
generates a new set of invaded sites in contact with the inlet, and a new set of
neighboring sites. The invasion process continues until the invading fluid reaches the
opposite side. Further injection will then not produce any further fluid displacement,
the fluid will flow through the system through the open path generated.
We then start to gradually increase the pressure in the fluid and allow the fluid to
invade from the left side of the lattice. Let us say we have increased the pressure to
the value p (.0 ≤ p ≤ 1). This would mean that all sites that have .pi ≤ p and that
are connected to the left side would be invaded.
This corresponds to a percolation problem. If we make a percolation system with
occupation probability p, then the fluid will have invaded all the clusters that are
connected to the left side. Thus, we have mapped the invasion percolation problem
onto a percolation problem. Let us implement this approach.
First, we generate a random lattice of critical pressures and an array of pressures
p that we will loop through:
import numpy as np
import [Link] as plt
from [Link] import measurements
L = 400
z = [Link](L,L) # Random distribution of thresholds
p = [Link](0.0,0.7,0.01)
We step gradually through this set of p-values, finding the clusters of connected
sites that have p-values smaller than p[npstep]
Then, we find the labels of all the clusters on the left side of the lattice. All the
clusters with these labels are connected to the left side and are part of the invasion
percolation cluster called cluster. We do this in two steps. First, we find a list of
unique labels that are on the left side. Then we find all the clusters with labels that
are in this list using the numpy-function isin:
leftside = lw[:,0]
leftnonzero = leftside[[Link](leftside>0)]
uniqueleftside = [Link](leftnonzero)
cluster = [Link](lw,uniqueleftside)
Then we make a matrix that stores at what time t (pressure value .p(t)) a particular
site was invaded. This is done by simply adding a 1 for all set sites at time t to a
matrix pcluster. The first clusters invaded will then have the highest value in the
pcluster matrix. We use the pcluster matrix to visualize the dynamics.
Finally, we check if the fluid has reached the right-hand side by comparing the labels
on the left-hand side with those on the right-hand side. If any labels are the same,
12.3 Invasion Percolation 201
there is a cluster connecting the two sides (a spanning cluster), and the fluid invasion
process stops:
Results for Fluid Displacement Process The resulting pattern of injected nodes is
illustrated in Fig. 12.3, where the colors indicate the pressure at which the injection
took place. It can be seen from the figure that the injection occurs in bursts. When
a site is injected, many new connected neighbors are available as possible sites to
invade. As the pressure approaches the pressure needed to percolate to the other
side, these newly appearing sites of the front will typically also be invaded, and
invasion will occur in gradually larger regions. These bursts have been characterized
by Furuberg et al. [12], and it can be argued that the distribution of burst sizes as
well as the time between bursts are power-law distributed.
Fig. 12.3 Illustration of the invasion percolation cluster. The color-scale indicates normalized
pressure at which the site was invaded
actually be invaded, since the invasion process requires a connected path from the
inlet to the site for a site to be filled. We will therefore expect that the width of the
invasion percolation front corresponds to the correlation length .ξ = ξ0 (pc − p)−ν
as p approaches the percolation threshold .pc , because this is the length at which
clusters are connected. That is, clusters that are a distance .ξ from the left side will
typically be connected to the left side, and therefore connected, whereas clusters that
are further away than .ξ will typically not be connected and therefore not invaded.
This shows that the critical pressure will correspond to .pc . This also shows that
when the fluid reaches the opposite side, the system is exactly at .pc , and we expect
the invasion percolation cluster to have the same scaling properties as the spanning
cluster at .p = pc . There will be small differences, because the invasion percolation
cluster also contains smaller clusters connected to the left side, but we do not expect
these to change the scaling behavior of the cluster. That is, we expect the fractal
12.3 Invasion Percolation 203
dimension of the invasion percolation cluster to be D. This implies that the density
of the displaced fluid decreases with system size.
Invasion Percolation with and Without Trapping The process outlined above
does, however, not contain all the essential physics of the fluid displacement process.
For displacement of an incompressible fluid, a region that is fully bounded by the
invading fluid cannot be invaded, since the displaced fluid does not have any place to
go. Instead, we should study the process called invasion percolation with trapping. It
has been found that when trapping is included, the fractal dimension of the invasion
percolation cluster is slightly smaller [10]. In two dimensions, the dimension is .D ≃
1.82.
This difference between the process with and without trapping disappears for
three-dimensional geometries because trapping become unlikely in dimensions
higher than 2. Indeed, direct numerical modeling shows that the fractal dimension
for both the ordinary percolation system and invasion percolation is .D ≃ 2.5 for
invasion percolation with and without trapping.
Gravity Stabilization
as p approaches .pc both the front position and the front width diverges, that is, both
the front position .h̄ and the width, w, are proportional to the system size L:
h̄ ∝ w ∝ L ,
. (12.13)
However, when the system size increases, we would expect other stabilizing effects
to become important. For a very small, but finite fluid injection velocity, the viscous
pressure drop will eventually become important and comparable to the capillary
pressure. Also, for any deviation from a completely flat system for a system with a
slight different in densities, the effect of the hydrostatic pressure term will eventually
become important. We will now demonstrate how we may address the effect of such
a stabilizing (or destabilizing) effect [6, 25].
Invasion Percolation in a Gravity Field Let us assume that the invasion perco-
lation occurs in the gravity field. This implies that the pressure needed to invade a
204 12 Dynamic Processes in Disordered Media
pore depends both on the capillary pressure, and on a hydrostatic term. The pressure
Pic needed to invade site i at vertical position .xi in the gravity field is:
.
Γ
Pic =
. + Δρgxi , (12.14)
ϵ
We can again normalize the pressures, resulting in
Δρg '
.piC = pi0 + x , (12.15)
Γ ϵ2 i
where the coordinates are measured in units of the pore size, .ϵ, which is the unit of
length in our system. The last term is called the Bond number:
Δρg
Bo =
. , (12.16)
Γ ϵ2
Here, we will include the effect of the bond number in a single number g, so that
the critical pressure at site i is:
g = 0.001
grad = g*[Link](range(L),range(L))[0]
z = z + grad
leftnonzero = leftside[[Link](leftside>0)]
uniqueleftside = [Link](leftnonzero)
cluster = [Link](lw,uniqueleftside)
pcluster = pcluster + 1.0*cluster
# Check if it has reached the right hand side
span = np.intersect1d(lw[:,1],lw[:,-1]) # Test perc
if (len(span[[Link](span > 0)])>0):
break
[Link]([Link](pcluster),origin="lower")
Front Width Analysis This problem is similar to the diffusion front problem. For
an applied pressure p the front will typically be connected up to an average distance
.xc given as
. p = p 0 + xc g . (12.18)
The front will also extend beyond the average front position. The occupation
probability at a distance a from the front is .p' = pc − ag, since fewer sites will be
set beyond the front due to the stabilizing term g. A site at a distance a is connected
to the front if this distance a is shorter to or equal to the correlation length for the
occupation probability .p' at this distance. The maximum distance a for which a site
is connected back to the front therefore occurs when
This gives
and
a ∝ g −ν/(1+ν) ,
. (12.21)
We leave it as an exercise to find the form of the position .h(p, g), and the width,
w(p, g), as a function of p and g. We observe that the width has a reasonable
.
g=0.0 g=0.0001
g=0.001 g=0.01
Fig. 12.4 Illustration of the gravity stabilized invasion percolation cluster for .g = 0, .g = 10−4 ,
−3 −2
.g = 10 , and .g = 10 . The color-scale indicates the normalized pressure at which the site was
invaded
to address the behavior of the system when the characteristic length diverges is also
a powerful method of both experimental and theoretical use.
Gravity Destabilization
The gravity destabilized invasion percolation process corresponds to the case when
a less dense fluid is injected at the bottom of a denser fluid. This is similar to the
process known as secondary migration, where the produced oil is migrating up
through the sediments filled with denser water. Examples of the destabilizing front
is shown in Fig. 12.5.
12.3 Invasion Percolation 207
g=0.0 g=0.0001
g=0.001 g=0.01
Fig. 12.5 Illustration of the gravity de-stabilized invasion percolation cluster for .g = 0, .g =
−10−4 , .g = −10−3 , and .g = −10−2 . The color-scale indicates normalized pressure at which the
site was invaded
We can make a similar argument for the case when .g < 0, but in this case the
front is destabilized, and the correlation length .ξ ∝ (−g)−ν/(1+ν) corresponds to the
width of the finger extending front the front. The extending finger can be modeled
as a sequence of blobs of size .ξ extending from the flat surface. This implies that the
region responsible for the transport of oil in secondary migration is essentially a one-
dimensional structure: lines with a finite width w. The amount of hydrocarbons left
in the sediments during this process is therefore negligible. However, there will be
other effects, such as the finite viscosity and the rate of production compared to the
rate of flow, which will induce more than one [Link] destabilized invasion
percolation is used as a modeling tool in studies of petroleum plays and commercial
software packages are available for such simulation.
208 12 Dynamic Processes in Disordered Media
Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0
International License ([Link] which permits use, sharing,
adaptation, distribution and reproduction in any medium or format, as long as you give appropriate
credit to the original author(s) and the source, provide a link to the Creative Commons license and
indicate if changes were made.
The images or other third party material in this chapter are included in the chapter’s Creative
Commons license, unless indicated otherwise in a credit line to the material. If material is not
included in the chapter’s Creative Commons license and your intended use is not permitted by
statutory regulation or exceeds the permitted use, you will need to obtain permission directly from
the copyright holder.
References
21. H. Kesten, Percolation Theory for Mathematicians. Progress in Probability (Birkhäuser, Basel,
1982)
22. P.R. King, M. Masihi, Percolation Theory in Reservoir Engineering (World Scientific (Europe),
London, 2018)
23. S. Kirkpatrick, Percolation and conduction. Rev. Mod. Phys. 45(4), 574–588 (1973)
24. B.J. Last, D.J. Thouless, Percolation theory and electrical conductivity. Phys. Rev. Lett. 27(25),
1719–1721 (1971)
25. P. Meakin, A. Birovljev, V. Frette, J. Feder, T. Jøssang, Gradient stabilized and destabilized
invasion percolation. Phys. A Stat. Mech. Appl. 191(1), 227–239 (1992)
26. Y. Meir, R. Blumenfeld, A. Aharony, A.B. Harris, Series analysis of randomly diluted nonlinear
resistor networks. Phys. Rev. B 34(5), 3424–3428 (1986)
27. S. Roux, Relation between elastic and scalar transport exponent in percolation. J. Numer. Math.
19(6), L351–L356 (1986)
28. M. Sahimi, Relation between the critical exponent of elastic percolation networks and the
conductivity and geometrical exponents. J. Phys. C Solid State Phys. 19(4), L79–L83 (1986)
29. M. Sahimi, J.D. Goddard, Elastic percolation models for cohesive mechanical failure in
heterogeneous systems. Phys. Rev. B Condens. Matter Mater. Phys. 33(11), 7848–7851 (1986)
30. M. Sahini, M. Sahimi, Applications Of Percolation Theory (Wiley Press, Hoboken, 2003)
31. B. Sapoval, M. Rosso, J.F. Gouyet, The fractal nature of a diffusion front and the relation to
percolation. J. Phys. Lett. 46(4), 149–156 (1985)
32. W. Sierpinski, Sur une courbe dont tout point est un point de ramification. Comp. Rend. Acad.
Sci. Paris 160, 302–305 (1915)
33. S. Solomon, G. Weisbuch, L. de Arcangelis, N. Jan, D. Stauffer, Social percolation models.
Phys. A Stat. Mech. Appl. 277(1), 239–247 (2000)
34. H.E. Stanley, Cluster shapes at the percolation threshold: and effective cluster dimensionality
and its connection with critical-point exponents. J. Phys. A Math. General 10(11), L211–L220
(1977)
35. H.E. Stanley, Introduction to Phase Transitions and Critical Phenomena, Reprint edn. (Oxford
University Press USA, New York, 1987).
36. H.E. Stanley, A. Coniglio, Fractal structure of the incipient infinite cluster in percolation, in
Percolation Structures and Processes, ed. by G. Deutsher, R. Zallen, J. Adler (Adam Hilger,
Bristol, 1983)
37. D. Stauffer, A. Aharony, Introduction To Percolation Theory, 2nd edn. (Wiley Press, Hoboken,
1992)
38. D. Wilkinson, J.F. Willemsen, Invasion percolation: a new form of percolation theory. J. Numer.
Math. 16(14), 3365–3376 (1983)
39. K.G. Wilson, Renormalization group and critical phenomena. I. Renormalization group and the
kadanoff scaling picture. Phys. Rev. B Condens. Matter Mater. Phys. 4(9), 3174–3183 (1971)
40. J.G. Zabolitzky, D.J. Bergman, D. Stauffer, Precision calculation of elasticity for percolation.
J. Stat. Phys. 44(1), 211–223 (1986)
Index
Symbols r, 172
D, 68, 79, 80, 90 mu, 149
D(p), 185 ν, 77
DB , 122 n(s, p), 21
Dw , 173 φ, 1
DSAP , 121 σ , 42
DSC , 120 τ , 57, 165
Dmax , 121 zeta, 148
Dmin , 121 ζ̃K , 166
E, 164 ξ , 86
G, 136, 138 c, 1, 2
G(ξ, ξ ), 148 dw , 191
G(p, L), 138 g, 136
Gi , 139 g(c), 12
Gi,j , 139 g(r), 29
Ii , 139 g(r; p), 71
Ii,j , 139 g(s, t), 47
L, 30 gs,t , 39, 48
Lmin , 121 ib , 156
MSAP , 121 k, 184
Mmax , 121 k , 191
Mmin , 121 kij , 164
Ns , 20 n(ib ), 156
P , 17, 28, 36 n(s, p), 17, 19, 21, 28, 29, 40, 48,
P (p, L), 10, 57, 87, 89 54
Q(p), 34 p, 2, 3, 5, 36
R(p), 103 pc , 3, 4, 34, 59, 97
Rs , 63, 183 r 2 (t), 178
RSC , 148 rn , 172
Rblob , 148 s, 17, 20
S, 39, 56 sx , 52
S(p, L), 91 sξ , 23, 42
Z, 3, 33 sn(s, p), 17
Λ, 106 t, 47
Lambda, 107 t0 , 184
Π , 8, 12
Π (p, L), 94
β, 38 A
γ , 39, 57, 91 Anomalous diffusion, 185
κ, 164 Average cluster radius, 69
© The Author(s) 2024 211
A. Malthe-Sørenssen, Percolation Theory Using Python, Lecture Notes
in Physics 1029, [Link]
212 Index
G Multifractal, 157
γ , 28 Multiplicity, 12
Gravity destabilization, 203, 206
N
H Nearest neighbor, 3
Histogram, 51 Nearest neighbor connectivity, 3
Hyper-scaling, 90 Next-nearest neighbor, 3, 23
Next-nearest neighbor connectivity, 3
Non-wetting fluid, 198
I Normalization, 21
Incompressible fluid, 203 Normalization of n(s, p), 19
Inf-dimensional percolation ν, 29, 69, 76
average cluster size, 38 Number of neighbors, 3
β, 38
characteristic cluster size, 42
cluster number density, 40 O
density of spanning cluster, 34 Occupation probability, 3, 34
γ , 39 Occupied, 2
n(s, p), 40 Ohm’s law, 138
occupation probability, 34 One-dimensional percolation, 17
p, 34, 36 average cluster size, 25
pc , 39 characteristic cluster size, 23
percolation threshold, 36 characteristic length, 29
S, 38 cluster number density, 17
scaling ansatz, 45, 55 correlation function, 29, 72
scaling relations, 46 correlation length, 29
σ , 42 finite size scaling, 30
spanning cluster, 34 g(r), 29
sξ , 42 g(r; p), 72
Interfactial surface tension, 198 γ , 28
Invasion percolation, 198 ν, 29
invperc, 201 percolation threshold, 17
invpercgrad, 205 renormalization, 108
Iterative fractal, 126 Rs , 64
scaling ansatz, 24, 55
σ , 23
L spanning cluster, 28
Lacunarity, 128 sxi , 23
Logaritmic binning, 51
P
M perccorrfunc, 72
Mandelbrot-Given curve, 126 Percolating cluster, 7
Marging fixpoint, 107 Percolation probability, 8
Mass, 87 Percolation threshold, 3, 4, 59, 94, 97
Mass of spanning cluster, 81 percwalk, 177, 181
Matrix, 2 Periodic boundary conditions, 180
Maximum path, 121 Permeability, 137
Minimal path, 121 Π , 17
Moments, 91 Pore, 1
Moments of a distribution, 129 Porosity, 2, 5
Moments of cluster number density, 91 Porous media, 1
214 Index
S U
S, 27 Uncorrelated media, 1
Saturated porous medium, 198 Unifractal, 129, 157
Scaling ansatz, 42, 45, 54, 89 Universal exponent, 69
Scaling argument, 69 Universal number, 69
Scaling exponents Unstable fixpoint, 105
bounds, 149
Scaling function, 54
Scaling relations, 46, 56, 57 V
Self-affine, 198 Viscosity, 137
Self-avoiding path (SAP), 121
Self-similar, 101, 198
Self-similar scaling, 80 W
Shear modulus, 164 Width of diffusion front, 197
Shuffling, 6
Sierpinski gasket, 80, 128
σ , 23 X
Singly connected, 120 ξ , 29, 69, 71
Singly connected site, 119
Site, 2
Site lattice, 109 Y
sitetobond, 140 Young’s modulus, 164
The scaling exponent for connection paths in a percolation cluster, such as self-avoiding paths (SAPs), is determined by examining the mass of the paths relative to the system size. For a minimal path, its mass scales as M_min ∝ L^D_min, where D_min is the scaling exponent. Similarly, maximal path scaling involves M_max ∝ L^D_max. The average mass of all SAPs scales with an exponent D_SAP, which lies between D_min and D_max .
The significance of the scaling behavior of cluster number density, n(s, p), lies in its ability to describe the distribution of cluster sizes as the system approaches the percolation threshold, p_c. Its scaling form n(s, p) = s^{-τ} F(s/s_ξ), where s_ξ ∝ |p − p_c|^{-1/σ}, captures how cluster sizes diverge near p_c. This form helps in understanding how physical properties change with system size and critical point proximity, providing insights into universal behavior across different percolation systems .
The backbone of the spanning cluster serves as the critical path that connects opposite sides of the network, supporting transport properties like electrical conductivity. Despite the mass of the backbone becoming a smaller fraction of the total mass of the cluster as the system size increases, it remains essential for maintaining the network's integrity and connectivity .
The resistivity of a spanning cluster, R, is determined by the resistivity of its components connected in series: the singly connected bonds, R_SC, and the blobs, R_blob. The total resistivity is the sum of these components, reflecting how each part of the cluster contributes to the overall resistance. Therefore, for any system, R > R_SC, emphasizing the influence of blobs on resistivity .
SAPs are relevant as they define unique non-overlapping paths across a percolation cluster, crucial for understanding connectivity and transport properties. The minimal SAP provides the shortest path, while the maximal SAP gives insight into possible redundant connections. These paths help analyze the geometric properties and critical exponents of large-scale connectivity within percolation theory .
Programs are used to simulate percolation systems to study emergent phenomena and critical behavior near the percolation threshold. In these simulations, parameters like the site occupancy probability p, system size L, and the number of samples N are varied. This allows for analysis of percolation probability Π(p, L) and density of spanning clusters P(p, L) as functions of these parameters .
Studying small finite systems provides insights into basic percolation mechanisms, boundary effects, and deviations from infinite system predictions. For example, in an L=1 system, percolation occurs with probability p, while L=2 requires developing more advanced strategies due to multiple potential cluster configurations. These studies help refine models and underscore the importance of finite size scaling .
Minimal self-avoiding paths are the shortest paths connecting two sides of the cluster, characterized by M_min ∝ L^D_min, with smaller mass compared to maximal paths. Maximal paths are the longest connecting paths, scale as M_max ∝ L^D_max, and encompass more cluster sites. The exponents reflect how path lengths diverge with system size, influencing connectivity and transport properties .
Finite size scaling improves measurement by systematically varying system size to interpolate results to infinite system sizes. This method addresses limitations of merely using the largest system possible, allowing for more accurate estimations of properties like the percolation probability Π(p, L) and density of spanning clusters P(p, L), which is crucial as real systems cannot be infinite .
At criticality, the cluster number density is expressed as n(s, pc) = gs A^s (1 − pc)^2, indicating its dependence on universal scaling behaviors. The function reflects the distribution of cluster sizes at the critical point and is determined by scaling relations: n(s, p) = n(s, pc) e^{-sB(p−pc)^2}, where parameters like g_s and B are derived from the percolation model specifics. This provides insights into universality and scaling within percolation systems .