Emergence and Complexity in Theoretical Models of Self-Organized Criticality

The concept of self-organized criticality (SOC) was introduced by Bak Tang and Wiesenfeld in 1987 [btw], to explain the abundant fractal structures in nature, e.g. mountain ranges, river networks, power law tails in the distribution of earthquake intensities etc. The SOC refers to the non-equilibrium steady state of slowly driven systems, which show irregular burst like relaxations with a wide distribution of event sizes. The power law correlations of different physical quantities, extending over a wide range of length and time scales is a signature of criticality. Usually, reaching a critical state requires fine tuning of some control parameters e.g. temperature and magnetic field for the Ising model with a given interaction strength. However in SOC the systems reach a critical state under their own dynamics, irrespective of the initial states and without any obvious fine tuning of parameters.

In the last two decades a large amount of study is focused on understanding the mechanism of SOC. The questions regarding the universality classes of the critical states has still not been completely settled. Many theoretical models have been studied to address these issues. Most of these are cellular automata models with discrete or continuous variables, evolving under deterministic or stochastic evolution rules (see [dharphysica06] for a review). Among them the Abelian Sandpile Model (ASM) is studied the most, mainly because of its analytical tractability using the Abelian property[dharprl].

The standard ASM first proposed in [btw] is defined on a lattice with height variables $z_{i}$ at each site $i$, which is equal to the number of sand grains at that site. There is a threshold value $z_{c}$ for each site, and any site with height $z_{i}\geq z_{c}$ is said to be unstable. The system in a stable configuration is driven by adding a single grain at a randomly chosen site. If this addition makes the system unstable, it relaxes by the following toppling rule: All the unstable sites at one time step transfers one grain each to all its nearest neighbors. Grains can move out of the lattice by toppling at the boundary. When the system reaches a stable configuration it is again driven by adding a grain and the process is repeated. The model reaches a steady state, in which the probability distribution of the size of events has power law tail.

The ASM defined on an infinite lattice when driven by adding grains only at a single site and relaxed, produces beautiful and complex patterns in height variables (see figure $1.1$). They are one of the examples where complexity arises from simple rules. In section 1.2 we analyze some of these patterns and develop a detailed and exact mathematical characterization of them.

There are different variants of the ASM that have been proposed. We study two well known models among them.

First is the Zhang model. This is similar to the ASM, but with continuous non-negative variables, usually referred as energy. At any time step all the unstable sites relax by equally distributing all its energy amongst its nearest neighbors, with their energy reducing to zero. Energy can also move out of the system by toppling at the boundary. The driving is done by adding energy to a randomly chosen site and the amount of the energy is chosen at random from a distribution.

The second is a stochastic variant of the ASM. The first stochastic sandpile model was proposed by Manna and it is known as Manna model[manna]. The model is non-Abelian, but one can construct stochastic relaxation rules with Abelian character. We consider one such stochastic Abelian sandpile model introduced in [dm]. The model is similar to the ASM with non-negative integer height variables $z_{i}$ and a threshold value $z_{c}$ defined at each site. The driving is also done by adding one sand grain at a randomly chosen site in a stable configuration. The difference is in the relaxation rules: On toppling $z_{c}$ number of grains are transfered, each grain moving independent of others to the nearest neighbors with equal probability and the height at the toppling site reduces by $z_{c}$. For the one dimensional model defined on a linear chain, with $z_{c}=2$, there are three possible events in a toppling at site $i$: Both the neighbors $\left(i-1\right)$ and $\left(i+1\right)$ gets one grain each. Probability of this event is $1/2$. Other two possibilities are that both the grains move either to the left or to the right neighbor, each event with probability $1/4$.

While real sand, poured at one point on a flat substrate, produces a rather simple conical pyramid shape, nontrivial patterns are generated this way in the ASM on an infinite lattice. One such pattern on a square lattice with threshold height $z_{c}=4$, produced by adding grains at the origin in an initial uniform distribution of heights $z=2$, is shown in the figure $1.1$.

The reason for interest in these patterns is two fold.

Firstly, these are analytically tractable examples of complex patterns that are obtained from simple deterministic evolution rules. Here complexity means that we have structures with variations, and a complete description of which is long. Thus, a living organism is complex because it has many different working parts, each formed by variations in the working out of the same, but relatively much simpler genetic coding.

Secondly, these patterns have the very interesting property of proportionate growth. This is a well-known feature of biological growth in animals, where different parts of the growing animal grow at roughly the same rate, keeping their shape almost the same. Our interest in studying the sandpile patterns comes from these being the simplest model of proportionate growth with non-trivial patterns. Compare the two patterns in figure $1.1$ produced on the same background but with different values of $N$. The pattern grows in size and finer features become discernible at the center, but the overall shape of the pattern remains same. Most of the other growth models studied in physics literature, such as the Eden model, the diffusion limited aggregation, or the surface deposition, do not show this property [eden, dla, barabasi]. In these models, the growth is confined to some active outer region. The inner structures, once formed are frozen in and do not evolve further in time.

The standard square lattice produces complicated patterns and it has not been possible to characterize them so far. We consider a pattern which is simpler but still complex. The pattern is produced on the F-lattice which is a variant of the square lattice with directed edges. The lattice with checker board distribution of heights is shown in figure 1.2. Each site has equal number of incoming and outgoing arrows. The threshold height is $2$ and any unstable site relaxes by giving away one particle each in the direction of the outgoing arrows. The pattern produced by centrally seeding grains on the checkerboard background is shown in the figure 1.3.

The Zhang model in one dimension has the remarkable property that in spite of the randomness in the amount of energy added during driving, the steady state energy per site has a very sharply peaked distribution in which the width of the peak is much less than the spread in the input amount. One such distribution is shown in figure 1.7. The threshold energy $E_{c}=1.5$ and the driving energy is chosen from a uniform distribution in the range $\left[0.76:1.24\right]$. The distribution has a spike at $E=0$ and a peak at $E=1.0$ with a standard deviation $\sigma=0.0135$. In general the width of the peak decreases with the increase of the system size.

This behavior was noticed in numerical simulations in both one and two dimensions, and it is called as “emergence of quasi-units”. It is argued that for large systems, the behavior would be the same as in the discrete model [zhang]. Fey et. al. [fey] have shown that for some choices of the distribution of input energy, in one dimension, the variance of energy does go to zero as length of the chain $L$ goes to infinity. However they did not show how fast the variance decreases with $L$.

We study the emergent behavior in the one-dimensional model by looking at how the added energy is redistributed among different sites in the relaxation process. Let the amount of energy for driving at time $t$ be chosen randomly from a uniform interval $\left[1-\epsilon,1+\epsilon\right]$. The time is counted in terms of the relaxation steps and at one time step all the unstable sites relax together. We write the amount of driving energy at time $t$ as

$$\Delta_{t}=1+\epsilon u_{t},$$

(1.8)

where $u_{t}$ is uniformly distributed in the interval $\left[-1,1\right]$. We decompose the energy variable at a site $x$ in a relaxation time step $t$ as

$$E\left(x,t\right)=Nint[E\left(x,t\right)]+\epsilon\eta\left(x,t\right)$$

(1.9)

where $Nint$ refers to the nearest integer value. The integer part evolves as the integer heights in the ASM. The function $\eta\left(x,t\right)$ is independent of $\epsilon$ and is a linear function of $u_{t}$. The precise function depends on the evolution history $H_{t}$ which is determined by the initial configuration and the sequence of addition sites $d_{t}$ up to time $t$. We assume that at starting time $t=0$, the variable $\eta\left(x,t=0\right)$ are zero for all $x$, and the initial configuration is a recurrent configuration of the ASM. We define $G\left(x,t|d_{t^{\prime}},t^{\prime},H_{t}\right)$ by

$$\eta\left(x,t|\left\{u_{t}\right\},H_{t}\right)=\sum_{t^{\prime}=1}^{t}G\left(x,t|d_{t^{\prime}},t^{\prime},H_{t}\right)u_{t^{\prime}}.$$

(1.10)

We show that this is equal to the probability $Prob\left(x,t|d_{t^{\prime}},t^{\prime},H_{t}\right)$ of a marked grain in the corresponding ASM, added at the site $d_{t^{\prime}}$ at time $t^{\prime}$ following the history $H_{t}$, to be found at site $x$ at time $t$

$$G\left(x,t|a_{t^{\prime}},t^{\prime},H_{t}\right)=Prob\left(x,t|a_{t^{\prime}},t^{\prime},H_{t}\right).$$

(1.11)

Then we show that the variance of energy in the steady state at a site $x$ can be written as

$$Var\left[E\left(x\right)\right]=L/(L+1)^{2}+\epsilon^{2}\Sigma^{2}\left(x\right),$$

(1.12)

where

$$\Sigma^{2}\left(x\right)=\lim_{t\rightarrow\infty}Var\left[\eta\left(x,t\right)\right]=\frac{1}{3}\sum_{\tau=0}^{\infty}\frac{1}{L}\lim_{t^{\prime}\rightarrow\infty}\sum_{x^{\prime}}\overline{G^{2}(x,t^{\prime}+\tau|x^{\prime},t^{\prime},H_{t})}.$$

(1.13)

The overbar denotes averaging over the evolution histories $H_{t}$.

The function $\overline{G}$ has been studied in [punya], but for $\overline{G^{2}}$ the calculation is much more difficult. However analyzing the behavior in different limits we introduce a phenomenological expression for it and show that in the large $L$ limit the variance of energy at a site $x$ has a scaling form $L^{-1}g\left(x/L\right)$. We determine an approximate form of the scaling function

$$g\left(\xi\right)=A\ln\left(1+\frac{1}{B\sqrt{\xi\left(1-\xi\right)}}\right),$$

(1.14)

where $A$ and $B$ are numerical constants. This expression agrees very well with the results of our numerical simulation.

The sandpile models with stochastic toppling rules are important subclass of SOC models. These models are able to describe the avalanche behavior seen experimentally in the piles of granular media much better than the deterministic models [osloexpt]. Also in the numerical studies, one gets better scaling collapse, and consequently, more reliable estimates for the values of the critical exponents, than for the models with deterministic toppling rules [stochnum]. There is good numerical evidence that these models constitute a universality class different from their deterministic counter parts [univ1, univ2, numbtw, numsand1, numsand2].

Unfortunately, at present, the theoretical understanding of the models with stochastic toppling rules is much less than the deterministic models, even the characterization of the steady state is not known for the one-dimensional Manna model.

We consider the Manna model with Abelian toppling rule on a linear chain of length $L$, described in section 1.1. Any stable state of the model is expressed as a linear combination of the stable height configurations with the coefficients being the probability of finding the system in that configuration. We define an addition operator $\mathbf{a}_{i}$ corresponding to the addition of grains at the site $i$, which acting on a stable state takes it to another stable state achieved following the stochastic relaxation rule.

The Abelian property of the operators is shown in [dm]. However, unlike the deterministic ASM, the inverse operator $\left\{\mathbf{a}_{i}^{-1}\right\}$ need not exist. This makes the determination of the matrix form of the operators difficult for this model. In a straightforward exact numerical calculation of the steady state one needs to invert a matrix of size $4^{L}\times 4^{L}$. We use the operator algebra of the addition operators to obtain an efficient method which requires inverting a matrix only of size $2^{L}\times 2^{L}$.

Using the conservation of sand grains during toppling we show that

$$\mathbf{a}_{i}^{2}=\frac{1}{4}\left(\mathbf{a}_{i-1}+\mathbf{a}_{i+1}\right)^{2}.$$

(1.15)

In general the operators $\mathbf{a}_{i}$ need not be diagonalizable. However, using the Abelian property we construct a common set of generalized eigenvectors for all the operators such that in this basis the matrices simultaneously reduce to the Jordan block form. Corresponding to each block there is at least one common eigenvector. Then from eq. (1.15), the eigenvalues satisfy

$$a_{i}^{2}=\frac{1}{4}\left(a_{i-1}+a_{i+1}\right)^{2},$$

(1.16)

with the boundary condition $a_{0}=a_{L+1}=1$. We reduce this set of coupled quadratic equations to a set of linear equations by taking square root

$$\eta_{i}a_{i}=\frac{1}{2}\left(a_{i-1}+a_{i+1}\right),$$

(1.17)

where $\eta_{i}=\pm 1$. There are $2^{L}$ different choices for the set of $L$ different $\eta$’s and for each such choice, we get a set of eigenvalues $\left\{a_{i}\right\}$. In general there will be degenerate sets of eigenvalues. We show that degeneracies are possible only for $L=3\left(\rm{mod~}4\right)$ and the set can at most be doubly degenerate. This implies that the largest dimension of a Jordan block is $2$. We determine the matrix elements inside each block in terms of the solution of a set of coupled linear equations.

We define a transformation matrix between this generalized eigenvector basis and the height configuration basis.

$$|\left\{z_{i}\right\}\rangle=\sum_{j}\mathbf{M}_{\left\{z_{i}\right\},j}|\psi_{j}\rangle,$$

(1.18)

where $|\left\{z_{i}\right\}\rangle$ is the basis vector corresponding to the height configuration $\left\{z_{i}\right\}$ and $|\psi_{j}\rangle$ is the $j$ th generalized eigenvector. Any height configuration can be generated by an appropriate sequence addition operators acting on the all empty configuration. As the action of the addition operators on the generalized eigenvectors are known, this determines the elements of the transformation matrix $\mathbf{M}$.

Given $\mathbf{M}$ we can get the eigenvectors of the addition operators in the configuration basis, in particular, the steady state vector, by the inverse transformation

$$|\psi_{j}\rangle=\mathbf{M}^{-1}|\left\{z_{i}\right\}\rangle.$$

(1.19)

We numerically calculate the inverse matrix $\mathbf{M}^{-1}$ and determine the exact steady state for systems of small sizes. This result is then extrapolated to determine the asymptotic density profile in the steady state. Our results suggest that the steady state density averaged over all sites approaches the asymptotic density as

$$\frac{1}{\rho_{L}}=\frac{1}{\rho_{\infty}}+\frac{B}{\left(L+\delta,\right)^{\nu}}$$

(1.20)

with $\rho_{\infty}=0.953$ which is close to the Monte Carlo estimate $0.949$. We also find that the ratio of probabilities of the most probable to the least probable configuration varies as $\exp{\left(0.94L\log L\right)}$. We show that the steady state is not a product measure state.

The method described is easily generalized to other stochastic Abelian sandpile models and we discuss some examples of them.

One of the most striking aspects of physics is the simplicity of its laws. Maxwell’s equations, Schrodinger’s equation, and Hamiltonian mechanics are simple and expressible in few lines. However every place we look, outside the textbook examples, we see a world of amazing complexity: huge mountain ranges, scale free coastlines, the delicate ridges on the surface of sand dunes, the interdependencies of financial markets, the diverse ecologies formed by living organisms are few examples. Each situation is highly organized and distinctive, but extremely complex. So why, if the basic laws are simple, is the world so complicated? The idea of Self Organized Criticality was born aiming to give an explanation for this ubiquitous complexity [jensen]. In this chapter the basic concepts related to SOC, that will be important for this thesis, are introduced.

The examples, cited above, share a common feature: a power-law tail of the correlations. Consider the two point correlation of a quantity $\Delta h\left(\mathbf{x}\right)=h\left(\mathbf{x}\right)-\bar{h}$, where $h\left(\mathbf{x}\right)$ is the height at a place $\mathbf{x}$ in a mountain range, and $\bar{h}$ is its average value. The function $\langle\Delta h(\mathbf{x+r})\Delta h(\mathbf{x})\rangle$ increases as $r^{\delta}$, with the exponent $\delta$ varying very little for different mountain ranges. Similar distribution with extended tails is observed in many other natural phenomena: Gutenberg-Richter law in earth quake [gutenberg], Levy distribution in stock market price variations [levy], Hacks law in River networks [hackslaw1, hackslaw2] etc. Such power-law distributions entail scale invariance — there are no macroscopic spatial scales other than the system size, in terms of which one can describe the system, making it complex.

Such features are familiar to physicist in equilibrium systems undergoing phase transition. In standard critical phenomena there are control parameters such as temperatures, magnetic field, which requires to be fine tuned by an external agent, to reach the critical point. This is unlikely to happen in naturally occurring processes such as formation of mountain ranges, earth quakes or even stock markets. These are non-equilibrium systems brought to their present states, by their intrinsic dynamics — and not by a delicate selection of temperature, pressure or similar control parameters. ¹¹1Per Bak, in his book [blind], puts this in an interesting comment—“The nature is operated by a ’blind watchmaker’ who is unable to make continuous fine adjustments”.

In the summer of 1987, Bak, Tang and Wiesenfeld(BTW) published a paper [btw] proposing an explanation to such ubiquitous scale invariance. They argued that the dynamic which gives rise to the robust power-law correlations seen in the non-equilibrium steady states in nature must not involve any fine tuning of parameters. It must be such that the systems under their natural evolution are driven to a state at the boundary between the stable and unstable states. Such a state then shows long range spatio-temporal fluctuations similar to those in equilibrium critical phenomena. The complex features appear spontaneously due to a cooperative behavior between the components of the system. They called this self-organized criticality as the system self-organizes to its critical steady state.

SOC nicely compliments the idea of chaos. In the latter, dynamical systems with a few degrees of freedom, say as little as three, can display extremely complicated behavior. However, a statistical description of this randomness is predictable in the sense that, the signals have a white noise spectrum, and not a power law tail. A Chaotic system has little memory of the past, and it is easy to give a statistical description of such behavior. In short, chaos does not explain complexity. On the other hand, in SOC, generally, we start with systems of many degrees of freedom, and find a few general features which are also statistically predictable, but has a power-law spectra leading to complex behavior. In certain dynamical systems, e.g., logistic maps, there are points (the Feigenbaum point [chaos]) in the parameter space, which separates states with a predictable periodic behavior and chaos. At this transition point there is complex behavior, with power-law correlations. SOC gives description of how systems, under their own dynamics, without external monitoring, reaches this very special point (“edge of chaos”), explaining the robust complex behavior in natural systems.

In the book “How nature works?”, Per Bak gives various kinds of natural examples of SOC, of which the canonical one is the sandpile. On slowly adding grains of sand to an empty table, a pile will grow until its slope becomes critical and avalanches start spilling over the sides. If the slope is too small, each grain just stays at the place where it lands or creates a small avalanche. One can understand the motion of each grain in terms the local properties, like place, the neighborhood around it etc. As the process continues, the slope of the pile become steeper and steeper. If the slope becomes too large, a large catastrophic avalanche is likely, and the slope will reduce. Eventually, the slope reaches a critical value where there are avalanches of all sizes. At this point, the system is far out of balance, and its behavior can no longer be understood in terms of the behavior of localized events. The system is invariably driven towards its critical state.

In order to have a mathematical formulation of SOC, BTW studied a so-called cellular automata known as the sandpile model [btw], which is discrete in space, time and in its dynamical variables. The model is defined on a two dimensional square lattice where each site $i$ has a state variable $z_{i}$ referred as height, which takes only positive integer values. This integer can be thought of as representing the amount of sand at that location or in another sense it represents the slope of the sandpile at that point. Neither of these analogies is fully accurate, the model has aspects of both. One should consider it as a mathematical model of SOC, rather than an accurate model of physical sand.

A set of local dynamical rules defines the evolution of the model: At each time step a site is picked randomly, and its height $z_{i}$ is increased by unity. In this thesis, this step will be referred as the driving. If the height now is greater than or equal to a threshold value $z_{c}=4$, the site is said to be unstable. It relaxes by toppling whereby four sand grains leave the site, and each of the four neighboring sites gets one grain. If there are any unstable sites remaining, they too are toppled, all in parallel. In case of toppling at a site at the boundary of the lattice, grains falling outside the lattice are removed from the system. This process continues until all sites are stable. This completes one time step. Then, another site is picked randomly, its height is increased by $1$, and so on.

The following example illustrates the dynamics. Let the lattice size be $3\times 3$ and suppose at some time step the following configuration is reached where all sites are stable.

2 3 2 3 3 0 1 2 3

We now add a grain of sand at randomly selected site: let us say the central site is chosen. Then the configuration becomes the following

2 3 2 3 4 0 1 2 3

The central site is not stable, and therefore it will topple and the configuration becomes

401133.\begin{tabular}[]{|c|c|c|}\hline\cr 2&{\color[rgb]{1,0,0}4}&2\\ \hline\cr{\color[rgb]{1,0,0}4}&0&1\\ \hline\cr 1&3&3\\ \hline\cr\end{tabular}.

		2	4	2

This configuration has two unstable sites, so both will topple in parallel. Since these are at the boundary, two grains will be lost, on toppling. The new result is

021233,\begin{tabular}[]{|c|c|c|}\hline\cr{\color[rgb]{1,0,0}4}&0&3\\ \hline\cr 0&2&1\\ \hline\cr 2&3&3\\ \hline\cr\end{tabular},

		4	0	3

and further toppling leads to

0 1 3 1 2 1 2 3 3

This is a configuration with all sites stable. One speaks in this case of an avalanche of size $s=4$, since there are four topplings. Another measure is the number of steps required for relaxation, which in this case is $t=3$. For large lattices, in the steady state, the distribution of avalanche sizes and durations display a long power-law tail, with an eventual cutoff determined by the finite size of the system.

Since the original sandpile model by BTW a large number of variations of the model have been studied (see [dharphysica06, jensen] for reviews). These are mostly extended systems with many components, which under steady drive reaches a steady state where there are irregular burst like relaxations and long ranged spatio temporal correlations. It is to be noted that in these models the complexity is not contained in the evolution rules itself, but rather emerges as a result of the repeated local interactions among different variables in the extended system.

In the rest of this chapter, I will introduce three of the most studied models of sandpile and the techniques used to analyze them.

Since the work by BTW, a large number of different models have been studied e.g. sandpile models with many variations of the BTW toppling matrix [kadanoff] or sand grain distribution rules [maslov], stochastic topplings [manna], with activity inhibition [manna2], continuous height models [zhang], loop erased random walk [loop], Takayasu aggregation model [takayasu], train model [train1, train2], non-abelian sandpile directed sandpile model [lee, pan, threshold, pruessner], forest-fire model [forest], Olami-Feder-Christensen model [ofc] etc (and many more could have been defined). Most of these models could only be studied numerically, and for a while it seemed that each new variation studied belong to a new universality class each with its own set of critical exponents. It is a fair question to ask, what is the generic behavior?

Although this question is not yet resolved completely, by now, there has been a fair amount of understanding of this problem. The universality classes with renormalization group flow in these models can be summarized in the Fig. 3.2.

There are sufficient numerical evidence that sandpile models with deterministic toppling rules (DASM) and stochastic toppling rules (Manna model) constitute different universality classes. There are also several other model which show critical exponents different from these two [sneppen1, sneppen2, grass]. They are related to the directed percolation (DP) [dp], which describes the active-absorbing state phase transition in a wide class of reaction-diffusion systems. The activity in avalanches in sandpile can grow, diffuse or die, and any stable configuration is an absorbing state. Thus one would expect that in general the sandpile should belong to the universality class of active-absorbing state transition with many absorbing states [path]. However, these models do not involve any conserved fields. In Manna and DASM-type models, it is this presence of conservation laws of sand that makes the critical behavior different from DP [vespignani].

Recently, the effect of non-conservation has been explicitly studied [mohanty1, mohanty2] by introducing stickiness in the toppling rules (i.e. there is small probability that the incoming particles to a site get stuck there, and do not cause any toppling until the next avalanche hits the site, thus in effect there is no conservation of grains within an avalanche). It has been argued that as long as the sand grains have non-zero stickiness, the distribution of avalanche sizes follows directed percolation exponents. The DASM, and the stochastic Manna-type models are unstable to this perturbation, and the renormalization group flows are directed away from them to the directed percolation fixed point, as schematically shown in the figure. This picture is exactly verified in directed sandpiles. However, the argument is less convincing for undirected models, and the issues is not settled [bonachela2].

Soon after the sandpile model was introduced, several experimental groups measured the size distribution of avalanches in granular materials. Unfortunately, real sandpile do not seem to behave as the the theoretical models. Experiments show large periodic avalanches separated by quiescent states with only limited activity. While for small piles one could try to fit the avalanche distribution with power-law over a limited range, the behavior would eventually cross over, on increasing the system size, to a state which is not scale-invariant [sandexpt1, sandexpt2]. It is later realized that inertia of rolling grains is the reason for non-criticality. A large avalanche propagating over a surface with slope $\theta_{c}$ scours the surface, and takes away materials from it. The final angle, after the avalanche has stopped, is bellow $\theta_{c}$. So if we want to see power-law avalanches, we have to minimize the effect of inertia of the grains. This is achieved in an experiment in Oslo by using rice grains. Because of the elongated shape of the rice grains (Fig. 3.3) frictional forces are stronger and these poured at very small rate gave rise to a convincing power-law avalanche distribution [osloexpt].

A similar power law distribution of avalanche sizes are obtained in motion of domain walls in ferro-magnets [feromag1, feromag2] and flux lines in type II superconductors [sup1, sup2]. A more recent experimental realization of SOC is obtained in suspensions of sedimenting non-brownian particles by slow periodic shear [corte].

It is worth mentioning that SOC has been invoked in several other situations in geophysics (atmospheric precipitation [presp], river pattern due to erosion [erotion], landslides [land]), biology (neural-network [nunet]), economics (stock-market crashes [smcr]) and many more. I have deliberately chosen the above experimental examples for which experimental observations are accurate and reproducible.

Originally, SOC was proposed with the aim of providing an explanation of the ubiquitous complexity in nature [btw]. The abundance of fractal structures in nature, temporal as well as spatial, was considered to be an effect of a generic tendency — pertinent to most many-body systems — to develop by themselves in to a critical scale-invariant state.

However, certainly not all systems that organize themselves into one specific state will, when gently driven, exhibit scale invariance in that self-organized state. The experiments of real sandpiles referred earlier are a prime example. Neither is all observed power law behavior are an effect of dynamical self-organization into a critical state. The work by Sethna and co-workers on Barkhausen noise [sethna] is an interesting example of this, what Didier Sornette has called “Power laws by sweeping of an instability” [sornet].

Since the introduction of the idea, a large amount of discussion went into understanding the minimum conditions for a model to be self-organized critical. Though a broad picture has emerged in last two decades, it is still not complete and controversial. In the rest of this section, I will discuss two of the most discussed properties, using both examples and counter examples.

• •

  Slow driving limit. There is a strong belief in the community that an essential ingredient of SOC is slow driving (driving and dynamics operating at two infinitely separated timescales, i.e. avalanches are instantaneous relative to the time scale of driving). This idea got widely accepted after an argument given by Dickman et. al..[path] They argued that the dynamics in the sandpile model implicitly involve tuning of the density of grains to a value where a phase transition takes place between an active state, where topplings take place, and a quiescent “absorbing” state. When the system is quiescent, addition of new particles increases the density. When the system is active, particles are lost to the sinks via toppling, decreasing the density. The slow driving ensures that these two density changing mechanisms balance one another out, driving the system to the threshold of instability.

  However in the Takayasu model of aggregation [takayasu] the driving is fast. A simple example, it can be defined on a linear chain on which particles are continuously injected, diffuse and coalesce. One can write the explicite rules as follows:

  • –

    At each time step, each particle in the system moves by a single step, to the left or to the right, taken with equal probability, independent of the choice at other sites.

  • –

    A single particle is added at every site at each time step.

  • –

    If there are more than one particle at one site, they coalesce and become a single particle whose mass is the sum of the masses of the coalescing parts. In all subsequent events, the composite particles acts as a single particle.

  The probability distribution of total mass at a randomly chosen site, has power law tail, with an upper cutoff that increases with time. This is a signature of criticality. The analogue of avalanches in this model is the event of coming together of large masses. In fact, it can be shown [dharphysica06] that the model is equivalent to a directed version of sandpile. In this example, it is clear that the driving is fast, and the rate is comparable to the local movements of the particles.

• •

  Conservation. Conservation of grains is also considered as a key property for the criticality to emerge in sandpile models. A simple intuitive argument goes as follows: the sand grains introduced in the pile can dissipate only by reaching the diffusive “boundary” of the lattice. Owing to this and because of the vanishing rate of sand addition, arbitrarily large avalanches (of all possible sizes) must exist for an arbitrarily large system size, yielding a power-law size distribution. In contrast, in the presence of non-vanishing bulk dissipation, grains disappear at some finite rate, and avalanches stop after some characteristic size determined by the dissipation rate. This clearly says that bulk dissipation is a relevant perturbation in the sandpile dynamics and breaks the criticality [bonachela1].

  There are also some other models of SOC like Forest fire [forest], OFC model [ofc] where it was shown, mostly numerically, that non-conservation in the dynamics leads to non-critical steady state.

  However, extrapolating these results and considering conservation as a neccesarry criteria for SOC, in general, is not correct. A model which is clearly non-conservative and still, when slowly driven displays power-law in the avalanche size distribution is discussed in [ice]. Another two non-conservative models of SOC are a sandpile model with threshold dissipation[thresdis], and Bak-Sneppen model of evolution [sneppen2].

Finally, one could ask: Has SOC, taught us anything about the world that we did not know prior to it? Jensen addresses this question very nicely in his book. The most important lesson is that, in a great variety of systems, particularly for slowly-driven-interaction- dominated-threshold systems, it is misleading to neglect fluctuations. In these systems, sometimes, the fluctuations are so large that the fate of a major part of the system can be determined by a single burst of activity. Dinosaurs may have become extinct simply as a result of an intrinsic fluctuation in a system consisting of a highly interconnected and interacting web of species; there may be no need for an explanation in terms of external bombardment by meteorites. Fluctuations are so large that the "atypical" events decides the future of the system. This new insight is sufficiently important to justify and inspire more theoretical, and experimental research in SOC.

The work in this thesis ranges from characterizing the spatial patterns in sandpile model, to quasi-units in the stationary distribution of Zhang’s model, and determining exact steady state of Manna model. The first three chapters in the following are about sandpile as a growth model, where we show how well-structured non-trivial patterns emerge at large length scales, due to local interactions in cellular lengths where the patterns are not obvious. In chapter 7 we discuss another emergent behavior in the Zhang model. Chapter 8 contains an operator algebraic analysis of the stochastic sandpile models. Here is a brief summary of these chapters.

• Chapter 4: While a considerable amount of research went into characterizing the universality classes of sandpile models and understanding the mechanism of SOC, very limited work is done about spatial patterns in sandpile models. Such patterns were noted around the time when sandpile was first introduced [liu]. Yet, very little is known about them.

  This chapter is devoted to the study of a class of such spatial patterns produced by adding sand at a single site on an infinite lattice with initial periodic distribution of grains and then relaxing using the DASM toppling rules. We present a complete quantitative characterization of one such patterns. We show that the spatial distances in the asymptotic (in the limit when large number of grains are added) patterns produced by adding a large number of grains, can be expressed in terms solution of discrete Laplace’s equation (discrete holomorphic functions [duffin, mercat, laszlo]) on a grid on two-sheeted Riemann surface.

  We also discuss the importance of these patterns as a paradigmatic model of growth where different parts of the structure grow in proportion to each other, keeping their shape the same. We call these kind of growth as proportionate growth. We discuss the importance of such growth in real world examples.

• Chapter 5: In this chapter we describe how the pattern changes in presence of absorbing sites, reaching which the grains get lost and no longer participate in the avalanches. We show that, again, the asymptotic pattern can be characterized in terms of discrete holomorphic functions, but on a different lattice. Similar effects of multiple sites of addition on the pattern are also calculated.

  The most interesting effect of the sink sites is the change in the rate of growth of the pattern. In absence of sink sites the diameter of the pattern, suitably defined, increases as $\sqrt{N}$ where $N$ is the number of sand grains added in the lattice. When the pattern grows with the sink sites inside, the growth rate of the diameter changes, in general, to $N^{\alpha}$, where the exponent $\alpha$ depends on the sink geometries. For example, $\alpha=1/3$, when the sink sites are along an infinite line adjacent to the site where grains are dropped. When the site of addition is inside a wedge of angle $\pi/2$ with the sink sites along the wedge boundary, this value of the exponent is $1/4$. We use an scaling argument and determine $\alpha$, for some simple sink-geometries.

• Chapter 6: The growth rate also depends on the arrangement of heights in the background, and this dependence is quite intriguing. When the initial heights are low enough at all sites, one gets patterns with $\alpha={1/d}$, in $d$-dimension. If sites with maximum stable height $\left(z^{c}-1\right)$ in the starting configuration form an infinite cluster, we get avalanches that do not stop, and the model is not well-defined. In this chapter, we study backgrounds in two dimensions. We describe our unexpected finding of an interesting class of backgrounds, that show an intermediate behavior: For any $N$, the avalanches are finite, but the diameter of the pattern increases as $N^{\alpha}$, for large $N$, with $1/2<\alpha\leq 1$, the exact value of $\alpha$ depending on the background. It still shows proportionate growth. We characterize the asymptotic pattern exactly for one illustrative example.

• Chapter 7: As mentioned, the Zhang model on one and two dimensional lattices displays a remarkable property: emergence of quasi-units, where the continuous heights, in spite of the randomness in the driving, are peaked around a few discrete non-random values. Fey et. al have shown that on a linear chain the width of the distribution vanishes in the infinite volume limit. However they did not show how it approaches zero.

  In this chapter, we show that, the sequence of toppling of the continuous height variables, when suitably discretized, have an one-to-one relation with that of integer heights in the corresponding DASM. We use this relation to show that the width of the distribution of heights decreases in inverse power of the length of the chain. We also determine how the variance of height at a site, changes with position of the site along the length of the system.

• Chapter 8: This chapter contains an algebraic approach of determining the steady state of a class of sandpile models with stochastic toppling rules. The original Manna model, as discussed in section 1.2.3, does not have the abelian property of its deterministic counterpart. However, a simple modification of the toppling rules makes the model abelian [dm]. A similar construction is possible for other stochastic toppling rules. However, analysis of these models are still difficult as the corresponding addition operators (see section 3.2.1), in general, does not have an inverse, and are not diagonalizable.

  We show that, in principle, the operators can be reduced to a Jordan block form, using the algebra satisfied by these. These are then used to determine the steady state of the models. We illustrate this procedure by explicitly determining the numerically exact steady for a stochastic model on a linear chain. Using the desktop computers at our disposal, we have been able to perform the calculation for systems of size $\leq 12$ and also studied the density profile in the steady state.

As discussed in details in chapter 3, the sandpile models were introduced in physics in the context of self-organized criticality, where the main interest has been the power-law tail in the distribution of avalanche sizes [btw]. In this chapter our interest is different. We study the pattern produced by adding large number of sand grains at a single site in a two dimensional DASM starting from a periodic background, and allowing the system to relax using the sandpile toppling rule (See chapter 1). For example, consider the ASM on an infinite square lattice with initial heights $z_{i}=2$, for all sites, and add large $N$ number of grains at the origin. The distribution of heights $z_{i}$ in the relaxed configuration, produces a very beautiful, but complex pattern (Fig. 4.1) ¹¹1A natural sandpile, formed by pouring sand grains at a constant rate on a flat table with boundaries, gives rise to singular structures like ridges, in the stationary state. This has been attracted much attention recently [hadeler, falcone]. In this chapter, we give a detailed quantitative characterization of similar patterns produced on two directed lattices, starting from a simple periodic background.

These patterns are examples of complex patterns obtained from simple deterministic evolution rules, and are also analytically tractable. Here complexity means that we have structures with variations, and a complete description of which is long. Thus, a living organism is complex because it has many different working parts, each formed by variations in the working out of the same, but relatively much simpler genetic coding. There are some earlier known theoretical examples of complex patterns obtained from local deterministic evolution rules. Most studied among them are the Conway’s game of life [earlierone], which is a cellular automata model with local update rules, and Turing patterns [earliertwo], which are reaction-diffusion models. In general, a detailed and exact mathematical characterization of such patterns has not been possible so far. In this aspect the sandpile patterns are important as they are analytically tractable. Understanding these should also help in studying the more general problem. The most important aspect of these patterns is that, these are the simplest examples that show nontrivial spatial features with proportionate growth (see Fig. 4.2), where these features grows in proportion, keeping the overall shape same. Examples of proportionate growth are abundant in the animal kingdom, where a young animal, typically, grows in size with time, with different parts of the body growing roughly at the same rate. Obviously, this kind of growth requires some coordination and communication between different parts. While there are many models of growing objects studied in physics literature so far, e.g. the Eden model, Diffusion-Limited Aggregation (DLA), invasion-percolation etc. [herrmann], none has this property. All of these are mainly models of aggregation, where growth occurs by accretion on the surface of the object, and inner parts do not evolve significantly (Fig. 4.3). It is worth mentioning, that modelling proportionate growth with simple structures is easy, e. g., growth of a water droplet as one injects more water into it. However, generating a complex pattern with large number of structures inside, all growing at same rate, maintaining their relative shape, is highly nontrivial. This is what happens in the patterns produced in the sandpile models.

Also, there is an astonishing qualitative similarity between the formation of these patterns and the way a fertilized egg develops into a well structured multicellular organism. The development starts with a single cell which divides into more cells, and then they divide in turn. At some stage of the development, few cells generate newer types of cells and form organs. This process continues and after a long time it forms a large complex, but highly coordinated cellular assembly. The same genetic code in each cell is responsible for the cell differentiation and structure formation. In the abelian sandpile when the sand grains are added at the same site, the grains redistribute themselves and a pattern emerges around the addition site, which grows as more and more sands are added. The same redistribution rule is used for all sites, and yet the pattern has large visibly distinguishable structures inside, which, one can think of as different organs in an animal. This simple mathematical model of growing sandpile captures these qualitative features, even though the actual mechanism of growth is much more complicated in the biological world.

The spatial patterns in sandpile models were first discussed by Liu et al. [liu]. The asymptotic shape of the boundaries of the sandpile patterns produced by adding grains at a single site in different periodic backgrounds was discussed in [dhar99]. Borgne et al. [borgne] obtained bounds on the rate of growth of these boundaries and later these bounds were improved by Fey et al. [redig] and Levine et al. [lionel]. The first detailed analysis of different periodic structures found in the patterns were carried out by Ostojic in [ostojic]. Other spatial configurations in the abelian sandpile models, like the identity [borgne, identity, caracciolo] or the stable state produced from special unstable states, also show complex internal self-similar structures [liu], which share common features with the patterns studied here. In particular, the identity configuration on the F-lattice has recently been shown to have spatial structure similar to what we study here [caracciolo]. There are other models, which are related to the abelian sandpile model, e.g. the Internal Diffusion-Limited Aggregation (IDLA), Eulerian walkers (also called the rotor-router model), and the infinitely-divisible sandpile, which also show similar structure. For the IDLA, Gravner and Quastel showed that the asymptotic shape of the growth pattern is related to the classical Stefan problem in hydrodynamics, and determined the exact radius of the pattern with a single point source [gravner]. Levine and Peres have studied patterns with multiple sources in these models recently, and proved the existence of a limit shape[levine_peres].

The pattern on a standard square lattice is rather complicated (Fig.4.1), and it has not been possible to characterize it so far. In this chapter, we consider two variations of the square lattice, assigning orientations to its edges, such that each site has two inward and two outward arrows, as shown in Fig.4.4a and Fig.4.4b. They are known as F-lattice and Manhattan lattice, respectively.

We define a position vector on the lattice, $\mathbf{R}\equiv\left(x,y\right)$. The DASM is defined on the lattice by a height variable $z\left(\mathbf{R}\right)$, at each site $\mathbf{R}$. In a stable configuration all sites have height $z\left(\mathbf{R}\right)<2$. The system is driven by adding grains at a single site and if this addition makes the system unstable it relaxes by the toppling rule: each unstable site transfers one grain each in the direction of its outward arrows. We start with an initial checkerboard configuration in which $z\left(\mathbf{R}\right)=1$ for sites with $(x+y)=$ even, and $0$ otherwise. Clearly, the average density of sand grains for the initial configuration is $1/2$ per site. For numerical purpose we use a lattice large enough so that none of the avalanches reaches the boundary. The result of adding $N=5\times 10^{4}$ grains at the origin is shown in Fig. 4.5 and Fig. 4.6, for the two lattices.

The asymptotic pattern in large $N$ limit, for the two lattices are indistinguishable from each other, at large scales, except that the thin lines of $1$’s forming two triangles outside the octagon are rotated by $45^{\circ}$ (Fig.4.6). Since the lattices are different, this is quite intriguing. Specially there is no obvious lattice symmetry, and it is easily checked that patterns produced for small $N$ are quite different (Fig.4.7).

This pattern is somewhat simpler than in Fig.4.1, which makes its study easier. We shall discuss here only the F-lattice, but the discussion is equally applicable to the Manhattan lattice. Taking some qualitative features of the observed pattern (e.g. only two types of patches are present, and they are all $3$- or $4$- sided polygons) as input, we show how one can get a complete and quantitative characterization of the pattern. We also show that the asymptotic pattern has an unexpected exact $8$-fold rotational symmetry, and determine the exact coordinates of all the boundaries in it. We will also discuss some other cases, where exactly the same asymptotic pattern is obtained.

We first describe a general method of formally characterizing a large number of sandpile patterns, not just the three patterns shown till now. We start by defining $2\Lambda\left(N\right)$ as the diameter of the pattern when $N$ grains have been added. The exact definition of $\Lambda$ is flexible, and the characterization does not depend on the choice. For the patterns in Fig. 4.1, 4.5, and 4.6 we choose $2\Lambda$ as the width of the smallest rectangle that encloses all sites that have toppled at least once.

As mentioned before, the patterns exhibit proportionate growth, i.e., all structures in the pattern grows at the same rate to the diameter. While there is as yet no rigorous proof of this important property, we assume this in the following discussion. Then, it is natural to describe the patterns in the reduced coordinates defined by $\xi=x/\Lambda$ and $\eta=y/\Lambda$. A position vector in these reduced coordinates is defined by $\mathbf{r}=\mathbf{R}/\Lambda\equiv\left(\xi,\eta\right)$. Then in the limit $\Lambda\rightarrow\infty$, the patterns can be characterized by a function $\rho(\mathbf{r})$ which gives the local density of sand grains in a small rectangle of size $\delta\xi\delta\eta$ about the point $\mathbf{r}$, with $1/\Lambda\ll\delta\xi$, $\delta\eta\ll 1$. We define $\Delta\rho\left(\mathbf{r}\right)$ as the change in density $\rho\left(\mathbf{r}\right)$ from its initial background value.

A large number of sandpile patterns, including the one in Fig. 4.1, 4.5, and 4.6, are made of a union of distinct regions, called “Patches”, inside which the heights are periodic in space. Then, inside each patch $\Delta\rho\left(\mathbf{r}\right)$ is constant. For example, in the pattern in Fig. 4.5, there are only two types of patches, and the change in density takes only two possible values, $1/2$ in a high-density patch (color yellow in Fig. 4.5) and $0$ in a low-density patch (color orange). There are few defect-lines, which move with $N$, and can also be seen in Fig.4.1 and Fig.4.5. But these can be ignored in discussing the asymptotic pattern.

Let $T_{\Lambda}\left(\mathbf{R}\right)$ be the number of topplings at site the $\mathbf{R}$ when the diameter reaches the value $2\Lambda$ for the first time. Define

$$\phi\left(\mathbf{r}\right)=\lim_{\Lambda\rightarrow\infty}\frac{1}{2\Lambda^{2}}T_{\Lambda}\left(\mathbf{R}^{\prime}\right),$$

(4.1)

where $\mathbf{R}^{\prime}\equiv\left(\lfloor\Lambda\xi\rfloor,\lfloor\Lambda\eta\rfloor\right)$, with $\lfloor x\rfloor$ being the floor function which gives the largest integer $\leq x$. From the conservation of sand-grains in the toppling process, it is easy to see that

$$\sum_{\mathbf{R^{\prime}}\in n.n.}T_{\Lambda}\left(\mathbf{R^{\prime}}\right)-\alpha T_{\Lambda}\left(\mathbf{R}\right)=\Delta z\left(\mathbf{R}\right)-N\delta_{\mathbf{R},0},$$

(4.2)

where the sum is over the sites nearest neighbors of $\mathbf{R}$, and $\alpha$ is the number of them. Then in the rescales coordinate, $\phi$ satisfies the Poisson equation

$$\nabla^{2}\phi\left(\mathbf{r}\right)=\Delta\rho\left(\mathbf{r}\right)-\frac{N}{\Lambda^{2}}\delta\left(\mathbf{r}\right).$$

(4.3)

In an electrostatic analogy, we can think of $\Delta\rho(\mathbf{r})$ as an areal charge density, and $\phi(\mathbf{r})$ as the corresponding electrostatic potential. A complete specification of $\phi\left(\mathbf{r}\right)$ determines the density function $\Delta\rho\left(\mathbf{r}\right)$ which in turn characterizes the asymptotic pattern.

The key observation that allows us to determine the asymptotic pattern is the following proposition.

A proof can be done in the following way. Within a patch, the function $\phi(\mathbf{r})$ is Taylor expandable around any point $\mathbf{r}_{O}\equiv\left(\xi_{o},\eta_{O}\right)$ inside the patch.

$$\phi\left(\mathbf{r}\right)=f+d\Delta\xi+e\Delta\eta+a\left(\Delta\xi\right)^{2}+2h\Delta\xi\Delta\eta+b\left(\Delta\eta\right)^{2}+\mathcal{O}\left(\Delta\xi^{3},\Delta\eta^{3}\right)\dots,$$

(4.4)

where $\Delta\xi=\xi-\xi_{o}$ and $\Delta\eta=\eta-\eta_{o}$. Consider any term of order $\geq 3$ in the expansion, for example the term $\sim(\Delta\xi)^{3}$. This can only arise due to a term $\sim(\Delta x)^{3}/\Lambda$ in $T(x,y)$. Then, considering the fact that $T\left(\mathbf{R}\right)$ is an integer function of the coordinates, it is easy to see that it will change discontinuously at intervals of $\Delta x\sim\mathcal{O}(\Lambda^{1/3})$. This leads to change in the periodicity of heights at such intervals inside each patch which themselves are of size $\sim\Lambda$. This would then result in an infinitely many defect-lines in the asymptotic pattern. However there are no such features in Fig.4.1 or Fig.4.5. Therefore inside each periodic patch of constant $\Delta\rho(\mathbf{r})$, $\phi(\mathbf{r})$ can at most be quadratic in $\xi$ and $\eta$.

The argument finally boils down to proving the two features of the pattern, i.e., there is proportionate growth, and that the pattern can be decomposed in terms of periodic patches.    

In each periodic patch the toppling function $T(x,y)$ is a sum of two terms: a part, that is a simple quadratic function of $x$ and $y$, and another is a periodic part. The periodic part averages to zero, and does not contribute to the coarse-grained function $\phi(\mathbf{r})$ ²²2In some patterns, with other backgrounds (not discussed here) there are regions that occupy finite fraction area of the full pattern, which show aperiodic height patterns. These cases are harder to analyse.. The quadratic part, when rescaled, can be written as

$$\phi\left(\mathbf{r}\right)=f+d\xi+e\eta+a\xi^{2}+2h\xi\eta+b\eta^{2},$$

(4.5)

where $a,h,b,d,e,f$ are constants inside a patch, and $a+b=\Delta\rho/2$, corresponding to the patch. Then each patch is characterized by the values of these parameters.

Now we will show that the continuity of $\phi$ and its first derivatives along the boundary between adjacent patches imposes linear relations among the corresponding parameters. Consider two neighboring periodic patches ${\bf P}$ and ${\bf P^{\prime}}$ with mean densities $\rho$ and $\rho^{\prime}$ respectively. Let the rescaled quadratic toppling function be $Q(\mathbf{r})$ and $Q^{\prime}(\mathbf{r})$ in these patches. Then the boundary between the patches is given by the equation $Q(\mathbf{r})=Q^{\prime}(\mathbf{r})$. As the derivatives of $\phi$ are also continuous across the boundary, the boundary between two periodic patches must be a straight line, and

$$Q^{\prime}(\mathbf{r})=Q(\mathbf{r})+\frac{1}{2}(\rho^{\prime}-\rho)l_{\perp}^{2},$$

(4.6)

where $l_{\perp}$ is the perpendicular distance of $(\mathbf{r})$ from the boundary. We can start with a periodic patch ${\bf P}$, and go to another patch ${\bf P^{\prime}}$ by more than one path. Since the final quadratic function at ${\bf P^{\prime}}$ should be the same whichever path we take, this imposes consistency conditions which restrict the allowed values of slopes of the boundaries. Consider a point $z_{0}$ where $n$ periodic patches meet, with $n>2$ (Fig.4.8). If the $j$th boundary at this point makes an angle $\theta_{j}$ with the $x$-axis, and the density of the patch in the wedge $\theta_{j}\leq\theta\leq\theta_{j+1}$ is $\rho_{j+1}$ (Fig.4.8), then the condition that the net change in the quadratic form is zero if we go around $z_{0}$ once, reduces to the following condition:

$$\sum_{j=1}^{n}(\rho_{j+1}-\rho_{j})e^{2i\theta_{j}}=0,$$

(4.7)

with $\rho_{n+1}=\rho_{1}$. For $n=3$, with $\rho_{1}\neq\rho_{2}\neq\rho_{3}$, this equation has only trivial solutions with $\theta_{j}$ equal to $0$ or $\pi$ for all $j$. Hence, only $n\geq 4$ are allowed.

These linear equations amongst the parameters corresponding to neighboring patches, can be solved and this will determine the complete potential function $\phi$, giving a quantitative characterization of the pattern.

We now apply this method, in the last section, to the F-lattice pattern in Fig. 4.5, and determine the exact potential function $\phi\left(\mathbf{r}\right)$. We note that in this pattern, there are no aperiodic patches, only two types of periodic patches, where $\rho(\mathbf{r})$ only takes values $1$ or $1/2$. Also, the slopes of the boundaries between patches only take values $0$, $\pm 1$ or $\infty$. The patches are typically dart shaped quadrilaterals, and some triangles. These simplifications, not present in Fig.4.1, make possible a full characterization of the pattern in Fig.4.5.

We start by determining the exact asymptotic size of the pattern. We note from Fig.4.5 that the boundary of the pattern is an octagon (we shall prove later that this is a regular octagon). In fact, there are four lines of $1$’s outside the octagon. But these have zero areal density in the limit $N\rightarrow\infty$, and do not contribute to $\rho(\mathbf{r})$. We will ignore these in the following discussion.

Let $B$ be the minimum boundary square containing all ($\mathbf{r}$) that have a non-zero charge density $\Delta\rho(\mathbf{r})$. We observe that $B$ can be considered as a union of disjoint smaller squares, each of which is divided by diagonal into two parts where $\Delta\rho(\mathbf{r})$ takes values $1/2$ and $0$ (Fig.4.9). This is seen to be true for the outer layer patches. Towards the center, the squares are not so well resolved. Assuming that this construction remains true all the way to the center, in the limit of large $N$, the mean density of negative charge in the bounding square $=1/4$. Given that the total amount of negative charge is $-1$, the area of the bounding square should be $4$. Hence, the boundaries of the minimum bounding square are

$$|\xi|=1,\textrm{ }|\eta|=1.$$

(4.8)

This means, with our choice of the diameter as the width of the box $B$, we have

$$\Lambda\left(N\right)=\sqrt{N}+\textrm{lower order terms}.$$

(4.9)

In Fig. 4.11, we have shown, the correction term appears to grow as $N^{1/4}$.

Most of the time the avalanches does not reach the boundary. They are often stopped by the defect-lines inside the patches, which breaks the periodicity of the heights. For example, there are lines of alternating $1$’s and $0$’s inside the dense (all $1$) patches. When an avalanche enters the patch, the defect line shifts its position, partially increasing the size of the patch. An example of such event is shown in the Fig.4.10. Because of this, the diameter increases in steps with the increase of $N$ (see Fig.4.12).

Let $N_{b}$ be the minimum number of particles that have to be added so that at least one site at $y=b$ topples. We find that for $b=10$, $50$, $100$, and $300$, $\sqrt{N_{b}}=10.770$, $49.436$, $98.894$ and $297.798$. This is consistent with Eq. (4.9).

Then, the Poisson equation, in Eq. (4.3), for this pattern becomes,

$$\nabla^{2}\phi\left(\mathbf{r}\right)=\Delta\rho\left(\mathbf{r}\right)-\delta\left(\mathbf{r}\right),$$

(4.10)

i.e., there is unit amount of point charge at the origin.

We now determine the parameters in the quadratic form of the potential function. In order to do that in a consistent way, we first look at the topological structure of the pattern. We note that the patches become smaller, and there are more of them in number, as we move towards the center. One can use a coordinate transformation $r^{\prime}=1/r^{2}$, $\theta^{\prime}=\theta$ to avoid this overcrowding (Fig.4.13). We can now draw the adjacency graph (Fig.4.15$a$) of the pattern, where each vertex denotes a patch, and a bond between the vertices is drawn if the vertices share a common boundary. It is convenient to think of the triangular patches in the pattern as degenerate quadrilaterals, with one side of length zero. Then we see that the adjacency graph is planar with each vertex of degree four, except a single vertex of coordination number eight corresponding to the exterior of the pattern. The graph has the structure of a square lattice wedge of wedge angle $4\pi$. The square lattice structure of the adjacency graph is seen more clearly, if rather than $1/r^{2}$ transformation, the transformation used is $z^{\prime}=1/z^{2}$ ( this has been used earlier in [ostojic]), where $z=\xi+i\eta$, and view it in the complex $z^{\prime}$-plane (see Fig.4.14). Thus, one can equivalently represent the graph as a square grid on a Riemann surface of two sheets (fig.4.15$b$).

We now use the qualitative information obtained from the adjacency matrix of the observed pattern, to obtain quantitative prediction of the exact coordinates of all the patches. Consider an arbitrary patch ${\bf P}$, having an excess density $1/2$. The potential function in this patch is a quadratic function of $(\xi,\eta)$ and we parametrize it as

	$\displaystyle\phi_{{}_{\bf P}}(\mathbf{r})$	$\displaystyle=$	$\displaystyle\frac{1}{8}(m_{{}_{\bf P}}+1)\xi^{2}+\frac{1}{4}n_{{}_{\bf P}}\xi\eta+\frac{1}{8}(1-m_{{}_{\bf P}})\eta^{2}$		(4.11)
			$\displaystyle+d_{{}_{\bf P}}\xi+e_{{}_{\bf P}}\eta+f_{{}_{\bf P}}.$

The potential function in another patch ${\bf P}$ having zero excess density is parametrized as

$$\phi_{{}_{\bf P}}(\mathbf{r})=\frac{1}{8}m_{{}_{\bf P}}(\xi^{2}-\eta^{2})+\frac{1}{4}n_{{}_{\bf P}}\xi\eta+d_{{}_{\bf P}}\xi+e_{{}_{\bf P}}\eta+f_{{}_{\bf P}}.$$

(4.12)

Now consider two neighboring patches ${\bf P}$ and ${\bf P^{\prime}}$ with excess densities $1/2$ and $0$ respectively. Then using the matching conditions (see Eq. (4.6)), it is easy to show that if the boundary between them is a horizontal line $\eta=\eta_{{{}_{\bf P}}}$, we must have

	$\displaystyle m_{{{}_{\bf P}}^{\prime}}$	$\displaystyle=$	$\displaystyle m_{{}_{\bf P}}+1\textrm{, }n_{{{}_{\bf P}}^{\prime}}=n_{{}_{\bf P}}\textrm{, }d_{{{}_{\bf P}}^{\prime}}=d_{{}_{\bf P}},$
	$\displaystyle e_{{{}_{\bf P}}^{\prime}}$	$\displaystyle=$	$\displaystyle e_{{}_{\bf P}}+\eta_{{}_{\bf P}}/2\textrm{, }f_{{{}_{\bf P}}^{\prime}}=f_{{}_{\bf P}}-{\eta_{{}_{\bf P}}}^{2}/4.$		(4.13)

There are similar conditions for other boundaries. These result a coupled set of linear equations for the coefficients $\{m_{{}_{\bf P}},n_{{}_{\bf P}},d_{{}_{\bf P}},e_{{}_{\bf P}},f_{{}_{\bf P}}\}$. The equations for $m_{{}_{\bf P}}$ and $n_{{}_{\bf P}}$ do not involve other variables. In the outermost patch, clearly $\phi(\mathbf{r})=0$, and for this patch both $m$ and $n$ are zero. It follows that $m_{{}_{\bf P}}$ and $n_{{}_{\bf P}}$ are integers, equal to the Cartesian coordinates of the vertex corresponding to the patch ${\bf P}$ in the discretized Riemann surface in Fig.4.15b. In the following, we denote a patch by integers $(m,n)$, and write the corresponding coefficients $d_{{}_{\bf P}}$, $e_{{}_{\bf P}}$, and $f_{{}_{\bf P}}$ as $d_{m,n}$, $e_{m,n}$ and $f_{m,n}$. With this convention, the matching conditions in Eq.(4.13) can be rewritten as

$$d_{m+1,n}=d_{m,n}\textrm{, }e_{m+1,n}-e_{m,n}=\eta_{m,n}/2\textrm{, $(m+n)$ odd}.$$

(4.14)

Using similar matching conditions for the boundary of patch $(m$, $n)$ with slope $\pm 1$, we get the conditions

	$\displaystyle d_{m,n+1}-d_{m,n}$	$\displaystyle=$	$\displaystyle e_{m,n}-e_{m,n+1},\textrm{ }(m+n)\rm{~odd},$
	$\displaystyle d_{m,n-1}-d_{m,n}$	$\displaystyle=$	$\displaystyle e_{m,n-1}-e_{m,n},\textrm{ }(m+n)\rm{~odd}.$		(4.15)

We can eliminate the variables $d_{m,n}$ and $e_{m,n}$ with $(m+n)$ even using Eq. (4.14) and Eq.(4.15). Then the equations become

	$\displaystyle e_{m+2,n}-e_{m,n}$	$\displaystyle=$	$\displaystyle\eta_{m,n}/2,$		(4.16)
	$\displaystyle d_{m-2,n}-d_{m,n}$	$\displaystyle=$	$\displaystyle\xi_{m,n}/2,$		(4.17)
	$\displaystyle d_{m-1,n-1}-d_{m,n}$	$\displaystyle=$	$\displaystyle e_{m+1,n-1}-e_{m,n},$		(4.18)
	$\displaystyle d_{m-1,n+1}-d_{m,n}$	$\displaystyle=$	$\displaystyle-[e_{m+1,n+1}-e_{m,n}].$		(4.19)

It is convenient to introduce the complex variables $z=\xi+i\eta$, $M=m+in$ and $D=d+ie$. In these variables we can write the potential function, in Eq. (4.11) and (4.12), as

$$\phi(z)=\frac{1}{8}z\bar{z}+\frac{1}{8}Re[z^{2}\bar{M}+\bar{D}z]+f,$$

(4.20)

where overbar denotes complex conjugation.

On the $(m,n)$ lattice, with $(m+n)$ odd, the natural basis vectors are $(1,1)$ and $(1,-1)$. Let us call these $\alpha$ and $\beta$. We define the finite difference operators $\Delta_{\pm\alpha}$ and $\Delta_{\pm\beta}$ by

	$\displaystyle\Delta_{\pm\alpha}f(z)=f(z\pm\alpha)-f(z),$
	$\displaystyle\Delta_{\pm\beta}f(z)=f(z\pm\beta)-f(z).$		(4.21)

Then the equations (4.16-4.19) can be written as

	$\displaystyle\Delta_{-\alpha}d=\Delta_{\beta}e,$
	$\displaystyle\Delta_{-\beta}d=-\Delta_{\alpha}e.$		(4.22)

These equations are the discrete analog of the familiar Cauchy-Riemann conditions connecting the partial derivatives of real and imaginary parts of an analytic function where the role of the analytic function is played by $D=d+ie$.

From Eq.$(\ref{a6})$ and Eq.(4.19), it is easy to deduce that $D$ satisfies the discrete Laplace’s equation

$$[\Delta_{\alpha}\Delta_{-\alpha}+\Delta_{\beta}\Delta_{-\beta}]D=0.$$

(4.23)

If $m$ and $n$ are large, the corresponding patch is near the origin ($|\xi|+|\eta|$ is small), and where the leading behavior of $\phi(\mathbf{r})$ is given by $\tilde{\phi}(\mathbf{r})\sim-\frac{1}{4\pi}\log(\xi^{2}+\eta^{2})$ (see Eq. 4.10). Consider a point $z_{0}$, such that at $z_{0}$

$$\partial^{2}\tilde{\phi}/\partial\xi^{2}\approx m/4;~~~~~\partial^{2}\tilde{\phi}/\partial\xi\partial\eta\approx n/4.$$

(4.24)

Then, $z_{0}$ would be expected to lie in the patch labeled by $(m,n)$. This gives $z_{0}\approx\pm(\pi\bar{M}/2)^{-1/2}$. Then, setting $\partial\tilde{\phi}/\partial z$ equal to $\bar{D}/2$ gives us

$$D_{m,n}\simeq\pm\frac{1}{\sqrt{2\pi}}\sqrt{m+in}.$$

(4.25)

The equation (4.23), subjected to the behavior at large $|m|+|n|$ given by Eq.(4.25) on the $4\pi$-wedge graph (for each value of $(m,n)$, $D_{m,n}$ has two values) has a unique solution. Clearly the solution has eight fold rotational symmetry about the origin in the $(m,n)$ space. This implies that

$$D_{-n,m}=i^{1/2}D_{m,n};\rm{~for~all~}(m,n).$$

(4.26)

Given $D_{m,n}$, its real and imaginary parts determine $d_{m,n}$ and $e_{m,n}$, and using Eq.(4.16, 4.17) we determine the exact positions of all the patch corners. The exact eight-fold rotational symmetry of the adjacency graph of the pattern, and the fact that $D$ satisfies Eq.(4.26) on the adjacency graph together imply the eight-fold rotational symmetry of all the distances in the pattern.