 

. u." I ‘6.
I: .tliix

.x. ‘
; .2. 5...
xi. a. um. i

 

 

 

 

r! a at. .
iw 14607.:
4343.? .u

Al.

my

 

55.4 C (1.11.4

TATE Ll

‘ mumlii‘ii'ﬁiiim m

3 1293 0182

 

 

W "

This is to certify that the
thesis entitled

Phase Transitions in a
Two-Component Site-Bond Percolation Model

presented by

Holger M. Harreis

has been accepted towards fulﬁllment
of the requirements for

M. S. degree in Physics

a' professor

Date August 20, 1999

 

0-7639 MS U is an Afﬁrmative Action/Equal Opportunity Institution

 

 

LIBRARY

Michigan State

University

 

 

PLACE IN RETURN BOX to remove this checkout from your record.
To AVOID FINE-S return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE

DATE DUE

DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

use WWW“

Phase Transitions in a Two-Component Site-Bond Percolation Model

By

Holger Martin Harreis

A THESIS

Submitted to
Michigan State University
in partial fulﬁllment of the requirements
for the degree of

MASTER OF SCIENCE

Department of Physics and Astronomy

1999

ABSTRACT

Phase Transitions in a Two-Component Site—Bond Percolation Model

By

Holger Martin Harreis

In standard percolation models either bond (randomly open edges) or site (ran-
domly occupied vertices) percolation is dealt with. Site-bond percolation combines
the two formulations, dealing with randomly occupied sites and randomly existing
bonds connecting these sites. However in this version of the model only one active
component exists, the other sites are considered unoccupied. A further generalization
is to consider several components, which was done for site percolation as well as bond

percolation previously and called polychromatic percolation.

In this work, a method to treat a N-component percolation model as effective one
component model is presented by introducing a scaled control variable p+. In Monte
Carlo simulations on 163, 323, 643 and 1283 simple cubic lattices the percolation
threshold in terms of p+ is determined for N = 2. Phase transitions are reported
in two limits for the bond existence probabilities p: and p¢. In the same limits,
empirical formulas for the percolation threshold pi as function of one component-
concentration, fb, are proposed. In the limit p: = O a new site percolation threshold,

ff 2 0.145, is reported.

ACKNOWLEDGEMENTS

First of all I wish to express my most grateful thanks to Professor Wolfgang Bauer
for giving me the opportunity to come to Michigan State University and work in his
research group at the NSCL as well as for his being my advisor during this wonderful

year.

I profoundly appreciated the very friendly atmosphere in the N SCL Nuclear Theory
Group in particular and at the NSCL in general. Furthermore I am indebted to Sen
Cheng, Jens Hofkens and Jens von Bergmann for numerous conversations. In particu-
lar it is a pleasure to acknowledge the fruitful discussions led with Jens von Bergmann

and the answers that Jens dekens provided to all computer related questions.

I appreciated the ﬁnancial support by the German National Scholarship Foundation,

the Studienstiftung des deutschen Volkes.

iii

Contents

LIST OF TABLES

LIST OF FIGURES

1 Introduction

2 The Foundations of Percolation Theory

2.1 Percolation Models ............................
2.1.1 Bond Percolation .........................
2.1.2 Site Percolation ..........................
2.1.3 Site—Bond Percolation ......................
2.1.4 Polychromatic Percolation ....................
2.1.5 Miscellaneous Percolation Models ................

2.2 Phase Transitions .............................
2.2.1 Deﬁnitions .............................
2.2.2 The Order Parameter .......................

2.2.3
2.2.4
2.2.5

The Correlation Length and Long Range Order ........
Critical Exponents and Universality ...............
The Thermodynamic Limit ...................

2.3 Monte Carlo Simulations .........................
2.4 Some Results ...............................

2.4.1
2.4.2
2.4.3
2.4.4
2.4.5
2.4.6

The Inﬁnite Cluster ........................
Cluster Size ............................
The Number of Open Clusters per Vertex ...........
Cluster Numbers .........................
Scaling Theory ..........................

Analogies of Percolative Phase Transitions to Thermal Phase
Transitions ............................

2.5 Finite Size Scaling ............................

vi

vii

(DOOQAOOOD

10
13
13
15
16
18
21
22
23
23
28

2.6

Percolation in Nuclear Physics ......................

Two-Component Site-Bond Percolation

3.1

3.2

3.3

3.4

3.5

3.6

3.7

3.8

General Case: N Components ......................
Two Components using the Natural Variables .............
3.2.1 General Considerations ......................
3.2.2 Simulation Results ........................
Simulation Techniques ..........................
3.3.1 General Features .........................
3.3.2 Cluster Analysis ..........................
3.3.3 Determination of the Critical Value ...............
Two Components using the Scaled Control Parameter .........
3.4.1 Entire Parameter Range .....................
The Percolation Probability .......................
3.5.1 19¢ —> 0 Limit ...........................
3.5.2 The ’Smearing-Out’ Effect ....................
3.5.3 19: —> 0 Limit ...........................
The Critical Value of the Scaled Control Parameter ..........
3.6.1 p¢ —> 0 Limit ...........................
3.6.2 p: —> 0 Limit ...........................
3.6.3 Concentration Dependence in the pat : 0 Limit ........
3.6.4 The Cut-Out Technique .....................
3.6.5 Concentration Dependence in the p: = 0 Limit ........
Finite Size Effects .............................
3.7.1 One-Dimensional Percolation Model ...............
3.7.2 Techniques in Higher Dimensions ................
3.7.3 A Proposal for a Novel Technique ................
Spin-Off Results ..............................
3.8.1 The Number of Open Clusters per Vertex ...........
3.8.2 The Number of Spanning Clusters per Vertex .........

3.8.3 Two-Dimensional Percolation ..................

Summary and Conclusion

Source Code
A.1 perc_prob-00.f ..............................
A.2 perc-ppintv0477.f ............................

47

50
50
52
53
54
56
56
58
59
61
62
69
69
77
78
82
82
84
85
90
92
96
96
97
101
103
103
104
106

109

112
112
122

List of Tables

2.1

2.2
2.3

Critical exponents in a ferromagnetic system (d is the dimension of the
system under observation) ........................ 19

Percolative functions and their analogous counterpart in thermal systems. 43

Percolative functions and the section in which they were ﬁrst intro-
duced or further discussed ......................... 43

vi

List of Figures

2.1

2.2

3.1

3.2

3.3

3.4

3.5

3.6

3.7

Example of a typical site percolation conﬁguration on a square lattice
in the subcritical phase. Black squares represent occupied vertices. . .

Example of a typical site percolation conﬁguration on a square lattice
in the subcritical phase. Black squares represent occupied vertices. . .

Probability to belong to the inﬁnite cluster, poo, in a simple cubic two—
component site-bond percolation lattice of size 1283 as a function of
the two parameters p: and p¢, calculated for a fraction of the blue
species of fb = 0.5. ............................

Probability to belong to the inﬁnite cluster, poo, in a simple cubic two-
component site-bond percolation lattice of size 1283 as a function of
the two parameters p: and p¢, calculated for a fraction of the blue
species of fb = 0.2. ............................

Probability to belong to the inﬁnite cluster, poo, in a simple cubic two-
component site-bond percolation lattice of size 1283 as a function of
the two parameters p: and 19,5, calculated for a fraction of the blue
species of fb = 0.0. ............................

Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter, p+, in a 1283 simple cubic lattice with fb =
0.5, for values of p:,p¢ 6 [0,1], ﬁtted with p00 or (p+ — pi)5 (solid
line). Here pi 2 0.251 for the ’upper’ branch, pi 2 0.280 for the
’lower’ branch and ﬂ = 0.41 in both cases. The inset shows the same
data in a double logarithmic representation. ..............

Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter, p+, in a 1283 simple cubic lattice with fb =
0.3, for values of p=,p¢ 6 [0,1], ﬁtted with poo o< (p+ — pi)ﬂ (solid
line). Here pi 2 0.251 for the ’upper’ branch, pi 2 0.218 for the
p¢ = 0 branch, pi 2 0.266 for the p: = 0 branch and ﬂ = 0.41 in all
cases .....................................

The same data as shown in Fig. 3.5 in a double logarithmic represen-
tation, with a SIOpe of all three curves of 6 = 0.41 ............

The normalized density of edges connecting two sites of equal flavor,
a:(fb) and the normalized density of edges connecting two sites of
different ﬂavor, a¢(fb). See also (3.14) and (3.15) ............

vii

54

55

55

62

66

67

68

3.8

3.9

3.10

3.11

3.12

3.13

3.14

3.15

3.16

3.17

3.18

3.19

Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter p+, in a 1283 simple cubic lattice with fb =
0.5, at ﬁxed p: = 1 for varied p¢. 50 independent simulation events
where taken into account and the average was taken over events in
which the number of br-bonds actually formed in the simulation, 17.1,,
was non-zero. ...............................

Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter p+, in a 1283 simple cubic lattice with fb =
0.5, at ﬁxed p: = 1 for varied p¢. 50 independent simulation events
where taken into account and the average was taken over all events. .

Percolation probability as a function of p¢, with a parametric depen-
dence on p+, shown for a 1283 simple cubic lattice with fb = 0.5. . . .

Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter p+, in a 1283 simple cubic lattice with fb =
0.5, at ﬁxed p¢ = 1 for varied pz. 50 independent simulation events
where taken into account and the average was taken over events in
which the number of br-bonds actually formed in the simulation, nbr,
was non-zero. ...............................

Critical value of the scaled control parameter, pi, in a 1283 simple
cubic lattice with fb = 0.5, plotted as a function of p¢. ........

Critical value of the scaled control parameter, pi, in a 1283 simple
cubic lattice with fb = 0.5, plotted as a function of p=. ........

Critical value of the scaled control parameter, pi, in a 1283 simple
cubic lattice with fb = 0.3, plotted as a function of p¢. ........

Critical value of the scaled control parameter, pi, plotted as a function
of the fraction of blue sites, fb, in the limit p9,; = 0 and obtained from
simulations on simple cubic lattices of size L = 163, L = 323, L = 643
and L = 1283 by a ﬁt to the scaling law Ipi(L) — pil oc L‘l/V. The

errors, estimated as described in [vdM97], are smaller than the symbol
sizes .....................................

Threshold p; versus concentration fb for two one-component site-bond
percolation models .............................

Threshold pi versus concentration fb for two one-component site-bond
percolation models. ............................

Percolation probability poo versus the probability that a site gets ’cut-
0116, pcut- .................................

Critical value of p75 and of the scaled control parameter, p; and pi,

plotted as a function of the fraction of blue sites fb in the limit p: = 0.
The data has been obtained by a ﬁt to |pi(L) — pil oc L'V". Again

the errors, estimated as described in [vdM97], are smaller than the
symbol sizes .................................

viii

78

82

83

85

86

87

89

90

91

93

3.20

3.21

3.22

3.23

3.24

3.25

Exact (solid line) and approximative (dashed line) effective threshold
in a one-dimensional percolation system as a function of the linear
lattice dimension L .............................

Finite size scaling for pi(L, p75 2 0) shown as a function of the linear
dimension of a simple cubic lattice at a concentration of the blue species
Of fb = 0.5. ................................

Probability d) that a spanning cluster exists, versus the lattice size L.
The percolation threshold pi is given by the common intersection point

1111mm: (pi) ..................................

The number of open clusters per vertex in a simple cubic 1283 lattice,
calculated for a concentration of the blue species fb = 0.5. Also shown
is the approximate theoretical expectation ................

Probability to belong to the inﬁnite cluster, poo(p+), and the probabil-
ity to belong to any percolating cluster, ppc(p+) in a simple cubic 1003
lattice at fb = 0.5 ..............................

Threshold p; versus concentration fb for two one-component site-bond
percolation models on the square lattice, compare Fig. 3.16. .....

ix

98

99

101

104

Chapter 1

Introduction

Critical phenomena are an everyday experience, not only in laboratories, but even
more in kitchens around the world do they have quite some importance in applied
experiments. Whoever has boiled an egg, whoever has made water boil (or freeze)
has let a phase transition happen. The theory of phase transitions has evolved to
a very important tool. A general feature of phase transitions is the fact that some
quality of the system under observation changes when a parameter which is freely
adjustable, is changed. The freezing of water or boiling of an egg are examples of
phase transitions in which the quality undergoing a transition changes abruptly as
a function of the varied external parameter. Here, this quality is the microscopic
structure, or, related to that, the density. Two distinct phases are present and are
separated, if one plots, for example, the density p versus the temperature T in a
phase diagram, by the phase transition line, on which, when ’crossed’, the system
goes from one phase to another. In the 1940’s, research on polymers and the process
of polymerization led to the question how the phase diagram for polymers might look
like and how one might ﬁnd a reasonably simple model for polymerization that yet
would be able to reproduce the qualitative characteristics of such a complex system
as a polymer solution. The process of polymerization is close to the ’hardening’ of
an egg while being boiled: Single molecules, or monomers, form larger and larger

1

macromolecules, large, but ﬁnite clusters of monomers, by means of an increasing
number of chemical bonds. Such a viscous solution of branched polymers is called sol
and it may result, at some point, in a very large network, of size on the order of the
system size, called a chemical gel (Physical gels in contrast are formed by means of
reversible bonds such as in the gelation of silica particles in water or NaCl solutions).
Examples of such sol-gel transitions in everyday life are pudding, gelatine or the milk-
to-cheese transformation. The ﬁrst model to describe the behavior of gelation was
essentially invented by Flory, [F1041], and Stockmayer, [St043], in 1941 and 1943. In
1957, Broadbent and Hammersley, [BH57], introduced the term ’percolation’ and a
more mathematical approach to the model. In subsequent papers, the percolation
model has been applied, in a broad diversity of techniques, to a very wide range of
problems, see [Sah94], such as dispersion and ﬂow in porous media, fracture and fault
patterns in rock, the morphological and transport properties of composite materials,
conductivity of (semi-) conductors, reaction kinetics, antigen-antibody reactions and

aggregations as well as the network formation on lymphocyte membranes.

In this work we will give, in Chapter 2, an introduction to some basic concepts
of percolation and statistical models and then present a multi-component site-bond
percolation model in Chapter 3. The results of Monte-Carlo simulations of this model
in the special case of two components will also be shown. In the same chapter, two
limits of the two component site-bond percolation model will be dealt with and their

characteristics discussed.

We shall furthermore abide by the following notation: Words deﬁned in the text
will be emphasized using italic font, words considered to be important for the partic-

ular section will be highlighted by means of bold font.

Chapter 2

The Foundations of Percolation
Theory

In percolation theory, in the most general case, one deals with a graph of dimension
d without a priori specifying a certain structure. In this discussion, however, we will
restrict ourselves to percolation theory on the graph associated with Zd and more

speciﬁcally, we will only discuss d = 3 in most cases.

2. 1 Percolation Models

Let us begin with yet another illustrating example. We imagine a very large porous
rock that we submerge in a pool. At its center we have installed a water sensor. What
is the probability that the water will ﬁnd its way to the center of the stone and wet
our sensor? To answer the question, we might try to model this random medium, the
stone, in the following way. Let p be a number that satisﬁes 0 g p g 1 and Z3 be the
simple cubic lattice. We then walk through the lattice and assign, at random, with
probability p, a bond to each edge (called open edge), whereas with probability (1 — p)
we leave the edge unoccupied (closed edge). After this we analyze the structure on
the lattice by deleting the closed edges, which leaves us with all open edges and the

set of their endpoints. All sites that are at least singly connected by bonds are said

to belong to one cluster. The question of the sensor-wetting is now related to the
existence of an inﬁnite cluster. Clearly, in a ﬁnite system, which we will always have
to work with in a simulation, it is not at all apparent what inﬁnite should mean,
we will come back to this question in the following sections, here it is sufﬁcient to
think of it as a cluster that spans the whole system. It is easy to imagine that with
increasing p also the size of the clusters increases. At some point pC however, we will
have an abrupt change in the structure: An inﬁnite cluster will exist in the lattice for
p > pc and the wetting of the center of the rock becomes possible. This value of p is
called the critical value or the percolation threshold. The model presented here has
bonds as the main ingredient and is thus called the bond percolation model. We will
introduce it in more depth in the following section, accompanied by an introduction
to other formulations of percolation models in the next sections. We will then present
results of percolation theory, which in most cases will generally hold for bond and

site percolation.

2.1 .1 Bond Percolation

We shall write Z = {. .. ,—1,0, 1,. . .} for the set of all integers and Z“ for the set
of all integral coordinate vectors 2 = (z1,22, . . . ,zd) where d is the dimension of the
system. The integral coordinate vectors 2: shall be called vertices. We call 151" the set
of all line segments < x, y >, (edges), between a pair of vertices x and y. Vertices
and edges together form the lattice, denoted by IL", which, in the case of d = 3 is the
simple cubic lattice. Let p and q satisfy 0 _<_ p S 1 and p+q = 1. We shall assign each
edge the condition open with probability p and the condition closed with probability
q. An equivalent formulation shall be to state that an edge is occupied by a bond
or unoccupied. Furthermore the assignments shall happen uncorrelated, each edge is

open or closed independently of all others. By Pp(X) we shall denote the probability

of some event X and we will write Ep(A) for the expectation value of an observable A.
A path of length m in the lattice ILd means a sequence zo, e0, 21, e1, . . . ,zm_1, em-1, zm
of vertices z,- and open edges e,- =< 23,-, 2,41 >. A circuit of length m + 1 shall mean
a path zo,e0, . .. ,zm satisfying em =< zm, 20 >. Path and Circuit are called open
if all of their edges are open, and closed if all of their edges are closed. A subgraph
of IL"1 is a set of edges ei and vertices 2:,— with e: 6 IE4 and z,- E Zd. Two subgraphs
are disjoint if they do not have any edge and vertex in common, they are edge-
disjoint if none of their edges is an element of both subgraphs. With the deﬁnition
of a subgraph we can ﬁnally put up the following deﬁnition: An open cluster or
just cluster is every subgraph that ﬁrst only consists of open edges and the vertices
connected to them and secondly is disjoint from all other such subgraphs. To specify
a particular cluster we will write C (2), indicating that it contains the vertex 2 and it
will mean the set of vertices connected to 2 through open paths. The size of a cluster
C (2) will be written as |C(z)| and is deﬁned as the number of vertices contained in
C (2) A cluster of size one would then be just one vertex without any open edge:
C (2) = {z}. If we furthermore introduce the notation z <—> y to indicate that there
exists an open path joining vertex z and vertex y, than a cluster may also be written
as C (2) = {y E Zd : z +——> y}. In the same spirit, z<~++ y means that there exists no
open path connecting z and y. Finally, the surface of a cluster C will be denoted by
BC and is deﬁned as the set of vertices in C that are adjacent to vertices not in C.
In the deﬁnitions made here, for example in deﬁning a cluster, we have assumed that
we only deal with nearest neighbor interaction. This is not a necessary condition,
percolation models with connections even longer than next-nearest neighbors have
been introduced, see 2.1.5. For a connection range that goes to inﬁnity, however, the

percolation threshold goes to zero.

2.1 .2 Site Percolation

Site percolation may be formulated in an analogous way to bond percolation. Again,
we consider Z3 and call the integral coordinate vector z E Zd a vertex. We also
add the set of edges IE“ to obtain the lattice IL“. As before, we need two random
variables, f and g which, again, shall satisfy 0 S f g 1 and f + g = 1. However
now, the objects of interest are not the edges, which are Openwith probability p
and closed with probability q in the bond percolation model, but rather the vertices
which in the site percolation model are assigned the quality open with probability
f and closed with probability g = (1 — f). Every edge is considered to be open,
closed vertices can be thought of as junctions that are blocked. With this notion we
can use the concept of a path and a circuit as introduced above. Also the deﬁnition
of a cluster as given in the previous section still holds. Phrthermore, no change is
necessary in the deﬁnitions of size and surface of a cluster. We show an example
of typical site percolation conﬁgurations on the square lattice in Figs. 2.1 and 2.2,
where occupied sites are represented by black squares, whereas unoccupied sites are
symbolized by white squares. An experiment that was actually carried out and for
which the shown site percolation problem on the square lattice can be seen as a model,
is to place conducting and non-conducting spheres in a bowl and check on percolation
as a function of the concentration of the conducting spheres. In Fig. 2.1 we show
the system with a relatively low concentration of conducting spheres and we can see
that the system is in a non-conducting state, there is no path of black squares that
leads through the entire lattice. The system is in a subcritical state. In Fig. 2.2
the system is shown at a higher concentration of black squares, that is at a higher
concentration of conducting spheres. As can be seen, a percolating path that spans
the entire lattice, exists. The system is found to be in a supercritical state. Obviously

some transition takes place in the lattice in between these two states. The question

6

 

 

 

 

 

 

Figure 2.1: Example of a typical site percolation conﬁguration on a square lattice in
the subcritical phase. Black squares represent occupied vertices.

is, how this transition happens in detail.

As already brieﬂy stated at the end of 2.1, many results hold for bond and site
percolation. A more strict way is to say that the bond percolation and the site
percolation model belong to the same universality class, [Sta79], meaning that they
are described by the same set of critical exponents, which will be introduced in Sec.
2.4. Still, although bond percolation was historically ﬁrst, site percolation is the more
general of the two models, since it is known that every bond percolation process on
a given lattice can be expressed in terms of a site model on a different lattice, see,
for example, [Fis61], [FE61] and [Ke582]. The converse statement, however, does not

hold generally.

 

Figure 2.2: Example of a typical site percolation conﬁguration on a square lattice in
the subcritical phase. Black squares represent occupied vertices.

2.1.3 Site-Bond Percolation

A very natural generalization of the two previous models is to simultaneously allow
open and closed edges as well as open and closed vertices. This means that no longer
only edges determine the ’passage’ of percolation as in the bond model, or only sites do
so, as in the site percolation model, but rather both the condition of edges and vertices
becomes important. We then have p,q 6 [0,1] with p + q = 1 and simultaneously
f, g 6 [0,1] satisfying f + g = 1. The concept of a cluster is still valid as in the two
models presented before. Both models from Secs. 2.1.1 and 2.1.2 represent certain
limits of the site-bond percolation model: Bond percolation is site-bond percolation
with f = 1 ﬁxed and site percolation is site-bond percolation in the limit p = 1.

The site—bond percolation model thus allows us to go continuously from bond to site

percolation and vice versa. As sited in [Sta79], some of the earlier work on site-bond

percolation can be found in [OCG+78, WP78, CSK79, ARRS79].

2.1.4 Polychromatic Percolation

If one does not think of each vertex to be associated with only two states, occupied
or unoccupied, one has, in the most general formulation, N states the vertex can be
in, of which one can think of as an occupation by different species. Zallen, [Zal77]
introduced such a generalization for bond and site percolation and also gave it the
name polychromatic percolation. He focused on the coexistence of percolating species
in highly connected lattices, giving a criterion for the occurrence of a panchromatic
regime in which all species percolate and showing percolation phase diagrams. In
such a generalized model one has several ’percolation systems’ on one lattice and if
the appropriate parameters are favorable, one can, in principle, have all ’systems’ per-
colating. One important parameter is the connectivity of the lattice. On the simple
cubic lattice, for example, with a coordination number Z = 6 (again only counting
nearest neighbors), percolation of both species is possible in a certain range of the
bond existence probability p or the fraction of sites f, but on the square lattice with
Z = 4, no such regime exists. Therefor, not only the structure, but also the dimen-
sion of the lattice is important in polychromatic percolation. Yet another version
of polychromatic percolation has been introduced by Halley and Holcomb, [HH78].
They considered a so called ’Three—Component Reactive Percolation Model’. The
essentials are as follows: Two types of atoms are considered which can furthermore
react and form dimer molecules. Each vertex in the lattice can be occupied by either
one type of the atoms or the molecule. The main focus lies on the resistivity, which

was shown to behave qualitatively the same as observed in experiments.

A further generalization of polychromatic percolation lies in applying it to the

9

site-bond percolation model, thus dealing with N different types of vertices and a
corresponding number of different edges, given by the number of combinations in
which neighboring vertices can occur. Work on this problem has been carried out in
the context of nuclear physics, see Sec. 2.6, [Baued], with two components, and in a
general approach for N components, where an approximate percolation criterion was
given, in[10395]. Further work has been done in the frame of this thesis and we will

deal with this problem in a later part of this work in further detail.

2.1.5 Miscellaneous Percolation Models

A variety of further different models have been proposed over the years, of which we

want to brieﬂy introduce some, without giving a detailed review of every model.

Long-Range Interaction

In our discussion a cluster has been deﬁned as the set of vertices and edges that are
connected by nearest-neighbor distances only. It is also possible to work with longer
range interactions. If one lets the coordination number z go to inﬁnity, however, the
percolation threshold goes to zero as 1/(z — 1). As stated in [Sta79], in work with
nearest and next-nearest neighbor interactions on bcc lattices, [HSB+79, QBH76],
as well as with higher order interactions, [HKM78], no differences from ordinary
percolation have been found, indicating that also taking longer range interactions

into account does not change the universality class of the model.

Percolation on the Cayley Tree

This model is also called percolation on the Bethe lattice. It presents a very special
case, which nonetheless has received some attention, as the percolation problem is

exactly solvable here. In the Bethe lattice, one starts with one point, the origin and

10

adds 2 neighbors to it, out of which z neighbors emanate again, one of which is the
connection to the origin, but 2 — 1 are new sites. Thus no closed loops are possible.
It turns out that in the high dimensionality limit, d ——> oo, percolation on hypercubic
lattices can be approximatively described by percolation on the Bethe lattice, since

in higher dimensions loops on hypercubic lattices become increasingly unimportant.

Trees on a Lattice

Analogously to the Bethe lattice, one can constrain percolation on other lattices to
be ’non-cyclic’, which causes a different behavior to emerge than the one observed in
ordinary percolation. As stated in [Sta79] no phase transition occurs in the case of

the square lattice.
Interacting Percolation

Until now, we have considered percolation to be an uncorrelated random process,
meaning that the state of each edge or vertex was completely independent of the state
of all of its neighbors. We may introduce some correlation by combining percolation
models with the Ising model of ferromagnets. One possible realization is, for example,
the so called Ising—correlated site-bond percolation discussed in [HS81]. In this model
the sites are occupied not randomly but correlated as the spins in an Ising model and
then bonds are formed with a certain bond existence probability, which gives three
free parameters: The site concentration f, the temperature T and the bond existence
probability p. Random site-bond percolation is still contained as T —> oo limit in this
model (And thus also ordinary bond and site percolation in the f ——> 1 and p —> 1

limits respectively).

11

Continuous Percolation

All discussed percolation models so far happened to be constrained to a lattice of
some structure (we focused on the lattices related to Z“), processes in nature however
take place in continuous systems. Therefore the question arises, if critical exponents
turn out to be different in percolation processes that take place in a random model
system with continuous variables. Numerical evidence seems to exist that, at least
in three dimensions, the lattice structure is not important for the critical exponents.

Also considerations from the renormalization approach are in favor of this statement.

Directed Percolation

In this version of percolation, one can exclude open edges in some directions, allowing
only bonds to exist in the other directions. This will make the percolation threshold

go up. A further discussion can be found in the literature, for example [Dua90].

Invasion Percolation

One can think of invasion percolation as a system to dynamically model the dis-
placement of one ﬂuid by another in a random medium as, for example porous rock.
The model uses a lattice where each site is assigned a number 1' E [0, 1] representing
the radius of the pore, which a vertex is considered to model. All edges are open
and represent the throats. According to the fact that capillary forces dominate, an
advancement of the invading ﬂuid to pores with small radius takes place. One thus
starts out with all sites occupied by the defending ﬂuid, the invading ﬂuid is injected
from one site of the lattice. In each step, the neighboring defending ﬂuid sites with
the smallest radius are occupied by the invading ﬂuid. Finite clusters of the defend-
ing ﬂuid that are completely disjoint from the rest of the ﬂuid can form. They give

rise to two different versions of the model: Either the defending ﬂuid is modeled as

12

incompressible (like oil, for example) and the invading ﬂuid can not advance in those
trapped ﬁnite clusters of the defending ﬂuid, or it is modeled as compressible and the
invading ﬂuid can further advance. In both cases the advancement stops as soon as
the invading cluster has reached the opposite site of the sample and percolates. One
then has the system at its critical point always. Systems exhibiting this feature are
called ’self-organized critical’, [BC89]. No ﬁnite clusters of the invading ﬂuid exist
and one does not have an equivalent of the occupation probability f. It turns out
that in the former case, that is with trapping, invasion percolation is different from
ordinary percolation, whereas for invasion percolation without trapping the results

are equivalent to ordinary percolation.

2.2 Phase Transitions

Before, in Sec. 2.4, we will arrive at a discussion of some important results of percola-
tion theory, we want to brieﬂy review a few aspects of phase transitions in general, as

well as, in Sec. 2.3, give a concise presentation of Monte-Carlo simulation techniques.

2.2.1 Deﬁnitions

A variety of physical phenomena can be seen as belonging to two categories. Systems
in the ﬁrst category are characterized in that their constituents do not interact with
each other or that at least the interaction is negligible. In this case a knowledge of the
energy level of all the constituents allows a direct calculation of the thermodynamic
functions of the entire system. Examples of phenomena belonging to this ﬁrst category
are the speciﬁc heat of gases and solids, where for the solids the interaction is a priori
not negligible but since for a very wide range of temperature the positions of the atoms
stay more or less ﬁxed in space, a transformation to normal coordinates is possible

so that the entity of atoms can essentially be treated as a system of (practically)

13

non-interacting harmonic oscillators. Another example is the spectral distribution
of black-body radiation. Systems belonging to the second category, however, do not
provide the basis for such a simple calculation of the thermodynamic functions, one
can not eliminate the interaction by means of any transformation. This interaction
leads, at some value of the temperature or another parameter, to a macroscopically
signiﬁcant form of cooperative behavior. As a result, the thermodynamic functions
of systems in this category exhibit analytic discontinuities or singularities which are
associated with the phenomenon of phase transitions. Prominent examples are the
melting of solids, the condensation of gases, the phenomena of ferromagnetism and
antiferromagnetism, order-disorder transitions as well as normal state to superﬂuid

and normal state to superconducting state transitions.

Fundamental is the notion of a phase, it means a possible state of a macroscopical
system in thermal equilibrium. In different phases, the same macroscopical observ-
ables of the same system can assume totally different values, depending on the values
of external parameters. Regions in the parameter space in which changes of these
parameters bring forth a change of the system from one phase into the other, are the

critical parameter ranges.

The earliest classiﬁcation of phase transitions was provided in 1933 by Ehrenfest,
in which the order of a phase transition is the means of classiﬁcation. An n-th or-
der phase transition is, following Ehrenfest, characterized in that the (n — 1)st ﬁrst
partial derivatives of the free enthalpy G after the natural variables of the system
(i.e. the temperature T and the pressure p in the ﬂuid system or the temperature
and magnetic ﬁeld in the case of the magnet) are continuous at the transition point,
whereas the nth derivative exhibits a discontinuity at this point. With increasing
order, however, the differences between two phases grow increasingly unimportant
and the question arises to which order a distinction between two phases is possible.

14

Furthermore it turns out that phase transitions who are not of ﬁrst order type are
characterized rather by singularities than by discontinuities as the Ehrenfest theory
would imply. There are other reasons, that go beyond the scope of this brief review,
that have led to the conclusion that the Ehrenfest scheme does not allow for a classiﬁ-
cation general enough. That is why the Ehrenfest classiﬁcation scheme has been more
or less abandoned, only two types of phase transitions remain: Discontinuous phase
transitions and continuous phase transitions, where the former correspond to the ﬁrst
order phase transitions in the Ehrenfest classiﬁcation and the latter essentially rep-
resent the second order phase transitions in a broader sense. The deﬁnition of the
discontinuous phase transition remains as for the ﬁrst order phase transitions. They
are characterized by discontinuities in the ﬁrst partial derivatives of some thermody-
namic potentials (i.e. discontinuities in the volume or in the spontaneous magnetic
moment). A continuous phase transition is present if there is non-analytical behavior
observed in the second partial derivatives at the critical point (e.g. experimentally in
the response functions like the heat capacity, the compressibility or the susceptibil-
ity). In most cases however, the notation remains according to the old classiﬁcation,

with ’ﬁrst order phase transition’ and ’second order phase transition’ used.

A phemenological approach for a uniﬁed description of all second order phase
transitions was given by Landau, we will just refer to the literature, for example

[Pat96].

2.2.2 The Order Parameter

We already cited the non-analytical behavior at a critical point as being characteristic
for a continuous phase transition. Aside this effect, a further typical characteristic
of a continuous phase transition is the order parameter. It signiﬁes a macroscopic

observable that only makes sense to deﬁne in one of the two phases in a phase tran-

15

sition. The name arises since phase transitions are often associated with a change
in the order of the system. In a multi-particle system there are always two ’com-
peting’ tendencies for minimizing the free energy F = U — T S: To minimize the
internal energy U, which can be achieved by a high order state, or to maximize the
entropy S, which goes along with a higher order state of the system. The effect of the
entropy however, scales linearly with temperature and thus the equilibrium will be
temperature dependent, which can lead to a phase transition, separating two phases
characterized by different order states of the system. Examples of order parameters

are I

o Ferromagnet: Spontaneous magnetization M, = m, / V.
> 0, if T < T 6

(m3: Magnetlc moment), With M, { = 0, if T > Tc

0 Gas—Liquid system: Density difference An = n; — n9.
For T < Tc there are two phases, gas and liquid with two different densities n f
and ng, whereas above TC there is just one homogeneous phase where it makes

no sense to deﬁne An.

2.2.3 The Correlation Length and Long Range Order

An important concept in the context of phase transitions is the correlation function

of an observable A, being deﬁned in the following way:

9(7‘3 7'“") = (a(7’)a(F')> — (0(0) (007%, (2-1)

where a0”) is the density of A:

A = / (13¢ am. (2.2)

The correlation function gives a measure of how much correlated the physical property

A is at f' and at r'. In a spatially homogeneous system we have that the correlation

16

function is invariant under translations and thus becomes

90“") = 9(IF- 7"'|)-

As an example, let us consider the density correlation

971057") = (”(7’)”(F'D - (”(77),

where n0?) means the particle density. In the special case of a homogeneous system
one usually encounters a damped, oscillatory behavior of g as a function of the distance

|F — F’ |, with vanishing correlations for |F — 77’] -+ 00:

gn(]r — F']) ————> 0.

|F—F’|—>oo

This happens, since for increasing distance the ﬁrst term in (2.1) factorizes:

(n(F)n(f")) ——> (E) 2.

IF—F’l—mo V

Generally, the behavior of the correlation function in the vicinity of a critical point

is conjectured to follow:

]F—F’]
-o 4! exP( £(T) )
g(T,T )O( W. (2.3)

This can be taken as a heuristic deﬁning equation for the correlation length 5 (T)
It obviously sets a length scale in the system, giving a measure for the range of the
correlation. A very important property of the correlation length is that it diverges at

the critical point:

This implies, that when approaching the critical point, the correlation length will

grow above the range of every interaction range in the system, spatial correlations

17

from microscopic constituents will assume macroscopic dimension. One speaks of
long-range order or also critical ﬂuctuations. The most important implication of this
behavior is that the physical properties of the system will, in the regime of critical
ﬂuctuations, no longer be determined by the microscopic form of the interaction or
structure, but rather by the value of the correlation length 6, so that different observ-
ables of entirely different systems can show remarkably qualitatively equal behavior
in the vicinity of their respective critical points. This phenomenon is known as uni-
versality. Furthermore the system is scale-invariant, it ’looks’ the same, no matter

on what level of magniﬁcation it is looked at.

2.2.4 Critical Exponents and Universality

Very often, one observes the following limiting behavior of some physical property
A in the vicinity of a critical point. A shows a power law dependence on a reduced
parameter 5. The reduced parameter 5 is obtained from the parameter 7' which leads
to the phase transition:

7 — ’rc
= , 2.5
e r. ( )

 

where rC denotes the critical value of the parameter r. 5 then determines the critical

behavior according to

A(5) 0: 5‘”. (215)

5—)0

This introduces a critical exponent a which determines the behavior in the direct
vicinity of the critical point. In the sense of the discussion in Sec. 2.2.3 and the
implications given by the divergence of the correlation length, a power law is to be
expected since at the critical point the system looses any characteristical length scale:
A power law does not introduce a characteristical length scale, as for example opposed

to an exponential behavior o< ex”, where a would set a length scale. As there are

18

some physical quantities, like, for example, the heat capacity in the Ising model, that
diverge logarithmically on approaching the critical point and not according to a power

law, a critical exponent is usually deﬁned in a more general way:

, ln|A(e)|
: l — .
w 5% Inc (2 7)
I _ . 1n|A(5)|
d) — gl/r‘rd lne ' (2.8)

In this deﬁnition the power law is contained as a special case. It is not a priori clear
that the same exponent 2p 2 1/1’ should hold for approaching the critical point from
both sides. In Table 2.1 we have compiled the critical exponents that usually appear
in a ferromagnetic system, see [N0198]. The Greek letters used here are convention.

For simplicity we have omitted the primed exponents. The exponent 7] seems, at ﬁrst

 

 

 

 

] Property Exponent Deﬁnition
Heat Capacity a CH oc 5‘“ e —> 0 H = 0
Order Parameter 6 M, or —53 e /‘ O H = 0
Susceptibility ’7 XT oc 5-7 e —> 0 H = 0
Critical Isotherm 6 M o< H1” 5 = 0 H —-) 0
Correlation Length 1/ 5 or 8‘” 5 —-> 0 H = 0
Correlation Function 77 g(r’, 7"") oc IF— F’I—d+2—" e = 0 H = 0

 

 

Table 2.1: Critical exponents in a ferromagnetic system (d is the dimension of the
system under observation)

sight, to have no real implication, since in (2.3) the behavior of g(r', 7'”) at the critical

point was already given. From there, 7) should, since 5 —0+ 00, simply be given by
5—)

17 = 3 — d. 77 however provides a means to quantify how much a real physical system

deviates from the simple conjecture in (2.3).

As already brieﬂy stated in the previous section, a most remarkable feature of
continuous phase transitions is the phenomenon of universality. Since in discontinuous

phase transitions the correlation length stays ﬁnite and does not diverge, critical

19

phenomena, along with universality are only observed in continuous phase transitions.
The concept of universality was introduced by Grifﬁths, [Gri70]. The universality
hypothesis states that the critical exponents are nearly universal, that is to say, they

are nearly the same for all thermodynamical systems. They only depend on

1. the dimension d of the system
2. the range of the interaction—potential in the system

3. the spindimensionality n.

For a classiﬁcation of the interaction-range the interaction potential is usually as-

sumed to follow a functional form as
V(r) oc r"(d+2+f),

where r is the inter-particle distance. For f > 0 one then speaks of a short-ranged
interaction, if f < (d/2) - 2, then the interaction is said to be long-ranged and
for (d/ 2 — 2) < f < 0 one speaks of medium-ranged interaction. In the ﬁrst case,
real universal behavior is expected, the interaction range is negligible compared to
the diverging correlation length. In the second case, the classical theories with a
special set of critical exponents become correct and in the third case on has to deal
with a rather complicated situation, where the critical exponents may depend on the
actual value of f. The spindimensionality gives the number of independent degrees
of freedom of the spins in the system. For example, n = 1 in the Ising model, n = 2

in the X Y model and n = 3 in the Heisenberg model.

The critical exponents are not independent, but rather show mutual dependencies,
which are manifested in the scaling laws. It turns out that there exist two independent

critical exponents, through which all others can be expressed. These are mostly results

20

of renormalization theory, to which we will brieﬂy come back at a later point in the

context of ﬁnite-size scaling in percolation theory.

2.2.5 The Thermodynamic Limit
Thermodynamic potentials of a macroscopic system are known to be extensive,
(oc N, OCV),

meaning that at constant temperature, the total energy of the system can be seen
as a sum of the energies of parts of the system. This however requires that the
interaction between particles in one part and particles in another be negligible, which

can rigorously be only true in the thermodynamic limit
N
N —> 00, V —> 00, such that V —> const. (2.9)

Also only in the thermodynamic limit do the predictions of the different statistical

descriptions, canonical, microcanonical and grandcanonical ensemble coincide.

Furthermore, and this is the most important result for a discussion of phase transi-
tions, it can be shown that in a rigorous sense, in a ﬁnite system, no phase transitions
are possible. As example, one can contemplate a classical N —particle system in the
volume V. Thermodynamic potentials involve the logarithm of the grandcanoni-
cal partition function A. Phase transitions should then be expected at the zeros
of A. The grandcanonical partition function A turns out to be a polynomial of or-
der N’ with one-particle canonical partition functions Z,- as parameters, which are
all positive deﬁnite (N ’ is the maximum number of particles to ﬁt in V). Thus for
N’ < oo <=$ V < 00 there can occur no phase transition in the system. A detailed

theory has been developed by Yang and Lee, see for example [Pat96, Nol98].

21

2.3 Monte Carlo Simulations

In investigating thermal systems one will almost never, aside from some special sit-
uations, encounter an exactly solvable problem. The Ising model, for example, is
exactly solvable, but only up to two dimensions and for some special lattices. In most
of the cases, thus, in order to calculate an observable, one will have to rely on some
other technique than direct calculation since the number of possible conﬁgurations
in a typical system forbids any practical analytical approach. Furthermore, such a
calculation would be rather inefficient, as a high number of conﬁgurations does essen-
tially not contribute to the expectation value of the observable since their probability
is negligible. That is, one will have to use numerical methods. In most cases then one
will not include the whole conﬁguration space S in the averaging, but just some ﬁnite
portion of it, depending on the precision one wants to achieve. The most desirable
way to sample the conﬁguration space, would be to include all those conﬁgurations
that contribute with a higher probability than others at the given temperature. This
procedure is called importance sampling. It is realized in an algorithm by deﬁning a
transition function P(d>|7r). It gives the conditional probability that one conﬁguration
(b E S is changed to some other conﬁguration 7r 6 S. This function has to fulﬁll the
following requirements that we only want to cite but not discuss in further detail, see
[Bin79, BH97] for a discussion. It shall be ergodic, provide for detailed balance and

reproduce the equilibrium probability distribution of the system.

A Monte-Carlo simulation essentially provides an implementation of such a tran-
sition function. For thermal systems one has to care about the last point from above,
that is one needs to have an equilibration process before a high number s of Monte-
Carlo Steps can be applied for the ’measurement’. The average is then taken over

the conﬁgurations that occurred in these 3 steps. In thermal systems, one possible

22

realization of the algorithm to implement the Monte-Carlo step is the Metropolis

Algorithm.

In the context of percolation theory however, we do not face most of the problems
that come up for thermal systems, namely there is no equilibrium conﬁguration,
there is no critical slowing down and also there is no need for such requirements as
ergodicity and detailed balance to be fulﬁlled. This implies that the implementation
of the Monte-Carlo step in percolation theory can be rather simple: We go through
the lattice, occupy each site with the given probability and do the same for the
edges. This essentially constitutes one Monte Carlo step and is repeated a multiple
number of times to produce a high number of independent lattice conﬁgurations. The
expectation value of the observable that is to be investigated is computed by taking

into account these conﬁgurations.

2.4 Some Results

In this section, after having reviewed some general features of phase transitions and
also having pointed out the differences one encounters in the simulation of a percola-
tive system as opposed to a thermal system, we we will present some major results of
percolation theory that are the basis for the further discussion of the present work.
We will mostly deal with ﬁndings that hold for bond and site percolation theory,

should the contrary be the case, we will give an explicit indication of this.

2.4.1 The Inﬁnite Cluster

We now combine the graphical concept and the probability part introduced in Secs.
2.1.1 and 2.1.2 and ask what happens in the lattice if the probability p for a vertex to

be occupied is changed. As already stated in the introductory example of Sec. 2.1.1,

23

it is intuitively clear that with an increase in the value of p the average cluster size
will increase, since a higher value of p generally leads to a higher number of Open
edges. In a simple picture then, at some critical value pc, a further small increase in p
will allow large clusters to agglomerate and form a cluster Cm“, that spans the entire
system, meaning that some vertex yb E Cm”. in the bottom layer and some vertex
zt E Cm” in the top layer are connected: yb <———> zb. Although in a ﬁnite system
there exists a ﬁnite probability that the spanning cluster is not the biggest cluster in
the system, this probability is very small and usually the percolating cluster can be
identiﬁed with the cluster of maximal size, Cm“. In the thermodynamic limit, where
the lattice size L ——> 00, this would yield ICmazl ——> 00. Moreover, in a ﬁnite system,
ﬁnite probabilities exist that more than one spanning cluster is present. In a ﬁnite

system, we deﬁne the inﬁnite cluster, Coo, by the following properties:
ICool = maXﬂCil}, (2-10)
and
vb, 2: E C... (2.11)

where yb and zt have to fulﬁll the following requirements in a lattice of linear dimension

L:
31:. = (231'.- . . ,0), (2.12)
z; = (k,l, . . .,L), (2.13)
if percolation is checked in the direction of the last coordinate.

The Percolation Probability

After having introduced the concept of an inﬁnite cluster, we can deﬁne one principal

quantity of interest: The percolation probability, also called strength of the inﬁnite
24

cluster, p0O (p) as the probability that a given vertex belongs to an inﬁnite cluster,

pooh?) = PMICI = 00)- (2.14)
This may also be written as
oo—l
10:1—23and=s, aim
3:1

since a given vertex either belongs to the inﬁnite cluster or to some ﬁnite cluster.
The notation ‘oo — 1’ as upper limit in the sum is meant to indicate that the inﬁnite
cluster is excluded in the summation. Pp can further be expressed as a sum over all
existing cluster conﬁgurations, which are also called lattice animals, in the following

way:

Pp(|C|— " 8) =;asmb 10'" q"- (2.16)

Here, asmb is the number of lattice animals which are characterized by their number

of vertices s, edges m and boundary edges b.

Intuitively it is clear that p00 (p) is a monotonous function in p and has the limiting
values poo(0) = 0 and 1000(1) = 1. Furthermore, according to the above considerations,
we expect a critical behavior of the percolation probability with a critical value, the

percolation threshold pc, yielding

_ :0: if p<pc
1a®f—{>O, ﬂ p>%. (2N)

The percolation threshold is formally deﬁned as
19.: — —Sup{p l 1900(1) 09—) — 0} (2-18)

Uniqueness of the Inﬁnite Cluster

As stated in the previous section, in a ﬁnite system, more than one percolating cluster
can exist. For an inﬁnite percolating system however, it can be proved that the inﬁnite

open cluster is unique, see [Gri99], where the following theorem is given:

25

Theorem 1 Iffor some p it is poo(p) > 0, then

Pp(3 exactly one inﬁnite open cluster) 2 1.

For the proof, we may refer to the literature.
Continuity of the Percolation Probability

One will notice that in (2.17) the case p = pC has been left out. As stated in Ref.
[Gri99], it is a known fact that for dimension d = 2 and for dimensions d 2 19
there exists no inﬁnite cluster right at the critical point: poo(pc) = 0. For dimensions
d 2 3 S 19 the situation is not clear. It is believed that poo(pc) = 0 for dimensions
at Z 3 also, but no rigorous results exist. On one hand this is, in a sense, required,
since, as will be discussed in more detail later, an important conjecture of percolation
theory is that poo(p) undergoes a second order phase transition at the percolation
threshold pc. It is thus required that p00 (p) be continuous. The only way to fulﬁll
this requirement on the whole interval [0,1] is to assign poo(pc) = 0. On the other
hand, the existence of an inﬁnite open cluster at the percolation threshold pc, called
incipient inﬁnite cluster is discussed in this context, see for example [SA94]. However
we will have to distinguish between our model, as an ideal, mathematical system that
can be inﬁnite in extent and the simulated physical system, which will always be
ﬁnite and subject to errors of all kinds. In such a system it will be impossible to
probe at exactly one value of the bond existence probability. Due to ﬂoating point
representation imprecision and the system discreteness, we will always end up with
an average taken over some small interval around the critical value pc. We will
thus deal with the subject in the following way: As long as we are considering ideal

characteristics of our system, we will consider the percolation probability to be zero

26

right at the percolation threshold:

1900(1).) = 0, (2.19)

which implies

pa = maXﬁD | 1900(1)) = 0}. (2-20)

In the context of physical systems however, we shall abide by the imprecise notation
of an inﬁnite open cluster at the threshold, meaning the incipient cluster if probed
in an interval around the critical value pc, with the width given by the limitations of

the numerical experiment.

The Critical Phenomenon

For percolation in one dimension, d = 1, there is no real critical phenomenon, as
in a one-dimensional inﬁnite system an inﬁnite chain of occupied vertices is needed,
which is only possible for a critical value p6:1 : 1, so that only one phase, the
non-percolating one is accessible on a set with non-zero measure. The existence of
a nontrivial critical phenomenon for dimensions d 2 2 was shown by Broadbent and

Hammersley, [BH57, Ham57, Ham59] as cited in [Gri89] with the following theoreml.
Theorem 2 Ifd 2 2 then 0 < pg < 1.

We can thus specify two phases of the process: A subcritical phase and a supercrit-
ical phase, where in the former we will almost ever only have ﬁnite open clusters,
whereas in the latter there will be almost always an inﬁnite open cluster present. The

probability that there exists an open inﬁnite cluster, 1b(p), follows

_ 0, if p<pc
w(p)—{ 1, ,, pm . (2.21)

 

1Again, we will not repeat the proof here, it can be found in the cited reference

27

This is only exactly true in the thermodynamic limit of an inﬁnite system: L —-> 00.
For ﬁnite L, ﬁnite-size-effects will ’smear out’ the step function to some differen-
tiable function. A more solid discussion of this effect is deferred to a later point, the

phenomenon of ﬁnite-size effects and ﬁnite-size scaling will be discussed in Sec. 2.5.

Site versus Bond Percolation

In the site percolation model, introduced in 2.1.2, the deﬁnition of the percolation
threshold is exactly the same as in the case of bond percolation, (2.18). However, if
we perform a simulation on exactly the same lattice, what will be the numerical value

of pf“, compared to p’c’m"? It has to be noted, that a bond percolation model on a

graph G may be reformulated as an equivalent site percolation model on a covering

graph Ge, see [Gri99], for which

szml = p2"‘€(G.). (222)
For both processes on the same graph, however, we have

p2°""(G) S 102“er (223)
On a high number of lattices, the strict inequality is true.

2.4.2 Cluster Size

Aside from the very important quantity, the percolation probability, another impor-

tant property of percolation is the mean cluster size

X(P) = Ep(|C|)- (2-24)

It can, just as the percolation probability, be expressed in terms of the cluster size

distribution, being the ﬁrst moment, whereas the percolation probability represented

28

the zeroth moment:

00

x<p> = :3 PpuCI = s), (225)

8:1

where again Pp(|C I = s) is the probability that a given vertex belongs to a cluster
of size 3. At this point a brief note concerning the difference between bond and site
percolation is necessary. Basically, in site percolation a calculation of the mean cluster
size is possible in exactly the same way: The question is, what the probability is, when
going through the lattice and hitting every vertex with the same probability, that this
arbitrarily chosen vertex is part of an s-cluster. Since in bond percolation every vertex
is occupied, a distinction of open and closed vertices is not necessary. This is different
in site percolation. When hitting every lattice site with equal probability, also closed
vertices are included, which introduces an implicit dependence of the mean cluster
size on the fraction of occupied sites, p. To avoid this inﬂuence, we have to introduce

a normalization, the fraction of open sites. We thus generalize (2.25) to

_ 2:13 Pp(]Cl = 5)

 

x(p) — so , (2.26)
28:1 Pp(|C])
where
°° 1 for bond percolation
Z Pp(|C|) = f . . (2.27)
3:, p or Site percolation

Also should we state that there is yet another way to deﬁne a mean cluster size. One
can take the average over all clusters, which would be equivalent to asking what the
mean cluster size is, if every cluster and not every site, as before, is hit with equal
probability. Given the actual number of clusters with size 3 to be NS, this amounts

to deﬁning the mean cluster size as

00 N3
S 2 "313%. (2.28)
3:1 3

We will subsequently stick to the deﬁnition of the mean cluster size as given in (2.26)

in the further discussion.

29

Although this is not a priori obvious, it is

x(p) < 00 if p < pc. (2.29)

For values of p above the percolation threshold the sum in (2.25) diverges because of

the existence of an inﬁnite cluster:

X(P) = 00 if p > pa. (2.30)

This can be remedied by introducing a truncated mean cluster size that does not take

into account the inﬁnite cluster:
X’(p) = Ep(|C|; ICI < 00)
oo—l
= Z s Pp(|C| = s). (2.31)
8:]

The symbolic notation ’oo — 1' means that the inﬁnite cluster is not included in the
sum. If we consider X! (p) as a continuation of x(p) on p 6 [pc, 1], then x(p) shows
a singularity at p : pC and is ﬁnite on the rest of the interval [0,1]. Whenever
poo(p) = 0, then x(p) = xi (p) naturally. The fact that both poo and X show a phase
transition at the same one point p.2 leads to calling percolation a model with a unique

phase transition.

2.4.3 The Number of Open Clusters per Vertex

For completeness, we will just give a short deﬁnition of this quantity and come back
to its postulated scaling behavior in Sec. 2.4.5. The Number of open clusters per

vertex is deﬁned in the following way:

Not—a

OO
:21); ,, (|C|— _ s) (2.32)

The ﬁrst terms can be calculated using elementary considerations for the proba-
bilities of clusters of a given size and then summing over the cluster size. The number

of open clusters per vertex is also referred to as f-function.

30

2.4.4 Cluster Numbers

In this section we will present a short discussion about cluster numbers in percolation
systems. At some points, however, when the possibility arises in a natural way, we
will also wander off to other quantities, always giving reference though to the section
where we ﬁrst introduced it or where we will give a more precise discussion. Also
do we have to clarify that by cluster numbers we do not mean the total number of
open clusters per vertex, but rather the number of open clusters of a given size 3 per
vertex, that is, the cluster size distribution. For the sake of simplicity we will mostly
argue on the basis of site percolation, which does not give a restriction on the general

results.

In the previous section we have already argued with the different moments of the
cluster size distribution Pp(|C| = s). It has to be noted, however, that this term is
somewhat imprecise. Pp(|C'| = 3) represents the probability that an arbitrary site
(every site in the lattice is hit with equal probability) belongs to a cluster with 3
sites. The real cluster numbers shall be denoted by n,, by which we mean the actual

number of clusters of size 3, N3, normalized by the number of sites in the lattice N:

n, = TV“ (2.33)

Then, the probability that a given vertex belongs to an s-cluster, is given by the

product of n, with the actual cluster size 3:
PP(]C] : 3) = 713(7)) 8° (2°34)

For a calculation of n, it is necessary to know the number of different conﬁgurations
clusters can be realized in with exactly 8 sites, which we denote by gst, compare (2.16).
In (2.16) we already introduced the corresponding quantity asmb for bond percolation.

Here, we have one index less, 3 being the number of vertices in the cluster and t being

31

the number of perimeter vertices, which is the number of all empty neighbor vertices,
including internal unoccupied, neighboring vertices. To get n, we then have to sum
over all of these conﬁgurations, multiplied with their respective statistical weight. For
a conﬁguration with 3 vertices and t perimeter sites, this is just p3(1 — p)‘, giving as

the number of size 3 clusters:

n. = 29.4230 — p)‘. (2.35)

One-dimensional exact Solution

In percolation we deal with statistically independent events. Thus, we have factorizing
probabilities always. Therefore, the probability to have 5 sites occupied is p‘. For
these occupied sites to be a cluster we furthermore need one unoccupied site at each
end of such a linear conﬁguration of occupied sites. The probability to have this
condition fulﬁlled is simply given by (1 — p)2. Neglecting chain end eﬂects, every site
has probability p3 (1 — p)2 of being, say, the right end of the linear cluster and with
L sites in the whole lattice the total number of 3 clusters is L p’ (1 — p)2, giving a

normalized number of s-clusters per vertex

1a®h=ﬁ(1-mf 93$

The percolation threshold in a one-dimensional system does not allow for a real phase
transition behavior: In a lattice of L sites, there will be
(1 — p) L ——> oo
L—ioo

empty sites in the linear chain, leading to a critical value of

2:1. (2M)

32

Deﬁning the probability for an occupied site chosen at random, to belong to an 3-

cluster as

nss

Zn, 3’

we can calculate the mean cluster size as, cf. (2.26)

X = 221283
= 2:57:23. (2.38)

 

’LUS:

 

This equation still represents a general result. The denominator is equal to p and
after evaluating the sum in the numerator, by turning it into derivatives of easier

sums, this gives

19
x(p) = ————1-9- for p < pc. (2.39)

With pC = 1, we thus encounter a divergence in the mean cluster size for p —+ pc,
which will also turn out to be the case in higher dimensions. In one dimension it is
also easy to calculate the correlation function and the correlation length, which we
introduced in the realm of the general discussion of phase transitions in (2.3). The
probability that a site at some position .v is in the same cluster as the site at the
origin requires the :1: connecting sites to be occupied, which happens with probability
p‘”. This is already the correlation function which we just rewrite in an exponential

form

"2’13

 

9(33) = 8244?), (2-40)
thus deﬁning the correlation length 5 as
t— - 1 <2 41)
ln(p)° '

For p 2 pc = 1 the logarithm can be expanded as ln(1 — 2:) 2 —a: for :1: << 1, which

gives

 

6: p. _ p. (2.42)

Thus also the correlation length diverges at the threshold with a critical exponent,

see Table 2.1,
1/ z —1. (2.43)

The Cayley Tree

The second exactly solvable lattice is the so called Bethe lattice or Cayley tree, which
we already introduced in Sec. 2.1.5. If we try to ﬁnd a way connecting the origin
to the surface, then we ﬁnd (z — 1) new sites on which we can continue at every
branching, p (z - 1) of which will be occupied on average. We need at least one site

occupied at every branching for a continuous path to exist and thus:

I
z—l'

 

pa = (2-44)

Since we can, without any argumentative change, replace site by bond in the above
considerations, (2.44) is valid for both site and bond percolation. If we now want to
calculate n,, we can use (2.35). It turns out that in the Bethe lattice the number
of perimeter sites is uniquely related to the number of sites 3 in the cluster, thus gst

reduces to 9,. The relation of perimeter sites t to 3 reads
t = 2 + s(z — 2), (2.45)

as a single site has 2 empty neighbors, two sites have (2 — 2) neighbors and every site
more adds another (z -— 2) empty neighbors. Without explicitly deriving 9, we can

then calculate the ratio ns (p) / n, (pc):

 

 

ns(p) (p)8(11:1:)2+s(z—2).

34

(2.46)

For the sake of simplicity we will set z = 3 in the further derivation and do the

following (exact) expansion:

f(p) : 53(3)

 

= id.) + f’(p.-)(p — p.) + great) — ps2

1—2pc 1 2
———— — ——P— c
pc(1"pc)( p)

where in the last equality we have made use again of the fact that we set 2 = 3 ﬁnally

arriving at

 

pc:

[\Dlr—t

Plugging this result back into (2.46) we have

 

 

7:1?) = (11:: )2(1 _ 4(p — pct): (2.47)
oc exp(—c s) (2.48)
where
c = — ln(1 — a(p — pc)2) (2.49)
: a(p — pc)2 + 0((17 — pc)3), (250)

which is true for 1 — a(p -— pc)2 << 1. Making use of the heuristic argument that the
mean cluster size X (second moment of the cluster size distribution) has to diverge at
the threshold pc, which, with (2.38), is not possible for an exponential decay of n,,

we conjecture the plausible form

”s(pc) (X S—Tv (251)

35

for large 3. This deﬁnes the so-called Fisher-exponent r. By calculation of X and

comparison with the general form for X (not derived here, see, e.g. [SA94]),

 

1
X 0C , (2.52)
p _ pc
we have
T = 5/2. (2.53)
Thus ﬁnally we arrive, with (2.46), at
n,(p) or 3-5/2 exp(—c s), (2.54)

valid for large 3.

To conclude this section, we want to point out why the Bethe lattice solution is
useful as an approximation of high-dimensional lattices. For ’normal’ lattices we have

for the relation of surface of some volume to the volume itself
81/ o< Vl-l/d. (2.55)

In the Bethe lattice with a branching factor 2, a cluster that has branched t times

will have

t—l

Nt 2 1+2 2(z-1)"

 

 

 

n=0
z(z — 1)‘ — 2
= 2.
z _ 2 ( 56)
sites in the cluster of which
St = z(z — 1)"1 (2.57)
sites are surface sites. This gives the ratio
33‘. __ z(z —1)‘“1(z — 2)
Nt _ z(z — 1)‘ — 2
’31 (Z _ 2’ (2.58)

I represents a ﬁnite limit. In normal lattices, this is only attained, as (2.55) shows, for
a dimension l/d = 0. Thus, percolation on a Cayley tree should approximatively de-
scribe percolation on non-lOOp-free lattices for higher dimensions d —) 00. In practice,

the results derived for the Bethe lattice hold for dimensions d > 6.

Conjectured Scaling Behavior of n,

For general lattices and dimensions, no exact solution for the cluster numbers is
known. Aside from the method of series expansion, where the cluster numbers for
clusters with a small number of vertices are calculated exactly and then extrapolation
methods for s —> 00 are employed, one will thus have to rely on conjectures and
checking them in numerical experiments. We might, with the derivations of the
previous section, try to postulate a general formula for the cluster numbers. It would

be very natural to assume

1-

exp(—c s). (2_59)

7150(3—

This is essentially the solution for the Cayley tree, but here, we will allow a general

exponent for c:
c oc |p — pCII/U for p —> pc. (2.60)

An obvious problem with this form is the fact that the one-dimensional solution is
not contained as special case, since, using exp(ln p) 2 exp(p — 1) for p 2 pc, the

one-dimensional result (2.36) can be cast into

n. = (p. — in)2 exp(-(pc — MS), (2-61)

which does not reproduce the power law in terms of s. A second problem arises with
respect to possible divergences in the cluster number. According to (2.35) the cluster

number is a ﬁnite polynomial in p, and hence n, and d;—::— for all n should stay ﬁnite.

37

However, with (2.60), this is not the case for non-integer values of 1/0. Thus, we take

1/0

z’ 2 cs o< |(p —pc)| s to the power ofa as new 2:

z = (p — pc) 3". (2.62)

In the postulate above, (2.47), the assumption for the ratio 72, (p) / 72, (pc) is an expo-

nential behavior,

"s (:0)
”3(196)

 

oc exp(z).

For a more general scaling assumption which should also reproduce the one-dimensional

result, the conjecture should just read

"30”) = f(z). (2.63)

 

We then have for the cluster numbers, with n,(pc) from (2.51),

ns(p) = s" f((p — pc)s") for p ——> pc, 3 —> oo. (2.64)

The exact form of f (2) has to be determined in simulations. A problem that occurred
for the ﬁrst scaling hypothesis (2.59) has now been remedied: With the assumption

that f (z) analytical everywhere, It, will show no singularities, as expected.

The one-dimensional case now is included with a = 1 and T = —2, compare (2.61):
n, = 3‘222 exp(z).
2.4.5 Scaling Theory

We already applied ideas of scaling theory in the previous section with the scaling
assumption (2.64) for the cluster numbers. Scaling theory makes predictions about
the behavior of macroscopic observables for p in the vicinity of the critical value pc.

The basis of scaling theory has not yet been proven in a strict mathematical sense,

38

[Gri99], but the conjectures have been veriﬁed to a very good extent in simulations.
We will, in the following, derive two scaling laws, but deprive us of a derivation of

the remaining scaling laws, but instead summarize the general idea of scaling theory.

Percolation Probability

With the new notation for cluster numbers introduced in the previous section, (2.34),
we have for the percolation probability in site percolation, in analogy to (2.15) for

bond percolation,

1300(1)) : p — :3 ”s(p)3- (265)

At the critical point, it is p00 2 0 and thus we can add and subtract pc in the following

form:

00—1 00—1

poo : (P _ pc) + Z "s(pc)5 — Z ”s(p)3

5:1 8:1

00—1

(1) — pa) + 2 (12.02.) - n. (19))8 (2-66)

3:1

where for the convergence of the ﬁrst sum, T > 2 is necessary. Although f (z) does
not have to follow an exponential decay according to the new scaling assumption
(2.63), it will in any case show a decay. This, together with the fact that for the
Bethe lattice we found an exponential decay, is a heuristic motivation that only large

cluster numbers essentially contribute in the sum, which allows for replacing the sum

39

by an integral,

p... = [d3 s“’(f(0)-f(2))

1 2—0—1’ _ Z
= /d2 W3 (f(O) f( l)

= up — paw—2V“ / d7 |zl““"‘2)/" (f(0) — m»

(X |(P — pc)lﬂ- (2°67)

This introduces the critical exponent ﬂ

 

s = . (2.68)

The limits of the integral are now case dependent: 0 to +00 for p > pC and —00 to
0 for p < pc, which is why the absolute values are introduced. We have furthermore

used

1
ds = ———sl“"dz , 2.69
a(p _ pa) ( )

s = z1/”(p—pC)—1/", (2.70)

in the derivation. For negative values of z, the factor of proportionality in the last

equation of (2.67) has to vanish, that is the integral should be zero:

f dz Izl‘H’ m0) — f(2)) = o.

This brings up a requirement to be fulﬁlled by f (z)

40

Mean Cluster Size

An analogous calculation for the second moment of the cluster size distribution X

yields another scaling law. It is, see (2.26),

X or Zszn,
0( [d3 sT—2f(z)

dz
T—3 a 3—7 0—1
IP—Pcl( V [7'3“ V f(Z)

oc lp — ml”, (2.71)

which introduces another critical exponent,

 

'y = a (2.72)

Equations like (2.68) and (2.72) are known as scaling laws and they put up rela-

tions between the different critical exponents for different macroscopic observables.
The Correlation Length

We can introduce a probability Tpf (y, 2) that for a given vertex y in a cluster C also
some other vertex 2. = (n, 0,0, . . . ), lying directly ’above’ y belongs to the same ﬁnite

open cluster. It is conjectured that Tpf follows

nQ—d-ﬂ, if p : pc

Tf Z 58 '
("3”) {em—map», if p¢o,p.,1 (273)

Here, (1 is the dimension, 77 is yet another critical exponent and 5 (p) is the correlation

length. For p near to pC the behavior of E (p) has the form

6(1)) % lp — pa”, (2-74)

41

a behavior which we already found with 1/ = 1 in the one dimensional example, see

(2.41). The notation ’z’ signiﬁes that 6 and |p — pcl satisfy

lim log—€01— = —1/. (2.75)
pxpc log lp - pal

Scaling Theory

In a brief summary, the basic ideas of scaling theory are as follows. Macroscopic
quantities behave as powers of |p — pcl in the vicinity of the critical point pc. Examples
are the correlation length E, the percolation probability pc,o and the mean cluster size
X- The correlation length gives a length scale for the system under observation, one
interpretation is to see it as the minimal length scale on which percolation at p can be
distinguished from percolation at pc. The correlation length diverges when the system
approaches the critical point. A further result of scaling theory are the scaling laws,
which the critical exponents are conjectured to satisfy and which reduce the set of
seven critical exponents to a set of two independent exponents. A further assumption
is that the critical exponents satisfy the hyperscaling relations, which involve the
critical exponents and the dimension d of the system. The hyperscaling relations are
only valid however for d S 6. .As already discussed for phase transitions in general,
the critical exponents are universal, depending only on the dimension of the system,
not the microscopic structure. A heuristic justiﬁcation takes as argument that at the
critical point, where the critical exponents determine the behavior of the macroscopic
observables, the correlation length diverges, which makes any information on the
microscopic structure of the system irrelevant. Related to that is the fact that at the
critical point, the inﬁnite cluster is fractal with a fractal dimension D = 91/48 2: 1.89
for site percolation on two-dimensional lattices, see [Fed88]. This so called incipient

inﬁnite cluster is self-similar on all length-scales.

42

2.4.6 Analogies of Percolative Phase Transitions to Thermal
Phase Transitions

Before, in the next section, we proceed to a discussion of Finite Size Scaling we want
to show the existing analogies between the quantities introduced above for percolation
theory and the thermal functions discussed in the realm of thermal phase transitions

in Sec. 2.2. We do so in the concise form of Table 2.2. In Table 2.3 we give a reference

 

 

 

Percolative Function Notation Thermal Function Notatio?!
Distance from critical point p — pc Distance from critical point T-Tc

F -Function a(p) Free Energy F
Percolation Probability poo (p) Spontaneous Magnetization M
Mean Cluster Size x(p) Susceptibility X
Correlation Length g CorrelatiOn Length 5

 

 

 

Table 2.2: Percolative functions and their analogous counterpart in thermal systems.

to the section in which we introduced or discussed the respective percolative function.

 

 

 

I Percolative Function Sec. ]
Number of Open Clusters per vertex 2.4.3
Percolation Probability 2.4.1
Mean Cluster Size 2.4.4, 2.4.5
Correlation length 2.2.3, 2.5, 2.4.5

 

 

 

Table 2.3: Percolative functions and the section in which they were ﬁrst introduced
or further discussed.

2.5 Finite Size Scaling

Let us consider a ﬁnite, d-dimensional site percolation system with occupation proba-

bility p. As soon as we have p > 0 in this ﬁnite system, there is a non-zero probability

43

that every site in the system will be occupied, given by
(pL)d = exp(Ld lnp). (2.76)

This would not be possible in an inﬁnite system and is an indication that the behavior
qualitatively changes in a ﬁnite system with respect to the inﬁnite one. In a ﬁnite
lattice the singularities associated with the critical behavior will be smoothed out.
In simulations, it is found that the conjectured power laws in |p— pcl still hold,
but with a shifted value of pc. How does the shifted value pc(L) approach the real
value p6 = pc(oo)? Let us consider the one-dimensional example again. We call the

probability that there exist a spanning cluster H. It is given by

11sz

= exp(—L/§), (2'77)

thus deﬁning the correlation length 5 as

1

{(29) = _E’ (2.78)

like in (2.41). In an inﬁnite system, we would expect H to show a step-function like
behavior, with H = 1 at p = pc 2 1 and H = 0 for p < pc. According to our above
consideration, there will be a smoothed behavior present for L < 00. We may then
deﬁne an effective threshold, at a p deﬁned by a more or less arbitrary value of II.

We choose to deﬁne the effective threshold in the following way:

_ 1
pear = H 1(2) (2.79)
With the expansion
1 1
—— z —, (2.80)
lnp 1 — p

44

valid for p z 1, we have for the effective threshold, as H é 1 / e, implying L = E,

19617: 1 - (2.81)

E.
Thus the effective threshold at a linear dimension L approaches the true threshold in

an inﬁnite system as
— = L'l/V 2 82
pc pet?" a ( ' )
with

1/ = 1. (2.83)

We now want to apply the same concept to higher dimensions. Although the
situation is not as easy as in one dimension, we conjecture an essentially similar, but
generalized behavior. With II again being the spanning probability, we postulate the

following scaling behavior for H in ﬁnite systems:
11(1), L) = 9((p — pc)L1/") for Ip - pal << 1, L >>1. (2-84)

Then the derivative of II with respect to p, which gives the probability that the system

starts percolating if p is increased by dp , reads

dH : diLl/v

315— dz (285)

Here now we take the effective threshold to be the average threshold at which, aver-

aged over p E [0, 1], the system starts percolating:

peﬂ' : pav
1 (ill
2 d —. 2.
f p p dp ( 86)
0

45

This can be rewritten in terms of z = (p — pc)L1/", yielding
1 d_g_
o

= [1) dz ———(zL 1/"+pc)

z(p= 0)
4/ z(p=1) dg 202:1) dg
= L " / dz zE+pC [dz 5. (2.87)
Z(10:0) z(p=0)
Since we have
2(sz dg 1 dn
/ dz 5 = fdp E
z(p=0) o
= H(1)—H(0)
= 1, (2.88)

we deduce that the effective threshold approaches the real threshold as

z(P=1)

d9
e _ c : L—l/u / d
1017 p z Zdz
z(p=0)
o< L_1/" for Ipeg— pcl << 1, L << 1. (2.89)

As this ﬁnite size scaling is only valid for large lattice sizes, in the application system-
atic errors will occur. In higher dimensions no estimate for these deviations exists, we
will come back to the question of the systematic errors at a later point, see Sec. 3.7.
In the deﬁnition of the effective threshold (2.86), it is not at all necessary to deﬁne

the effective threshold as the average threshold, as we have done. It would be equally

46

possible to choose the inﬂexion point of H, or, as we have done in the discussion for
one dimension, to set peg = II'1(p’), where p’ can be some value in the smoothed
transition region of II. The result for how this effective threshold approaches the real

one, (2.89), does still hold.

2.6 Percolation in Nuclear Physics

Percolation has been introduced in nuclear physics by Bauer to propose a model
of fragmentation reactions, see [BDMP85, BPDM86, BBDG87, Bau88, Baued]. We
will brieﬂy present the essential features of this model, called the nuclear lattice
model. Although this is, by no means, a restriction of the model, we just consider

fragmentation reactions of the type
p+AT=Ap+X. (2.90)

As underlying lattice the simple cubic lattice is chosen and sites on this lattice are
populated in an approximately spherical way to model the target nucleus. Deformed
nuclei could also be considered, but the results are not changed essentially. Thus the
lattice sites model the nucleons of the target nucleus and at the, initial state each
nucleon is connected by bonds to its six nearest neighbors. An impact parameter b is
chosen and according to the so-called ﬁreball geometry the nucleons in the cylindrical
ﬁreball channel, that is, the participants, are removed from the lattice. However, as
stated in [BPDM86], since mostly only 6 - 8 nucleons are removed, the effect is not
very important and basically the same results are obtained for a simulation that does
not take care of the ﬁreball scenario. This has the further implication that also the

impact parameter b does not have a big inﬂuence on the results.
In this model, the breaking probability

Pb = (1 - p) (2.91)

47

is chosen as input parameter, where p, as before, denotes what we called the bond
existence probability. Using Monte Carlo techniques, a computer simulation is per-
formed in which nearest neighbor bonds are broken according to pb. In the basic
version of this model, spin and isospin degeneracy are assumed, so that it sufﬁces
to have one bond breaking probability pb. In the next step, cluster counting is per-
formed and the cluster size is recorded. In the simulation it is then integrated over
the impact parameter b, with the mass and multiplicity distribution for one value
of pb as result. Assumptions made for the realization of this model are as follows.
The motion of the nucleons in the target nucleus can be neglected, which amounts to
saying that the projectile sees a frozen image of the target, which is approximately
valid in the region treated with the model, where the proton energy lies far above
the Fermi energy of the target nucleons. As a further point, no energy conservation
is present in the model. It seems however, that, since the proton energy is sufficient
to completely fragment the target, energy conservation is not an important factor
in the determination of the mass spectra. A physical interpretation of the breaking
probability is possible through an analogy with the Ising model of ferromagnets. It
is possible, as cited in [BPDM86], to show that the breaking probability is related to

temperature in the following way:
Pb = 9XP(—Eb(Tl/T), (2.92)
where for real nuclei at low temperatures

Eb(t) = 8 — $172. (2.93)

The simulation data reproduces the mass yield curves and multiplicity distributions

found in nuclear experiment of proton induced Multi-Fragmentation very well.

The fact that in the basic version of the model only one bond breaking probability

was used is no inherent restriction of the model and extensions of it have been inves-

48

tigated, [Baued], using site-bond percolation with two components and two different
bond breaking probabilities. This system, in which one can model varied isospin
asymmetries, showed the same phase transition behavior as a one-component bond

percolation model with isospin symmetry. Further details about a two-component

percolation model will be discussed in Chapter 3.

49

Chapter 3

Two-Component Site-Bond
Percolation

After reviewing the underlying ideas of percolation and phase transitions in general in
Chapter 2, we want to focus on one special percolation model, that is percolation with
two components on a simple cubic lattice using the combined site-bond percolation
approach. Some aspects of this model were investigated by us in the framework of

this thesis’ research and we will present the results thereof in the following sections.

3.1 General Case: N Components

Let us begin with a description of the model we have used in the most general formu-
lation. Although later on, we will consider a graph of dimension d = 3 in the special
case of the structure of a simple cubic lattice, we don’t make any assumptions about
the lattice dimension or the lattice structure in this general section. The approach

that is presented here is completely free of all such restrictions.

Thus we have a lattice of dimension d, as described in 2.1.1, and we have N

different components with which the vertices are populated. Then there exist

Abp = (N +22 — 1) (3.1)

50

different bond existence probabilities, one for each possible combination of sites to be
connected. The bonds are assumed to be directionless, meaning that the probability
of an edge to be open does only depend on the state of the two vertices it is connecting,
not on their respective localization. There is another set of free parameters associated
with the concentrations of the different species that populate the lattice. These are

the component fractions

N.

f.- = —T'—.
Zi=1Ni

(3.2)

where the N,’s are the actual number of particles of component i. Naturally, the
component concentrations f,- have to satisfy
N
2 f, = 1. (3.3)
i=1
Because of the normalization condition, the component concentrations constitute
(N — 1) free parameters. Overall we then have

Apar=(N—1)+(N+2_1)

2 (3.4)

free parameters a,- to vary.

The main question we are pursuing is at which parameter regime {aznf} in the AW,-
dimensional parameter space an inﬁnite network Coo of connected vertices is present
in the lattice, that is, at which region {bjnf} the percolation probability poo({a,-}) is

IlOIlZCfOI

{b:?"’} = {az’ : p..(a.~) aé 0}. (3.5)

and, more speciﬁc, where and in which form, the percolation probability changes
from zero to a ﬁnite value. The particular type of bond shall be irrelevant for it to
belong to the inﬁnite cluster. Although there exists no principle argument against

this approach, it is quite impractical for a system with N 2 3. It would be highly
51

desirable to be able to reduce the dimensionality of the parameter space. Ideally one
would hope for a reduction of the dependence of the order parameter, pc,o in this case,
on one variable at ﬁxed component concentrations. We propose one such variable in
following the deﬁnitions of [Baued], generalizing them to N components. In analogy
to the bond existence probability in ordinary bond percolation processes, which gives
the density of bonds present in the system under observation, we introduce the scaled
control parameter p+,
N
P+ = Z aijpij- (3.6)
i27=1

Here, the pij’s denote the bond existence probability on an edge that has endvertices
occupied by species i and j. The aij’s give the probability to ﬁnd a nearest neighbor

edge that connects two sites of ﬂavors i and j:

0113‘ = 2fifj

1’ aji- (3.7)

Also the a,j’s are normalized and thus subjected to the constraint

2 an“ = 1. (3.8)

This should provide the essential tool for a simulation aimed at answering the question

related to (3.5) and we will present results in Sec. 3.2.

3.2 Two Components using the Natural Variables

We do not want to go into the detail of the simulation techniques used, at this
point, but merely present the results, deferring a deeper discussion of the simulation

techniques to Sec. 3.3. In Appendix A we show the source code of the programs

52

used in the simulations. We have performed simulations on 163, 323, 643 and 1283
simple cubic lattices and set N = 2, the latter having as only argument the ease of
visualization of the results, which would hardly be possible for N > 2. However, as
there exists no fundamental constraints by setting N = 2, the results shown in this

section should also hold for N > 2.

3.2.1 General Considerations

For simplicity we shall call one species blue, the other red. With N = 2 we have as

number of free parameters
Apar : 4a (39)

which are the fraction of one component, say of the blue sites, fb, and three bond
existence probabilities: pbb for bonds connecting two blue sites, p,.r for bonds with

two red endvertices and p¢ for b-r-bonds. We will set

I): E pbb : prra (3'10)

thus introducing a symmetry in the system. This can be motivated by, for exam-
ple, considerations like isospin symmetry, where the e+e+ and e‘e‘ interactions are

identical. We furthermore have

Pye = pbr (3.11)

= (1 — pz). (3.12)
Given these conditions, (3.6) now reads

19+ = 0:19: + 094% (3-13)

53

where (3.7) has become

01¢

2fb(1— fb)
%’ 2(fb - %)2
ft2 + (1 — fb)2
é‘f 2(fb — ":92-

(3.14)

(3.15)

Here we have replaced the double index by a more intuitive notation for just two

components. Due to the normalization condition (3.8) we have

= = (1—a¢).

(3.16)

We carry out simulations according to these principles and with the techniques as

described in Sec. 3.3

3.2.2 Simulation Results

In Fig. 3.1 we show the results of the simulations, namely p00 as a function of the

two control parameters p: and p¢.

 

We see that the percolation probability poo

 

 

Figure 3.1: Probability to belong to the inﬁnite cluster, poo, in a simple cubic two-
component site-bond percolation lattice of size 1283 as a function of the two param-

eters p: and p¢, calculated for a fraction of the blue species of fb = 0.5.

54

changes from 0 (front corner) to 1 (back corner) for p= —) 1 and p¢ —> 1, with a
critical line of a second order phase transition in the (p.__,, p¢)-plane. We expect some
parametric dependence on the concentration of one of the two species. We chose to
take fb, also possible would be the concentration of the red species in the lattice, f,.
In Fig. 3.1 we have set fb = 0.5. We show the effect of varied fb in Figs. 3.2 and

3.3, where fb : 0.2 and fb = 0.0 respectively. It can be seen that, according to

 

 

 

Figure 3.2: Probability to belong to the inﬁnite cluster, poo, in a simple cubic two-
component site-bond percolation lattice of size 1283 as a function of the two param-
eters p: and p92, calculated for a fraction of the blue species of fb = 0.2.

 
   
     
 

o
u
mo

(WW

9
990000909

9

’9

l
l

to
o)
o

lo
Mm

    
      
    
 

     

   
 

..
o
o O O O - O - -
6%??330 :z::::5:..
' 00....--ooooo-
’ O, '0 O 0.....0.9......’
1 0' aggao’.:.:.:.:.:.,o /
~ " '0 00%?”ot'f/ ’
‘~-. II ”5%?60303“
~ I I’ '0 0.4’

    

 

 

 

‘1'
\i‘
5
§
~
'V
s
6‘

'o
8
lo
0’th
NM‘
N”
N“
H
OH
N
0e
W
M ‘
’~"’l
—A
S

o)
o
,9
'1
Q

 
    
 
 
 
  
   

. o
A

o

,0

9

:v
i
. 3

O

A

o
(o
o
o
(,3

O

i)
.'

   
  
 
  
 

/A
v

Figure 3.3: Probability to belong to the inﬁnite cluster, poo, in a simple cubic two-
component site-bond percolation lattice of size 1283 as a function of the two param-
eters p: and p¢, calculated for a fraction of the blue species of fb = 0.0.

55

the expectations, for only one component populating the lattice, fb = 0 or fb = 1, we
regain a two-dimensional curve which seems to scale, at least on an eye-analysis basis,
as predicted by scaling theory, compare Sec. 2.4.5 and in particular (2.67 ) Before, in
Sec. 3.4, we will come to a realization of the ideas presented in our general discussion
in Sec. 3.2.1, we will explain in more depth the simulation techniques used, in Sec.

3.3.

3.3 Simulation Techniques

The simulation programs are coded in FORTRAN and the source code can be found in
Appendix A. Several programs have been written to accommodate different needs.
First, there is a general program, which determines the percolation probability poo,
the probability that a spanning cluster exists, 1/2, and the number of open clusters per
vertex. Three other programs have been written, aimed at more speciﬁc tasks, two to
determine the critical value of pr through the critical values of p:, p¢, respectively,
and one to ﬁnd the threshold of fb. In the following sections, we will, in a concise

manner, describe the main features of the simulation as implemented in the programs.

3.3. 1 General Features

We aim at simulating a simple cubic lattice of dimension n3, with a maximum lattice
size maxn3. The implementation has been realized as a three dimensional logical array
neutron (maxn,maxn,maxn) in FORTRAN. Using a (pseudo) random number generator
taken from [PTVF92], modiﬁed to work with double precision, we sweep through
the lattice in order to determine which site is occupied by which of the two possible
species. In the source code, the notation is neutron(i , j ,k), meaning in a more neu-
tral formulation the blue components. As already stated in Sec. 2.3, a Monte Carlo

step in the case of percolation just consists of generating one particular conﬁguration

56

of the lattice. We do so by generating a pseudo random number ran at every vertex
3] = (i, j, k) and then comparing it to the probability of having a blue site, fb. The
array neutron(i , j ,k) is then assigned the following values:

TRUE, if ran < fb

neutron(1,j,k) = { FALSE, if ran > fb .

(3.17)

In the next step we have to decide about the existing bonds in the lattice. The bond

status is implemented in the form of a four-dimensional logical array,
connected(maxn,maxn,maxn,3),

which gives the existing bonds at every vertex in the three spatial directions. For
populating the edges with bonds, we essentially proceed in the same way as in the
case of the sites, by sweeping through the lattice and calling the random number
generator at each edge. However there is one difference. Using neutron(i , j ,k) , we
ﬁrst have to decide whether we have an edge with different endvertices or one that

connects two sites of the same ﬂavor. The program then assigns

. . __ TRUE, if ran < p
connected(1,j,k,x) — { FALSE, if ran >p , (3.18)
where, in the case of a: = 3, which is the positive z-direction,
_ p2, if neutron(i,j,k) = neutron(i,j,k+1) (319)
p ’ p¢, if neutron(i,j ,k) 7e neutron(i,j ,k+1) '

We use free boundary conditions, so that all values of connected(i, j ,k,x) on

the boundaries are initialized as FALSE:

connected(n, i ,j , 1) = FALSE
connected(i ,n,j ,2) = FALSE
connected(i,j ,n,3) = FALSE, i,j 6 [Ln] (3.20)

Furthermore a third array, an integer value array this time, is used to save the infor-

mation on the cluster status of every vertex. The integer value assigned to every

57

entry of clusterNumber(O:maxn1,0zmaxn1,0:maxn1) is the number correspond-
ing to the place in row at which the cluster ﬁnd algorithm, see Sec. 3.3.2, found
the cluster to which the present vertex belongs. Here maxnl 2 man + 1. The
value clusterNumber(i, j ,k) = 0 is meant to indicate that (i, j, k) has not yet been
checked by the algorithm. The zeroth and n + 1st layer vertices (in each direction)
are assigned cluster numbers hcln equal to the highest number of clusters occurring

if the lattice is fully fragmented, plus one:
hcln = n3 +1. (3.21)

This provides for a ‘surface cluster‘ that seems already checked, which lets the cluster
algorithm stop at the natural lattice boundaries. In the next step the lattice con-
ﬁguration is analyzed with the cluster ﬁnd algorithm and the appropriate values are
recorded. The whole Monte Carlo step is repeated new“ times, using independent
random number seeds each time, which are obtained from the system clock, with

new“ on the order of 100. Then the averages of the observables are calculated.

3.3.2 Cluster Analysis

We use a cluster algorithm ﬁrst described in [BPDM86]. We abide by the same
notation as used in the cited reference for the description of the algorithm. The set
of all vertices in the lattice is denoted by A. The set of vertices that have not yet
been checked by the cluster ﬁnd algorithm and are thus still uncounted, is written
as Aml’", with Am“ C A. The task is to ﬁnd all clusters 0,, and record their size
|C,-|. The algorithm accomplishes this task by scanning the lattice for the ﬁrst vertex
2", belonging to the subset of unchecked vertices: 2'; E Aml’". Now all other sites

2",- 6 Am’", connected by open edges to 2",, 2} +—> 2", have to be detected. Am“

58

denotes the next generation of unchecked sites,
A(1),i : A(0),i _ 0(0).i.

0(0)" 2 {2}} is called the zeroth generation of 0,, where C,- = {2", . . . 2",} is the cluster
with |C’.-| = n which we want to detect. In the next step all previously unvisited
neighbors of 2", are checked if they are connected to the site from which the search
is originating. All neighboring sites that are connected to 2", form the so—called ﬁrst

generation 0(1)" of 0,. At some n we will have

A(n+l),i ___ A(n),i _ 0(n—1),i : A(n),i,

as C("l’" = (0, which is the termination condition for the algorithm. The cluster C,- is

then given by

a=UwW ﬁx)

The algorithm restarts screening the lattice until a new starting vertex 2",“, belonging
to Am)“1 2 Aml'i — 0,, is found. All clusters have been found if Awl’N = (l) and N is

then the number of clusters in the present conﬁguration.

3.3.3 Determination of the Critical Value

The second type of program written is aimed at the question of determining the
critical value in a more precise manner than from extrapolation methods based on
phase transition plots, of, for example, poo. Three different versions of this type of
program have been written, which determine the threshold for the onset of percolation
in terms of pz, p¢ and fb. In the description here, we will only describe the one that
gives pi through determining pi. The ﬁrst steps in all three programs are the same

as described in Sec. 3.3.1, in that sense as that the procedure to populate the lattice

59

and the edges are the same. They stand, however in a different framework, which is
described in the following. Every Monte-Carlo step now comprises an nested interval
method approach to the critical value. We repeat each Monte Carlo step nevent times,
with nevent again on the order of 100 and using different random seeds for each Mont
Carlo step. For each value of fb the program is given the values of p¢ that are to be
probed as well as a start interval for p: E [a, b]. The actual value of p: that is to be

probed in the ﬁrst sweep of the lattice is determined according to
p==c=/\b+(1—/\)a. (3,23)
A typical value that we used for /\ was
A = 0.9.

With the random number assigned for this Monte Carlo step the lattice is populated
and the edges determined to be open or closed. Then, the cluster ﬁnd algorithm,
as described in Sec. 3.3.2, goes through the lattice. This time, however, as the aim
is to determine the critical value, not the whole cluster information is needed, the
only interesting question is whether there exists an inﬁnite cluster as deﬁned in Sec.
2.4.1 with (2.10) and (2.11). Therefore, the cluster ﬁnd algorithm is constrained
to one boundary layer of the lattice in one direction, which in the program is the
(i, j, n), i, j E [1, n] layer. From there, the algorithm follows the branching of every
cluster but stops immediately if it reaches the bottom of the lattice, that is the
(i, j, 1), i, j E [1, n] layer. We then only have two possibilities: Either the algorithm
stops after checking all clusters originating in the (i, j, n), i, j 6 [Ln] layer and
not ﬁnding any spanning cluster, where we assign 1b = 0 or it stops after ﬁnding a

spanning cluster and 11) = 1 is assigned. Then the new pZ-interval is

l

(3.24)

 

a,c, if 112:0
c,b, if $21'

60

This process is repeated a desired number of times, depending on the size of the ﬁnal

interval one wants to attain. Usually, we set as ﬁnal interval size
|d — cl 2 104,

where

depending on the outcome of the last step. The value of Cﬁnal from the ﬁnal step then

gives the critical value

1952‘ = cﬁnaz. (3.25)

After completion of the new", Monte-Carlo steps averaging over all obtained p? yields

the threshold:

1 " .
(i = — E'. 3.26
p- n gt). ( )

The critical value of the scaled control parameter p+ is then simply given as

p: = 0:1): + a¢p¢. (3.27)

A calculation of the statistical errors is performed as well. However, to obtain the
total error, the systematic error, due to ﬁnite size effects, has to be added to the
Statistical errors. We will come back to a discussion of the systematic errors and

Ways to estimate them in Sec. 3.7.

3.4 Two Components using the Scaled Control Pa—
rameter

After giving a review of the techniques used in the simulations of which we present

the results, we now want to come back to the discussion led in Sec. 3.1. In (3.6) we

61

introduced a scaled control parameter with the intention of being more capable to
extract information about a possible phase transition and to be able to describe the
phase transition behavior in terms of just one variable. In Sec. 3.2 we then discussed
our model, in the case of two components, still in terms of two variables p: and p¢.

Now we want to realize a description in terms of the variable p+.

3.4.1 Entire Parameter Range

First, we will aim at a discussion of the data from simulations where, as before,

fb = 0.5 and then come to present the results for fb 79 0.5.

Blue Species Concentration fb = 0.5

We take the same data as shown in Fig. 3.1, that is, we have p._., p¢ E [0, 1] and analyze

it in terms of the scaled control parameter p+. The result is shown in Fig. 3.4. Two

 

 

 

 

 

 

 

 

 

 

1 . ' ' .............. ..................................... ...................................... ................................. .
1 .-
0.8 - i
8
Q.
0.6 r
8
a .................. . 0.1
0.4 ~ l
‘, " 0*(p+'pS)B
' poo(p==0. 12¢) A
0-2 " f; Poo(P=1P¢=O) ° 0.01....“ ....m. ...“. ...-.m‘
;; p..(p=¢0.p¢¢0) o 0.001 0.01 60.1 1
o , s... : 1 I 1 p+.p1+
0.3 0.4 0.5 0.6 0.7 0.8

17+

Figure 3.4: Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter, p+, in a 1283 simple cubic lattice with fb = 0.5, for values
of p=,p¢ 6 [0,1], ﬁtted with pc,o oc (p+ — pi)” (solid line). Here pi 2 0.251 for the
’upper’ branch, pi 2 0.280 for the ’lower’ branch and ﬂ = 0.41 in both cases. The
inset shows the same data in a double logarithmic representation.

62

distinct branches can be seen, both in the shape of a second order phase transition.
Fig. 3.4 ﬁrst suggests that the scaled control parameter p+ is a good variable, as
we only have two universal scaled curves. However, it suggests too, that some other
phase transition happens in the system as well. On a ﬁrst eye-analysis it seems that
we see the system in two different states, which both show, as a function of p+, a
similar behavior in the percolation probability poo. Before we come to an analysis
about how the two phases, of which the two universal scaled curves are an indicator,

come about, we will devote some time to an analysis of the curves themselves.

In the main part of Fig. 3.4 we show the data on linearly scaled axes. Both
are in the form of a second order phase transition and we ﬁt both curves with the

expectation from ordinary one component percolation theory, see (2.67),

poo = a (19+ - Pi)ﬂ- (328)
With the ﬁtting procedure, we ﬁnd for the critical value of the ’upper branch’

pi = 0.251 :1: 0.002. (3.29)

Taking ﬁnite size effects into account, for which we did not correct at this point, this
stands in reasonable agreement with other numerical results for one component bond

percolation theory, see for example, [SA94], which state
p: = 0.2488, (3.30)

or [LZ98], where the threshold is given very precisely as
pi = 0.2488126 :l: 0.0000005. (3.31)

The value of the lower branch, however, is shifted to a higher threshold. With the

same ﬁtting procedure we ﬁnd for the critical value of the data belonging to the ’lower

63

branch’
pi = 0.280 :t 0.002. (3.32)

Another important quantity that should elucidate the situation is the critical expo-
nent 6. To test upon it, we print the same data as shown in the main part of Fig.
3.4 in a double logarithmic plot, the result of which is shown in the inset of Fig. 3.4.
We can see that in the double logarithmic display both branches are well-ﬁtted by
linear functions with the same slope B, that is, on a linear scale with the same critical

exponent, which we ﬁnd to be
(3 = 0.41, (3.33)

for both branches. This is the same critical exponent as found by other authors for
one-component bond percolation theory and as given, for example, in [SA94]. We
should note again that the critical value of p+ depends on the dimension and the
lattice structure of the system, whereas the critical exponents just depend on the

dimensionality.

We now come back to the question which two states we see manifest in form of
the two poo curves. A further analysis of the data unveils that the ’upper branch’
constitutes of points with both p: and p¢ non-zero’, whereas all points with (p: =
0, psi 75 0) and (p: 96 0, p¢ = 0) fall on the ’lower’ branch. In a shorter notation we

can summarize this behavior in the following way

1?: 11¢ # 0, ’upper branch’ (3.34)

1): 13¢ = 0, ’lower branch’ (3.35)

This leads to the conclusion that some fundamental change has to take place in the

 

1We have to note that the term ’non-zero’ is somewhat imprecise, we will further elaborate this
point in Sec. 3.5.2

64

system, when the system is brought to either one of the two limits
I): —> 0.

or
p¢ —~> 0.

Together with the results stated above we can now conclude that, as long as both
bond types are active, which corresponds to the (p: 96 0,p¢ 7t 0)-regime, the two-
component site-bond percolation system under observation here and the ordinary
one-component bond percolation system show an identical phase transition behavior,
which is in accordance with the ﬁndings presented in [Baued]. Obviously this means
further that the actual values of the individual bond existence probabilities p: and p¢
do not have a direct inﬂuence on the percolation probability, but rather that the actual
bond density p, is the direct equivalent of the control parameter in one-component
percolation and thus fulﬁlls the expectations under which it was introduced. As soon
as one of the bond types is disabled, however, the system fundamentally changes its
behavior, going over to a new critical value of p+ for the order parameter poo. The
basic features are retained, a second order phase transition of the order parameter,
poo, is still observed and furthermore with the same critical exponent 3 as in the

non-zero regime of p: and p¢.
Blue Species Concentration fb 75 0.5

Probing at another value of fb yields essentially the same effect. As an example, in
Fig. 3.5, we show data for fb 2: 0.3. Here, the curve that we called ’upper branch’,
that is the data for values p= ;£ 0 and p¢ 71$ 0 can be shown to be exactly the same

as in the previous case, where fb = 0.5, we retain precisely the same critical value

pi r: 0.251
65

 

 

0.8 - _

90°
00°
00°.
,0
,0
o
o
o

   

0.6 r /
08.
04 >- ‘ """"""""""" .,
p..(p.=0. 9..) A
p00(pa’ p¢=0) 0
0.2 ~ p0,,(p.¢0, p¢¢0 o

 

 

a*(p+-p‘.’)

 

 

l L

0 'I .. 1 1 1 1 1
0.2 0.25 0.3 0.35 0.4 F(3.45 0.5 0.55 0.6 0.65
4.

Figure 3.5: Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter, p+, in a 1283 simple cubic lattice with fb = 0.3, for values
of pz, p¢ 6 [0,1], ﬁtted with poo oc (p+ — pi)ﬂ (solid line). Here pi 2 0.251 for the
’upper’ branch, pi 2 0.218 for the p9,; = 0 branch, pi 2 0.266 for the p= = 0 branch
and [i = 0.41 in all cases.

and the same critical exponent 3. For the data that constitutes of points where either
p: = 0 or p¢ = 0, the situation is changed. First, the data for p= = 0 and p¢ = 0 do
not fall on the same universal scaled curve, but rather show different limiting values

for p+ —> 1 and different critical valuesz. We ﬁnd

pi(p= = 0,p¢) = 0.266 :1: 0.002. (3.36)
for the ’p: = 0’-branch and

pi(pz, 19¢ = 0) = 0.218 :1: 0.002. (3.37)

for the ’p¢ = 0’-branch. In Fig. 3.6 we show the same data3 as depicted in Fig. 3.5

in a double-logarithmic plot. As can be seen, all three curves are very well ﬁt with

 

2When we say p+ -> 1 limit here, we mean to indicate the direction in which we let p+ go, being
well aware that 104. = 1 can not be attained for most values of fb whilst p= = 0 or p; = 0.

3 We have to note that, since the data shown in Fig 3.5 is lower quality data then the one depicted
in Fig. 3.4, we have, in Fig. 3.6, corrected for some of the stronger ﬂuctuations seen in Fig. 3.6 and
left out some data points

66

 

 

p...

0.1 r P049990. p¢¢0) o
: poo(p==o! p¢)

: poo(p:! p¢=og <>
' a*(p+'pE-)

 

 

 

 

0.01 . LAM-1...-..2 -
0.01 c 0.1 1
p+'p+

Figure 3.6: The same data as shown in Fig. 3.5 in a double logarithmic representation,
with a slope of all three curves of 6 = 0.41.

linear functions of the same slope 6, that is the same critical exponent H in the linear

scale. As was already the case for fb = 0.5, we again have
B = 0.41. (3.38)

for all three branches. A different behavior in the limiting values of p+ and the critical
value pi for the ’p= = 0’-branch and the ’p¢ = 0’-branch has to be expected, since in
the system, we have a symmetry in the concentration with respect to fb = 0.5. This

symmetry is made possible through us setting

P: : Pbb : prri

which leaves edges connecting two sites of the same species and edges connecting two
sites of different species as only possibilities. We then have a ’binary’ system with
a: and cry; as the two quantities that show a symmetry with respect to f), = 0.5 and

furthermore have

axn=oo=aan=99.

67

see (3.14), (3.15) and Fig. 3.7, which leads to the coinciding of the (p= = 0,p¢ ¢ 0)

and (p: # 0, p¢ = 0) curves as shown in Fig. 3.4. The fact that the limiting values

 

 

 

 

1 \\‘ I T
\‘ a=(fb) ........
0 8 09(0)) —
=3" 0.6
3+ ,,,,,,,,,,,,,,,,,,,,,,,
if:
\Tl 0.4
d
0.2
0 1 1 1 1
0 0.2 0.4 fb 0.6 0.8 1

Figure 3.7: The normalized density of edges connecting two sites of equal ﬂavor,
a:(fb) and the normalized density of edges connecting two sites of diﬂerent ﬂavor,
a¢(fb). See also (3.14) and (3.15).

poo(p+ —) 1) of both curves differ can be understood in the following way. Having set
p: = 0 and letting p+ —> 1 requires letting p75 —> 1. Then the maximum value of p00

is determined by fb:

1900(1): = 0.792 = 1) = min(fb. (1 - ft»). (3.39)

with a symmetry again with respect to fb = 0.5. The same line of thought applies to

the case where p7; = 0 and p: = 1, however due to the behavior of a: we have

poo(p= :1,p¢ = 0) = max(fb,(1— .00) (3.40)

It has to be noted that the equality in (3.39) and (3.40) is not strict, since geometrical
effects play an important role. By geometrical effects we mean the fact that due to
the geometric position of the particles on the lattice some of the particles will be

inaccessible by either, depending on the limit, bb—, rr- or br-bonds and thus a part

68

of the whole fraction fb or fr = (1 — fb) will be excluded. We want to elaborate on
the footnote (2) made above to clarify the notation p+ —) 1. Having set one of the
bond existence probabilities equal to 0, the maximum value of p+ is given by a: or
a¢ respectively. Thus, unless fb = 0 or fb = 1 the real limit p+ = 1 can never be
reached once either p: = 0 or p¢ = 0, rather will it be impossible to probe above a

maximum value of p+, which is given as

max

p+ : a¢1

OI'

p+ : (1:1

depending on whether (p: = 0,p¢ : 1) or (p: = 0,p¢ = 1) was set.
3.5 The Percolation Probability

In the afore—held discussion on the behavior of the percolation probability as a function
of the scaled control parameter p+ it became clear that in a wide range of values of the
’contributing’ parameters p: and p¢ the percolation probability exhibits features sim-
ilar to the ones observed for the percolation probability in an ordinary one-component
bond percolation system, but that when setting one of the contributing parameters
to zero, the system would behave different. In this section we want to present results
that show how this different behavior comes about. We will more or less focus on the
limit p¢ = 0, as the same line of arguments applies in the p: = 0 limit, but we will

still present the results from an investigation in the second limit too.
3.5.1 p95 —> 0 Limit

Going back to Fig. 3.4 we can choose an arbitrary value of p+ > pi(’lower branch’),

so that the corresponding percolation probabilities will be non-zero and compare the

69

percolation probability that we have on the upper branch with the one that one

obtains through the lower branch. One realizes that there exists a ﬁnite difference:

2980(1): ,4 0,292 79 0) - poo(p=.p¢ = 0) > 0- (3-41)

The question arises, how the percolation probability behaves when the system is let

into the p¢ = 0 limit.

For simplicity, we ﬁrst choose
p: = 1.
We will furthermore, for the sake of ease of argument, set
fb = 0.5

in this ﬁrst discussion and then later on state the differences that occur for other
values of fb. Then setting p¢ = 0 corresponds to a lattice in which all available
bb—edges and all available rr-edges are open, but all br-edges are closed. With a
very high probability there will be an inﬁnite network present. This comes about
in the following way. Effectively with this choice of parameters, we are modeling a
two-component site percolation model with two disjunct types of clusters. These are
clusters of b-type and clusters of r-type, which are entirely disconnected as p¢ = 0 was
chosen. Then, for site percolation, for the given fb, we are above the site percolation

threshold in the case of a simple cubic lattice, which is, compare [SA94],

C = 0.3116. (3.42)

site

Since the second species has a concentration f, = (1 - f b) scite above the threshold

as well, it is highly probable that both species will have a spanning cluster. According

to the deﬁnition we set forth for an inﬁnite cluster (2.10) and (2.11), there will be

only one inﬁnite cluster Coo, however, which has to fulﬁll the additional requirement

70

of being the biggest cluster. Let us just assume that the biggest, spanning cluster,
which we identify with the inﬁnite cluster, is of blue color. Following (2.66) and

(2.65), we can express the percolation probability as

1W(¢=0)=fb-:n2,s (3.43)

where, like before, the symbolic notation ’oo — 1’ for a ﬁnite system is meant to
indicate that in the sum the inﬁnite cluster is excluded. From a simple argument,
we can approximately evaluate the sum, without exactly calculating every term in
the cluster size distribution n2 .Given the fact that p: = 1 and p¢- — 0, according
to the argument presented in the preceding section, see (3.40), we will encounter a

maximum value of poo,
= f, = 0.5. . (3-44)

It has to be noted that this value is an upper limit, typically the percolation proba-
bility will come to lie slightly below this value, as conﬁrmed in simulations, where we

ﬁnd
p00 2 0.48. (3.45)

The reason are geometrical effects. By geometrical effects it is meant that some blue
sites will be inaccessible from other blue sites, thus leading to the fact that not the
whole fraction of fb blue sites is accessible. In (3.43) 71.2 denotes the number of b-
clusters of size 8. As already stated, there is a high probability that also the second
color has a spanning cluster. If not so, then at least the biggest r-cluster, 05m,

will have a size on the order of the size of the inﬁnite cluster, since both colors are

symmetric,

|~ ~ |le- (3.46)

lCr max

71

Analogously to (3.43) we can formulate the probability that a given lattice site is part

of the biggest r-cluster, C' as

man: 7

max— 1

192.1002 = 0) = fr — Z n; s, (3.47)

3:1

and with the same line of thought as for the blue species we ﬁnd

r,maz

pm... = f.— = 0.5. (3.48)
Again, the true value will be slightly below that, where, again,
pin“ 2 0.48, (3.49)

is a typical value. Starting from this conﬁguration where poo S, 0.5, and pin“ S, 0.5,

we now set

p95 2 E (3.50)

This will introduce a very low br-bond density in the lattice. The inﬁnite cluster in
this new conﬁguration of parameter values can then consist of both 0- and r-vertices

and we can write the percolation probability p00 as

oo—l
poo(p;é = e) = 1 — Z n. s. (3.51)
3:1

This yields a difference in the percolation probability at p91 : 0 and the percolation

probability at p¢ : 5:

600(5) = 1900022 = 8) - 90000.4 = 0) _
: 1—ZnSs—(fb—ans). (3.52)

By rewriting the concentration of the blue species, fb, as

fb = (1_ fr)
: 1 — Zn; 8, (3.53)
8:21

72

we obtain

mas: oo—l

600(5) : (00:16; 3 + 2n; 3) — Z n, s. (3.54)

3:1

The difference (500 has the property that, if we let 5 —+ 0, for all e > 0, the limit will

be non-vanishing:
lirré 600(5 > 0) > 0. (3.55)

This comes about, since, with a very high probability, the bonds making up for this
low density will connect the two biggest clusters in the lattice, C00 and C' as

man: 3

according to our consideration above, we have,
9.... + 192...: § 1- (3-56)

Thus almost surely do we have (3.55) fulﬁlled.

How is the situation in a simulation on a ﬁnite lattice? In this case the situation
is slightly different. Setting the bond probability p¢ to some ﬁxed value p¢ > 0 does
yield a corresponding bond density, but just on average. Events where no bonds at
all are introduced in the lattice do have ﬁnite probabilities. What is the implication?
Let us consider a case where we set

P¢=€

= 1 an)

1
N edges -

 

which introduces on average one bond < 5,1" > in the lattice. For each event where
such a br-bond does occur, almost surely < 5,? > will connect the inﬁnite cluster

Coo and the biggest cluster of the red colored species, C" We rely on the same

man: '

argument as above, which states that we have probabilities to belong to the inﬁnite

cluster and to belong to the biggest red cluster, both close to 0.5,

Poo 1:93.151-
m

Then, the ’new’ inﬁnite cluster will be made up out of the ’former’ Cc,o and 0;,”
together. This will yield, as in the case of the inﬁnite lattice, a higher percolation
probability compared to the p¢ : 0 system. The rare events in which the br-bond
< 5,7" > does not connect C00 and Cfnax or where it connects C00 to some ﬁnite
r-cluster, CE", as well as the events in which no br-bond at all is formed, will exhibit

a smaller increase in the percolation probability or none at all. In a simulation an

average over all events will yield an eﬂective value of (50°, denoted by 6:017 , with
0 < 655,7 < 600. (3.58)

The arguments presented, with (3.55) for the inﬁnite lattice and (3.58) for the ﬁnite
lattice, lead to the conclusion that in the limit p¢ ——> 0 a ﬁrst order phase transition
takes place in poo.

The simulation results support this conjecture, as shown in Fig. 3.8. It has to be

1

 

 

? l ? \i} {i} T r
0.9 -
pi (1251.19...) :9:
0.8 ~
3
5 0.7 ~
0.6 ~ .
0.5 ’ e

 

 

 

0 1 2 3 7 4 5
A=(p+-0.5)10
Figure 3.8: Probability to belong to the inﬁnite cluster, poo, as a function of the scaled
control parameter p+, in a 1283 simple cubic lattice with fb : 0.5, at ﬁxed p: : 1
for varied p91. 50 independent simulation events where taken into account and the

average was taken over events in which the number of br-bonds actually formed in
the simulation, nbr, was non-zero.

noted that for the data shown here in the analysis the average was taken only over

74

events in which the number of br-bonds, nb,, actually produced in the simulation,

was non-ﬁnite:

Tim-7&0.

Another possibility is to average over all events, which gives a smaller 637, as one also
includes events in which the lattice is completely unchanged and the same percola—
tion probability is found. Data from the same simulation, analyzed with the second
averaging technique, where all events were taken into account, is shown in Fig. 3.9.

Naturally the question arises, what the behavior is like for p: 31$ 1. Still, a ﬁrst order

 

1. . . ..
0.9- 61¢f‘l’

0.8 . j (2 p.23 <p==1.p..) 1—6—1

i

9

poo

0.7 ~
0.6 r

0.5 " e

 

 

 

2 3 7
A=(p+-O.5) 10

Figure 3.9: Probability to belong to the inﬁnite cluster, poo, as a function of the scaled
control parameter p+, in a 1283 simple cubic lattice with fb : 0.5, at ﬁxed p= : 1
for varied p¢. 50 independent simulation events where taken into account and the
average was taken over all events.

phase transition will be observed, since essentially nothing is changed, compared to
the case where p: : 1. The effect will be somewhat less pronounced, that is a smaller

500 will be seen. This comes about since

ICoo(P= <1)| < lCoo(P= =1)|.

75

and

ICIMAI): < 1)| < ICEaAp: =1)|-

It is clear, however that

 

p_ > 191(1): 74 0,192 75 0)

_ a=(fb) (3.59)

is necessary to obtain this phase transition, as otherwise
poo : 0-
Equation (3.59) represents a condition p: and fb combined are subjected to. For
101(1): at 0.1% at 0) pi(p=,p¢ = 0)

a=(fb) < p: < a=(fb) (3’60)

 

 

the above argument, which relies on the existence of an inﬁnite cluster in the p¢ : 0
limit does not hold anymore in exactly the same way. However we merely have to
replace the term ’inﬁnite cluster 000’ by ’biggest b—cluster’ and then rely on the same
conduct of argument again. There will be no percolation being realized with p¢ : 0,
but connecting C3,“ (instead of Coo, as before) with 0;,” at p¢ : 5 creates an
inﬁnite br-cluster. In the case of p: satisfying (3.60), the effect will be diminished ,

as already explained afore, as the sums in (3.54) turn out smaller.

What is the situation like for fb # 0.5? In order to have the original argument
hold for p: : 1, the system has to stay above the site percolation threshold as given
in (3.42). Below, the same consideration as above applies, there will only be ﬁnite
clusters present and again, the term inﬁnite cluster has to be replaced by biggest
ﬁnite cluster to follow the same line of argument. The ﬁrst order phase transition

exhibits a smaller discontinuity with decreasing fb for

fb < 0.5,
76

or with increasing fb for
fb > 0.5,

and the effect continuously vanishes for fb —> 0 and fb —> 1.
3.5.2 The ’Smearing—Out’ Effect

Some arguments of the discussion from the preceding section can be applied in the
question of an effect that we call ’smearing—out’ effect, which we already brieﬂy com-
mented on in a footnote in Sec. 3.4.1. The observable we are interested in is the
percolation probability poo. We set p¢ and p= to some ﬁxed value p=, p¢ E [0, 1], and
determine the percolation probability. We then show the percolation probability p00
as a function of p¢ and take p: as parameter, the result of which is shown Fig. 3.104.
It is noteworthy that we transformed the real parameter dependence on p: that was
present in the simulation to a parametric dependence on p+ since then the known
values of thresholds are more easily compared. What leads to the functional forms
depicted in Fig. 3.10? At p+ : 0.25 there exist only ﬁnite b- and r-clusters, since the
system is below the threshold, pi 2 0.28. Then, increasing p71 only connects ﬁnite
clusters, leading to a small 600, leading to the nearly continuous looking function. For
a slightly higher value of pr : 0.26 the eﬂ'ect is quite more important and yet more
astounding is the change from the curve where p+ : 0.26 to the one for p+ : 0.28,
where the claim of a ﬁrst order phase transition seems more obvious than in the cases
cited before. For completeness we also show curves for p+ : 0.30 and p+ : 0.5, where
it becomes clear that the (50° observed keeps becoming bigger with increased p+. At
the same time, the width of the transition region, where effects discussed above in

the context of 637 are of importance, becomes increasingly smaller. This transition

 

4We would like to note that the small ’hole’ in the ’right hand side’ of the simulation data is no
inherent effect, but rather due to a lack of simulation in this region, which however is symmetric to
the ’left hand side’ of the curves.

77

 

 

 

 

 

 

 

 

 

0.9 - Doom... 9. = 025) o .

0.8 ~ 9.49... 9+ = 0.26) o .

0.7 - "w p..(p,., p+ = 0.28) . _

4” 1%., 9.49... 9+ = 0-30) .

0'6 ' 1w: ““““““““““““ 1__ . °p.,.,(p,., p+ = 0.50) . ‘
0’5 0.5 137‘ 1

0.4 4 .

W

03 ”- 000 0° o

0.2 - o o

0.1 ° o

0 “. 1 . 1, ‘ . .
0 0.2 0.4 0.6 0.8 1

pat-

Figure 3.10: Percolation probability as a function of p96, with a parametric dependence
on p+, shown for a 1283 simple cubic lattice with f1, : 0.5.

region and its width have an effect when the percolation probability is shown as a
function of p+. It causes the two ’branches’ shown in Fig. 3.4 not to be clear, but
’smeared out’ around the most probable ’plateau’ value observed in Fig. 3.10. In
Fig. 3.4 we corrected for this smearing out by only taking into account data points
lying on the ’plateau’ seen in Fig. 3.10. Thus we can now precise or notation of ’p:
and p¢ non-zero’, clarifying that we mean values of p: and p¢ chosen in such a way
that the effect of the transition region, due to the ﬁrst order phase transition having
smaller discontinuities for smaller p: (smaller p+ in the ﬁgure), can not be seen, that

is, values chosen in such a way that p00 comes to lie in the plateau region.

3.5.3 p: ——> 0 Limit

In this section we want to brieﬂy show that the same results hold for p: —-> 0 as
in the preceding section for p¢ —+ 0. The backbone of the argumentation, which
is the presence of two clusters of essentially the same size in the lattice at p: : 0,

is retained, however the argumentation how they come about has to be changed.

78

Again, for ease of argument, we ﬁrst set 19,5 2 1 ﬁxed and vary 1): from p: = 0
to p= = 8. For (p: = 0,p¢ = 1) the situation will be as follows: There are only
br-bonds allowed in the lattice, leading, with this combination of parameter values,

to an inﬁnite br-cluster Coo. The percolation probability, poo, is given again, compare

(3.39) and (3.44), by
poo é fb = 0.5. (3.61)

Although, contrary to the p¢ = 0 limit, no second species is left to provide a second
cluster with a size on the order of the size of the inﬁnite cluster Coo, there is a second
br-cluster present in the system that we denote by Cm” and that, as was in the 19¢ = 0
limit the biggest red cluster to the biggest blue cluster (assumed to be the inﬁnite
cluster), is symmetric to the inﬁnite cluster. Thus, with a very high probability, Cm“

spans the entire system as well. For 000 and Cm“ we have

lCmaxl z ICool, (3.62)
yielding, as before,

poo + pmaz g 1. (3.63)

How does Cm” come about? This feature is unique to lattices where, from a given
vertex 2 it is not possible to reach every other vertex 2:, on the lattice with a path
2 <—+ z,- of length m odd and a path of length n even. This is the case for the simple
cubic lattice. Starting from some blue vertex 2? = (k, l, m), which coordinate integers

satisfy
k+l+m=(2n+1),n€N, (3.64)
it will be impossible to reach another blue vertex 2? = (r, s, t) with

r+s+t=2n,n€N, (3.65)

79

since in this conﬁguration, with all bb- and rr-bonds disabled and only br-edges open,
always an even number of bonds is needed to get from a vertex with one color to
another vertex of the same color. Thus in this conﬁguration clusters incorporating
’odd’ vertices of one color and clusters incorporating ’even’ vertices of the same color
will be entirely disjunct, providing for the two systems that both have, with a very

high probability, spanning clusters. If we now let

> 0, (3.66)

a low density of bb- and rr-bonds is introduced in the lattice. The bond(s) introduced
by p: = a open, with a high probability, an edge connecting Cm” and Coo, leading

to a non-vanishing 600,
lim 600(5 > 0) > O. (3.67)
e—>0

Another way of arguing, which is restricted to the special case where f), = 0.5 involves
the system as discussed in Sec. 3.5.1 and a symmetry transformation. Before we pro-
ceed to give the transformation, let us introduce some short notation. All operations

are done in Z / Z2. The color of a vertex is denoted by

_ 0, if (k,l,m) blue
c(k,l,m) _ { 1, if (k,l,m) red ’ (3'68)

and we deﬁne the parity of a vertex as

0, if (k+l+m)=2n,nEN

p(k,l,m)={ 1, if (k+l+m)=(2n+1), nEN' (3'69)

The symmetry transformation is described by the following transformation rules:

19’: = 10¢ (3.70)
I); = p: (3-71)
c(k,l,m)' = c(k,l,m)+p(k,l,m). (3.72)

80

This changes the color of every ’odd’ blue and of every ’even’ red vertex. This
means that on average half of the blue species are changed to red and half of the red
species are changed to blue. Thus the concentrations are preserved and obviously, by
conserving the species concentrations f), = 0.5 and f, = 0.5, the densities of accessible
edges, a:(fb) and a¢(fb) are conserved too. Every br-edge gets transformed into a
bb— or rr-edge and every bb- and rr-edge becomes a br-edge. Since, however, the bond
probabilities p: and p¢ are transformed according to (3.70), these edges get occupied
with the same probability as the br-edges before. Thus, as already stated, (3.70)

represents a symmetry transformation for the system where f), = 0.5.

Applying this transformation to the special case (1): =0, 19¢ =1) yields a transfor-
mation to the system discussed in Sec. (3.5.1), where (p: = 1, 13,5 = 0) or vice-versa.
Thus the simplest explanation of the behavior encountered here is to refer to the sys-
tem discussed in Sec. 3.5.1 and to the arguments presented there and then transform

the system with (3.70) to the system interested in in this section.

For simulations on ﬁnite lattices the slight changes in the line of argument as
discussed in Sec. 3.5.1 apply, to which we merely want to refer at this point. We
support the claim of conjecture (3.67) with Fig. 3.11, which reproduces the behavior
shown in Fig. 3.8. The two methods of averaging discussed above are also applicable
here, but, since Fig. 3.8 already exhibits the same features as Fig. 3.11 we restrict
ourselves to this one averaging technique. In this limit here, (19: = 0, p95 2 1), the
system has to be above a different percolation threshold, ff, see [HB99], which we

will discuss in Sec. 3.6.5.

81

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

)r 5 . G 9 .1 .
f w (l) (I) (l) ‘1’ (P (P
0.9 ~ I a
poo (p¢=1ip:) h—Gﬁ
0.8 -
08.
0.7 ~
0.6 -
0.5 " e
0 1 2 3 4 5

A = (p+ - 0.6)107

Figure 3.11: Probability to belong to the inﬁnite cluster, poo, as a function of the
scaled control parameter p+, in a 1283 simple cubic lattice with f), = 0.5, at ﬁxed
p¢ :2 1 for varied 19:. 50 independent simulation events where taken into account and
the average was taken over events in which the number of br-bonds actually formed
in the simulation, 71.5,, was non-zero.

3.6 The Critical Value of the Scaled Control Pa-
rameter

Obviously the question arises, what the behavior of the critical value of the scaled
control parameter, pi, is like. In Sec. 3.4.1 the results presented pointed at a concen-
tration dependence of the threshold. We will again, as already done in the previous
section, ﬁrst focus on the limit p¢ ——> 0 and then turn our attention to the case where

p: —) 0.
3.6.1 p75 —-> 0 Limit

From Fig. 3.4 in Sec. 3.4.1 it was already clear that when setting the system into
a limit in which only one type of edges can be open, the percolation threshold, or
to be exact, the critical value of the scaled control parameter, pi, would take a

different value than in the case where all types of edges could be open. What Fig. 3.4

82

could not show was the way this transition came about. How is the critical value pi
transformed? In Sec. 3.5.1 we found that in the limit we focus on here, p¢ —> 0, the
percolation probability undergoes a ﬁrst order phase transition. Is the same true for
the percolation threshold? In a simulation specially aimed at answering this question,

we set p¢ to ﬁxed values,

19¢ 6 [0,1]

and then apply the algorithm described in Sec. 3.3.3, which uses the nested interval
method, to determine the corresponding critical value of pz, pi. With p¢ and p;
the percolation threshold pi can be calculated. We show the results in Fig. 3.12.

It can be seen that the transition is continuous. Unlike the percolation probability,

 

0.28 I f

4?-

0.27 ~

0.26 ~

 

0.25 -

024 _ Pi (Qt) ”—‘H
1/(u+v*p¢)+t

 

0.23 r

0.22 -

 

 

0.21

 

0 0.02 0.04 0.06 0.08
Pat

Figure 3.12: Critical value of the scaled control parameter, pi, in a 1283 simple cubic
lattice with f), = 0.5, plotted as a function of p¢.

the percolation threshold pi does not exhibit a discontinuity when p¢ ——> 0. In the
Figure, we only show the interval p¢ E [0.0, 0.1], as, for p¢ = 0.1 the system is at the
critical value that is encountered for the regime in which both p._. and p¢ are non-zero

in the sense of Sec. 3.5.2, (p: 7é 0,p¢ 75 0): pi z 0.2488. In Fig. 3.12 we also show

83

a ﬁt to the data, which is of the form

1

C 2, —> 0 = —— t. .
The ﬁt parameters were found to be
u = 34.6 :1: 0.3 (3.74)
v = 1823 :l: 72 (3.75)
t = 0.2462 :1: 0.0003. (3.76)

It has to be noted that this formula only stands on empirical grounds, ﬁts by expo-

nential functions may also be useful.

3.6.2 p: —> 0 Limit

Again, for completeness, we also want to discuss the second possible limit where one
of the bond types is disabled. We thus set p: to some ﬁxed value p: E [0, 1] and then
subject p¢ to a variation with the nested interval method algorithm in order to ﬁnd
its critical value pi. Again, having determined pi, we calculate the corresponding
pi and then show pi versus the bond existence probability p=. This has been done
in Fig. 3.13. The curve does not look any different from the one displayed in Fig.
3.13 and a ﬁtting procedure with the function given in (3.73) reveals the following ﬁt

parameters, which are in very good accordance with (3.74).

u = 3374:03 (3.77)
v = 1744i73 (3.78)
t = 0.2461i0.0003 (3.79)

84

 

 

0.28 4 r
0.27 ~
0.26 1

0.25 - *

0.24 » 950%) “‘4‘“
1/(u+v*p=)+t

 

0.23 r

0.22 -

 

 

0.21 1 % . .
0 0.02 0.04 0.06 0.08
p=

 

Figure 3.13: Critical value of the scaled control parameter, pi, in a 1283 simple cubic
lattice with f), = 0.5, plotted as a function of pz.

3.6.3 Concentration Dependence in the p¢ = 0 Limit

In the preceding section we showed the behavior of the system for either one of the
bond probabilities going to zero. We observed that the critical value of the scaled
control variable, p+ under goes a continuous phase transition when letting p: —> 0
or p¢ —> 0. Until now, in this question of the limiting behavior we have only been
looking at a system with a ﬁxed concentration of blue species of f), = 0.5. We can
presume that for values of the blue species concentration fb aé 0.5 the general behavior
of a continuous transition will be conserved. This indeed is true, as simulation results
indicate, as shown in Fig. 3.14. Also shown is a ﬁt to the data, as in (3.73) of the

form

1

85

0.26 4

 

0.25 r

0.24 P

0.23 -

0.22 -

 

0.21 - pi (pat) r——l—«

1 /(u+v*p¢)+t

 

0.2 -

 

 

0.19

 

0 0.02 0.04 0.06 0.08
P¢

Figure 3.14: Critical value of the scaled control parameter, pi, in a 1283 simple cubic
lattice with f), = 0.3, plotted as a function of p¢.

where the ﬁt parameters were found to be

u = —17.640.7 (3.81)
v = —129.0414.6 (3.82)
t = 027440.003 (3.83)

in this case. It can be seen that the limiting value of pi (p¢ = 0) is not the same as
for f), = 0.5 and we suspect a concentration dependence of the threshold in the limit
13¢ —+ 0. We undertook simulations, again of simple cubic lattices of size L = 163,
L = 323, L = 643 and L = 1283, at different concentrations of the blue species,
determining the critical value of the bond density for p¢ = 0. The results shown in
Fig. 3.15 are ﬁts of the data obtained from the simulations on the four different

lattice sizes to the scaling law given in Sec. 2.5 in (2.89):

[pi(L) — pil oc L—l/V, (3.84)

86

 

0.28 T . . , >

0.27 - 9366.13.40) <>
a,(fb)*1/(h+m*fb) ——
0.26 -

0.25<
0.24 -
0.23 -

0.22 -

 

 

0.21 ‘ L ‘ 1
0 0.1 0.2 0.3 0.4 0.5

lb

 

Figure 3.15: Critical value of the scaled control parameter, pi, plotted as a function
of the fraction of blue sites, fb, in the limit p¢ = 0 and obtained from simulations on
simple cubic lattices of size L = 163, L = 323, L = 643 and L = 1283 by a ﬁt to the
scaling law |pi(L) — pi| or L‘l/V. The errors, estimated as described in [vdM97], are
smaller than the symbol sizes.

where by pi we mean to indicate the threshold in a lattice of inﬁnite extent. The

critical exponent V was kept ﬁxed at
1/ = 0.88. (3.85)

We already alluded to the problem of the estimation of systematic errors in the
context of ﬁnite size scaling. We will come back to a discussion of different estimation
techniques in Sec. 3.7 and for the moment content ourselves with the statement that
the estimates we performed led to errorbars smaller than the symbol size in Fig. 3.15.

Also shown in Fig. 3.15 is a ﬁt of the data with

2 _l2 l
Pi(fb,P¢=0)= (f;+:7:f:2-

 

(3.86)

This is Eq. (3.13) for the critical values with the purely empirical assumption of
a hyperbola for pi(fb), which is in agreement with the results of Heermann and

Stauffer for a one component site—bond model [HS81]. Fitting the parameters to our

87

simulation data results in

h = 400740.002 (3.87)

m = —4.42840.005 (3.88)

A comparison with the formula given by Heermann and Stauffer yields

 

 

1
h = C
pbond
= 4.019 (3.89)
1 _ pgond
m _
pﬁmd(f:ite - 1)
= —4.386, (3.90)
where pgmd = 0.2488 and site 2 0.3116 are the percolation thresholds for one com-

ponent bond- and site- percolation on a three dimensional simple cubic lattice, re-
spectively. However, in contrast to the one component model, due to the symmetry
introduced in the system, notably that p: E pm, 2 p", see (3.10) in Sec. 3.2.1 we
only have f), E [0.0, 0.5] as independent regime here, with the interval f), E [0.5,1.0]
being symmetric to the one shown here, with only the roles of blue and red sites
being switched. This is a manifestation of the two components behaving like two su-
perposed, non-interfering one component site-bond percolation systems. We further
illustrate this interpretation in the following. We presume that in the limit we are
discussing right now, p¢ = 0, there is no cross inﬂuence from the two species. We
check this presumption by taking the empirical formula as shown above or as given
in [H881] for a one component site-bond percolation model and plotting them for a
concentration f), and a concentration f, = (1 - f0). This corresponds to an isolated
system of blue sites with concentration f), and a second isolated system of red sites
with concentration f, = (1 - f0) and yields Fig. 3.16. For ease of comparison, we
express both concentrations in terms of one and thus for both curves do we use the

88

 

 

 

 

 

0.8 - 1
1/(h+m*fb) __
1/(h+m*(1-fb)) ----------- 4
0 II
D.
0 1 i 1 t
0 0.2 0.4 0.6 0.8 1

'0

Figure 3.16: Threshold p": versus concentration f), for two one-component site-bond
percolation models.

same variable fb. Following the presumption that both behave independently, we
just consider the two curves as superposed and obtain a phase diagram for the two
component system, by regarding the left axis as showing p: and not its critical value
pi. We ﬁrst realize that indeed we do reproduce a picture with a symmetry with

respect to f), = 0.5 and secondly that there are three phases:

1. Below both curves: No percolation.
2. Above one of the curves, below the second: One component percolates.

3. Above both curves: Both components percolate.

We can also transform the natural variable that we have used so far in the analysis,
p; to the scaled control variable pi. We show the phase diagram in terms of pi in
Fig. 3.17. The curves now are in the same shape as our simulation results. We can
summarize two facts: First the fact that indeed in the limit p¢ = 0 the two—component
site-bond percolation system may be described by two superposed, independent one

component site-bond percolation systems and second that the two-component site-

89

 

‘ a=(fb)*1/(h+m*fb)
013(1 -fb)'1/(h+m*(1 -1b)) - ........

    

0.8 -
0.6
0.4 -

0.2 P ............................ _

 

 

 

0 0.2 0.4 0.6 0.8 1
lb

Figure 3.17: Threshold pi versus concentration f), for two one-component site-bond
percolation models.

bond percolation system essentially behaves as a one-component bond-percolation
system when both bond types are active, see Sec. 3.4.1. These two facts let us
conclude that the phase transition in pi observed for p¢ —+ 0, see Sec. 3.6.1, may be
interpreted as an effective transition from a one-component bond percolation model

to a one—component site-bond percolation model.

3.6.4 The Cut-Out Technique

In addition to the conﬁrmation of the ﬁnding (3.86) through the results of Heermann
and Stauffer, we also performed another type of simulation to check on the decreasing
behavior of pi in the regime for small fb. It has to be stated that this method does
not represent a thorough veriﬁcation, however it does give a good conﬁrmation of the
results. We simulated a one-component bond percolation system on a simple cubic
lattice at some ﬁxed p: > pi, let the cluster structure evolve and then determined
the percolation probability. In a second step an algorithm swept through the whole

lattice, cutting out sites at random with a given probability pm), which was varied,

90

pm 6 [0,1], and determined the percolation probability again. In Fig. 3.18 we show
the percolation probability poo as a function of the probability that an arbitrary site

gets cut out, pm. A linear, decreasing dependence on the cut-probability is observed

0.38 1 r

 

  
   

pco(pcut) '—'A_‘

qw‘pcut

poo
.o u
0)

 

 

 

0 0.005 0.01 0.015 0.02
pcut

Figure 3.18: Percolation probability p00 versus the probability that a site gets ’cut-
01“,) pout-

in the ﬁrst regime, for pout ,2 0.012 a stronger decrease sets in. For completeness, we

state the values of the ﬁt parameters for the linear ﬁt

2900mm) = q + 0 pm, (3.91)

of the data also shown in Fig 3.18, which are

q = 036140.001 (3.92)

0 = —5.840.1. (3.93)

We can relate pm) to f), and through the observation of a decreasing p00 can make
plausible a linear rise of pi with pout or. In the framework of the two-component model
in the p7; = 0 limit, pout corresponds to fb. In the same regime however, we have
a quadratically decreasing function a:(fb), see Fig. 3.7. To arrive at the threshold

pi, we have to multiply pi with 01:, which then will be decreasing, dominated by

91

 

the quadratic decrease in 01:. It has to be noted that the argument is only valid
for relatively small values of the cutout probability pcut since at higher values the
system shows a rapid decrease to a vanishing percolation probability. It is at least an
intuitive, heuristic argument for explaining the seemingly counter-intuitive effect of a
decreasing threshold with an increasing number of sites of another species, which in
the limit p¢ act as geometrical obstructions for the percolating major concentration
species. We have to emphasize again that this is by far no thorough argument but
rather one of mere plausibility, giving however a good intuitive explanation for the

decrease of pi for small fb.

3.6.5 Concentration Dependence in the p: = 0 Limit

We now come to the question of the concentration dependence in the case where
all edges connecting sites occupied by the same species are closed always. In the
preceding sections we found that the percolation probability does not show a different
behavior whether p94 ——> 0 or p: —> 0. Here however, it will become clear, that contrary
to the percolation probability, the critical value of the bond density, pi, does show a
different behavior in the two different limits. The simulation is essentially the same
as described in Sec. 3.6.3, p: 4. 0 is set and with the nested interval method the
critical value of p73 is determined, from which pi is calculated. Again we undertook
simulations on simple cubic lattices of size L = 163, L = 323, L = 643 and L = 1283
and then ﬁtted the results with the scaling relation (3.84), providing us with the
’true’ threshold in the limit of the lattice size L —+ 00. The same discussion as led in
Sec. 3.6.5, concerning the systematic errors, applies and again we refer to Sec. 3.7
for a more complete discussion. We present the simulation results in Fig. 3.19. In
the upper curve of Fig. 3.19 we go back to the bond existence probability p¢. Its

critical value, as a function of the concentration, pc f), , is well reproduced with an
,4

92

 

1

 

p; (fb, p==0) <>

 

 

0.9 .
(a*exp(-d*fb)+6)
0.8 1
0.7 -
“at“
91+ 0.6 ‘
Q.
0.5 -

pi(fb.p.=0) 0
0'4 ’ Mb)*(a*exp(-d*fb)+c) """""

 

 

 

0'3 ' _G_0-6-9_9-9-9-e--ou43-49--o~o--o--o--o--e--€>
0.2 1 1 1 1 1 1 1
0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5

fb

Figure 3.19: Critical value of p¢ and of the scaled control parameter, pi and pi,
plotted as a function of the fraction of blue sites f), in the limit p= = 0. The data
has been obtained by a ﬁt to |pi(L) — pi] oc L‘l/V. Again the errors, estimated as
described in [vdM97], are smaller than the symbol sizes.

exponential ﬁt:

80 = P;(fb)

= a exp(—d fb) +c (3.94)

Fitting this empirical formula to the simulation data gives

a = 2.121: 0.07 (3.95)
d = 10.9 :1: 0.2 (3.96)
c = 0.547 :1: 0.002 (3.97)
We deﬁne
ff E w"(p¢=1) (3-98)
and numerically get
f; 2 0.14 (3.99)

93

11']

Si

from (3.94). We then also conduct an independent simulation, making use this time,
as described in Sec. 3.3, of the nested interval method to determine the critical value

of fb. We ﬁnd
ff 2 0.145 i 0.001, (3.100)

which again is a result from a ﬁt with data obtained on simple cubic lattices of size
L = 163, L = 323, L = 643 and L = 1283 to the scaling relation (3.84). What is the
signiﬁcance of the parameter ff? Let us ﬁrst consider the inverse case, where p74 = 0
and p: = 1. In this limit we effectively deal with a one-component site percolation
system, since one bond type is entirely ’switched off’ and the other edges are open
with probability one. This implies that we can decide to call one component the
inactive one, say the red species, and focus on the second species. If two blue sites
have one edge in common, it is open with probability one. This however, is exactly the
rule in a one-component site percolation system where the site percolation threshold

is known to be = 0.3116, (3.42), for a simple cubic lattice. In the exactly inverse

:88
case, (p: = 0, p¢ = 1), if two vertices occupied by different species have one edge in
common, it is open always. We thus conclude that ff can be regarded as a new site
percolation threshold in a simple site percolation system in which, unlike the ordinary
site percolation model, nearest neighbors only belong to the same cluster if they are
of opposite ﬂavor. In the lower curve of Fig. 3.19 we show the critical value of the

bond density, pi. It is simply given as the product of pi(fb) with the edge density
a¢(fb), see (3.13) and (3.14):

pi(fb,p==0) = a¢(fb)p;(fb)

= ($46—92) (aw—41.) +6 (3101)

Here, ff determines the critical point of the phase transition line pi(fb, _ = 0).

An argument against a novel behavior could be the statement that in this case, the

94

 

concentration of one species is no longer the appropriate parameter to characterize
the state of the system. To show that this argument does not hold for this purpose, we
propose one such new parameter that should be universal for the percolation systems
dealt with in this work. This is the density of accessible edges, which was already
introduced in (3.7), or, more speciﬁcally, for a two-component model in (3.14) and
(3.15). Although usually egdes are not included as ingredients in a site percolation
system, this quantity makes sense for a site percolation system too, as edges are
readily included in the line of argument of a site percolation system by setting them
open with probability one. For the model considered here the density of accessible

edges is given by a¢(fb) and at the new critical site percolation threshold ff we have

1

4:10;) = —,— — 201: — g)? (3.102)

For the ordinary site percolation system in which neighbors of the same ﬂavor form

clusters, we have the density of accessible edges given as

am = (fa-“)2, (3.103)
and the critical value

0:143 = (fine)? (3104)

Then, if both systems were equal with respect to this variable, the critical density of

accessible edges in both systems would be expected to be the same:
ai E. agm. (3.105)
This requires
C 1 1
fb : 5 ¥ §\/1 — 2 Oxide, (3.106)
giving two solutions, symmetric with respect to f), = 0.5,

f; = 0.051, (3.107)

95

T1

O-I'M-J. ..

 

 

an

3!

ll

and

f1? = 0949. (3.108)

This stands in obvious contradiction to our ﬁndings. We conclude that this is another

argument in favor of a qualitatively new behavior.

3.7 Finite Size Effects

In the preceding sections we alluded several times to the subject of error estimates.
The statistical errors are not the ones that are troublesome, how a reasonably precise
estimate of the systematic errors resulting from ﬁnite size effects can be achieved, is
a more complicated subject. We will ﬁrst focus on the one-dimensional case again

and then discuss brieﬂy techniques that we used.

3.7 .1 One-Dimensional Percolation Model

In Sec. 2.5 we derived the ﬁnite size scaling law for one dimensional percolation, see

(2.82) which read
pc — peﬁ oc L‘W (3.109)
In the derivation we stated that it would be necessary to fulﬁll the condition
I? — pol << 1-

for the scaling law in the approximative form, as derived, to be valid. We now want
to know, how big the systematic errors are when we apply the scaling law to lattices
of small linear dimension. Why are systematic errors to be expected at all? The
reason lies exactly in the fact that in the derivation in Sec. 2.5 the approximation

Ip - pol << 1 was used. This led to the relatively simple expression of

approx: _ 1 _

peﬂ

b1] 14

96

Had we not used this approximation, we would have obtained
pi?“ = exp(—1/L), (3.110)
resulting from
“(134217) = 1/6,

see (2.79). By relying on the approximative formula, however, the necessary condition
for its validity can not be satisﬁed any more, since for small L the system is not close

to the threshold. We therefore implicitly also have the condition

<< 1. (3.111)

ﬁll“

The systematic error entering the ﬁnite size scaling can be readily calculated,

exp(—1/L) — 1+ %

= f: ‘Lm. (3.112)

1
"2:2 71.

CTsys

 

 

We illustrate this effect in Fig. 3.20. In the next section we will go to a discussion
of the same phenomenon in higher dimensions, where it will become clear that the
analytical treatment presented here in the case of one dimension will not be possible

any more and one has to rely on estimation techniques.

3.7.2 Techniques in Higher Dimensions

In Sec. 2.5 we had conjectured a certain scaling behavior, (2.84) and from that derived

the scaling relation (2.89). In both equations, it was a requirement that
lpe6(L) - pol << 1 and L >> 1.

It is the second condition that one has to worry about. It is not a priori clear where

the linear dimension L of a lattice has to lie in order to be able to neglect the errors

97

 

 

 

 

Inc-931‘?“ = 1-exp(-1/L) —
-1
8 0'8 Pc-péﬁpm" = L -----------
$5, 0.6
0
D.
E: 0.4
00
9-
a:
0.2 ]
0

 

 

 

0 2 4 6 8 10 12 14 16
L

Figure 3.20: Exact (solid line) and approximative (dashed line) effective threshold in
a one-dimensional percolation system as a function of the linear lattice dimension L.

that are due to applying the ﬁnite size scaling equation in a regime where it only
holds in a very approximative way. We want to illustrate the situation that one faces
when applying the ﬁnite size scaling relation, (2.89), in Fig. 3.21, in which we show
the critical value of the bond density, pi, as a function of the linear lattice dimension,

as well as a ﬁt with
|pi(L) — pg| = a. L—l/V. (3.113)

What is the ’true’ value of pi? One way to proceed is to ﬁt the data with (3.113).
The true threshold is obtained as a ﬁt parameter, together with an error estimate.
Another possibility constitutes of plotting the data in a double logarithmic plot and

then ﬁtting with pi as ﬁt parameter such that the expectation of a straight line
C C 1
lnlp+(L) —p+| =lna— ZlnL, (3.114)

where a = const, is fulﬁlled. This would also provide an error estimate for the ﬁt
parameters. We applied the ﬁrst method, as shown in Fig. 3.21. This led to errorbars

smaller than the symbol size in Fig. 3.15. We have to note that this method does not
98

911-1

 

0.226

 

0.225 ~ p9, (L, tb=0.5, p¢=0) 1—1~
0.224 - "I'M” pi —
0.223 -
0.222 -
0.221 1

0.2 ~

 

 

 

 

0.219 4 . . . . .
0 20 40 60 80 100 120 140

L

Figure 3.21: Finite size scaling for pi(L,p¢ = 0) shown as a function of the linear
dimension of a simple cubic lattice at a concentration of the blue species of f), = 0.5.

explicitly provide for an estimate of the systematic errors themselves. This is, as the
only weighting of the data points used in the ﬁtting procedure comes in through the
statistical errors attached to each data point, no weighting that takes into account
the ’distance’ from the validity region of the scaling law is included. On the other
hand, the statistical errors do show a ﬁnite size scaling: The errorbars for smaller
lattice sizes turn out to be bigger than those for larger lattices, see Fig. 3.21, even
though the same number of events has been taken into account for calculating the
averages and their errors. Thus we already include a systematic ﬁnite size error in
the procedure. Yet this should not account for the second kind of systematic errors,
the ones due to the application of (3.113), which is only valid for L —-> 00. Is there
a way to estimate the errors solely due to the ﬁnite size scaling law? One way that
we can think of, is to try to ﬁnd a sort of a ﬁnite size weighting factor scaling for
the error inﬂuence of each data point in a ﬁnite size scaling procedure. This scaling

function could, for example be in the form of an exponential.

A search of the literature did not reveal any thorough discussion on an estimate

99

of the systematic errors in a determination of the ’true’ threshold pi(oo) due to ﬁnite
size scaling. We found two ways of dealing with the systematic error, where one
method tries to estimate them and the second tries to avoid them in the calculation

of the true threshold. We will present both methods.

The ﬁrst method was proposed in [vdM97]. It relies essentially on ﬁtting the data
with (3.113). The thresholds are calculated for different lattice sizes. Two ﬁts are
done, one including the last three data points, giving pi(Lmax, Lma$_1, Lmax_2) the
other including the last two data points resulting in pi(L,,m, Lmax_1). The difference

between the two ﬁt parameters is then taken as estimate for the systematic error,

0sys : |pi(Lma21Lmax—11Lmaz—2) — pi(Lmazn Lmaz—l)] - (3115)

We also performed calculations following this method to estimate the systematic
errors. Again, the errorbars obtained were smaller than the symbol size in Fig. 3.15.
A way to implement essentially this same method in what might be an improved way
would be to do the ﬁnite size scaling for all combinations of available data points,
not only ﬁtting the last two and the last three, and then taking the largest difference

between the resulting ﬁts for pi as estimate for the error am.

The second method is described in [BH97] In this case, not the percolation
probability p00 but rather the probability that a cluster spans the entire system, 1/2,

is taken to be the observable. For the inﬁnite lattice, we have

_ 0. if psi):
14—{1’ if P>Pi . (3.116)

As a consequence, all 2/2L(p+) have a common intersection point dam,“ (pi), which
we illustrate in Fig. 3.22. Extrapolation to Lmax —> 00 yields an estimate for the
percolation threshold pi. The advantage of this method is that it does not involve the

choice of a critical exponent, as 1/ in (3.113). Furthermore it does not involve a scaling

100

 

 

 

 

 

1 I I T :T"==ﬂ‘r:ﬁ:
(VI-rt "‘i"
‘p(L=20)I———+—1 f I ’1 i135
i
0.8 ~ ‘I‘(L=40) ----- x ..... g ' x I, 4
a ’ 3 if
\II(L=60)’ ----- l- ----- : a j: I};
:I 0.6 ~ \y(L=1oo) .W9-_, 35” 11* {if ‘
'+ J I
9’ l” I}
5" 0.4 - g?! 4
{1
i .
0.2 ~ ink?
iii
x33"
0:"" .. .- . . . .
0.23 0.235 0.24 0.245 0.25 0.255 0.26 0.265 0.27
9.

Figure 3.22: Probability w that a spanning cluster exists, versus the lattice size L.
The percolation threshold pi is given by the common intersection point 19L...“ (pi).

relation. However, still an extrapolation to L —) 00 is needed and the determination

of the percolation threshold itself is not without error.

3.7 .3 A Proposal for a Novel Technique

In Sec. 3.3.3 we described the approach that we took to determine, in a simulation, the
critical values pi, pi and ff. In summary, the simulation ﬁrst carried out the nested
interval method and thus determined the threshold for one particular random number
sequence, then the algorithm was started with another random number sequence and
after a sufﬁcient number of independent simulations, the average gave the critical
value. Another way to proceed, however, would be to just do the inverse: First
carry out many independent simulations, determine the percolation probability or the
spanning probability and then, decide, according to an appropriately chosen value a,
whether the system percolates or not. The parameter a could be chosen such as to

give an effective threshold of l/e, for example, as in (2.79), or, in general,

peﬂ: 212—101). (3.117)

101

 

 

Based on this averaged information, whether poo < a or p00 > a, the decision about
the nested interval method would be made and the new ’test’ value for the parameter
1):, p¢ or f), chosen, at which again a number of independent simulations would
yield the decision about the further direction of the parameter value. How does one
obtain a way to estimate the systematic errors when applying this method? The idea
is to apply the same procedure to several values of a and furthermore at different
lattice sizes L, where three different lattice sizes would be a minimum requirement.
Another requirement is the knowledge of an interval Api for which, from theoretical

considerations, we know deﬁnitely

pi denotes the ’real’ threshold in a system of inﬁnite extent. Through 1p this interval
Api gets mapped onto an interval of values 0, A0,. We would then get three different
curves: pi(Aa'i,L1), pi(Aa§,L2) and pi(Aa§,L3). In the following the ﬁnite size

scaling law should be applied to all combinations of points

{131(1’1,L1)1Pi(b2,L2)1Pi(bs,L3)} (3.119)
with the condition that
01 6 A01, 02 E Aug, 03 6 A0,?” (3.120)

where the case of a higher number of lattice sizes to take into consideration, works in

an analogous way. This eventually leads to a hyperplane
25.01.14.113. L —> 00). (3.121)

Now the minimum pi(bl,bg,b3,L —+ oo)‘,,,,-,, and maximum pi(b11b21b31L “* 00);.“

have to be determined, giving an estimate of the error,

Usys = |pi(bl,bg,b3, L —> 00y — pi(bl’ 02, 03, L —') my I, (3.122)

min max

102

 

 

and, as average taken over the hyperplane, the critical value pi. The interval

Apt-:1 : [pi-(b1) b21b39 L _) 00)::71z'n1pi-(b11b2) b3) L —) 00):. l) (3123)

man:

may be taken as new input in the procedure and in successive iterations the critical
value will be known with an error estimate. This method is based on ideas of von
Bergmann [vB99] and has evolved in several discussions with the author. It has to be
noted however that its practical aspects could not yet be tested and thus it remains
unclear whether it can live up to the expectations when applied to the determination

of the threshold and an estimate of its error.

3.8 Spin-Off Results

In this section, we will present some topics and related results that were not the main
interest in the line of research, but yet, in the pursuit of the main objectives, were
investigated in the line of work. We will not give a broad discussion of each of the

results but rather summarize in a concise manner the main aspects.

3.8.1 The Number of Open Clusters per Vertex

We already introduced the number of open clusters per vertex in Sec. 2.4.3. The
programs used determined this quantity that has also been called the f-functz'on and
we show a result for two different values of the concentration of the blue species in
Fig. 3.23. Also shown in the ﬁgure is the theoretical expectation, which we derive
in the following way, see [Bau]. The total number of open clusters per vertex has
to be calculated as given in (2.32). The ﬁrst terms of this sum can still be given

analytically, we consider the ﬁrst two:

N
“(19+) : N(1— 13+)z + 3310+“ " 19+)2Z—2 + - ~-
Z
= N(1— 5p.) + 002.2.)- (3.124)

103

 

 

 

 

 

1 I
K(p+) 1—4—4
0.8 - (1-(2/2)*p+) —
0.6 ~
’1
O.
37
0.4 r- J
0.2 ~ 4
o 1 1 - .
0 0.2 0.4 0.6 0.8 1
D.

Figure 3.23: The number of open clusters per vertex in a simple cubic 1283 lattice,
calculated for a concentration of the blue species f), = 0.5. Also shown is the approx-
imate theoretical expectation.

Here, N is the total number of vertices in the lattice,
N 2 L3,

and z is the coordination number of the lattice type, thus 8 = 6 in our case of the
simple cubic lattice. In Fig. 3.23 we see that indeed this approximation, which,
because of the expansion made in the derivation is valid for p+ << 1, does hold in this
regime. Finite size effects are not taken into account here. Thus we ﬁnd exactly what
is to be expected: For low values of pr the lattice is highly fragmented, showing a
high number of small clusters, whereas for high values of p+ gradually larger clusters

take over.

3.8.2 The Number of Spanning Clusters per Vertex

For the inﬁnite lattice, according to Theorem 1, there can almost always be only
one inﬁnite cluster present in the lattice and this will be the largest cluster. For

ﬁnite lattices however, a ﬁnite probability exists that more than one cluster spans

104

 

the entire system. Furthermore, also a ﬁnite probability exists for a smaller cluster
than the largest one to percolate. The question is, wether this affects the percolation
threshold in a ﬁnite lattice, and if it does, in which way. In an attempt to estimate
the inﬂuence of this ﬁnite size effect on the percolation probability poo, we introduced
a new quantity, which we called ppc and which represents the probability that an
arbitrary vertex of the lattice is part of any of the percolating clusters. If there is
only one percolating cluster, then we have poo and ppc equal. Thus, as we have the
theorem on the uniqueness of the inﬁnite cluster in the limit of an inﬁnite system, we

know that

lim poo = [iim ppc. (3.125)

L—)oo

We tested the new order parameter for several lattice sizes and focused on the question
of a difference in the percolation thresholds in terms of both variables. Within the
limits of the errorbars we could detect no different behavior, indicating that the
deﬁnition of the inﬁnite cluster in a ﬁnite system and the checking mechanisms for it,
see (2.10) and (2.11) are reasonable. We show one such comparison of poo and ppc in
Fig. 3.24. Analogously to ppc we can also deﬁne a ’new’ spanning probability, given
as the probability that any cluster, not necessarily the one that complies with the
requirements for being the inﬁnite cluster, Coo, percolates. We denote this quantity by
10,6. Similarly to the results shown above, no difference in the percolation thresholds
derived from 1,0 or 1pm, could be detected. With the theorem of the uniqueness of the

inﬁnite cluster we again have
232.3 = 332.1% (3126)

With these two results it is conﬁrmed that the percolation probability as used by us
in the framework of the research presented here in order to determine the percolation

threshold, is a sensible observable.

105

 

 

 

0.35 - T ' ' «)
ppC(p+) ’—.—‘ I

0.3 7
poo(p+) """ x """

0.25 -

.o
N
I

0.15 -

ppcm.). p-10.)

.0
—L

0.05 r 1

 

 

 

o 4: .4 l U 1 1
0.244 0.248 0.252 0.256 0.26
9.

Figure 3.24: Probability to belong to the inﬁnite cluster, poo(p+), and the probability
to belong to any percolating cluster, ppc(p+) in a simple cubic 1003 lattice at f), = 0.5.

3.8.3 Two-Dimensional Percolation

In Sec. 2.2.4 we discussed that in general critical exponents are expected to be
universal, meaning that they only depend on the dimension of the system, but not
its microscopic structure. The critical exponents for two dimensional ordinary bond
percolation problems are well known, see for example, [SA94]. Obviously the question
remains to be answered whether the results presented in this work for two-component
site-bond percolation on a three dimensional simple cubic lattice do in principle hold
on a two dimensional lattice too. We investigated this question in some aspects on
a square lattice. More speciﬁcally we focused on the question of the concentration
dependence of the threshold in the limits of one of the two bond existence probabilities
set to zero: p: = 0 or p¢ = 0, see 3.6.3 and 3.6.5. Before we proceed to state
the results we have to note that there is a fundamental difference concerning two-
component site-bond percolation on a simple cubic lattice and two-component site-
bond percolation on a square lattice. This difference lies in the fact that on the simple

cubic lattice a panchromatic regime does exist, which is not the case for the square

106

 

Ezra”- .

lattice. This comes about as the site percolation threshold on the square lattice is
noticeably higher than on the simple cubic lattice, leading to non-overlapping curves.
We illustrate our statement with Fig. 3.25 and also want to refer to Sec. 3.6.3, where
we led the corresponding discussion for the simple cubic lattice. Fig. 3.25 essentially

is Fig. 3.16 for the square lattice. We take again the assumption of a hyperbola, see

 

 

 

 

1 .
0.8 ~
0.6 ~
1/(h+m*fb) —
0-4 * 1/(h+m*(1-fb)) ----------- ‘
0.2 ~ 1
o 1 1 1 1
0 0.2 0.4 0.6 0.8 1

f0

Figure 3.25: Threshold pi versus concentration f), for two one-component site-bond
percolation models on the square lattice, compare Fig. 3.16.

(3.86),

1

PiUmpys = 0) = m

(3.127)

This is justiﬁed by simulations performed by us, results of which could very well be
reproduced by a hyperbola of the given form with the parameters as calculated in

(3.128). We calculate the parameters as in (3.89), with pgond = 0.5 (exact, see [SE64])

107

 

and c = 0.5927 for the square lattice yielding

site

 

 

1
h = C
pbond
= 2 (3.128)
m = c 1_fbond
pbond( site - 1)
= —-2.455. (3.129)

In contrast to the curves for the simple cubic lattice it becomes clear in Fig. 3.25
that for the square lattice, due to the higher site percolation threshold site, there are
only two phases present, no polychromatic regime is encountered where both species
percolate at the same time. Although this is no principle constraint in the sense
that yet mixed clusters are possible and thus also the limiting behavior of the system
when one of the bond existence probabilities is let to zero, should be investigable, we
did not further involve in research in two dimensions, aside from also checking the
cut-out technique, see 3.6.4, which yielded qualitatively the same results as on the
simple cubic lattice. The problem could be remedied by modeling the system on a

lattice of higher connectivity in two dimensions, where the site percolation thresholds

are lower, see [SA94] or [Gri99].

108

Chapter 4

Summary and Conclusion

Since its ﬁrst introduction in the context of polymer science, percolation theory has
evolved to an important ﬁeld, providing possibilities to study phase transitions in a
framework relatively easy to model. In percolation theory a somewhat unusual phase
transition is dealt with, since it does not involve any concept of temperature. Rather
all interesting functions are purely geometric properties, like for example, the perco-
lation probability or the spanning probability, and not thermal properties like speciﬁc
heats as functions of temperatures as dealt with in the realm of thermal phase tran-
sitions. Nonetheless there is quite a deep similarity existent between thermal phase
transitions and percolative phase transitions, which we already alluded to in Table
2.2. It also becomes clear when looking at the scaling laws that describe the asymp-
totic behavior of the respective systems close to their phase transition point. As in
thermal phase transition simple power laws govern the behavior in the vicinity of
the critical point and also do the critical exponents show the property of universal-
ity. This principle of universality suggests that one should be able to ﬁnd the same
critical exponents in real-life experiments as in computer experiments even though in
the real experiments impurities and defects in the lattice will alter and modify the
microscopic structure. As stated in [DZA83], metal-insulator ﬁlms do indeed show
the same behavior as found in computer experiments. Other possible applications of

109

percolation theory can be found again in [DZA83] or in [Sah94]. In addition, with
these correspondences between thermal and percolative phase transitions, including
universality of the critical exponents, an understanding of the percolative phase tran-
sitions, which are in some aspects easier to simulate and investigate, can help better

the understanding of more complicated thermal phase transitions.

In summarizing, this work presented, after a general introduction to percolation
theory and phase transitions in Chapter 2, research on a two-component site-bond
percolation system. We ﬁrst described a novel method to treat N-component per-
colation as described in [Baued] and [HB99]. This approach was then applied in
Monte-Carlo simulations to a percolation model on a simple cubic lattice of various
sizes for N = 2 components. We introduced a scaled control variable, p+, giving the
bond density in the system under observation, being a function of the concentration
of one of the species in the lattice, fb, and the two bond existence probabilities p___
and p91. Here the intuitive notation p: means to indicate the probability for an edge
to be open if the endvertices are of the same ’color’ whereas the notation p¢ stands
for the probability of an edge with endvertices of different color to be open. In the
simulations we determined the percolation probability p00 as a function of p+. We
found that for a regime where either one of the bond existence probabilities is set
to zero, p: = 0 or p¢ = 0, the system shows a different phase transition behavior
in the order parameter poo, see 3.5. We could furthermore establish that there is a
discontinuous phase transition taking place in the percolation probability poo whence
either p._. —-+ 0 or p¢ —> 0. As to the behavior of the percolation threshold pi in the
system we found that when both bond existence probabilities are non-zero, p: 76 0
and p91 ¢ 0, the critical value of the scaled control variable, pi, coincides with the
known percolation threshold in an ordinary one component bond-percolation system,

pgond x 0.2488, independent of the composition of the system characterized by fb. In

110

 

any one of the limits p: = 0 or 19;; = 0 however, a different percolation threshold pi
is found which furthermore depends on the concentration of one of the components:
pi(fb). As described in Sec. 3.6, the functional form turns out to be different in
both limits, we thus have pi(fb, p: = 0) and pi(fb, p¢ = 0). For the latter case we
could report that an empirical formula given by Heermann and Stauffer [H881] for a
one-component site-bond percolation system applies in the two-component site-bond
percolation model too. In the ﬁrst case, however, a different behavior was found, for
which we proposed a novel empirical formula. In a brief notation this is summarized
by stating that for p: p¢ ¢ 0 we ﬁnd a behavior similar to one—component bond
percolation, whereas for p: p¢ = 0 a different behavior is encountered. Aspects of

the simulation techniques and related results were also discussed.

The ﬁeld for future work in this area seems vast. As a ﬁrst generalization one might
try to apply the same method as described and applied here, to multi-component sys-
tems on lattices of higher dimension and/or higher connectivity. It should also be
interesting to check on this percolation model with a ’bond discrimination’ scheme
built in, differentiating between different types of bonds for the percolation process.
This approach should also ﬁnd a broad range of possible applications. One might
think of special networks, or gelation phenomena where several components are in-
volved, which interact in special ways, as well as wetting phenomena. It should be
possible to ﬁnd applications where either one of the limits discussed above or a model
with all types of bonds active, is applicable. Furthermore an application to stock—
market simulations seems possible and is being undertaken by the authors of [HB99].
In concluding, the present work and the works it is based on show that, although per-
colation theory has grown to a ﬁeld that much attention is devoted to, open questions,
aside from still lacking proofs of widely accepted conjectures, remain, the answers of

which have possibly useful applications to other ﬁelds of science.

111

 

 

Appendix A

Source Code

Following are the source codes of two of the programs used in the simulations. We
display these two as being representative of the two general techniques used, as de-
scribed in Sec. 3.3. All programs were written in FORTRAN and were compiled on
different machines with different FORTRAN-compilers. The random number generator
function listed at the end of the respective programs was taken from Ref. [PTVF92],
as the standard random number generators provided on the platforms used showed
correlations, especially for linear lattice dimensions L = 128. For an explanation of

the algorithms and concepts used in the program we refer to Sec. 3.3.

A.1 perc_prob_00 . f

The program displayed in the following is, as described in the section on simulation
techniques used, 3.3, the general purpose program aimed at determining the percola-
tion probability poo, the probability that a spanning cluster exists, 1,1), and the number
of open clusters per vertex.

program perc,prob00

t uses tortran 77 conventions
* uses a random number generator taken from numerical recipes in
t fortren 2nd ed.: ren1(iseed), modified to use double precision

implicit double precision (A-H,0-Z)

112

 

 

call readpar(n. nevent, prob_nHin, dprob-n, nprob_n)

 

"{1}

Loop over prob_n:

 

prob_n I prob_nHin
do iprob-n I O, nprob_n
call init(n, nevent, prob_n. iprob_n,
peqSt. dpeq, npeq, pneqSt, dpneq, npneq)
pneq I pneqSt
do ipneq I O. npneq
M ' paqSt
do ipeq I O, npeq
call Perc(n, nevent, prob-n, peq, pneq)
poq - M + dp-q
end do
pneq I pneq + dpneq
end do
prob-n I prob,n + dprob_n
end do

 

if.“

End Loop over prob-n:

 

1"

end

6”!

 

 

SUBROUTINE readpar(n, nevent. prob_nHin, dprob-n, nprob-n)
implicit double precision (A-H.0-Z)

parameter (maxn I 128)

input parameters:

open(unitI12.fileI’parall.dat’,statusI’old’)
read (12,*) n
I number of lattice sites in one direction
if (n .gt. maxn) then
print 4, ’ > Fatal Error: n must be less or equal maxn’
end it
read I, iseedO
iseed I iseedO
I seed for random number generator
read (12,.) nevent
I number of simulated events
read (12,!) prob_nMin, dprob_n. nprob_n
I probability that a given site is a neutron

END

SUBROUTINE init(n. nevent, prob-n. iprob-n,
peqSt, dpeq. npeq, pneqSt, dpneq. npneq)

implicit double precision (A-H.0-Z)

character lengthintot1, lintotl, nintoIB. probninfot4.
numinrot2, informationt40, pneqStintoI4

read (12.’(a)’) blankl
read (12.’(a)’) b1ank2
read (12.’(88,18.0)’) prob_n_info. prob_n_num

read (12,!) peqSt, dpeq. npeq
read (12.¢) pneqSt. dpneq, npneq

open(unitI13,fileI’runnumber.dat’,form=’tormatted’,status=’old’)
read(13,’(a)’) lengthinfo

read(13.’(i2.2)') icurrentnumber

close(13)

113

 

iprobn I prob_n I 1000

ipneqStI pneqSt I 1000
write(nin10.’(i3.3)’) n
write(probninto.’(i4.4)’) iprobn
write(pneqStinto,’(i4.4)’) ipneqSt
write(numinto,’(i2.2)’) icurrentnumber
write(linto,’(a)’) lengthinfo

open(unitI14,fileI’zlog_’//ninfo//’_’//probninfo//’-’
& //pneqStinto//'_’//linto//"//numinto//’.dat’,
& formI’formatted’,statusI’unknown’)
open(unitI16,file=’zpint_’//ninfo//’,’//probninfo//’_’
& //pneqStinfo//’_’//linfo//"//numinfo//’.dat’,
& formI’tormatted’,statusI’unknown’)
open(unitI17,fileI’zpperc-’//ninto//’_’llprobninfo/l’_’
& l/pneqStinfo//’_’//1into//"//numinto//’.dat’,
& formI’formatted’,statusI’unknown’)
open(unitI18,fileI’szunc-’l/ninfol/’-’//probninfo//’_’
& //pneqStinto//’-’//1in10//"//numinfo//’.dat’.
& tormI'rormatted’,statusI’unknown’)

if ((nevent-l) .eq. 0) then

write (14,I)’SHARNING:standard deviation of average gives
& division by 0’

write (16.*)’QHARNING:standard deviation of average gives
& division by 0’

write (17,I)’#HARNING:standard deviation of average gives
& division by 0’

write (18,I)’8HARNING:standard deviation of average gives
& division by 0’

endif

information I ’II calc. with perc-prob00 on It’
write (14,I) information

write (14,’(a20,i$.4)’) ’n I’,n
write (14,’(a20.e15.8)’)’prob_n I’.prob_n
write (14,’(a20,e15.8)’)’peqStart =’,peqSt
write (14.’(a20.e15.8)’)’dpeq I’.dpeq
write (14,’(a20.15.4)’) ’npeq I’,npeq
write (14.’(a20,e15.8)’)’pneqStart I’,pneqSt
write (14,’(a20,e15.8)’)’dpneq I’,dpneq
write (14.’(a20,i5.4)’) ’npneq =’,npneq
write (14,*)

write (14,I)

END

 

SUBRDUTINE Perc(n, nevent, prob_n. peq, pneq)

implicit double precision (A-H.0-Z)

parameter (maxn I 128)
parameter (maxnl I 129)
parameter (maxn2 I 16384)
parameter (maxn3 I 2097162)

logical neutron(maxn,maxn,maxn), maxperc

integer

c1usterNumber(0:maxn1,0:maxn1,0:maxn1)

double precision nn, np

do irun
call
call
call
call

call

I 1,nevent

rndinit(iseed)

lattice(n. prob-n, iseed, neutron, nn, np)
rndinit(iseed)

PercChk(n, prob_n. peq, pneq, irun, iseed. neutron.
ppbonds, pnbonds. marClusSize, multip, maxperc)
Statis(n. nevent. prob-n. peq, pneq, irun. nn, np,
ppbonds. pnbonds. maxClusSize. multip, maxperc)

114

 

 

and do

END

 

SUBROUTINE rndinit(iseed)
implicit double precision (A-B,0-Z)
integer iseed. iseedO, iseed1

iseedO I O

iseed I 0‘

call system-clock(iseed0,iclock2,iclock3)
iseed I - (iseedO)

iseed1 I iseed

zinitret I ran1(iseed1)
write(14,’(a20,i13.12)’)’randomseed I',iseed

write(14.I)
write(14,t)

END

 

 

SUBRDUTINE lattice(n. prob_n, iseed1. neutron,
& theneutrons, theprotons)

implicit double precision (A-H,0-Z)

parameter (maxn I 128)
logical neutron(maxn.maxn.maxn). trueFalse

Determine which sites are neutrons:
theneutrons I 0.00
theprotons I 0.00

do is I 1.n
do iy I 1.n
do ix I 1.n
trueFalse I .TRUE.
it (ran1(iseed1) .lt. prob-n) then
theneutrons I theneutrons + 1.00
else
trueFalse I .FALSE.
theprotons I theprotons + 1.00
end if
neutron(ir.iy,iz) I trueFalse
and do
and do
and do
write(14,’(a20,e15.8)’)’theneutrons I’,theneutrons
write(14,’(a20,e16.8)’)’theprotons I’,theprotons
END

 

SUBROUTINE PercChk(n, prob-n. peq, pneq, irun, iseed,
& neutron, ppbonds, pnbonds, maxClusSize, multip, maxbotIn)

implicit double precision (A-H,0-Z)
parameter (maxn I 128)

parameter (maxnl I 129)

parameter (maxn2 I 16384)
parameter (maxn3 I 2097162)

logical connected(maxn.maxn.maxn.3),
t neutron(maxn,marn,maxn).

115

--.JI

 

& trueFalse, maxbotIn, maxtopIn

integer nevPoint(maxn3,3).
t clusterNumber(0:maxn1.0zmaxn1,0zmaxnl)
integer iseed. iseed1. hcln
I Set random seed back to the initial value (ct Staufferbook p 73)

iseed1 I iseed
zinitret I ran1(iseed1)

n2 I n**2
n3 I nIt3
hcln I (n3 + 1)
I One more than highest possible clusterNumber if the whole lattice
I is fragmented
I Determine existing bonds
I a: x-bonds -----------

ppbonds I 0.00
pnbonds I 0.00
do iz I 1.n
do iy I 1.n
do ix I 1.n-1
trueFalse I .FALSE.
if (neutron(ix,iy,iz) .eqv. neutron(ix+1.iy,iz)) then
if (ran1(iseed1) .lt. peq) then
trueFelse I .TRUE.
ppbonds I ppbonds + 1.00
end if
else
if (ran1(iseed1) .lt. pneq) then
trueFalse I .TRUE.
pnbonds I pnbonds + 1.00

end it
end if
connected(ix,iy.iz,1) I trueFalse
end do
end do
end do
I b: y-bonds -----------

do is I 1.n
do iy I 1.n-1
do ix I 1.n
trueFalse I .FALSE.
if (neutron(ix,iy.i2) .eqv. neutron(ix.iy+1,iz)) then
if (ran1(iseed1) .lt. peq) then
trueFalse I .TRUE.
ppbonds I ppbonds + 1.00
end if
else
it (ran1(iseed1) .lt. pneq) then
trueFalse I .TRUE.
pnbonds I pnbonds + 1.00
end if
end if
connected(ix,iy,iz,2) I trueFalse
end do
end do
end do
I c: z-bonds -----------
do iz I 1.n-1
do iy I 1.n
do ix I 1.n
trueFalse I .FALSE.
it (neutron(ix.iy,iz) .eqv. neutron(ix,iy,iz+1)) then
it (ran1(iseed1) .lt. peq) then
trueFalse I .TRUE.
ppbonds I ppbonds + 1.00

116

 

ﬁlil’i’

end if
else
if (ran1(iseed1) .lt. pneq) then
trueFalse I .TRUE.
pnbonds I pnbonds + 1.00
end if
end if
connected(ix,iy,iz,3) I trueFalse
end do
end do
end do
Initialize connected values on boundaries as false
a: x-bonds -----------
do is I 1.n

do iy I 1.n
connected(n,iy,iz.1) I .FALSE.
enddo
enddo
b: y-bonds -----------

do is I 1.n
do ix I 1.n
connected(ix,n,iz.2) I .FALSE.

enddo
enddo
c: z-bonds -----------
do iy I 1.n
do ix I 1.n
connected(ix,iy,n,3) I .FALSE.
enddo
enddo

Assign the highest possible clusterNumber (hcln) to the surface
lattice points, so that they are treated as already belonging to a
cluster and are not dealt vith in the cluster find algorithm

do iy I 0,n+1
do ix I 0.n+1
c1usterNumber(ix.iy.0) I hcln
end do
end do
do iy I 0,n+1
do ix I 0.n+1
clusterNumber(ix,iy,n+1) I hcln
end do
end do

do iz I 0,n+1
do ix I 0,n+1
c1usterNumber(ix,0,iz) I hcln
end do
end do
do is I 0.n+1
do ix I 0.n+1
c1usterNumber(ix.n+1,iz) I hcln
and do
end do

do iz I 0,n+1
do iy I 0,n+1
clusterNumber(0.iy.iz) I hcln
end do
end do
do is I 0,n+1
do iy I 0.n+1
c1usterNumber(n*1,iy,iz) I hcln
end do
end do

117

 

 

ﬁ

toot-new cluster, start

i

1400

Assign the clusterNumber 0 (not part of any cluster) to all points
in the lattice that don’t belong to the surface

do is I 1.n
do iy I 1.n
do ix I 1.n
clusterNumber(ix.iy,iz) I 0
end do
end do
end do
CLUSTER FIND SUBROUTINE
multip I 0
maxClusSize I O

 

do iZI 1.n
do iy I 1.n
do ix I 1.n
IF(c1usterNumber(ix,iy,iz) .eq. 0) then

 

clusterNumber(ix.iy,iz) I multip + 1
newPoint(1,1) I ix

nevPoint(1,2) I iy

nevPoint(1,3) I iz

IHINO I 1

IMAXO I 1

continue

IMIN I IMAXO + 1

IMAX I IMIN

do I I IMINO,IHAXO
IIX I nevPoint(I.1)
IIY I nevPoint(I,2)
112 I nevPoint(I,3)

IF(clusterNumber(IIX-1.IIY.IIZ) .eq. 0) THEN
IF( connected(IIX-1,IIY.IIZ,1)) THEN

nevPoint(IMAX,1) I IIX-l
nevPoint(IHAX.2) I II?
neuPoint(IHAX.3) I IIZ
IMAX I IMAX + 1
clusterNumber(IIx-1,IIY,IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX+1,IIY,IIZ) .eq. 0) THEN
IF( connected(IIX.IIY,IIZ.1)) THEN

nevPoint(IMAX.1) I IIX+1
nevPoint(IHAx.2) I IIY
newPoint(IMAX,3) I III
IMAX I IMAX + 1
clusterNumber(IIX+1,IIY.IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX,IIY-1,IIZ) .eq. 0) THEN
IF( connected(IIX,IIY-1,IIZ,2)) THEN

nevPoint(INAX.1) I IIX
nevPoint(IMAX,2) I IIY-l
newPoint(IHAX,3) I IIZ
IMAX I IMAX 4* 1
clusterNumber(IIX,IIY-1,IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX,IIY+1,IIZ) .eq. 0) THEN
IF( connected(IIX,IIY,IIZ,2)) THEN

nevPoint(IHAX,1) I IIX
nevPoint(IMAX,2) I IIY+1
nevPoint(IMAX,3) I IIZ

118

 

 

IHAX I IHAX + 1
clusterNumber(IIX,IIY+1,IIZ) I multip+1
END IF
END IF
IF(clusterNumber(IIX,IIY,IIZ-1) .eq. 0) THEN
IF( connected(IIX,IIY.IIZ-1.3)) THEN

newPoint(IHAX,1) I IIX
newPoint(IHAX,2) I II?
newPoint(IHAx,3) I IIZ-1
IMAX I IHAX + 1
clusterNumber<IIX.IIY,IIZ-l) I multip+1
END IF
END IF

IF(clusterNumber(IIX,IIY,IIZ+1) .eq. 0) THEN
IF( connected(IIX.IIY,IIZ,3)) THEN

newPoint(IMAX,1) I IIX
newPoint(IHAx,2) I IIY
newPoint(IHAX,3) I II2+1
IMAX I IHAX + 1
clusterNumber(IIX,IIY,IIZ+1) I multip+1
END IF
END IF

end do

IF(IMAX .NE. ININ) THEN
IMINO I IHIN
IHAXO I IMAX - 1
GOTO 1400
ENDIF
; increase multiplicity
multip I multip + 1
; IHAXO contains cluster size
if (imaxo .gt. maxClusSize) then
maxClusSize I imaxo
indexHax I multip
end if

END IF
end do
end do
end do
END CLUSTER FIND SUBROUTINE
Check on Percolation of infinite cluster
maxtopIn I .FALSE.
maxbotIn I .FALSE.
do iy I 1.n
do ix I 1.n
if (clusterNumber(ix.iy,n) .eq. indexMax)
& maxtopIn I .TRUE.
end do
and do
if (maxtopIn) then
do iy I 1.n
do ix I 1.n
if (clusterNumber(ix.iy,1) .eq. indexHax)
& maxbotIn I .TRUE.
and do
end do
end if

pinf I dble(maxClusSize) / (dble(n)II3)

write (14,’(a20,15.4)’) ’irun I’.irun

write (14,’(a20.e15.8)’) ’peq I’,peq

write (14.’(a20,e15.8)’) ’pneq I’,pneq

write (14,’(a20.e15.8)’) ’ppbonds I’,ppbonds
write (14,’(a20,e15.8)’) ’pnbonds I’,pnbonds
write (14,’(a20,i13.12)’)’maxClSizeI’,maxClusSize

119

 

§§§§

write (14,’(a20,e15.8)’) ’pinf I’,pinf
write (14,’(a20.i13.12)’)'multip I’.multip
write (14,’(a20,12)’) ’maxbotIn I’.maxbotIn
write (14.I)

write (14,I)

END

 

SUBROUTINE Statis(n, nevent. prob_n. peq, pneq. irun, nn, np,
& ppbonds, pnbonds. maxClusSize, multip, maxperc)

implicit double precision (A-H.0-Z)

parameter (maxn I 128)
parameter (maxnl I 129)
parameter (maxn2 I 16384)
parameter (maxn3 I 2097152)

logical maxperc
double precision nn. np, anu. invanu

SAVE sumP, sumN, sumPPBonds. sumPNBonds. sumInfClus. sm2InfClus.
t sumPercP, sm2PercP, sumMult. sm2Hu1t

if (irun .eq. 1) then

cum? I 0.00
sumN I 0.00
ave? I 0.00
aveN I 0.00
sumPPBonds I 0.00
sumPNBonds I 0.00
sumInfClus I 0.00
sm2InfCLus I 0.00
aveInfCLus I 0.00
av2InfClus I 0.00
stdInfClus I 0.00
sumPercP I 0.00
sm2PercP I 0.00
avePercP I 0.00
av2PercP I 0.00
sthercP I 0.00
sumHult I 0.00
sm2Mult I 0.00
aveHult I 0.00
av2Mult I 0.00
sthult I 0.00

endif

Summation

Particles

sun? I sum? + up

sumN I sumN + nn

Bonds

sumPPBonds I sumPPBonds + ppbonds

sumPNBonds I sumPNBonds + pnbonds

Percolation

if (maxperc) sumInfClus I sumInfClus + dble(maxClusSize)
if (maxperc) sm2InfClus I sm2InfClus + (dble(maxClusSize))II2
if (maxperc) sumPercP I sumPercP + 1.00

if (maxperc) sm2PercP I sm2PercP + 1.00

sumMult I sumHult + dble(multip)

sm2Hult I sm2Hult + (dble(multip))II2

Averaging

if (irun .eq. nevent) then

120

 

if ((nevent-1) .eq. 0) then
srt I 0.10-20
else
srt I sqrt(db1e(nevent-1))
endif
I Particles
aveP I sumP / dble(nevent)
aveN I sumN / dble(nevent)
I Bonds
avePPBonds I sumPPBonds / dble(nevent)
avePNBonds I sumPNBonds / dble(nevent)
I Percolation
aveInfClus I (sumInfClus / (dble(n)**3)) / dble(nevent)
av2InfC1us I sm2InfClus / ((dble(n)II3)II2)
stdInfClus I sqrt(abs(av2InfClus/dble(nevent)-aveInfClus*I2))
& /srt
avePercP I sumPercP / dble(nevent)
av2PercP I sm2PercP
sthercP I sqrt(abs(av2PercP / dble(nevent) - avePercPII2))
& /srt
aveNult I sumHult / (dble(n)II3) / dble(nevent)
av2Hult I sm2Hult / ((dble(n)II3)II2)
sthult I sqrt(abs(av2ﬂult / dble(nevent) - aveHultII2))
l / srt

anu I 2.00Iprob_nI(1.DO-prob_n)
invanu I (1.00 - anu)
pplus I invanqueq + anquneq

§

Output
I LogFile
write(14.’(a20.e15.8)’)’anu I’,anu
write(14,’(a20,e15.8)’)’invanuI’,inwanu
write(14,’(a20,e16.8)’)’pplus I’.pp1us
write(14,’(a20.2e15.8)’)’pinf. stdpinf I’,aveInfCIus.stdInfClus
write(14.’(a20.2e15.8)’)’pperc. stdppercI’,avePercP. sthercP
write(14.’(a20.2e15.8)’)’ffunc, stdffunCI’,aveHult, sthult
I Files
write (16,’(5e15.8)’) peq, pneq, pplus.
& aveInfClus, stdInfClus
write (17,’(5e15.8)’) peq, pneq, pplus.
& avePercP. sthercP
write (18.’(5e15.8)’) peq, pneq, pplus,
t aveHult, sthult
endif
END

FUNCTION ran1(idum)
IMPLICIT DOUBLE PRECISION (A-H,U-Z)

INTEGER idum. IA, IN, IO. IR, NTAB, NDIV
DOUBLE PRECISION ran1. AH, EPS, RNMX
PARAMETER (IAI16807, IHI2147483647, AMI1./IH, IQI127773, IRI2836,
& NTABI32, NDIVI1+(IN-1)/NTAB, EPSI1.2d-7. RNHXI1.d0-EPS)
INTEGER j, k. iv(NTAB). iy
SAVE iv. iy
DATA iv INTABIO/. iy /0/
if (idum .le. 0 .or. iy .eq. 0) then
idumImax(-idum,1)
do jINTAB+8, 1, -1
indum/IQ
idumIIAI(idum-kIIO)-IRIk
if (idum .lt. 0) idumIidum+IH
if (j .le. NTAB) iv(j)Iidum
enddo
iniv(1)

121

 

 

endif

indum/IO
idumIIAI(idum-RIIQ)-IRIk

if (idum .lt. O) idumIidum+IM
jI1+iy/NDIV

iniv(j)

iv(j)Iidum
ran1Imin(AHIiy,RNMX)

return

END

A.2 perc_ppintv0477 . f

The second program, perc_ppintv0477 . f was written to accommodate the more
speciﬁc need of determining the critical value of the bond density, pi, through deter-
mination of pi. Similar versions aiming at the critical value of 19,5 and fb respectively
have also been used. Again, for a more thorough discussion of the programs features,
we refer to Sec. 3.3.

program perc-ppintv0477

uses fortran 77 conventions

varies prob_n, pp at fixed pm and determines ppc, pplusc

uses a random number generator taken from numerical recipes in
fortran 2nd ed.: ran1(iseed), modified to use double precision
Lattice population with neutrons has own subroutine: lattice

{iii}

implicit double precision (A-H.0-Z)

call readpar(n, nevent, prob-nHin, dprob_n. nprob_n, delta,
t zlambda)

 

Loop over prob_n:

*Iilr

 

prob_n I prob_nﬂin
do iprob_n I 0. nprob_n
call init(n, nevent. prob_n, iprob_n. delta. zlambda. apeqSt,
& bpeqSt, pneqSt, dpneq, npneq)
pneq I pneqSt
do ipneq I 0, npneq
call iter(n. nevent, prob-n, apeqSt, bpeqSt, pneq,
& delta. zlambda)
pneq I pneq + dpneq
end do
prob_n I prob-n + dprob_n
end do

if

 

End Loop over prob_n:

ﬂi

 

end

 

SUBROUTINE readpar(n. nevent, prob_nHin, dprob_n. nprob_n, delta,
& zlambda)

implicit double precision (A-H,O-Z)
parameter (maxn I 128)
I input parameters:

open(unitI12,fileI’parall.dat’,statusI’old’)

122

 

read (12,I) n
I number of lattice sites in one direction
if (n .gt. maxn) then
print I, ’ > Fatal Error: n must be less or equal maxn’
end if
read I. iseedO
iseed I iseedO
I seed for random number generator
read (12.I) nevent
I number of simulated events
read (12.I) prob_nHin, dprob_n. nprob_n
I probability that a given site is a neutron
read (12.I) delta
precision to be reached in the final interval step
read (12.I) zlambda
interval division

END

 

SUBROUTINE init(n, nevent, prob_n. iprob_n. delta. zlambda,
l apeqSt, bpeqSt. pneqSt. dpneq, npneq)

implicit double precision (A-H.O-Z)

character lengthinfoI1. linfoIl. ninfoI3, probninfoI4,
& numinfoI2, informationI40, pneqStinfoI4

 

read (12,’(a)’) blankl
read (12.’(a)’) blank2
read (12,’(8a.f8.0)’) prob_n_info, prob_n-num

read (12.I) apeqSt. bpeqSt
read (12.I) pneqSt, dpneq, npneq

open(unitI13,fileI’runnumber.dat’,formI’formatted’,statusI’old’)
read(13,’(a)’) lengthinfo

read(13,’(i2.2)’) icurrentnumber

close(13)

iprobn I prob_n I 1000

ipneqStI pneqSt I 1000
write(ninfo,’(i3.3)’) n
write(probninfo,’(i4.4)’) iprobn
write(pneqStinfo,’(i4.4)’) ipneqSt
write(numinfo,’(i2.2)’) icurrentnumber
write(linfo,’(a)’) lengthinfo

open(unitI14,fileI’zlog_’llninfo/l’_’//probninfo//’_’

t //pneqStinfo//’_’//linfo//"//numinfo//’.dat’,

t formI’formatted’,statusI’unknown’)
open(unitI16,fileI’zpplusclog_’//ninfo//’_’//probninfo//’-’

& //pneqStinfo//’_’//linfo//"l/numinfo/l’.dat’,

& formI’formatted’.statusI’unknown’)

if (iprob_n .eq. 0) then
open(unitI17,fileI’zpeqc_’llninfo/l’_’//probninfo//’_’

& //pneqStinfo//’_’//linfo//"llnuminfo/l’.dat’,
& formI’formatted’,statusI’unknown’)
open(unitI18,fileI’zpplusc-’//ninfo//’-’//probninfo//’_’
& //pneqStinfo//’-’//linfo//"l/numinfo/I’.dat’.
& formI’formatted’.statusI’unknown’)
endif

if ((nevent-1) .eq. 0) then
write (14,I)’8HARNING:standard deviation of average gives
& division by 0'
write (16.I)’8HARNING:standard deviation of average gives
t division by 0’
write (17,I)’8HARNING:standard deviation of average gives

123

 

t division by 0’
write (18,I)’8HARNING:standard deviation of average gives
& division by 0’
endif

information I ’II calc. with perc_ppintv0477 on II’
write (14.I) information

write (14,’(a20,i4.4)’) ’n I’,n

write (14,’(a20.e15.8)’)’delta I’,delta
write (14,'(a20.e15.8)’)’zlambda I’,zlambda
write (14,’(a20,e15.8)’)’prob_n I’,prob_n
write (14,’(a20.e15.8)’)’apeqStartI’,apeqSt
write (14.’(a20,e15.8)’)’bpeqStartI’,bpeqSt
write (14,’(a20.e16.8)’)’pneqStartI’,pneqSt
write (14,’(a20,e15.8)’)’dpneq I’.dpneq
write (14,’(a20,i4.4)’) ’npneq I’.npneq
write (14,I)

write (14,I)

END

 

SUBROUTINE iter(n, nevent, prob_n, apeqSt, bpeqSt, pneq,
& delta. zlambda)

implicit double precision (A-H,O-Z)

 

parameter (maxn I 128)
logical neutron(maxn.maxn,maxn)
integer perc, percount

double precision invanu

sumpeqc I 0.00
sm2peqc I 0.00
sumpplusc I 0.00
sm2pplusc I 0.00

anu I 2.00Iprob-nI(1.DO-prob_n)
invanu I (1.00 - anu)

write (14.’(a20.e15.8)’)’anu I’.anu

do irun I 1.nevent

call rndinit(iseed)

call lattice(n, prob_n, iseed, neutron)

call rndinit(iseed)

a I apeqSt

b I bpeqSt

percount I 0

do while ((b - a) .gt. delta)
c I zlambdaIb + (1-zlambda)Ia
call PercChk(n, prob_n. c, pneq. irun, iseed, neutron. perc)
if (perc .E0. 1) then

b I c
percount I percount + 1
else
a I c
endif
end do
peqc I c

write(14,’(a20,e15.8)’)’anu’,anu
write(14,’(a20,e15.8)’)’invanu’.invanu
write(14.’(a20,e15.8)’)’peqc’,peqc

pplusc I invanqueqc I anquneq
write(14,’(a20.e16.8)’)’pplusc’.pplusc

write (16,’(i4.4,3e15.8)’) irun.pneq,peqc.pplusc
sumpeqc I sumpeqc + peqc

sm2peqc I sm2peqc + (peqc)II2

sumpplusc I sumpplusc + pplusc

124

 

sm2pplusc I sm2pplusc + (pplusc)II2
end do

srt I sqrt(db1e(nevent-1))
if ((nevent-l) .eq. 0) then
srt I 0.10-20
endif
avpeqc I sumpeqc / dble(nevent)
stdpeqcI sqrt(abs(sm2peqc/dble(nevent)-avpeqcII2))/srt
avpplusc I sumpplusc / dble(nevent)
stdppluscI sqrt(abs(sm2pplusc/dble(nevent)-avppluscII2))/srt
deltapplus I invanuIdelta

write (17.'(i4.4,5e15.8)’) n,prob-n, pneq, avpeqc, delta. stdpeqc
write (18,’(i4.4.5e16.8)’) n,prob_n, pneq, avpplusc. deltapplus,
& stdpplusc

END

 

if

SUBROUTINE rndinit(iseed)
implicit double precision (A-H.0-Z)
integer iseed, iseedO. iseed1

iseedO I 0

iseed I 0

call system_clock(iseed0,iclock2.iclock3)
iseed I - (iseedO)

iseed1 I iseed

zinitret I ran1(iseed1)
write(14.’(a20,i13.12)’)’randomseed I’,iseed

Brito(14,*)
write(14.I)

END

 

SUBROUTINE lattice(n, prob_n. iseed1, neutron)
implicit double precision (A-H.O-Z)

parameter (maxn I 128)
logical neutron(maxn,maxn,maxn), trueFalse

Determine which sites are neutrons:
theneutrons I 0.00

theprotons I 0.00

sumNeutron I 0.00

sumProton I 0.00

do iz I 1.n
do iy I 1.n
do ix I 1.n
trueFalse I .TRUE.
if (ran1(iseed1) .lt. prob-n) then
theneutrons I theneutrons + 1.00
else
trueFalse I .FALSE.
theprotons I theprotons + 1.00
end if
neutron(ix.iy.iz) I trueFalse
end do
end do
end do
sumNeutron I sumNeutron + theneutrons

125

 

 

sumProton I sumProton + theprotons
write(14,’(a20.e15.8)’)’sumNeutron I ’,sumNeutron
write(14,’(a20,e15.8)’)’sumProton I ’,sumProton
END

 

SUBROUTINE PercChk(n, prob_n, peq, pneq, irun. iseed,
t neutron, perc)

implicit double precision (A-H.O-Z)

parameter (maxn I 128)
parameter (maxn1 I 129)
parameter (maxn2 I 16384)
parameter (maxn3 I 2097152)

logical connected(maxn.maxn.maxn,3),
& neutron(maxn,maxn,maxn).
& trueFalse
integer newPoint(maxn3.3),
& clusterNumber(0:maxn1.0:maxn1,0:maxn1)
integer iseed, perc, hcln
I Set random seed back to the initial value (cf Staufferbook p 73)

iseed1 I iseed
zinitret I ran1(iseed1)

n2 I nII2
n3 I nII3
hcln I (n3 + 1)
I One more than highest possible clusterNumber if the whole lattice
I is fragmented
I Determine existing bonds
I a: x-bonds -----------

ppbonds I 0.00
pnbonds I 0.00
do is I 1.n
do iy I 1.n
do ix I 1.n-1
trueFalse I .FALSE.
if (neutron(ix.iy.iz) .eqv. neutron(ix+1.iy.iz)) then
if (ran1(iseed1) .lt. peq) then
trueFalse I .TRUE.
ppbonds I ppbonds + 1.00
end if
else
if (ran1(iseed1) .lt. pneq) then
trueFalse I .TRUE.
pnbonds I pnbonds + 1.00
end if
end if
connected(ix,iy.iz,1) I trueFalse
end do
end do
end do
I b: y-bonds -----------
do 12 I 1.n
do iy I 1.n-1
do ix I 1.n
trueFalse I .FALSE.
if (neutron(ix,iy,iz) .eqv. neutron(ix.iy+1,iz)) then
if (ran1(iseed1) .lt. peq) then
truePalse I .TRUE.
ppbonds I ppbonds + 1.00
end if
else

126

if (ran1(iseed1) .lt. pneq) then
trueFalse I .TRUE.
pnbonds I pnbonds + 1.00
end if
end if
connected(ix,iy,iz,2) I trueFalse
end do
end do
end do
I c: z-bonds -----------
do iz I 1.n-1
do iy I 1.n
do ix I 1.n
trueFalse I .FALSE.
if (neutron(ix.iy.iz) .eqv. neutron(ix,iy.iz+1)) then
if (ran1(iseed1) .lt. peq) then
truePalse I .TRUE.
ppbonds I ppbonds + 1.00
end if
else
if (ran1(iseed1) .lt. pneq) then
truePalse I .TRUE.
pnbonds I pnbonds + 1.00
end if
end if
connected(ix,iy.iz.3) I trueFalse
end do
end do
end do
I Initialize connected values on boundaries as false
I a: x-bonds -----------
do is I 1.n
do iy I 1.n
connected(n.iy,iz,1) I .FALSE.
enddo
enddo
I b: y-bonds -----------
do is I 1.n
do ix I 1.n
connected(ix,n.iz,2) I .FALSE.
enddo
enddo
I c: z-bonds -----------
do iy I 1.n
do ix I 1.n
connected(ix.iy,n,3) I .FALSE.
enddo
enddo

Assign the highest possible clusterNumber (hcln) to the surface lattice
points, so that they are treated as already belonging to a
cluster and are not dealt with in the cluster find algorithm

§§§§§

do iy I 0,n+1
do ix I 0,n+1
clusterNumber(ix,iy,0) I hcln
end do
end do
do iy I 0,n+1
do ix I 0,n+1
clusterNumber(ix.iy.n+1) I hcln
end do
end do

do iz I 0,n+1
do ix I 0,n+1
clusterNumber(ix.0,iz) I hcln
end do

127

 

Giff

#

t¢¢¢-new cluster, start

*

1400

end do
do iz I 0,n+1
do ix I 0,n+1
clusterNumber(ix,n+1.iz) I hcln
end do
end do

do iz I 0,n+1
do iy I 0,n+1
clusterNumber(0,iy,iz) I hcln
end do
end do
do iz I 0,n+1
do iy I 0,n+1
clusterNumber(n+1,iy.iz) I hcln
end do
end do

Assign the clusterNumber 0 (not part of any cluster) to all points in
the lattice that don’t belong to the surface

do iz I 1.n
do iy I 1.n
do ix I 1.n
clusterNumber(ix.iy,iz) I 0

 

end do
end do
and do
CLUSTER FIND SUBROUTINE
perc I 0
multip I 0
ipclsize I 0

maxClusSize I 0

izI1
do iy I 1.n
do ix I 1.n
IP(clusterNumber(ix,iy.iz) .eq. 0) then

 

clusterNumber(ix,iy,iz) I multip + 1
newPoint(1,1) I ix

newPoint(1,2) I iy

newPoint(1.3) I is

IHINO I 1

IMAXO I 1

continue

IMIN I IMAXO + 1

IHAX I IMIN

do I I IMINO,IHAXO
IIX I newPoint(I.1)
IIY I newPoint(I,2)
112 I newPoint(I.3)

IF(clusterNumber(IIX-1,IIY,IIZ) .eq. 0) THEN
IP( connected(IIX-1.IIY,IIZ,1)) THEN

newPoint(IMAx,1) I Ill-1
newPoint(IHAX.2) I IIY
newPoint(IHAX,3) I IIZ
IMAX I IMAX + 1
clusterNumber(IIX-1,IIY,IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX+1,IIY,IIZ) .eq. 0) THEN
IF( connected(IIX,IIY.IIZ,1)) THEN
newPoint(IHAx.1) I IIX+1

128

 

newPoint(IHAX.2) I II?

newPoint(IHAx,3) I 112
IHAX I IHAX + 1
clusterNumber(IIX+1,IIY.IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX.IIY-l.IIZ) .eq. 0) THEN
IF( connected(IIx.IIY-1,IIZ.2)) THEN

newPoint(IHAX,1) I IIX
newPoint(IHAX.2) I IIY-1
newPoint(IHAX.3) I 112
IMAX I IHAX + 1
c1uaterNumber(IIX,IIY-1,IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX,IIY+1,IIZ) .eq. 0) THEN
IF( connected(IIX,IIY,IIZ,2)) THEN

newPoint(IHAX,1) I IIX
newPoint(IHAX,2) I IIY+1
newPoint(IHAX,3) I IIZ
IHAX I IHAX I 1
clusterNumber(IIX,IIY+1,IIZ) I multip+1
END IF
END IF

IF(clusterNumber(IIX,IIY,IIZ-l) .eq. 0) THEN
IF( connected(IIX,IIY,112-1.3)) THEN

newPoint(IHAX,1) I III
newPoint(IHAX,2) I IIY
newPoint(IHAx,3) I IIZ-l
IHAX I IHAX + 1
clusterNumber(IIX,IIY,IIZ-1) I multip+1
END IF
END IF

IF(clusterNumber(IIX,IIY.IIZ+1) .eq. 0) THEN
IF( connected(IIX,IIY,IIZ,3)) THEN

newPoint(IHAX,1) I IIX
newPoint(IHAX,2) I IIY
'newPoint(IHAX,3) I IIZ+1
IHAX I IMAX + 1

clusterNumber(IIX.IIY,IIZ+1) I multip+1

IF(IIZ+1 .EO. n) THEN
perc I 1
multip I multip+1
ipclsize I imaxo
write (14,’(a20.15.4)’) ’irun I’,irun
write (14,’(a20,e15.8)’) ’peq I’,peq
write (14,’(a20,e15.8)’) ’pneq I’,pneq
write (14.’(a20,e16.8)’) ’ppbonds I',ppbonds
write (14,’(a20,e15.8)’) ’pnbonds I’,pnbonds
write (14,’(a20.i13.12)’)’multip I’.multip
write (14,’(a20,i13.12)’)’ipclsize I’,

ipclsize

write (14,’(a20.i2.1)’) ’perc I’,perc
write (14.I)
write (14.I)
RETURN

END IF

END IF
END IF

end do

IF(IHAX .NE. ININ) THEN
ININO I IHIN
IHAXO I IHAX - 1
GOTO 1400
ENDIF
; increase multiplicity

129

multip I multip + 1
; IMAXO contains cluster size
if (imaxo .gt. maxClusSize) then
maxClusSize I imaxo
indexHax I multip

end if
END IF
end do
end do
END CLUSTER FIND SUBROUTINE
write (14,’(a20,i5.4)’) ’irun I’,irun
write (14.’(a20,e15.8)’) ’peq I’,peq

write (14.’(a20,e15.8)’) ’pneq I’,pneq
write (14,’(a20,e15.8)’) ’ppbonds I’,ppbonds
write (14.’(a20,e15.8)’) ’pnbonds I’,pnbonds
write (14.’(a20,i13.12)’)’multip I’,multip
write (14,’(a20,i13.12)’)’ipclsize I’,ipclsize
write (14,’(a20.i2.1)’) ’perc I’,perc
write (14,I)

write (14,*)

END

 

FUNCTION ran1(idum)
IMPLICIT DOUBLE PRECISION (A-H,O-Z)

INTEGER idum, IA, IM, IO, IR. NTAB. NDIV
DOUBLE PRECISION ran1, AH. EPS, RNHX
PARAMETER (IA-16807, IHI2147483647. AHI1./IH, IOI127773, IRI2836.
& NTABI32, NDIVI1+(IH-1)/NTAB, EPSI1.2d-7, RNHXI1.d0-EPS)
INTEGER j. k, iv(NTAB), iy
SAVE iv, iy
DATA iv /NTABIO/, iy /0/
if (idum .le. 0 .or. iy .eq. 0) then
idumImax(-idum,1)
do jINTAB+8, 1, -1
indum/IO
idumIIAI(idum-kIIO)-IRIk
if (idum .lt. 0) idumIidum+IH
if (j .le. NTAB) iv(j)Iidum
enddo
iniv(1)
endif
indum/IQ
idumIIAI(idum-kIIQ)-IRIk
if (idum .lt. 0) idumIidum+IM
jI1+iy/NDIV
iniv(j)
iv(j)=idum
ran1Imin(AHIiy,RNMX)
return

END

1130

Bibliography

[ARRS79]

[Bab97]
[Bau]

[Bau88]
[Baued]

[3395]
[BBDG87]

[B089]
[BDMP85]

[BF 94]
[BF95]
[BH57]
[BH97]
[Bin79]
[BPDM86]

[01397]
[0099]

P. Agrawal, S. Redner, P. J. Reynolds, and H. E. Stanley, J. Phys. C 12
(1979).

F. Babalievski, cond-mat/9711304 (1997).
W. Bauer, (private communication).
W. Bauer, Phys. Rev. C 38 (1988), 1297.

W. Bauer, Proceedings of the Workshop on Nuclear Matter in Different
Phases and Transitions, (Les Houches, France), 1999, to be published.

W. Bauer and A. Botvina, Phys. Rev. C 52 (1995), R1760.

W. Bauer, G. F. Bertsch, and S. Das Gupta, Phys. Rev. Lett. 58 (1987),
863.

P. Bak and K. Chen, Physica D 38 (1989), 5.

W. Bauer, D. R. Dean, U. Mosel, and U. Post, Phys. Lett. 150 B (1985),
53.

W. Bauer and W. A. Friedman, ??/nucl-th/9411012 ?? (1994), ??
W. Bauer and W. A. Friedman, Phys. Rev. Lett. 75 (1995), 767.

S. R. Broadbent and J. M. Hammersley, Percolation processes i. crystals
and mazes, Proceedings of the Cambridge Philosophical Society, vol. 53,
1957, pp. 629—641.

K. Binder and D. W. Heermann (eds), Monte carlo simulation in statis-
tical physics, an introduction, 3 rd ed., Springer Verlag, Berlin, 1997.

K. Binder (ed.), Monte carlo methods in statistical physics, Springer Ver-
lag, Heidelberg, 1979.

W. Bauer, U. Post, D. R. Dean, and U. Mosel, Nucl. Phys. A 452 (1986),
699.

R. Cont and J. P. Bouchaud, cond-mat/9712318 (1997).
J. Cardy and F. Colaiori, Phys. Rev. Lett. 82 (1999), 2232.
131

 

[osmg]

[Deb97]
[DG72a]

[DG72b]

[Dua90]
[Dup99]
[Dur88]

[DZA83]

[F E61]
[F ed88]
[FisGl]
[F1041]
[Gar95]

[G BP98]
[GGH92]
[G M96]
[GM97]
[Gri70]
[Gri89]
[Gri99]
[Ham57]

[Ham59]

[H899]

A. Coniglio, H. E. Stanley, and W. Klein, Phys. Rev. Lett. 42 (1979),
518.

J.-M. Debierre, Phys. Rev. Lett. 78 (1997), 3145.

C. Domb and M. S. Green (eds), Phase transitions and critical phenom-
ena, vol. 2, Academic Press, 1972.

C. Domb and M. S. Green (eds), Phase transitions and critical phenom-
ena, vol. 5a, Academic Press, 1972.

J. A. M. S. Duarte, Z. Phys. 80 (1990), 299.
B. Duplantier, cond-mat/ 9901008 (1999).

R. Durrett, Lecture notes on particle systems and percolation, Wadsworth
and Brooks/Cole, 1988.

G. Deutscher, R. Zallen, and J. Adler (eds), Percolation structures and
processes, Adam Hilger, Bristol, 1983.

M. E. Fisher and J. W. Essam, J. Math. Phys. 2 (1961), 609—619.
J. Feder, Fractals, Plenum Press, 1988.

M. E. Fisher, J. Math. Phys. 2 (1961), 620—627.

P. J. Flory, J. Am. Chem. Soc. 63 (1941), 3083.

C. Garrod, Statistical mechanics and thermodynamics, Oxford Univesity
Press, 1995.

T. Gharib, W. Bauer, and S. Pratt, Phys. Lett. B 444 (1998), 231.
S. Das Gupta, C. Gale, and K. Haglin, nucl-th/9210020 (1992).

S. Galam and A. Manger, Phys. Rev. E (1996).

S. Galam and A. Manger, Phys. Rev. E 56 (1997), 322.

R. B. Grifﬁths, Phys. Rev. Lett. 24 (1970), 1949.

G. Grimmett, Percolation, Springer Verlag, New York, 1989.

G. Grimmett, Percolation, 2nd ed., Springer Verlag, Berlin, 1999.

J. M. Hammersley, Percolation processes. lower bounds for the critical
probability, Annals of the Mathematical Statistics, vol. 28, 1957, pp. 790—
795.

J. M. Hammersley, Barnes supérieures de la probabilité critique dans
un processus de ﬁltration, Le Calcul des Probabilité et ses Applications,
CNRS, Paris, 1959, pp. 17—37.

H. M. Harreis and W. Bauer, cond-mat/9907330 (1999).
132

[HBM97]
[HH78]
[HK76]
[HKM78]

[HS81]
[HSB+79]

[Hua87]
[10395]
[KADS83]

[KBK97]

[Kes82]
[Kes87]

[Kre91]

[L298]
[N0193]

[N0198]
[oco+78]
[page]
[PTVF92]

[QBH76]
[SA94]

J. Hoshen, M. W. Berry, and K. S. Minser, Phys. Rev. E 56 (1997), 1455.
J. W. Halley and W. K. Holcomb, Phys. Rev. Lett. 40 (1978), 1670.
J. Hoshen and R. Kopelman, Phys. Rev. B 14 (1976), 3438.

J. Hoshen, R. Kopelman, and E. M. Monberg, J. Statist. Phys. 19 (1978),
219.

D. W. Heermann and D. Stauffer, Z. Phys. 44 (1981), 339.

J. Hoshen, D. Stauffer, G. H. Bishop, R. J. Harrison, and G. P. Quinn, J.
Phys. A 12 (1979).

K. Huang, Statistical mechanics, 2nd ed., John Wiley and Sons, 1987.
A. S. Ioselevich, Phys. Rev. Lett. 74 (1995), 1411.

A. Kapiltunik, A. Aharony, G. Deutscher, and D. Stauffer, J. Phys. A 16
(1983), L269—L274.

G. Kortemeyer, W. Bauer, and G. J. Kunde, Phys. Rev. C 55 (1997),
2730.

H. Kesten, Percolation theory for mathematicians, Birkhauser, 1982.

H. Kesten (ed.), Percolation theory and ergodic theory of inﬁnite particle
systems, Springer-Verlag, 1987.

U. Krengel, Einfu'hrung in die Wahrscheinlichkeitstheorie und Statistik,
3rd revised ed., Vieweg, 1991.

C. D. Lorenz and R. M. Ziﬂ', Phys. Rev. E 57 (1998), 230.

W. Nolting, Grundkurs T heoretische Physik, .4 Thermodynamik, 2 nd re-
vised ed., Verlag Zimmermann-Neufang, 1993.

W. Nolting, Grundkurs Theoretische Physik, 6 Statistische Physik, 3 rd
ed., Vieweg Verlag, 1998.

H. Ottavi, J. Clerc, G. Giraud, J. Roussenq, E. Guyon, and C. D. Mitescu,
J. Phys. C 11 (1978), 1311.

R. K. Pathria, Statistical mechanics, 2nd ed., Butterworth-Heinemann,
1996.

W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery,
Numerical recipes in fortran, 2nd ed., Cambridge University Press, 1992.

G. P. Quinn, G. H. Bishop, and R. J. Harrison, J. Phys. A 9 (1976), L 9.

D. Stauffer and A. Aharony, Introduction to percolation theory, revised
2nd ed., Taylor and Francis, London, 1994.

133

[Sah94]

[SE64]
[Sin81]

[Sta71]

[Sta79]
[Sti89]
[Sto43]
[vB99]
[vdM97]
[vdM98]
[WP78]
[Ye092]

[Zal77]
[Za183]
[Zif92]

M. Sahimi, Applications of percolation theory, Taylor and Francis, Lon-
don,1994.

M. F. Sykes and J. W. Essam, J. Math. Phys. 5 (1964), 1117.

YA. G. Sinai, Theory of phase transitions: Rigorous results, Pergamon
Press, 1981.

H. E. Stanley, Introduction to phase transitions and critical phenomena,
Oxford University Press, 1971.

D. Stauffer, Phys. Rep. 54 (1979), 3.

K. Stierstadt, Physik der Materie, VCH Weinheim, 1989.

W. H. Stockmayer, J. Chem. Phys. 11 (1943), 45.

J. von Bergmann, (private communication), 1999.

S. C. van der Marck, Phys. Rev. E 55 (1997), 1514.

S. C. van der Marck, Intl. J. Mod. Phys. C 9 (1998), 529.

H. J. Wintle and P. H. Puhach, J. Statist. Phys. 18 (1978), 557.

J. M. Yeomans, Statistical mechanics of phase transitions, Clarendon
Press, Oxford, 1992.

R. Zallen, Phys. Rev. B 16 (1977), 1426.
R. Zallen, The physics of amorphous solids, John Wiley, New York, 1983.
R. M. Ziff, Phys. Rev. Lett. 69 (1992), 2670.

134

"Illll‘lllllllllllll