Some Properties of Sandpile Models as Prototype of Self-Organized Critical Systems
SSome Properties of Sandpile Models as Prototype of Self-Organized Critical Systems
M. N. Najafi, S. Tizdast, and J. Cheraghalizadeh Department of Physics, University of Mohaghegh Ardabili, P.O. Box 179, Ardabil, Iran ∗ This paper is devoted to the recent advances in self-organized criticality (SOC), and the concepts.The paper contains three parts; in the first part we present some examples of SOC systems, in thesecond part we add some comments concerning its relation to logarithmic conformal field theory,and in the third part we report on the application of SOC concepts to various systems ranging fromcumulus clouds to 2D electron gases.
PACS numbers: 05., 05.20.-y, 05.10.Ln, 05.45.DfKeywords: sandpile model, invasion, fluid dynamics, critical exponents
I. INTRODUCTION
Since the fascinating work of Bak, Tang and Wiesen-feld (BTW) coined the self-organized criticality (SOC)in 1987 [1] a huge number of papers appeared toexplore various aspects of this term. It forms now alarge class of critical phenomena. These systems showcritical properties without tunning of any externalparameter, for which the BTW sandpile model wasthe first prototype. Dhar discovered for the first timethe Abelian structure of sandpiles so that we call themabelian sandpile models (ASM) [2]. Despite its simpledynamics, ASM has various interesting features andnumerous works, analytical and computational, havebeen done on this model [3–11]. Thanks to conformalfield theory (CFT), it is known that the BTW modelis described with c = − c being the centralcharge [4]. Additionally the geometrical aspects ofthis model are understood in terms of its relation toloop-erased random walks (LERW) [12], which itself isrelated to Schramm-Loewner evolution (SLE) with thediffusivity parameter κ = 2 [13, 14].In this paper, after introducing the original work ofBTW, we explore some properties of the model, withan emphasis on the dynamics of the avalanches in sand-piles. Its relation to Logarithmic Conformal Field The-ory (LCFT) is explored. In the last part of the paper,we review various aspects of SOC in various systems, in-cluding the SOC in the fluid propagation in porous me-dia, in cumulus clouds, in excitable complex networks,and in imperfect supports. We also describe the SOCtechniques to explain the 1 /f noise and metal insulatortransition in two-dimensional electron gas (2DEG). Weintroduce also some generalizations of BTW, which arethe invasion BTW model, and the BTW model on vi-brating systems.The paper is organized as follows: in the following sec-tion, we introduce shortly the SOC in Nature and explainbriefly the examples and avalanche dynamics and the ba-sic ingredients. Section III is devoted to the definition of ∗ [email protected] the BTW model and other variants. We will present rela-tion to the logarithmic conformal field theory, includingthe ghost-free fields and W-Algebra in the Sec. IV. Vari-ous applications of SOC concepts to natural processes ispresented in Sec. V. II. SOC IN NATURE
In nature, there are circumstances that unlike thethermal critical systems, no external parameter needs tobe adjusted in order to reach and maintain in criticality,which are called self-organized critical (SOC) systems.Such a system automatically reaches and organizes itselfin a critical state [15]. These systems, which are usuallyopen and absorb and dissipate energy, change with timebut their general properties are almost unchanged onobserved time scales. They need external energy tocompensate the dissipation.The original aim of the BTW model (to be explainedin the next section) was to explain the ”1 /f noise” phe-nomena that is seen in many natural systems, like rainfall [16], sun flares [17, 18], real piles of rice and otherobjects [19, 20], earthquake [21–28], forest fire [29], andclouds [30–34]. The aim of this section is to introducethese natural phenomena as a motivation for analyzingthe models of SOC. A. Examples
1. SOC in Earthquake
One of the most popular examples of the SOC systemsis earthquake for which the frequency of earthquakes ( N )with energy ( E ) follows Gutenberg-Richter’s power-lawrelation as follows [35, 36] N = aE − b (1) a r X i v : . [ c ond - m a t . s t a t - m ec h ] S e p where N is the number of earthquakes with energy E , a is the constant number that is a measure of the size andthe amount of vibrational activity in the area, and b isa critical exponent being often between 0 . . A asfollows N = cA − d (2)where d ≈ .
40 [29] is another exponent. SOC inthe earthquake is reported in many cases [21]. ManySOC models have been developed to capture thephysics of earthquakes, like block-spring models [37, 38],sandpile based models [21, 26, 27], and sandpile onearthquake network [39]. In these models, the dy-namics are predicted to be avalanche-like, based ona local stimulation (by increasing the local stress andtension), and the spread of stress throughout the system.Generally, two strategies are often taken for explainingthe observations of earthquakes: the quenched-disorderbased models ascribing the observations of the seismicactivities to the geometric and material irregularitiesin the earth, and the dynamical-instability modelsattributing the complexities to the stochastic forcingarising from the dynamic nonuniformities [25]. In theformer, the power-laws observed in an earthquake isrelated to geometric features of the fault structure [40].Whether the earth is operating according to one ofthese schemes or in a hybrid one remains an open andfundamental problem. The application of SOC ideasfor earthquake is very efficient and provides realisticresults [21, 26, 27]. The Olami-Feder-Christensen earth-quake model model [26] is a two-dimensional coupledmap lattice model which is known as a simplified versionof the Burridge-Knopoff spring-block model [41] forearthquakes. This model is famous and attracted muchattention for it serves as a paradigm for nonconservativeSOC systems, and also reproducing the most importantstatistical property of real earthquakes [26, 42], and alsoOmori’s law [43, 44], and the statistics of foreshocks andaftershocks [45].In a recent study, the ideas of SOC models were appliedto the Rigan earthquake [39], in which a close relationshipwas observed between the dynamics of the SOC modeland the real data of the earthquake. This relation can beunderstood by focusing on the stress propagation due tothe tectonic motion of the continental plates, which is aslow steady process (like the other “slowly driven” SOCsystems), but the release of stress occurs sporadically inbursts of various sizes.
2. SOC in Forest Fire
Forest fire is another natural phenomenon thatshows SOC behaviors, like power-law and scaling be-haviors [46, 47]. Various surveys of firefighting data in different parts of the United States and Australia haveshown a range of the size critical exponents between1 . .
5, depending on the area [29, 47]. The SOCstructure of forest fire was discovered by Malamud et. al. [47]. Many models have emerged in order tocapture the physics of this phenomenon. Among them,an important one is the Drossel-Schwabl model [48, 49]which in some limit gives acceptable exponents.As a model for forest fire, let us consider the Drossel-Schwabl model [48–51], defined on a d -dimensional latticewith lattice length L . In each time, each site of thesystem is empty or occupied by a green tree or a burningtree. At the initial time t = 0 we suppose that thelattice sites are either occupied by green, or empty. Thelattice state is updated at any time by the following rules: . The burning site will be vacated in the nexttime, . The site where the green tree is located will catch fireat the next time if at least one of its nearest neighbor isburning, otherwise, it fires spontaneously with lightningprobability f , . In an empty site, a tree grows with a p probability.Starting from an initial tree configuration, one Themodel can become critical only in the limit p → f /p →
0. In this limit, the length of the correlation isdivergent [52]. The latter provides the conditions un-der which the time scale of tree growth and burning outthe forest clusters are well-separated. The scale of theaverage size of clusters is controlled by the growth rate θ ≡ pf . In the simulations, at each time step θ sites arerandomly chosen for being occupied (if it is already oc-cupied nothing happens, otherwise it turns to occupied).Then a randomly chosen site is ignited so that all sitesin the connected cluster to which the start site belongsburn. This model becomes critical in the limit θ → ∞ ,although there has been much discussion whether thismodel is SOC or not [50].
3. SOC in Sun Flares
After the observation of R. C. Carrington and R.Hodgson in 1859 on the solar flares in white light, muchattention has been paid to this problem. The cause ofsolar activities is the presence of a solar magnetic fieldin the hot plasma around the outer layer of the sun orthe convective zone created by the collision of particles.The convective zone is the highest inner layer of the sunthat extends from the radiative zone to the surface ofthe sun. This area is made up of convective effervescentcells. The magnetic flux forms active regions whichinclude sunspot [18, 53].Price et al. [54] criticized the hypothesis of chaos in so-lar activity and found no evidence for a low-dimensionaldeterministic nonlinear process using analysis of thesunspot number time series. Continuing this criticaltrend, the theory of Self- Organized Criticality (SOC)was proposed as the basic mechanism for explaining solaractivity [55, 56]. According to the concept of SOC the-ory, the solar corona operates in a self-organized criticalstate, while the solar flares constitute random avalancheevents with a power-law profile like the earthquake pro-cess. L.P. Karakatsanis and G.P. Pavlos [54] presentednew results that strongly support the concept of low di-mensional chaos and SOC theory. In this study, theyfound a co-existence between the self-organized criticalstate according to the soc theory and low dimensionalchaotic dynamics underlying to the solar activity. Theseresults are obtained by nonlinear analysis of the sunspotindex. For the original signal, the largest Lyapunov expo-nents were found to be zero, while the correlation integralslope profile was similar to the alternative data slope pro-file, showing a high-dimensional stochastic process and acritical state based on SOC theory [57].
4. SOC in Rain-Falls
Another natural example of the SOC phenomenon isthe rainfall, which is witnessed by the power-law behaviorfor the number of rain events versus size and the numberof droughts versus duration [16, 58]. In this case, sta-ble stimulation is provided by the heat of the sun, whichcauses the oceans to evaporate, for which rain relaxationwith the continuous and uninterrupted event, rain. Pe-ters et. al showed that the accumulated water columndisplays scale-less fluctuations and also the number den-sity of rain events per year N ( M ) versus event size M be-haves like power-law, with exponent 1 .
36, and the num-ber density of droughts per year N ( D ) versus droughtduration D with an exponent 1 .
42. To understand theother quantity that was shown to be in power-law, let usdefine the rain-fall rate q ( u ) ≡ (cid:80) i n i V i u i , where n i is thethe density number of droplets with a volume of V i thatreaches the earth at the speed of u i . Then it was shownthat R ( τ ) /S ( τ ) ∼ τ H , (3)where R ( τ ) ≡ max ≤ t ≤ τ X ( t, τ ) − min ≤ t ≤ τ X ( t, τ ) (4)and X ( t, τ ) ≡ (cid:80) tu =1 ( q ( u ) − (cid:104) q (cid:105) τ ), and (cid:104) q (cid:105) τ ≡ τ (cid:80) τt =1 q ( t )∆ t [16]. Self-organized criticality in rainfallswas observed in many other studies [59–63], which stim-ulated many theoretical studies on the subject [64–68].
5. SOC in Clouds
Self-affinity and scaling properties in clouds havebeen found from satellite images [69], and in particular in cumulus clouds [70] on several scales. Variousobservables were shown to exhibit scaling behavior,like the area-perimeter relation [69–75], the nearestneighbor spacing [76], the rainfall time series [77], clouddroplets [78], and the distribution function of geometricalquantities [79–82]. After these observations, and consid-ering the multi-fractality of clouds [69, 70, 73, 83–86],attempts for classifying clouds into universality classeswere carried out based on cloud field statistics [87–89]and cloud morphology [90]. The self-organized crit-icality in atmospher was first detected by Peters byanalysing the precipitation [91], and developed furtherfor atmospheric convective organization [92]. The fractaldimension of perimeter of self-organized vorties shown tobe near using the quasi geostrophic vorticity equation.The areaperimeter relation with the D = 1 . ± .
02 ofcirrus, and D = 1 . ± .
05 for cumulunimbus tropicalclouds. These fractal dimensions are in agreementwith the relative turbulent diffusion model, predicting1 .
35, which is also confirmed by means of some otherobservations.As a main building block of the atmosphere dynam-ics, turbulence seems to be essential in the dynamics andformation of clouds. Analysis of the images from thefair weather cumulus clouds reveals that they addition-ally exhibit self-organized criticality degrees of freedom,leading us to use the term SOC turbulent state. Observa-tions (in our submitted paper) support the fact that thissystem, when projected to 2D, demonstrates conformalsymmetry compatible with c = − c = 0 conformalfield theory. Using a mix of turbulence and cellular au-tomata, namely, the coupled map lattice model [93], oneobtains the same exponents as the observations. Also,in a separate (unpublished yet) work we developed a 2Dmonte carlo based stochastic model including the compe-tition between avalanche dynamics and cohesive energybetween water droplets that generates the same proper-ties. The fractal geometry of clouds was seen in manyreal observations, like the multi-fractality structure ofclouds, universality classes of cloud fields, analysis rain-fall time series, nearest-neighbor spacing statistics, cu-mulus cloud morphology, ”variable” and ”steady” cloudyregions for warm continental cumulus cloud, fractal anal-ysis of high-resolution cloud droplet measurements, thefractal dimension of noctilucent clouds, the fractal di-mension of convective clouds around Delhi, scaling prop-erties of clouds, self-similarity of clouds in the intertrop-ical convergence zone, depending on the equivalent blackbody temperature.
6. SOC in real piles
SOC has been observed in labs for real piles, like pile ofbeads [94], rice pile [95–99], and other granular piles withvarious aspect ratio [100]. It has been observed that fordry sandpiles, macroscopic behaviors can be determinedfrom the angle θ c [4]. This angel is called the thresholdangle, which depends on the structural details of its con-stituent grains. A sandpile with a local slope less than θ c anywhere is stable, and adding a small amount of sandwill cause a small reaction, but if this operation resultsin a slope larger than θ c , then an avalanche is formed,which sometimes is of the system size. On a pile with anaverage slope a bit less than θ c , the pile’s response to theaddition of sand is not very predictable. One possibil-ity is that there will be no relaxation, or it may cause amedium-size avalanche or a catastrophic avalanche thatwill affect the entire system, called the critical state. Bak et. al. [57] observed that if we build a pile by slowly pour-ing sand on a flat circular plate, we would get a conicalpile with a slope equal to θ c . This system is constantlymoving towards its critical state, this is an indication ofSOC. The steady states of this process are described bythe following properties:Sand is being added to the system at a constant smallrate, but it leaves the system in a very irregular manner,with long periods of apparent inactivity interspersed byevents that may vary in size and which occur at unpre-dictable intervals. The interesting thing about these sys-tems is that in steady-state, where the average amountof input and output energy is equal, the system exhibitscritical properties. For example, power-law are foundfor various observations in the system, and the lengthof the correlation is infinite. Bak, Tang, and Wiesenfeldprovided the BTW model that is an automatic cellularmodel for sandpile. This model is defined on the finitelattice. There is a positive integer variable at each siteof the lattice, called the height of sandpile at the site( z i ) and we have the threshold height for our lattice thatcalled critical height ( z c , it is often 2 d that d is latticedimension.) At each time step a site is picked randomly,and it’s height z i is increased by unity ( z i → z i + 1).If z i > z c , this site is unstable. It relaxes by topplingwhereby four sand grains leave the site, and each of thefour neighboring sites gets one grain. If there is any un-stable site remaining, it is toppled too. This process isan avalanche. The system in the BTW model reachesa stable state after the transition from transient states,which are called recurrent state because they are likelyto be repeated during the casual process.The rice pile model is another example of self-organizedcriticality phenomena. The mechanism, in this case, isthat we consider two pieces of glass with a specific diam-eter of 5 mm and a distance of 16 mm from each other.We pour rice with a certain length from the middle ofthese two pieces of glass unit the system is stable, thenwe add colored grains similar to rice to the system andfollow their movement so that we record the time of theirentry and exit into the system. By obtaining the differ-ence between the two terms, we measure the distributionof the length of time, the grains have been in the system.The following equation is obtained by plotting the time distribution in terms of lattice length ( L ): P ( T, L ) = L − β F ( TL ν ) (5)That ν and β are the critical exponents. The statisti-cal distribution of quantities in self-organized criticalityphenomena is the power distribution.There are other examples of SOC systems, like riverbasins [101, 102], in air pollution [103], in climatechange [104], in brain plastisity [105], in stock mar-kets [106, 107], in magnetosphere [108], in midlati-tude geomagnetic activity [109], in magnetohydrodynam-ics [110], in Kardar-Parizi-Zhang growth model [111] inBean state in YBCO thin films [112], in granular sys-tems [113], and in much more systems [114] which areout of scope of this paper. B. Avalanche Dynamics and the Basic Ingredients
It is a common belief that avalanche dynamics are theunderlying mechanism that is responsible for SOC behav-iors. Although normal diffusive transport and branchingare believed to be very basic ingredients of SOC, it wasshown in [115] that the avalanches in SOC can be linearas branchless random walks. In this case, their scalingbehavior is different from that of branched avalanches.Sandpile models were introduced by Bak, Tang, andWiesenfeld [1] (BTW) as a prototypical example for aclass of models that show self-organized criticality. Thesemodels show critical behavior without fine-tuning of anyexternal parameter. BTW model includes avalanche-based dynamics in which the system is slowly stimulated,i.e. subjected to small external perturbations. Largeevents in these systems, which are the result of thesesmall stimuli, occur less frequently and on a larger scale,i.e. the energy is gradually absorbed and is excreted outon a larger scale. An interesting feature in these systemsis that in the steady-state, where the average amount ofenergy input and output is equal, the system exhibitscritical properties, which is a SOC state. The Abelianstructure of the sandpile model was first discovered byDhar so that it was thereafter named as Abelian sand-pile model (ASM) [25]. Despite its simplicity, ASM hasvarious interesting features and numerous works, analyt-ical and computational, have been done on this model.Among them one can mention different height and clusterprobabilities [40], the and avalanche distribution [116],and also its the connection to the other models like thespanning trees [117], the ghost model [118, 119], and the q -state Potts model [120]. For a good review see [4, 121].Moreover, some of these results are analyzed in light ofconformal field theory description with central charge c = − κ = 2 [119]. III. THE BTW MODEL
Let us consider the BTW on a two-dimensional square d − dimensional lattice hypercubic lattice with a linearsize L and coordination number z = 2 d . To each site aheight variable h i is assigned which takes its value fromthe set 1 , , ..., z . This height variable shows the numberof sand grains in the underling site. The dynamics ofthis model is as follows: in each step, one grain of sandis added to a randomly chosen site i , i.e., h i → h i + 1. Ifthe resulting height becomes more than z (i.e. becomesunstable), the site topples and looses 2 d grains of sand,each of which is transferred to one of 2 d neighbors ofthe toppled site. As a result, the neighboring sites maybecome unstable and topple, and in this way, a chainof topplings may happen in the system until the sys-tem reaches a state with no unstable site. The chain oftopplings then is called an avalanche. If a boundary sitetopples, one or two grains of sand (for the sites in the cor-ners of the lattice two grains, and for the other boundarysites one grain) will leave the system. After reaching astable configuration (the avalanche is finished), the pro-cess is repeated starting from another random site forgrain injection. The toppling rule for the site i can alsobe written in the form h j → h j + ∆ ij , where∆ ij = − z if j = i +1 if j and i are neighbors0 otherwise (6)which is discrete Laplacian operatore. A. Transient v.s. Recurrent Configurations
Let us suppose that we start from a random heightconfiguration. Then during the system evolution, manyconfigurations come about, some of which are transient,meaning that they do not occur again, and some of whichare recurrent. In fact, the primitive configurations aretransient, during which the average height grows withtime (let us define the time as the number of injections).This linear growth cannot definitely last infinitely, andthe system saturates at some stage (becomes stationary),after which the average height becomes nearly constant,meaning that on average energy input and output arethe same. Recurrent states live in this regime. The totalnumber of recurrent configurations is det ∆. One of theimportant questions is how we can identify a configura-tion to be recurrent or transient. Fortunately, there aretests that do this for us. One of our tests is the lack offorbidden subconfigurations (FSC) which does not existin a recurrent configuration. FSCs can be identified sim-ply by the requirement that they can never be created bythe addition of sand and relaxation, if not already presentin the initial state [4]. We introduce two following testfor checking transient or recurrent configurations:We use the instability test of boundary site for the given configuration as follows: we add sand to any boundarysite of configuration and we allow the system to evolveuntil achieving stable configuration. If the new config-uration is the same as the first configuration, then theoriginal configuration has been recurrent. The Secondtest is called the burning test. This test is a dynamicone in which we burn sites one by one. First, we assumethat all the sites are unburned. In the next step, we burnall the sites whose heights are higher than the number oftheir unburned neighbors. If all the sites are burned atthe end of the burning process, the desired configurationis recurrent [4 ? ]. B. Other Sandpiles
1. Manna Model
The BTW model is a representative member of theBTW universality class. A relevant question is how onecan change the details of this model to change its uni-versality class. Stochasticity is one candidate to do this,which was tested for the first time by Manna by introduc-ing a two-state SOC system, known also as the Mannamodel [122]. It includes randomness in the toppling rule,i.e. in a two-dimensional system if a site has more thanone sand, it is unstable and topples. During a toppling,a direction is chosen randomly (each direction is chosenwith the probability of 0 .
5) and all of the grains are dis-tributed in this direction (no sand grain is transferred tothe other direction). This model is interpreted as a re-alization od a system with particles experiencing a localinfinite repulsive force between each other.Much attention has been paid to identify whether theManna model belongs to the BTW universality class ornot [123–131], which is still open. The most challengingproblem to this end is the accurate determination of ex-ponents (controlling the finite-size effects, controlling thenoise etc.).
2. Zhang Model
The Zhang model [132, 133] indicates the continuousstate of the BTW model, and the height of each site isconsidered to be a real number. The dynamics govern-ing this model is such that at any moment, a continuousand random value between zero and one is added to thesite. If its height exceeds a critical value, then that siteis unstable, and its value is distributed to nearby sites,and its grain content becomes zero. This continues un-til all sites become stable. This model does not havean Abelian property, because the amount transferred ineach toppling to neighboring sites depends on the initialvalue [4].
3. The General Abelian Sandpile Model
In sandpiles, the Abelian property is that the order oftopplings in an avalanche dose not matter, both the in-terchanged topplings reaching to the same configuration.Let us describe a general set up for the Abelian sandpiles.We consider the model on a graph with N sites labeledby integers i = 1 , , ..., N . We express the height of eachsite i with z i , which is a positive and integer number,and assign a threshold height z ci for each sites. At eachtime step, a site is chosen randomly and one sand grain isadded to it. We are also given an integer N × N topplingmatrix ∆, and a set of N integers { z i,c } , i = 1 , , ..., N .If for any site i , z i > z ci , then the site is unstable and ittopples. If z i > z ci then z j → z j − ∆ ij , for every j . Wemay choose z ci = ∆ ii , for which the allowed values of z i in a stable configuration are 1 , , ..., ∆ ii . Evidently thematrix ∆ has to satisfy some conditions to ensure thatthe model is well behaved [4]. . ∆ ii > i (otherwise toppling is never termi-nated). . For every pair i (cid:54) = j , ∆ ij ≤
0. This condition is re-quired to establish the Abelian property. . (cid:80) j ∆ ij ≥ j (this condition states that sandis not generated in the toppling process). . There is at least one site i such that (cid:80) j ∆ ij >
0, calleddissipative sites.
4. Oriented/Directional Abelian Sandpile Model
If we define an Abelian sandpile model on a directionallattice (e.g. the tilted square lattice), then we have ori-ented the Abelian sandpile model. In this model, themovement of the sand is in a certain direction and thethreshold height is commonly considered to be unity. Byadding one sand in a random site, after which the numberof sand grains reaches to two (becomes unstable), thensand moves randomly to one of the bottom sites [4]. Theexponents here depend on the direction of the propaga-tion, i.e. the time direction (top to bottom), or spacedirection (left to right).
IV. RELATION TO THE LOGARITHMICCONFORMAL FIELD THEORY
Using the number of recurrent states in sandpiles(which is det ∆ as stated above), a connection is es-tablished with the free ghost field. From the proper-ties of the Grassmann algebra, one can easily shown thatdet ∆ = (cid:82) (cid:81) Li =1 d θ i d¯ θ i exp (cid:2)(cid:82) d z∂θ ( z ) ¯ ∂ ¯ θ (¯ z ) (cid:3) , where θ i and ¯ θ i are independent Grassmann variables. Thereforethe connection to the free ghost field is established de- fined by the following action S = (cid:90) d z∂θ ( z ) ¯ ∂ ¯ θ (¯ z ) = 12 π (cid:90) ε αβ ∂θ α ¯ ∂θ β , (7)where the pair of free grassmanian scalar fields are de-fined as θ α = ( θ, ¯ θ ), and (cid:15) αβ is the cononical symplecticform, ε = +1, and ε αβ = − ε αβ . Using the q -state Pottsmodel in the limit q →
0, Majumdar and Dhar showedthat ASM is equivalent to c = − q -state Potts model decays with distance r like r − x T , where [135] x T = 1 + y − y (8)where y ≡ π cos − (cid:0) √ q (cid:1) . For q → x T = 2,which is a known result for height-height correlation issandpiles [134]. Other correlations and probabilities insandpiles can be expressed in terms of the grassmannfields and theirs derivatives [3]. A. Ghost Free Fields and W -Algebra In this section, we turn to the continuum limit of thesandpile model, which is a nonunitary field theory. It isshown in [134] that it is a c = − logarithmic is used here due to appearing logarithmic correlations inCFT, i.e. LCFT. In these theories, each primary filed hasa logarithmic partner, which enters in the ordinary op-erator product expansions (OPE) [136]. It can be shownthat this is equivalent to adding a nilpotent exponent tothe conformal dimension of primary fields [137, 138]. Inthis way, much of what we have seen about a commonconformal field theory can be easily generalized to LCFT.The existence of ghostly fields in a field theory means theexistence of a state with a negative magnitude, or in theother words, the theory in question is nonunitary. Inordinary CFT, the OPE of the energy-momentum ten-sor and a primary field gives singular terms, generatingthe corresponding conformal family for which one usesthe Virasoro algebra [139]. In LCFTs, in addition to theprimary field φ , we have a logarithmic partner ψ , whichsatisfies the following OPE with the energy-momentumtensor T [136] T ( z ) φ ( ω ) = hφ ( ω )( z − ω ) + ∂φ ( ω ) z − ω + ...T ( z ) ψ ( ω ) = hψ ( ω )+ φ ( ω )( z − ω ) + ∂ψ ( ω ) z − ω + .... (9)where h is the conformal dimension of φ . The action ofthe Virasoro operator at the zero levels L on the pair φ and ψ has a jordan form, giving rise to logarithmicterms in the corresponding correlation functions. Theinfinitesimal transformation of the pair is δ ε φ ( z ) = ( h∂ z ε + ε∂ z ) φ ( z ) δ ε ψ ( z ) = ( h∂ z ε + ε∂ z ) ψ ( z ) + ∂ z εφ ( z ) , (10)The calculations become simpler if we use the mixed op-ertor φ ( z, λ ) = φ ( z ) + λψ ( z ) where λ is a nilpotent num-ber ( λ = 0). Then the above equations cast to thefollowing abbrivated form δ ε φ ( z, λ ) = (( h + λ ) ∂ z ε + ε∂ z ) φ ( z, λ ) . (11)using of which one obtains the finite transformation ( z → ω ( z )) φ ( z, λ ) = (cid:18) ∂ω∂z (cid:19) h + λ φ ( ω, λ ) . (12)One can derive the two point function of the mixed fieldsthat is invariant under the translation, scale, rotationand special conformal transformations (cid:104) φ ( z , λ ) φ ( z , λ ) (cid:105) = a ( λ , λ )( z − ω ) h + λ + λ , (13)where a ( λ , λ ) = a ( λ + λ ) + a λ λ , resulting to (cid:104) φ ( z ) φ ( z ) (cid:105) = 0. Similar results can be derived for oth-eer two point correlations, and also the three point func-tions (cid:104) φ ( z , λ ) φ ( z , λ ) φ ( z , λ ) (cid:105) = f ( λ , λ , λ ) z a z a z a , (14)where a ij = h i + h j − h k + λ i + λ j − λ k f ( λ , λ , λ ) = (cid:88) i =1 c i λ i + (cid:88) ≤ i 0) and the OPE T ( z ) ψ α ( ω ) = θ α ( ω )2( z − ω ) + φ α ( ω ) + ψ α ( ω )( z − ω ) + ∂ψ α ( ω ) z − ω + ... (23)which is anomalous OPE, having an extra singular term(the first term in the right hand side). This make R representation more complex than R , and its opera-tor content is larger. In fact W -algebras come about inthis representation which works with operators at levelthree [141] W + = ∂ θ∂θW = ( ∂ θ∂ ¯ θ + ∂ ¯ θ∂θ ) W − = ∂ ¯ θ∂ ¯ θ. (24)Note that this set is isospin one with spins +1, 0 and − , , (1 , 1) and 0 , W i ( z ) W j ( ω ) = g ij ( z − ω ) − T ( ω )( z − ω ) − ∂T ( ω )( z − ω ) + ∂T ( ω )( z − ω ) + 4 T ( ω ) z − ω + ∂T ( ω ) z − ω − ∂T ( ω ) z − ω ) − f ijk ( W k ( ω )( z − ω ) + ∂W k ( ω )( z − ω ) + ∂ W k ( ω ) z − ω + 125 ( T W k )( ω ) z − ω )(25)Where g ij is the metric on the isospin one representa-tion, g + − = g − + = 2 and g = − 1, and f ijk are thestructure constants of SL(2). One can write the W alge-bra, by using the above OPEs. Following Gaberdial andKausch [141], we have: (cid:2) L m , W in (cid:3) = (2 m − n ) W im + n (cid:2) W im , W jn (cid:3) = g ij (2( m − n )Λ m + n + ( m − n )(2 m + 2 n − mn − L m + n − m ( m − m − δ m + n )+ f ijk ( (2 m + 2 n − mn − W km + n + V km + n )(26)Where Λ =: T : − ∂ T and V a =: T W a : − ∂ W a are quasiprimary normal ordered fields. This W alge-bra is different from Zamolodchikov’s W algebra [143],because f ijk is different. By using the above W algebraand Gaberdiel and Kausch [141] found the null vectors,from which, using the zero mode the following equationfor highest weight field φ were found [141, 144]: L (8 L + 1)(8 L − L − φ = 0 . (27)implying that h must be from the set (cid:8) , − , , (cid:9) whichare represented by V h . For V we have L φ = 0 which issatisfied for I and ˜ I that is the only logarithmic highestweight representation of c = − 2, and the other loga-rithmic representations are not highest weight (see R representation for example). Two of the three remaininghighest weight representations are related to the twistedsector and the other is nontwisted. For a more completereference see [145]. V. APPLICATION OF SANDPILES TONATURAL PROCESSES In this section, we turn to the application of SOC con-cepts to real systems. All subjects considered here aresimulations (except a part for SOC in clouds). Mostparts of this section are carried out by the authors ofthis paper. A. SOC in Fluid Propagation in Porous Media The fluid propagation in porous media (involvingavalanche-type dynamics) is a complex procedure, result-ing in various interesting patterns. The avalanches arisefrom a non-linearity in the laws governing the dynamics of fluid in this system, namely the critical water satu-ration . The fluid in a typical region of the reservoir isstatic (it does not macroscopically transfer to the neigh-boring regions), until the accumulated fluid saturation inthat region exceeds a certain saturation known as criti-cal saturation S C above which the fluid overflows freely(governed by the pressure gradient) to the neighboring re-gions. Therefore, for a region (consisting of many pores),the total water saturation is small, no fluid transfers tothe other regions. Once the amount of fluid in that re-gion reaches S C , then some small droplets of fluid in thepores aggregate so that the fluid acquires the ability tomove to the neighboring regions. In terms of the porousmedia parameters, one common choice for the relationbetween the relative permeability of phase α (shown by k rα ) and the saturation of that phase ( S α ) is k rα = (cid:26) S α − S C S α ≥ S C S α < S C (28)which shows the above-mensioned dynamics. In [120], itwas shown that the set of Darcy equations (known as thereservoir flow or RF model) for two-phase propagationin porous media is very similar to the BTW dynamics,except that the former is directional (the fluid movesin the direction of the pressure gradient), whereas thelatter is not. In this work the ordinary BTW model hasbeen used defined on a percolation lattice which realizesa porous media (uncorrelated) in which some pointsare occupied with a probability of p and the othersare unoccupied (with the probability 1 − p ). Fig.(1)shows the toppling rule for a typical site that has threeoccupied neighbors. To find the connection betweenthis model and the Darcy model, the Schramm-Loewnerevolution theory (SLE) theory was used and some otherstandard statistical analyses were performed such asfractal dimension on their domain-walls to investigateits behavior in terms of p . The results of the SLEtheory are depicted in Fig.(2) where ξ t , known asthe driving function in the SLE theory is a continuousreal-valued function that is shown to be proportionalto the one-dimensional Brownian motion ( ξ t = √ κB t , κ being the diffusivity parameter) for the conformalinvariant curves [146]. Two critical models with thesame diffusivity parameter are believed to be in a sameuniversality class. This graph shows that the RF and theSOC models are equivalent on the percolation threshold p = p c , being compatible with the Ising universality class.Evidently the permeable pores of the porous mediaare not completely independent in the natural systems.As an example, consider the sedimentation process ofthe reservoir rock in which some parts of the reservoirrock become impermeable to flow. Therefore, the corre-lation of the unoccupied sites depends on the dynamics ofthe sedimentation process. In [147] it was assumed thatthe correlation of the occupied sites is given by a zero-magnetic field, ferromagnetic Ising model. The strengthof the correlations is controlled by the Ising coupling con-FIG. 1: SEC. V:(a) The toppling rule for a typical site which has three occupied neighborsFIG. 2: SEC. V :(a) The plot of (cid:104) log( l ) (cid:105) versus (cid:104) log( r ) (cid:105) (fractal dimension) for the reservoir flow model and SOCmodel. (b) (cid:104) ξ (cid:105) − (cid:104) ξ (cid:105) versus t for the case p = p c stant and the artificial temperature T . This model showsthat at the Ising critical temperature T c , the model iscompatible with the universality class of two-dimensional(2D) self-avoiding walk (SAW). The mixing of two con-formal symmetric models is also of interest to the theoret-ical side. Mathematically this problem can be tracked interms of the Zamolodchkivs c -theorem, in which knowingthe scaling perturbing filed, one can obtain the change ofcentral charge of the conformal field theory δ c = c IR − c UV [148]. The effect of the second CFT model can presum-ably be codded in a scaling field from the operator con-tent of the original (first) CFT model.In another work [149] Najafi et al. concentrated onthe SOC dynamics on 3D correlated porous media. Itwas found evidence for a new nonequilibrium universal-ity class that is reached by changing the geometry of theunderlying graph upon which the model is defined. This might be applicable to experiments with spatial flow pat-terns of transport in heterogeneous porous media [150]. B. SOC in Cumulus Clouds As stated in the previous sections, there are manypapers reporting on the fractal structure of the clouds,some of which are based on SOC. Convective clouds(such as cumulus clouds as a member of the cumuliformclouds) develop in unstable air due to the buoyancyforce, resulting from water vapor, supercooled waterdroplets, or ice crystals, depending upon the ambienttemperature. Cumulus clouds have flat bases and areoften described as cotton-like low-level clouds, less than2 km in altitude unless they are more vertical (cumuluscongestus form). These clouds which may appear in0FIG. 3: SEC. V B:(a) The log-log plot of trace lengths l in terms of L (the box linear size). The dashed line is alinear fit with slope D af = 1 . ± . R in terms of N , and the lower insetis a semi-log plot of the loop green function in terms of r , with the exponent ν = 0 . ± . 01. (b) Ensemble averageof log l in terms of log r ( (cid:104)(cid:105) means the ensemble average) with slope D bf = 1 . ± . 02. The log-log plot of thedistribution function of r and l are shown in the upper and lower insets , with exponents τ r = 2 . ± . 03 and τ l = 2 . ± . 02 respectively.lines or in clusters (producing little or no precipitation)are often precursors of other types of clouds, suchas cumulonimbus when influenced by weather factorssuch as instability, moisture, and temperature gradient.When they grow into the congestus or cumulonimbusclouds, they are more probable to precipitate. Theheight of the cloud (from its bottom to its top) dependson the temperature profile of the atmosphere.Cumulus clouds form via atmospheric convection asair warmed by the surface begins to rise, resulting in thetemperature decrease and humidity rise. At a threshold,named as lowest condensation level (LCL), in which therelative humidity reaches 100%, then condensation to thewet phase (known as wet-adiabatic phase) starts. Thereleased latent heat (due to condensation) warms up theair parcel, resulting to further convection. At LCL, thenucleation process starts on various nuclei present in theair. The process of formation of raindrops and rainfallhas been explained successfully by Langmuir [116].Although the liquid water density within a cumuluscloud changes with height above the cloud base [121](for the non-precipitating clouds the concentration ofdroplets ranges from 23 to 1300 droplets per cubiccentimeter [151]), the density can be thought of asbeing approximately constant throughout the cloud theheight of the cumulus clouds depends on the amountof moisture in the thermal that forms the cloud, andhumid air will generally result in a lower cloud base.In stable air conditions in which their vertical growthis not that high, they are considered to be effectivelytwo-dimensional.In places, cumulus clouds can have holes where thereare no water droplets [151]. This fact causes to createthe fractal structures that it provides powerful tools to classify them in terms of the circumstances in which theyform. The fractal structure of clouds has been reportedin some previous works [152–154].The self-organized criticality in the atmosphere andclouds was first detected by Peters et al. by analyzingthe precipitation [91]. They used satellite data and de-fine a critical value of water vapor as a tuning param-eter, and precipitation as the order parameter shows anon-equilibrium continuous phase transition to a regimeof strong atmospheric convection and precipitation.In an unpolished work, we uncovered the SOC state ofclouds directly by analyzing some earth to sky images ofcumulus clouds under fair air conditions. The analysisof the level lines of two-dimensional cloud fields stronglysuggests that these are in a SOC state, which is alsoconfirmed by a Schramm-Loewner evolution (SLE) anal-ysis, finding that these belong to the sandpile univer-sality class, i.e. c = − L is the box linear size in the box counting method, l ∼ L Df a is the trace length, R is end-to-end distance, N is the number of steps along the trace, (cid:112) (cid:104) R (cid:105) ∼ N ν ,and (cid:104) log l (cid:105) = D bf (cid:104) log r (cid:105) in which r is the gyration radiusof the loop. Their analysis show D af = 1 . ± . D bf = 1 . ± . ν = 0 . ± . 01 where ν = D f . Alsothis analysis show that the two-dimensional images of theclouds are very close to the loop-erased random walkers(LERW) traces with D af = 5 / C. Vibrating Piles A very important question concerning the SOC modelis its stability against external manipulations, like vibra-tions. As stated above, some experiments have also been1FIG. 4: SEC: V D: Scheme showing the invasion natureof the two-species sandpile model. In (1) the left (blue)site becomes unstable since h r > h th . In (2) howeverboth grains are lower than h th for the right (white) site,but h r + h b > H . In this situation r or b is randomlychosen for toppling, here b is chosen. Then in (3) the b content of the upper (dark yellow) site is increased.Therefore, effectively r has invaded b and pushed ittowards another site.done to test the critical properties of real piles in thepresence of external vibrations [155] which were mod-eled and simulated [156]. Recently we have simulatedthe BTW model under vibration conditions, affecting thetoppling rules in one direction, namely x -direction. Tothis end, we manipulated the toppling rules, which de-pend on time. The toppling matrix was considered to be∆ ( i,j ) , ( i (cid:48) ,j (cid:48) ) = n i = i (cid:48) , j = j (cid:48) − n i = i (cid:48) , j = j (cid:48) ± − n (1 + (cid:15) sin ωt ) i = i (cid:48) + 1 , j = j (cid:48) − n (1 − (cid:15) sin ωt ) i = i (cid:48) − , j = j (cid:48) ω = 2 π/T is the angular frequency, T is the timeperiod, and (cid:15) is the vibration strength parameter. Thistoppling rule states that the system is vibrating in the x direction making the model anisotropic. The proper-ties of the model were investigated in terms of T and (cid:15) .We uncovered that the exponents run with ω and (cid:15) . Im-portantly increasing the strength of vibrations makes theavalanches more smooth with exponents that depend onthe (time and space) directions. D. Invasion Sandpile Model In the previous section, we provided evidence that theDarcy model has similarities with the BTW model atthe critical occupation. Invasion percolation (IP) is an-other natural choice, in which one phase invades the otherphase towards the production well. In fact, IP [157] is a standard model to study the dynamics of two immisciblephases (commonly denoted by wet and non-wet phases)in a porous medium [158, 159]. During this process, thewet phase invades the non-wet phase, and the front sepa-rating the two fluids advances by invading the pore throatat the front with the lowest threshold [159]. This modelsuffers the lack of the notion of critical saturation intro-duced in the previous sections. In a parallel line of think-ing, one may try to make the BTW model two-phase withthe extra tool of invasion.Let us consider a L × L square lattice and initiallyassign two random integers, h r and h b , to each site,uniformly from the interval { , , , ..., h th } , h th beingthe threshold for a one species. h r and h b are thenumber of red and blue sand grains in our two-speciessandpile model, representing the two (wet and non-wet)phases in the reservoir. The reported results are inde-pendent of the value of h th , so we set it to 20. A site i is considered stable if three conditions are fulfilledsimultaneously. The first two are the standard ones forthe one-species sandpile model, namely, h r ( i ) ≤ h th and h b ( i ) ≤ h th , where h th represents the CFS. The thirdone is h r ( i ) + h b ( i ) ≤ H , where H < h th is the secondthreshold. This additional condition is motivated by thefact that in the non-linear Darcy equations, there is anauxiliary equation expressing that the sum of two-phasesaturations S w + S o is a constant that depends onthe capillary pressure. Thus, a site i is unstable andtopples if at least one of the following conditions are met:C1: h r ( i ) > h th ,C2: h b ( i ) > h th ,C3: h r ( i ) + h b ( i ) > H The dynamic goes as follows. Initially, all h r and h b are chosen randomly from a uniform distribution,such that no site is unstable. Then, iteratively, we firstchoose a species (either r or b , with equal probability)and a site i at random to add a particle of that species,i.e. h x ( i ) → h x ( i ) + 1 where x is the selected type.If that site becomes unstable, it topples, accordingto the following rule: If condition C1 is met, then h r ( i ) → h r ( i ) − h r ( j ) → h r ( j ) + 1 where j isthe neighbor of i with the lowest red-grain content.If condition C2 is met, then h b ( i ) → h b ( i ) − h b ( j ) → h b ( j ) + 1 where j is the neighbor of i with thelowest blue-grain content. If condition C3 is met, then h x ( i ) → h x ( i ) − h x ( j ) → h x ( j ) + 1 where j is theneighbor of i with the lowest x -grain content, and x israndomly chosen to be r (red) or b (blue). As a resultof the relaxation of the original sites, the neighboringsites may become unstable and also topple. Therefore,the toppling process is repeated iteratively until all sitesare stable again. This collective relaxation is calledan avalanche . The sand grains can leave the samplefrom the boundaries, just like in the ordinary BTWmodel [57]. Note that, with two species, an avalanche2FIG. 5: SEC: V D: (a) (cid:104) log l (cid:105) in terms of (cid:104) log r (cid:105) for which the relation gives the fractal dimension D f . Upper insetshows r ∗ , the separator of the small and large scale regimes, in term of L − . (b) Distribution function of l r (redgrains). Bottom inset shows l ∗ ,the separator of the small and large scale regimes, in term of L − . Along with thedefinition of the crossover points by plotting R + R in terms of the tentative crossover point bottom and upperinset for (a) and (b), respectively.TABLE I: SEC: V D: The β and ν exponents from thefinite-size analysis, τ , and τ for m , s , l , and r corresponding to the red avalanches. The last rowcontains the exponents for the 2D BTW model for thesake of comparison with τ ( L → ∞ ) [119, 161]. quantity m s l rτ ( L → ∞ ) 1 . ± . 04 0 . ± . 03 2 . ± . 03 3 . ± . τ ( L → ∞ ) 1 . ± . 02 1 . ± . 04 1 . ± . 03 1 . ± . β − − . ± . 05 1 . ± . ν − − . ± . 03 0 . ± . τ 2D BTW . ± . 01 1 . ± . 01 1 . ± . 03 1 . ± . of one species might trigger an avalanche of the otherone, see example in Fig. 4 and details in the caption.The reason that we call this invasion is that here onespecies pushes the other one due to the finite capacityof the pore, i.e. the total volume of the particles cannotexceeds a threshold (see C3), as in real situations. Inthe Darcy reservoir model, C3 is an auxiliary equation,where H plays the role of the maximum finite saturationthat is possible in a pore [120, 160] and is the source ofthe invasion in the invasion percolation model [157].This model has two types of avalanches: one-speciesand two-species avalanches. The first one involvesonly the redistribution of grains of one species. Inthe second, there is mass transport of the two species.The symmetrical one-species avalanches have the sameresults for blue and red avalanches. Also, it has twodifferent regimes: for large avalanche sizes the fractaldimension D (2) f is consistent with 5 / D (1) f = 1 . ± . 02 that it is shown in Fig.(5). For TABLE II: SEC: V D: The exponents β , ν , τ , and τ for m , s , l , and r corresponding to the two-speciesavalanches. quantity m s l rτ ( L → ∞ ) 0 . ± . 05 0 . ± . −− . ± . τ ( L → ∞ ) 1 . ± . 02 1 . ± . 03 1 . ± . 03 2 . ± . β − − . ± . 05 1 . ± . ν − − . ± . 03 0 . ± . extract two regimes the crossover point was used the R test [163] ( R = R + R ) that it is shown in bottominset of Fig.(5a), and the upper inset shows crossoverpoint of fractal dimension as r ∗ in term of 1 /L where L is system size. For all measures (gyration radius( r ), avalanche size ( s ), and avalanche mass ( m )) isseen this behavior. Table.(I) shows all exponent forthis measures in one-species avalanches regime. Forexample in Fig.(5b) is seen distribution functions ofloop length for red grains( l r ) same Fig.(5a). In the twospecies avalanches, the regime fractal dimension hastwo-state. The fractal dimension of large avalanches, D f = 1 . ± . 01, consistent with 2D BTW model.However, the fractal dimension for small avalanches D f = 1 . ± . 01, which is different from the exponentfound for one-species avalanches. Table.(II) shows allexponent for this measures in two-species avalanchesregime. For further information see [164]. E. SOC in Exitable Complex Networks Probably the most important reason for considering of3FIG. 6: SEC: V E: (a) The plot of M in terms of R and the corresponding exponents γ M R . Upper inset: γ UV M R and γ IR M R . Lower inset: the cross-over radius R ∗ in terms of α with the exponent 0 . ± . 03. (b) The same as (a)for various lattice sizes L . The finite size dependent (UV and IR) slopes γ M R have been shown in the inset. M M (0) R R (0) n toppling n toppling (0) τ ( α = 0) 1 . . τ . . . . . . τ ( α = 1) 3 . . . . γ τ . . . . α = 1) 2843 4601 11 . . γ cut . . . . . . TABLE III: SEC. V E: The asymptotic values of the exponents. For each quantity there is a ”cut” value in which across over between small scale behavior (which are consistent with regular 3D BTW model) and large non-universalbehavior occur. It has been found these cut-values scale with α in a power-law fashion. For example M ∗ ≡ M cut3 = m cut ( α = 1) α − γ cut which have been shown by ”cut” in the table. In contrast to τ , τ runs with α forall quantities; τ = τ ( α = 1) α + γ τ which have been shown separately in the table. After [165].the BTW dynamics on top of the complex networks hasbeen the important observation of Beggs et. al. [166] inwhich it was shown that the propagation of spontaneousactivity in cortical networks is self-organized criticalphenomena governed by avalanches much similar to theBTW model. The focus of researches in this area hadbeen on the structural and functional properties of ran-dom lattice models. One of the most challenges in thesesystems is finding the circumstances under which thesystem shows the critical behaviors [ ? ]. In order to mon-itor critical behaviors, different time series are usuallyanalyzed in which power-law behavior is expected [105].Therefore one may be encouraged to investigate timeseries of topplings in the SOC model on various randomlink lattices, like scale free networks [167, 168], multiplexnetworks [167], optimized scale-free network on Eu-clidean space [169], Watts-Strogatz small-worlds [170],directed small-world networks [171], and on Scale FreeNetworks with preferential sand distribution [172].An example of the network that is embedded in theEuclidean space is a graph with random links with finite range interaction (RLFRI), i.e. two nodes are connectedif their distance is smaller than a control parameter R .Therefore the topology of this graph is tuned by n and R , where n is the number of links per site. Then dynamicis defined via the following toppling matrix∆ i,j = z i if i = j − i and j are connected0 other (30)where z i is the degree of node i .The other system of interest is the small world net-works that are interpolation between regular and randomnetworks. In addition to the regular links between neigh-bors, there are some long-range links between random-chosen sites. For this system, we have two dependentfurther random fields in addition to the height field. Theconnection matrix L( i, j ) is unity if sites i and j (notneighbors) are connected by a long-range link and zerootherwise. The distribution of lengths and the degree ofnodes are chosen to be uniform in the interval of allowed4values (naturally the lengths are restricted to the linearsize of the system). The other one is z c ( i ) = 6+ (cid:80) j L( i, j )which accounts for the number of total links in the node i . In this language if the height of a node exceeds z i ittopples according to the rule h ( i ) → h ( i ) − ∆ i,j in which:∆ i,j = − i and j are neighbors or L( i, j ) (cid:54) = 0 z i i = j α (the percent of long range links) is defined as follows: α ≡ × / (cid:80) i,j L( i, j ) is to prevent double counting.Fig. 6(a) shows the plot of (cid:104) log( M ) (cid:105) in terms of (cid:104) log( R ) (cid:105) whose slope is γ M R ≡ D M F which is the 3D mass frac-tal dimension for L = 300 and various α ’s. We notethat D M F ( α = 0) (cid:39) . ± . 02 [173]. Interestinglyit is seen that the graphs smoothly cross over to thelarge scale regions in which the slope (fractal dimension)( m IR ≡ γ IR M R ) is different from the slope in the small-scale region with the slope m UV ≡ γ UV M R . We namethe small scales as UV limit and the large scales as IRlimit . The point of this behavior change depends on α .This point can easily be calculated using the linear fitof the graphs in each individual region. The transitionpoint ( R ∗ ) is simply the point in which the fits meet eachother. The fact that m UV ( α ) is nearly α -independent and m IR ( α ) runs crucially by varying α can be seen in the up-per inset of Fig. 6(a). We interpret R ∗ as the point atwhich the cross-over takes place to the large scale prop-erties since for r (cid:46) R ∗ the results are very close to theregular BTW model, whereas for r (cid:38) R ∗ the behavior isdifferent and not universal (presumably mixed with thefinite-size effects). More interestingly we have observedthat R ∗ is a decreasing function of α , i.e. R ∗ ∼ α − ζ in which ζ = 0 . ± . 05. Since α can be interpretedas the measure of how directly a randomly chosen siteis connected to a boundary site at which dissipation oc-curs, we can say that effectively (on average) a fraction ofgrains are dissipated in a bulk toppling depending on theamount of α . The other effect is the sink role of connec-tion sites in micro-avalanches. In fact when α increases,the probability that a micro-avalanche involves a site thathas a long-range link to the other micro-avalanches in-creases. Since roughly speaking, such sites play the roleof sink points, one may expect that the effective modelfor micro-avalanches is a dissipative one. It is known thatthe dissipative BTW model is equivalent to the massiveghost action S = (cid:90) d z (cid:18) ∂θ ¯ ∂ ¯ θ + m θ ¯ θ (cid:19) (33)where θ and ¯ θ are complex Grassmann variables and m is the number of sand grains dissipated in each toppling ( m can be fractional). On the other hand it isknown that R ∗ ∼ m − [174]. From these two points,one concludes that effectively our model is equivalentto the dissipative BTW model with m ∼ α . Thiscorrespondence is acceptable only for r (cid:46) R ∗ and showsthat the large scale regime is directly affected by thedissipations in the boundary sites and the finite-sizeeffects. This result is reasonable since the amount ofgrain dissipation in a single component of an avalanche(the number of sand grains that are transferred out ofthat area) is proportional to the number of nodes withlong-range links in that area. The problem descriptionis not however as simple as stated above since thereare surely some other links that return energy to theoriginal micro-avalanche which partly compensatesthe dissipation effects. It is worth noting a commentconcerning the numerical value of ξ which is claimedto be 1 /d ≈ . 33 in three dimensions [175] which istrue for the total avalanches. For micro-avalancheshowever, the statistics is different and it is acceptablethat R ∗ should be proportional to m − ∝ α − / which isrepresentative of the grain dissipation towards the othermicro-avalanches and is a well-known property of thedissipative sand-pile models. We have also consideredthe finite size effect of the results which have been shownin Figs. 6(b). The constant trend of m UV is seen in theinset of this figure, in which it is seen that its numericalvalue (for α = 1) is nearly robust against varying latticesize L , whereas m IR changes considerably by lattice size.This reflects the universal behavior of m UV . F. SOC in Imperfect Supports There are many phenomena in nature in which thedynamics are defined on the percolation lattice whichprovides motivation for studying this field. Some ex-amples are the recent experiments in which the voidsof percolating clusters were filled by (commonly mag-netite) nano-particles of ferromagnetic fluids [176–179],the loop-erased random walk on the percolation lattices(which is related to watersheds [180]) and fluid propaga-tion in porous media [157, 181–183]. Mixing two statisti-cal models (one as the dynamical model and the other asthe host for the first one) may be interpreted as the in-terplay between two statistical models from which somenew non-trivial critical behaviors can emerge. One ofthe important questions in the contex of self-organizedcriticality is the effect of stochasticity on its critical be-haviors, and on the universality class of the determinis-tic model. This stochasticity can be annealed (like theManna model), quenched where the disorder is pinned tolattice points. For the latter case, consider a situationin which some points of the lattice is inaccessible for thesand grains to pass through, namely impermeable sites.One may distribute the impermeable cores throughout5FIG. 7: SEC: V F: (a) Fractal dimension of the curves for p = p c and L = 1024. Inset shows finite size effect for thefractal dimension, i.e. D p = p c f ( L ) = 1 . − . L ). (b) Finite size effect for the fractal dimension for p > p c . Insetshows the dependence of γ ( L ) on L the lattice using a single parameter p : each site is per-meable with the probability p , and impermeable withthe probability 1 − p . Then the host system is simplydescribed by the percolation theory. In this case, let usdefine z i as the number of active neighbors of the site i . Then the sand grains are governed by the followingtoppling matrix:∆ i,j = z i i = j − i, j are neighbors, and j is permeable0 other (34)This toppling role has been shown in Fig.(1). This modelwas shown to have the same properties as the Darcymodel in 2D (square) lattices at p = p c , whereas they aredifferent for other p values [120, 184]. Note that the sandgrains are allowed to move on the spanning percolationcluster, that connects two opposite boundaries. In twodimensions, we first construct the percolation lattice asstated in the previous section for various rates of p andthen study the BTW model on the percolation clusters.We start from random h configurations for each amountof p and add sand grains randomly throughout thesample. In this case is seen that after a number of sandgrain injections that it depends on the amount of p , thesystem reaches the steady-state, e.g. the average heightin the percolation lattice (¯ h ) becomes constant.Each avalanche has an exterior boundary that formsa loop. To identify these loops, we label each siteof the system in each toppling process. it is whiteif that site is unoccupied or untoppled and black ifthat site is occupied and toppled. Then, the loops arewell-defined as the separators of the black and whitesites. For these avalanches, we define gyration radius, loop length, cluster mass, and the number of topplingsin each avalanche, represented here by r, l, m , and n t respectively.The samples for p ≥ p c are self-similar and show criticalbehaviors, e.g. the fractal dimension is well-defined forthem. The dependence of the fractal dimension on p hasbeen shown in figure 7 in which the finite-size scalingargument has been presented. We also directly observedthat the largest system size, the fastest BTW-like( p = 1) behaviors starts. This shows that in the thermo-dynamic limit the behavior of BTW-like ( p = 1) is thedominant behavior for all p c < p ≤ 1. This hypothesishas also been confirmed by further analysis of e.g. thedistribution function of r, l and m . We find that thefractal dimension of avalanche exterior boundary (loop)is compatible with the fractal dimension of the exteriorboundaries of the geometrical spin clusters of the 2DIsing model,i.e D p = p c f ( L = ∞ ) (cid:39) D Isingf = .In this part along with the three-dimensionalanalysis, we study the energy propagation throughtwo-dimensional slices, i.e. cross-sections of the three-dimensional system. The problem of two-dimensionalpropagation of sand grains (energy) in three dimen-sional systems seems to be very important from boththeoretical and experimental sides. More precisely theimportant question in the theoretical physics is howthe information in d + 1dimensions would be reflectedin its d dimensional sub-system. For this purpose oneshould map the original d + 1 dimensional model to ad-dimensional one and measure how some informationis lost and how the degrees of freedom in the subtracteddimension affect the d − dimensional model, i.e. whichmodel lives in the lower dimensional system. We hadtwo separate studies: critical p = p c and off-critical6FIG. 8: SEC. V F: (a) The order parameter ζ ( T ) in terms of T for various lattice sizes. Top inset: data collapse of ζ with exponents β = 0 . ± . 02 and ν = 0 . ± . 05 ( (cid:15) ≡ T − T c T c ). Lower-right inset: χ ( T ) ≡ ∂ζ/∂T in terms of T ,showing a peak at T c . Lower-left inset: average coordination number ¯ Z of the Ising clusters as a function of T forvarious systems sizes. (b) The critical exponent of the distribution function of gyration radius τ r in terms of T .Upper inset: the finite size dependence (linear in terms of L − ) of τ r at T c , revealing that τ T c r ( L → ∞ ) = 1 . ± . 03. Lower inset: the same finite size extrapolation of τ r for all temperatures. The exponentundergoes a clear jump at T c . p c < p pp c violating thehyper- scaling relations obtained for the critical case.For the 2D induced model in the off-critical regime, weshowed that there is a p-value ( p ∈ (0 . , . p Dc ≤ p < p Dc in the cross-sections does not have a thermodynamiclimit, whereas for p ≥ p Dc the system is identical to the p = 1 system. For details see [173, 185]The above analysis was for the case where the imperfec-tions are uncorrelated. Turning on the correlations forthe spatial configuration of the imperfections over thelattice makes the results different.Based on the determination of the fractal dimensionof the external perimeter of the avalanches, and theSchramm-Loewner evolution, it was suggested in [186]that the BTW on the critical percolation, results tocritical Ising universality class, and also the BTW modelon the Ising-correlated percolation lattice results toself-avoiding walk universality class.More interesting is the BTW model on the three- dimensional Ising correlated lattice, for which themagnetic phase transition is not accompanied by a per-colation transition, allowing us to measure the propertiesof the model across the transition point. In this case,we implement the dynamics of the Bak-Tang-Wiesenfeldmodel (BTW) [187], also known as the Abelian sandpilemodel, on a diluted cubic lattice. This lattice is com-prised of sites that are either active (through which sandgrains can pass) or inactive (completely impermeable tosand grains), which are labeled by the quenched variable s (called spin ) that is +1 for the active case, and − T ) to obtain the spin configuration, whichis expressed by the Hamiltonian H = − J (cid:80) s i s j ,where s i is the spin on site i , J > (cid:104) i, j (cid:105) means that i and j arenearest neighbors. The 3D Ising model undergoes amagnetic phase transition at T = T c ≈ . 51 for the cubiclattice. Since the spin clusters (two sites belong to thesame cluster if they are nearest-neighbors and have thesame spin) of the 3D Ising model on the cubic lattice percolate at any temperature , no percolation transitiontakes place at T c . This can be understood by notingthat the critical site percolation threshold for the cubiclattice is around 0 . < . T → ∞ of the Ising model). After constructing an Isingconfiguration at a given temperature using Monte Carlo,we implement the BTW dynamics on top of the spanning(majority) spin cluster (SSC), i.e. a cluster comprised ofspins with the same orientation connecting two oppositeboundaries of the lattice. Free boundary conditionsare imposed in all directions. In the BTW dynamics,we consider on each site i a height h i (the number ofsand grains) taking initially randomly (independently7and uncorrelated) with the same probability one integerfrom { , ..., Z i } , in which Z i is the number of activeneighbors of the i th site. Then we add a sand grain at arandom site i , so that h i → h i + 1. If this site becomesunstable ( h i > Z i ), then a toppling process starts,during which h j → h j − ∆ i,j , where ∆ i,j = − i and j are neighbors, ∆ i,j = Z i if i = j , and is zero otherwise.After a site topples, it may cause some neighbors tobecome unstable and topple, and so on, continuing untilno site is unstable anymore. Then another random siteis chosen and so on. The average height grows withtime until it reaches a stationary state after which thenumber of grains that leave the system through theboundary is statistically equal to the number of addedones. The dynamics can be implemented with eithersequential or parallel updating. Criticality in threedimensions also induces two-dimensional (2D) criticalproperties, which enables us to apply 2D techniques likeconformal loop ensemble theory [165, 173, 185, 188].Here we consider three-dimensional (3D) avalanches, aswell as their two-dimensional (2D) projections on thehorizontal plane.To characterize more precisely the two SOC phases, weanalyze the average number of topplings per site (top-pling density) in avalanches. We may define as the orderparameter ζ ( T ) ≡ f perc − f ( T ), where f ( T ) = m ( T ) N ( T ) , m ( T ) being the number of topplings, N ( T ) the num-ber of sites in the spanning (majority) spin cluster, and f perc ≡ f ( T = ∞ ). Fig. 8(a) reveals that ζ ( T ) is zerofor T > T c , and starts to grow continuously in a power-law fashion when T is decreased below T c signaling aphase transition at T = T c , at which χ ( T ) ≡ ∂ζ∂T showsa distinct peak. The finite size scaling relation for ζ is ζ L ( T ) = L − β/ν G ζ (cid:0) (cid:15)L /ν (cid:1) (see upper inset of Fig. 8(a))in which (cid:15) ≡ T − T c T c , G ζ ( x ) is a scaling function with G ζ ( x ) | x →∞ → x β , and β = 0 . ± . 02 and ν = 0 . ± . T → ∞ corresponds to a site percolation cluster with occupationprobability p = . For T > T c , therefore we will call thisphase SOC p = . The behavior in the T < T c region, how-ever, is dominated by T = 0 i.e. the regular lattice, andtherefore we call this phase BTW. At the transition pointall of these exponents show a sharp change. For exam-ple, τ r in the lower inset of Fig. 8(b) abruptly changesits value from BTW to SOC at T = T c , its value be-ing completely different for T < T c and T > T c . It isstep-like in the limit L → ∞ which is obtained by linearextrapolation in terms of 1 /L . We observed that for alltemperatures T < T c , τ r extrapolates to 1 . ± . 04, andfor T > T c , it is 1 . ± . 04. At T = T c , this exponent is1 . ± . 03 which is different from both values. G. Diffusive Sandpiles let us introduce the local smoothings. A local smooth-ing is defined as the action in which a site is chosenrandomly and is checked for a more stable configuration.To do this, the height of its neighbors is checked. Supposethat the selected site is i with nearest neighbors i , i , i and i . Among these neighbors, label i max as the sitein which E ( i max ) = Max { E i } i =1 , and i min as the site inwhich E ( i min ) = Min { E i } i =1 . A local smoothing is com-posed of two updates: δE ≡ int [( E ( i max ) − E ( i )) / i max to i , and then δE ≡ int [( E ( i ) − E ( i min )) / 2] grains (if positive) flowfrom the site i to i min , in which int[ x ] is the integerpart of x . In case of more-than-one sites having thesame (maximum or minimum) height, the site from/intowhich the grains flow is chosen randomly. If the site i is locally maximum (minimum), automatically no grainflows into (from) the site to the neighbors. No localsmoothing is applied to the unstable sites. Therefore,we have two kinds of relaxations: the sites which areunstable topple and the sites which are chosen for localsmoothings moderate their local height gradient. Wecall the first procedure as the toppling and the secondone as the local smoothing . This problem can also becalled a diffusive sandpile model , in which the grains arelubricated such that they have the chance to slip to theneighboring sites. more detail in [189]Let us consider the problem in the mean-field (MF)level. Consider a square lattice with N = L sites and4 L boundary sites. Suppose that the avalanche mass(number of toppled sites in an avalanche) at time T hadbeen A ( T ), and one energy unit is added at T + 1, andalso the average height at time T ( T th injection) is con-sidered to be ¯ E ( T ). When one energy unit is added tothe system, then the average energy is raised by 1 /N .But some energy is dissipated from the boundaries. Inthe mean-field level, we should first calculate the prob-ability that a boundary site had been unstable in theprevious avalanche. This probability is simply the num-ber of boundary sites (4 L ) times the probability that arandomly chosen site is involved in the avalanche. Thelatter is equal to A ( T ) /N . All in all, we reach the fol-lowing result for the time-dependent average energy:¯ E ( T + 1) = ¯ E ( T ) + 1 N − L A ( T ) N . (35)The above analysis reveals that in the steady-state inwhich ¯ E ( T + 1) = ¯ E ( T ), we have ¯ A ( T ) = L . Nowconsider the effect of local smoothing. We suppose thatits effect is decreasing the range of the correspondingavalanches A (cid:48) ( T ), and also suppose that A (cid:48) ( T ) is pro-portional to A ( T ) (the same avalanche in the absence oflocal smoothing). The proportionality constant is surely ζ -dependent, i.e. A (cid:48) ( T ) = f ( ζ ) A ( T ). In this case, underthe conditions that the zero-smoothing ζ = 0 system isin the steady-state, we have:¯ E ( T + 1) = ¯ E ( T ) + 1 N (1 − f ( ζ )) . (36)8FIG. 9: SEC. V G: (a) The (shifted) height average in terms of T (the number of injections) for various rates of ζ .Inset: ¯ E ¯ E in terms of ζ . (b) Mass as a function of time T for various rates of ζ and L .This means that ¯ E grows linearly with time (with theproportionality constant N (1 − f ( ζ )). This growth takesplace up to a time that ¯ E → E th at which the avalanchemass reaches the system size. In this case, A ( T ) ≈ N = L . In this case the first equation of this document showsthat ¯ E ( T + 1) = ¯ E ( T ) − (4 L − . (37)In this case, the average energy decreases abruptly by4 L − 1. To continue, we introduce the probability of¯ E ( T + 1) = z , conditioned to ¯ E ( T ) = M , which is foundto be: P ( ¯ E T +1 = z | ¯ E T = M ) = (cid:40) δ z,M + N (1 − f ( ζ )) M < E th δ z,M − L +1 M ≈ E th (38)One of the important quantities in the analysis of thedynamical systems is the branching ratio defined by b ( M ) ≡ E (cid:104) ¯ E T +1 M | ¯ E T = M (cid:105) . For a given M , if b ( M ) > E grows, and if b ( M ) < b ( M ) = (cid:40) − f ( ζ ) NM M < E th − L − M M ≈ E th (39)in which E [ ] is the expectation value. Note that when b ( M ) > b ( M ) < M , the aver-age number of grains of the system will increase (herelinearly) (decrease, here abruptly) with T . This rela-tion predicts that a bifurcation takes place at a non-zero ζ , above which some oscillations occur. For the firstbranch, the mean height increases linearly with T up tothe time at which ¯ E ≈ E th . At this point, the averageheight drops abruptly by δ ¯ E ≈ − L (the lower branch).This is accompanied with some large avalanches, whichare named as spanning avalanches (SA) . This droppingshould be independent of ζ . To test these predictions,we have calculated and plotted ¯ E in terms of T for var-ious rates of ζ in Fig. 9(a). Two separate regimes are distinguishable in this figure: In the primitive times itincreases linearly, and for large enough times it enters anew regime, e.g., for ζ = 0 it is nearly constant. How-ever, for nonzero ζ we see that some oscillations arise inwhich the average grain number drops abruptly after alinear part in accordance with the prediction of the MFapproach. In the inset of this figure, we have plottedthe difference between these two limits (among which ¯ E oscillates) ¯ E − ¯ E in terms of ζ , which quantifies these os-cillations. This figure characterizes the bifurcation pointat which the transition to the oscillatory regime takesplace. Actually ¯ E − ¯ E starts from zero in small enough ζ and at some L -dependent bifurcation point ( ζ ∗ ) growsrapidly, and then saturates immediately.Fig. 9(b) visualizes this event, in which a new charac-teristic reference point appears for large enough ζ . Inthis figure, we have shown the mass of the avalanches (the number of distinct toppled sites in the avalanche) as afunction of time for various rates of ζ and L . Consider forexample L = 64 in the regime (cid:38) 8, for which the massesof some avalanches reach the system size, i.e., the toppoints in the figure whose mass is almost (64) = 4096.These avalanches are the mentioned SAs and are absentin small ζ . The SAs and the abrupt drop of averageheight occur simultaneously and therefore have the sameorigin (both belong to the lower branch of Eq.(39)). Theavalanches that belong to the first branch, whose meansizes grow linearly with the injections are called deformedavalanches (DAs). The mean size of DAs depends on E th ¯ E . It is notable that the microstates which grow withtime in the observed quasistationary state are transient.The existence of transient states in the quasi steady statemay lead to new studies, and new insights may come upfor the phenomena of SOC as a whole. H. propagation of electrons in 2D electron gas /f is a well-observed phenomenon in condensed mat-ter systems, especially in two-dimensional electron gas(2DEG). Whether its origin is in avalanche-like dynam-ics or not need a detailed analysis of the system. Herewe report on this possibility, by introducing an avalanchebase model for two-dimensional electron gas.In [190] the author presents a possibility based on whichthe percolation of the electrons with avalanche dynamicscan be a source for the MIT of the two-dimensional elec-tron gas in zero magnetic fields. They call it the semi-classical localization of electrons, which corresponds to percolative - non-percolative phase transition, although itis different from the conventional percolation theory inessence. Although we don’t have enough reasons to callit SOC, since the dynamics are much similar to BTW-likeavalanches, we brought it in this section. The percolative phase has the property dd T σ < σ ≡ the conductivity)which is the characteristics of the metallic phase. Inter-estingly this MIT occurs in the diffusion regime of 2DEGand therefore has nothing to do with the Anderson lo-calization. In this model, we consider a two-dimensionalelectron gas in contact with some electronic reservoirs.The dynamics of the electrons are divided in two cate-gories according to the phase relaxation time τ φ associ-ated with inelastic or spin-flip scattering up to which theelectrons retain their coherence. Corresponding to this,we divide the spatial dynamics of the electrons to twoscales: l (cid:28) r (cid:28) l φ and r (cid:29) l φ in which l is the meanfree path due to the electron-electron or the electron-phonon interactions, l φ ≡ (cid:112) Dτ φ is the phase relaxationlength, r is the length scale of the electron dynamics intime t which can be estimated classically as r ∼ √ Dt and D is the diffusion coefficient. In the first scale theelectrons retain their quantum phase, whereas, for thelatter case, the picture can be semi-classical, since quan-tum fluctuations in this scale do not play a vital roleand one can use the classical Boltzmann transport equa-tion [191]. This approach has been proved to be useful inmany situations and physical interpretation of some phe-nomena, like the interpretation of finite-size power-lawconductivity of 2DEG [192], the self-averaging [193], andthe percolation prescription of 2DEG [194] each of whichconsiders the linear size ∆ L ∼ l φ as an important spatialscale. We have treated the electron gas inside these cellspurely quantum mechanically, but for the transport ofthe particles to the neighboring cells some semi-classicalrules have been developed. To be most symmetric, thecells have been chosen to be hexagonal. This model isessentially different from conventional percolation theory(used for example in Ref. [194] and Ref. [195]). It is mostsuitable to be called a cellular automaton model in whichsome electrons can propagate throughout the system ac-cording to some local dynamical rules. In this model,the electrons propagate through the system according tothe (temperature-dependent) energy content and also thechemical potentials of the cells. In some cases, some elec-trons can reach from one side to the opposite boundary,which is called percolated . 1. general set up The energy of the electron gas and the chemical poten-tial inside each cell is calculated by means of the Thomas-Fermi-Dirac (TFD) approach. The average energy of the i th cell, inside which the charge is supposed to be uni-form, is (cid:104) E i (cid:105) = K ( T, ˜ N i ) + V ee ( T, ˜ N i ) + E imp ( T, ˜ N i ) inwhich the terms are finite temperature averages of the ki-netic, the electron-electron interaction and the impurityenergies respectively and ˜ N i is the number of electrons inthe cell. The total energy of a cell is shown to be [190]: E T = L (cid:88) i =1 (cid:104) − αT Li (cid:16) − e N i /T (cid:17) + β N i − γ i N i (cid:105) (40)in which α = 2 m (cid:16) πk B ζ φ ( T ) (cid:126) (cid:17) , β = √ (cid:15) ζ φ ( T ) (cid:16) eαk B (cid:17) , γ i = sinh − (1) αe π(cid:15) k B ζ φ ( T ) Z i , N i = k B α ˜ N i and L is the to-tal number of cells. Apparently this calculation containssome simplification which takes a part some complexi-ties that are unnecessary for the physics of the proposedMIT.The chemical potential of a cell ( µ i = ∂A i /∂N | V,T inwhich A i is the Helmholtz free energy of the i th cell)as the main building block of the transition rules ofelectrons between cells is obtained using the relation A ˜ N ( V, T ) − T (cid:0) ∂A∂T (cid:1) ˜ N,V = (cid:104) E (cid:105) . By considering the factthat µ ( T → → ζ φ ( T ) = aT − / for two dimen-sional electron gas [191] ( a is a proportionality constant),one finds that (see APPENDIX A Ref.[196]): µ i = k B T ln (cid:0) e h i − (cid:1) + U T h i − IZ i T (41)in which U = k B m √ Dae π (cid:15) (cid:126) , I = sinh − (1) e π(cid:15) √ Da , h i = N i T and i stands for the i th cell. The effect of randomnessof Z i ’s (that are supposed to be random noise with anuniform probability measure), which captures the on-site(diagonal) disorder is investigated. The probability ofadding a particle to the i th cell of the system is shownto be proportional to exp [ − β ˜ µ i ] in which ˜ µ i ≡ µ i − µ and µ is the average chemical potential of the system,whereas the probability of the transition between twosites (say cell 1 → cell 2) is obtained byrelative probability = e − β ( µ − µ ) (42)for which the following relation is used: µ − µ = k B T ln (cid:18) e h − e h − (cid:19) + U T ( h − h ) − IT ( Z − Z ) . (43)these relations are of especial importance in the follow-ing sections. This is the base of the model proposed.Consider figure 10a in which an electron system has0 (a) (b) FIG. 10: SEC. V H: (a) A schematic graph of dividing the 2D lattice into many hexagons. Inside the hexagons wehave a pure quantum electron gas. The transfer between the cells occurs semi-classically. The green points showelectrons and the red ones show the impurities. (b) A schematic set up of a 2D electron gas surrounded by chargereservoirs. The same partitioning has been carried out in this case. The electron can enter and exit the 2D systemat any random point (with some energy considerations), e.g. from the boundaries. (a) (b) FIG. 11: SEC. V H: (a) The electron transition to the neighboring cells. (b) An schematic movement pattern of anelectron in the virtual lattice. The gray cells are boundary sites (outer sites with µ = w φ ) from which the electronscan leave the system. The black circle is the site at which the electron has been injected and the gray circlesrepresent the sites which have been unstable and relaxed through a chain of charge transfersbeen divided into some hexagonal cells. We have showna cell and its neighbors in fig. 11a (see also Fig. 11b)each of which has its own local chemical potential µ i , i = 0 , , , ..., 6. Let us consider for a moment that theorder of potentials is µ i < µ j for i < j . According to Eq.(A.2). Ref. [196] the site 0 is said to be unstable (hasthe potential to give an electron to its neighbors) if itschemical potential µ (0) exceeds the chemical potentialof the bulk µ . If the mentioned site is unstable, it hasthe potential to release electrons to the neighbors, andthe first candidate for this charge transfer is the neighborwith the smallest chemical potential µ , i.e. µ here. Af-ter this process (whether the charge transfer to the firstneighbor has taken place or not) the next candidate for the electron transfer is the site with the nearest µ to µ ,i.e. µ , etc. The Metropolis Monte Carlo method [197]is employed for these charge transfers, i.e. the electrontransport from 0 to any site i is occurred with the prob-ability: P → i ∼ (cid:40) Θ( µ (0) − µ ) × Max (cid:8) , e − β ( µ ( i ) − µ (0)) (cid:9) i areneighbors, and the second line is for the other cases, andΘ( x ) is the step function.Using this probability, Najafi showed that the system1undergoes a percolation transition inline in the T − ∆phase space, realizing the observed metal-insulator tran-sitions in 2D electron gases (2DEG) [198]. A separatestudy on the 1 /f noise in 2DEG is currently being donebased on the same physics explained above. VI. CONLUSION In this paper, we reviewed the SOC concepts in varioussystems. First, we presented some examples, includingthe systems that show SOC, like earthquake, rain falling etc . In the second part, we presented the evidence show-ing that the BTW sandpile model is c = − W -algebras. The simulation results for SOC invarious systems were presented in the last part. There weconsidered the SOC in fluid propagation in porous me-dia, in cumulus clouds, in an excitable random system,in imperfect supports, and in 2DEG. We also consideredvibrating ASM, invasion sandpile model, and diffusivesandpiles. [1] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. A ,364 (1988).[2] D. Dhar, Physical Review Letters , 1613 (1990).[3] S. N. Majumdar and D. Dhar, Journal of Physics A:Mathematical and General , L357 (1991).[4] D. Dhar, Physica A: Statistical Mechanics and its Ap-plications , 29 (2006).[5] E. V. Ivashkevich, D. V. Ktitarev, and V. B. Priezzhev,Physica A: Statistical Mechanics and its Applications , 347 (1994).[6] D. Dhar and S. S. Manna, Phys. Rev. E , 2684 (1994).[7] D. V. Ktitarev and V. B. Priezzhev, Phys. Rev. E ,2883 (1998).[8] S. Majumdar and D. Dhar, Physica A: Statistical Me-chanics and its Applications , 129 (1992).[9] S. Mahieu and P. Ruelle, Phys. Rev. E , 066130(2001).[10] H. Saleur and B. Duplantier, Physical review letters ,2325 (1987).[11] A. Coniglio, Phys. Rev. Lett. , 3054 (1989).[12] S. N. Majumdar, Phys. Rev. Lett. , 2329 (1992).[13] O. Schramm, Israel Journal of Mathematics , 221(2000).[14] M. N. Najafi, S. Moghimi-Araghi, and S. Rouhani,Phys. Rev. E , 051104 (2012).[15] D. Markovi´c and C. Gros, Physics Reports , 41(2014).[16] O. Peters, C. Hertlein, and K. Christensen, Physicalreview letters , 018701 (2001).[17] P. Charbonneau, S. W. McIntosh, H.-L. Liu, and T. J.Bogdan, Solar Physics , 321 (2001).[18] L. P. Karakatsanis, G. Pavlos, and D. Sfiris, Interna-tional Journal of Bifurcation and Chaos , 1250209(2012).[19] M. A. Munoz, R. Dickman, R. Pastor-Satorras,A. Vespignani, and S. Zapperi, in AIP Conference Pro-ceedings , Vol. 574 (American Institute of Physics, 2001)pp. 102–110.[20] R. Dickman, M. Alava, M. A. Munoz, J. Peltola,A. Vespignani, and S. Zapperi, Physical Review E ,056104 (2001).[21] A. Sornette and D. Sornette, EPL (Europhysics Letters) , 197 (1989). [22] L. Telesca, V. Cuomo, V. Lapenna, and M. Macchiato,Geophysical research letters , 3765 (2001).[23] B. Gutenberg and C. F. Richter, Science , 183 (1936).[24] A. Davis, A. Marshak, W. Wiscombe, and R. Cahalan,Journal of Geophysical Research: Atmospheres , 8055(1994).[25] J. M. Carlson, J. S. Langer, and B. E. Shaw, Reviewsof Modern Physics , 657 (1994).[26] Z. Olami, H. J. S. Feder, and K. Christensen, PhysicalReview Letters , 1244 (1992).[27] P. Bak and C. Tang, Journal of Geophysical Research:Solid Earth , 15635 (1989).[28] Y. Huang, H. Saleur, C. Sammis, and D. Sornette, EPL(Europhysics Letters) , 43 (1998).[29] D. L. Turcotte and B. D. Malamud, Physica A: Statis-tical Mechanics and its Applications , 580 (2004).[30] U. Lohmann, F. L¨u¨ond, and F. Mahrt, An introductionto clouds: From the microscale to climate (CambridgeUniversity Press, 2016).[31] S. Lovejoy and D. Schertzer, Journal of Geophysical Re-search: Atmospheres , 2021 (1990).[32] H. Hentschel and I. Procaccia, Physical Review A ,1461 (1984).[33] R. F. Cahalan and J. H. Joseph, Monthly weather re-view , 261 (1989).[34] J. H. Joseph and R. F. Cahalan, Journal of AppliedMeteorology , 793 (1990).[35] B. Gutenberg and C. F. Richter, Bulletin of the Seismo-logical society of America , 163 (1942).[36] B. D. Malamud, D. L. Turcotte, F. Guzzetti, and P. Re-ichenbach, Earth Surface Processes and Landforms ,687 (2004).[37] D. Sornette, Journal de Physique I , 2089 (1992).[38] J. Nussbaum and A. Ruina, pure and applied geophysics , 629 (1987).[39] M. Najafi, M. Rahimi-Majd, and T. Shirzad, EPL (Eu-rophysics Letters) , 20001 (2020).[40] Y. Y. Kagan and L. Knopoff, Science , 1563 (1987).[41] R. Burridge and L. Knopoff, Bulletin of the seismologi-cal society of america , 341 (1967).[42] Z. Olami and K. Christensen, Physical Review A ,R1720 (1992).[43] F. Omori, On the after-shocks of earthquakes , Vol. 7(The University, 1894). [44] S. Hergarten and H. J. Neugebauer, Physical review let-ters , 238501 (2002).[45] A. Helmstetter, S. Hergarten, and D. Sornette, PhysicalReview E , 046120 (2004).[46] S. Clar, B. Drossel, and F. Schwabl, Journal of Physics:Condensed Matter , 6803 (1996).[47] B. D. Malamud, G. Morein, and D. L. Turcotte, Science , 1840 (1998).[48] B. Drossel and F. Schwabl, Physical review letters ,1629 (1992).[49] B. Drossel and F. Schwabl, Physica A: Statistical Me-chanics and Its Applications , 183 (1993).[50] P. Grassberger, New Journal of Physics , 17 (2002).[51] P. Grassberger, Journal of Physics A: Mathematical andGeneral , 2081 (1993).[52] R. Krenn, Natural hazards and self-organized criticality (na, 2012).[53] A. Ruzmaikin, in Symposium-International Astronomi-cal Union , Vol. 138 (Cambridge University Press, 1990)pp. 343–353.[54] L. Karakatsanis and G. Pavlos, Nonlinear Phenomenain Complex Systems , 280 (2008).[55] P. Bak, “How nature works: the science of self-organizedcriticality. 1996, new york: Copernicus,”.[56] L. Vlahos, SP-505 , 186 (2002).[57] P. Bak, C. Tang, and K. Wiesenfeld, Physical reviewletters , 381 (1987).[58] H. E. Hurst, Nature , 494 (1957).[59] R. F. S. Andrade, H. Schellnhuber, and M. Claussen,Physica A: Statistical Mechanics and its Applications , 557 (1998).[60] Z. Wang and C. Huang, Advances in Meteorology (2012).[61] A. Sarkar and P. Barat, Fractals , 289 (2006).[62] C. Nnaji, Journal of Science and Technology (Ghana) (2011).[63] R. Bove, V. Pelino, and L. De Leonibus, Communi-cations in Nonlinear Science and Numerical Simulation , 678 (2006).[64] A. Garc´ıa-Mar´ın, F. Jim´enez-Hornero, and J. Ayuso,Hydrological Processes: An International Journal ,295 (2008).[65] A. Deluca, N. R. Moloney, and ´A. Corral, Physicalreview E , 052808 (2015).[66] S. T. Pinho and R. F. Andrade, Physica A: StatisticalMechanics and its Applications , 483 (1998).[67] R. F. Andrade, S. T. Pinho, S. C. Fraga, and A. P.Tanajura, Physica A: Statistical Mechanics and its Ap-plications , 405 (2002).[68] R. F. S. Andrade, Brazilian journal of physics , 437(2003).[69] S. Lovejoy, Science , 185 (1982).[70] P. Austin, M. Baker, A. Blyth, and J. Jensen, Journalof the atmospheric sciences , 1123 (1985).[71] R. Chatterjee, K. Ali, and P. Prakash, (1994).[72] C. von Savigny, L. A. Brinkhoff, S. M. Bailey, C. E.Randall, and J. M. Russell III, Geophysical ResearchLetters (2011).[73] S. P. Malinowski and I. Zawadzki, Journal of the atmo-spheric sciences , 5 (1993).[74] A. Batista-Tom´as, O. D´ıaz, A. Batista-Leyva, andE. Altshuler, Quarterly Journal of the Royal Meteoro-logical Society , 983 (2016). [75] F. S. Rys and A. Waldvogel, Physical review letters ,784 (1986).[76] J. H. Joseph and R. F. Cahalan, Journal of AppliedMeteorology , 793 (1990).[77] J. Olsson, J. Niemczynowicz, and R. Berndtsson, Jour-nal of Geophysical Research: Atmospheres , 23265(1993).[78] S. P. Malinowski, M. Y. Leclerc, and D. G. Baumgard-ner, Journal of the atmospheric sciences , 397 (1994).[79] T. C. Benner and J. A. Curry, Journal of GeophysicalResearch: Atmospheres , 28753 (1998).[80] S. M. Rodts, P. G. Duynkerke, and H. J. Jonker, Jour-nal of the atmospheric sciences , 1895 (2003).[81] J.-I. Yano and Y. Takeuchi, Journal of the Meteorolog-ical Society of Japan. Ser. II , 661 (1987).[82] K. Gotoh and Y. Fujii, Journal of applied meteorology , 1283 (1998).[83] S. Lovejoy and D. Schertzer, in Non-Linear Variabilityin Geophysics (Springer, 1991) pp. 111–144.[84] S. Lovejoy, D. Schertzer, and A. Tsonis, Science ,1036 (1987).[85] R. F. Cahalan and J. H. Joseph, Monthly weather re-view , 261 (1989).[86] P. Gabriel, S. Lovejoy, D. Schertzer, and G. Austin,Geophysical research letters , 1373 (1988).[87] S. Lovejoy and D. Schertzer, Journal of Geophysical Re-search: Atmospheres , 2021 (1990).[88] Y. Tessier, S. Lovejoy, and D. Schertzer, Journal ofApplied Meteorology , 223 (1993).[89] J. D. Pelletier, Physical review letters , 2672 (1997).[90] S. Sengupta, R. Welch, M. Navar, T. Berendes, andD. Chen, Journal of Applied Meteorology , 1245(1990).[91] O. Peters and J. D. Neelin, Nature physics , 393 (2006).[92] J.-I. Yano, C. Liu, and M. W. Moncrieff, Journal of theatmospheric sciences , 3449 (2012).[93] M. R., Y. S., D. Y., and N. T., Proceedings NinthPacific Conference on Computer Graphics and Applica-tions. Pacific Graphics, IEEE , 363 (2001).[94] O. Ramos, E. Altshuler, and K. M˚aløy, Physical reviewletters , 078701 (2009).[95] L. A. N. Amaral and K. B. Lauritsen, Physical reviewE , R4512 (1996).[96] C. Aegerter, R. G¨unther, and R. Wijngaarden, PhysicalReview E , 051306 (2003).[97] V. Frette, K. Christensen, A. Malthe-Sørenssen,J. Feder, T. Jøssang, and P. Meakin, Nature , 49(1996).[98] K. Christensen, ´A. Corral, V. Frette, J. Feder, andT. Jøssang, Physical review letters , 107 (1996).[99] L. A. N. Amaral and K. B. Lauritsen, Physica A: Sta-tistical Mechanics and its Applications , 608 (1996).[100] D. Denisov, Y. Villanueva, K. L˝orincz, S. May, andR. Wijngaarden, Physical Review E , 051309 (2012).[101] T. J. Coulthard and M. J. Van De Wiel, Geomorphology , 216 (2007).[102] M. J. Van De Wiel and T. J. Coulthard, Geology ,87 (2010).[103] K. Shi and C.-Q. Liu, Atmospheric Environment ,3301 (2009).[104] Z. Liu, J. Xu, and K. Shi, Theoretical and appliedclimatology , 685 (2014). [105] L. de Arcangelis, C. Perrone-Capano, and H. J. Her-rmann, Physical review letters , 028107 (2006).[106] M. Bartolozzi, D. Leinweber, and A. Thomas, PhysicaA: Statistical Mechanics and its Applications , 451(2005).[107] D. Stauffer and D. Sornette, Physica A: Statistical Me-chanics and its Applications , 496 (1999).[108] J. Valdivia, J. Rogan, V. Munoz, L. Gomberoff, A. Kli-mas, D. Vassiliadis, V. Uritsky, S. Sharma, B. Toledo,and L. Wastavino, Advances in Space Research , 961(2005).[109] J. Wanliss and V. Uritsky, Journal of Geophysical Re-search: Space Physics (2010).[110] V. M. Uritsky, A. Pouquet, D. Rosenberg, P. D.Mininni, and E. Donovan, Physical Review E ,056326 (2010).[111] G. Szab´o, M. Alava, and J. Kert´esz, EPL (EurophysicsLetters) , 665 (2002).[112] C. Aegerter, M. Welling, and R. Wijngaarden, EPL(Europhysics Letters) , 753 (2004).[113] M. Bogu˜n´a and ´A. Corral, Physical review letters ,4950 (1997).[114] M. Paczuski, S. Maslov, and P. Bak, Physical ReviewE , 414 (1996).[115] M. S.S. and S. A.L., Physica A: Statistical Mechanicsand its Applications , 135 (2002).[116] I. Langmuir, Journal of meteorology , 175 (1948).[117] S. L¨ubeck and K. D. Usadel, Phys. Rev. E , 5138(1997).[118] M. Najafi, S. Moghimi-Araghi, and S. Rouhani, Journalof Physics A: Mathematical and Theoretical , 095001(2012).[119] M. Najafi, S. Moghimi-Araghi, and S. Rouhani, Physi-cal Review E , 051104 (2012).[120] M. Najafi, M. Ghaedi, and S. Moghimi-Araghi, PhysicaA: Statistical Mechanics and its Applications , 102(2016).[121] H. Stommel, Journal of Meteorology , 91 (1947).[122] S. Manna, Journal of Physics A: Mathematical and Gen-eral , L363 (1991).[123] A. H., M.-A. S., and N. M.N., Physica A: StatisticalMechanics and its Applications , 196 (2015).[124] L. P. Kadanoff, S. R. Nagel, L. Wu, and S.-m. Zhou,Physical Review A , 6524 (1989).[125] S. Manna, Physica A: Statistical Mechanics and its Ap-plications , 249 (1991).[126] Y.-C. Zhang, Physical Review Letters , 470 (1989).[127] D. Dhar and R. Ramaswamy, Physical Review Letters , 1659 (1989).[128] A. Ben-Hur and O. Biham, Physical Review E ,R1317 (1996).[129] K. Christensen and Z. Olami, Physical Review E ,3361 (1993).[130] A. Chessa, H. E. Stanley, A. Vespignani, and S. Zap-peri, Physical Review E , R12 (1999).[131] R. Dickman and J. Campelo, Physical Review E ,066111 (2003).[132] Z. S.D., Physical Review E , 259 (1999).[133] L. S., Physical Review E , 1590 (1997).[134] S. N. Majumdar and D. Dhar, Physica A: StatisticalMechanics and its Applications , 129 (1992).[135] W. Janke and A. M. Schakel, Nuclear Physics B ,385 (2004). [136] V. Gurarie, Nuclear Physics B , 535 (1993).[137] S. Moghimi-Araghi, S. Rouhani, and M. Saadat, Nu-clear Physics B , 531 (2001).[138] S. Moghimi-Araghi, S. Rouhani, and M. Saadat, Inter-national Journal of Modern Physics A , 4747 (2003).[139] P. Francesco, P. Mathieu, and D. S´en´echal, Conformalfield theory (Springer, 1996).[140] K. Kyt¨ol¨a and D. Ridout, Journal of mathematicalphysics , 123503 (2009).[141] M. R. Gaberdiel, International Journal of ModernPhysics A , 4593 (2003).[142] F. M. K. A. N. W. R. A. Blumenhagen, R. and R. Varn-hagen, Nuclear Physics B , 255 (1991).[143] A. B. Zamolodchikov, in W-Symmetry (World Scientific,1995) pp. 221–229.[144] M. R. Gaberdiel and A. Neitzke, Communications inmathematical physics , 305 (2003).[145] M. Rajabpour, S. Rouhani, and A. Saberi, Fortschritteder Physik: Progress of Physics , 1289 (2007).[146] J. Cardy, Annals of Physics , 81 (2005).[147] J. Cheraghalizadeh, M. Najafi, H. Dashti-Naserabadi,and H. Mohammadzadeh, Physical Review E , 052127(2017).[148] O. Schramm, in Selected Works of Oded Schramm (Springer, 2011) pp. 791–858.[149] M. N. Najafi, J. Cheraghalizadeh, M. Lukovi´c, and H. J.Herrmann, Phys. Rev. E , 032116 (2020).[150] S. Oswald, W. Kinzelbach, A. Greiner, and G. Brix,Geoderma , 417 (1997).[151] J. Warner, Journal of the Atmospheric Sciences , 1049(1969).[152] S. Lovejoy, Science , 185 (1982).[153] R. Chatterjee, K. Ali, and P. Prakash, (1994).[154] K. Madhushani and D. Sonnadara, (2012).[155] H. Jaeger, C.-h. Liu, and S. R. Nagel, Physical ReviewLetters , 40 (1989).[156] A. Mehta and G. Barker, Physical review letters , 394(1991).[157] D. Wilkinson and J. F. Willemsen, Journal of PhysicsA: Mathematical and General , 3365 (1983).[158] R. Glass and L. Yarrington, Geoderma , 231 (1996).[159] A. P. Sheppard, M. A. Knackstedt, W. V. Pinczewski,and M. Sahimi, Journal of Physics A: Mathematical andGeneral , L521 (1999).[160] M. Najafi, Physics Letters A , 2008 (2014).[161] S. L¨ubeck and K. D. Usadel, Phys. Rev. E , 4095(1997).[162] S. L¨ubeck and K. Usadel, Physical Review E , 5138(1997).[163] S. A. Glantz, B. K. Slinker, and T. B. Neilands, Primerof applied regression and analysis of variance , Vol. 309(McGraw-Hill New York, 1990).[164] M. Najafi, Z. Moghaddam, M. Samadpour, and N. A.Ara´ujo, arXiv preprint arXiv:2003.02482 (2020).[165] M. Najafi and H. Dashti-Naserabadi, Physical ReviewE , 032108 (2018).[166] J. M. Beggs and D. Plenz, Journal of neuroscience ,11167 (2003).[167] K.-M. Lee, K.-I. Goh, and I.-M. Kim, Journal of theKorean Physical Society , 641 (2012).[168] D.-S. Lee, K.-I. Goh, B. Kahng, and D. Kim, PhysicaA: Statistical Mechanics and its Applications , 84(2004). [169] R. Karmakar and S. Manna, Journal of Physics A:Mathematical and General , L87 (2005).[170] J. Lahtinen, J. Kert´esz, and K. Kaski, Physica A: Sta-tistical Mechanics and its Applications , 535 (2005).[171] G.-J. Pan, D.-M. Zhang, Y.-P. Yin, and M.-H. He,Physica A: Statistical Mechanics and its Applications , 435 (2007).[172] H. Bhaumik and S. Santra, arXiv preprintarXiv:1705.10646 (2017).[173] H. Dashti-Naserabadi and M. Najafi, Physical ReviewE , 052145 (2015).[174] M. N. Najafi, S. Moghimi-Araghi, and S. Rouhani,Phys. Rev. E , 051104 (2012).[175] H. Bhaumik and S. Santra, Physical Review E ,062817 (2013).[176] A. R. Kose, B. Fischer, L. Mao, and H. Koser, Proceed-ings of the National Academy of Sciences , 21478(2009).[177] H. Kikura, J. Matsushita, M. Matsuzaki, Y. Kobayashi,and M. Aritomi, Science and Technology of AdvancedMaterials , 703 (2004).[178] M. Matsuzaki, H. Kikura, J. Matsushita, M. Aritomi,and H. Akatsuka, Science and Technology of AdvancedMaterials , 667 (2004).[179] H. Kikura, J. Matsushita, N. Kakuta, M. Aritomi, andY. Kobayashi, Journal of materials processing technol-ogy , 93 (2007).[180] E. Daryaei, N. Ara´ujo, K. Schrenk, S. Rouhani, andH. Herrmann, Physical review letters , 218701(2012).[181] M. Najafi and M. Ghaedi, Physica A: Statistical Me-chanics and its Applications , 82 (2015).[182] C. Oliveira, A. Ara´ujo, L. Lucena, M. Almeida, andJ. Andrade, Physica A: Statistical Mechanics and its Applications , 3219 (2012).[183] W. Li, J. L. Jensen, W. B. Ayers, S. M. Hubbard, andM. R. Heidari, Journal of Petroleum Science and Engi-neering , 180 (2009).[184] M. Najafi, Journal of Physics A: Mathematical and The-oretical , 335003 (2016).[185] M. Najafi and H. Dashti-Naserabadi, Journal of Statis-tical Mechanics: Theory and Experiment , 023211(2018).[186] M. Najafi, arXiv preprint arXiv:1801.08978 (2018).[187] P. Bak, C. Tang, and K. Wiesenfeld, Phys. Rev. Lett. , 381 (1987).[188] H. Dashti-Naserabadi and M. Najafi, Physical ReviewE , 042115 (2017).[189] M. N. Najafi and Z. Moghadam, Phys. Rev. E ,042120 (2019).[190] M. Najafi, Solid State Communications , 84 (2018).[191] B. L. Altshuler and A. G. Aronov, in Modern Problemsin condensed matter sciences , Vol. 10 (Elsevier, 1985)pp. 1–153.[192] D. Backes, R. Hall, M. Pepper, H. Beere, D. Ritchie,and V. Narayan, Physical Review B , 235427 (2015).[193] H. Bruus and K. Flensberg, Many-body quantum theoryin condensed matter physics: an introduction (Oxforduniversity press, 2004).[194] Y. Meir, Physical review letters , 3506 (1999).[195] S. D. Sarma, M. Lilly, E. Hwang, L. Pfeiffer, K. West,and J. Reno, Physical review letters , 136401 (2005).[196] M. N. Najafi, The European Physical Journal B , 172(2019).[197] H. Gould and J. Tobochnik, Computer SimulationMethods (Addison-Wesley Reading, 1996).[198] S. V. Kravchenko, G. Kravchenko, J. Furneaux, V. M.Pudalov, and M. dIorio, Physical Review B50