Dynamic polarization random walk model and fishbone-like instability for self-organized critical systems

We study the phenomenon of self-organized criticality (SOC) as a transport problem for electrically charged particles. A model for SOC based on the idea of a dynamic polarization response with random walks of the charge carriers gives critical exponents consistent with the results of numerical simulations of the traditional ‘sandpile’ SOC models, and stability properties, associated with the scaling of the control parameter versus distance to criticality. Relaxations of a supercritical system to SOC are stretched-exponential similar to the typically observed properties of non-Debye relaxation in disordered amorphous dielectrics. Overdriving the system near self-organized criticality is shown to have a destabilizing effect on the SOC state. This instability of the critical state constitutes a fascinating nonlinear system in which SOC and nonlocal properties can appear on an equal footing. The instability cycle is qualitatively similar to the internal kink (‘fishbone’) mode in a magnetically confined toroidal plasma where beams of energetic particles are injected at high power, and has serious implications for the functioning of complex systems. Theoretical analyses, presented here, are the basis for addressing the various patterns of self-organized critical behavior in connection with the strength of the driving. The results of this work also suggest a type of mixed behavior in which the typical multi-scale features due to SOC can coexist along with the global or coherent features as a consequence of the instability present. An example of this coexistence is speculated for the solar wind–magnetosphere interaction.


Introduction
In 1982, the challenge to understand 1/ f 'noise' inspired Montroll and Shlesinger [1] to suggest that distributions with long tails should be a consequence of a simple generic stochastic process in much the same way as Gaussian, or normal, distribution in the theory of probability is a consequence of the central limit theorem. It took five more years before Bak, Tang and Wiesenfeld (BTW) [2] introduced, in 1987, the paradigmatic concept of self-organized criticality, or SOC, operating on systems natural attraction to fractal structures. The claim was that irreversible dynamics of systems with many coupled degrees of freedom ('complex' systems) would intrinsically generate self-organization into a critical state without fine tuning of external parameters. This conjecture initiated a burst of research activity in experiment, theory and simulation. Diverse formulations of SOC have been proposed, emphasizing the strengths and weaknesses of particular approaches [3]. The phenomenon was demonstrated on a number of automated lattice models, or 'sandpiles', displaying avalanche dynamics and scale invariance [4]- [6]. An impressive list of publications 2 has been produced in the attempt to prove or disprove the SOC hypothesis for the various systems. The notion of SOC can be thought of as belonging to the nascent 'science of complexity', which addresses the commonalities between apparently dissimilar natural, technological and socio-economic phenomena, e.g. market crashes and climate disruptions [7]. Despite its promising performance, 3 the SOC hypothesis is a subject of strong debate in the literature, and many issues related to it remain controversial or demand further investigation.
In this paper, we analyze the SOC problem as a transport problem for electric charges. We provide analytical methods to calculate the various critical exponents by employing the formalism of dielectric-relaxation theory. Our theoretical findings are in good agreement with computer simulations of the traditional sandpiles [4,5]. The lattice dynamics are described in terms of random walk hopping of charged particles on a self-adjusting percolation system. We show that the control parameter must vanish faster than a certain scaling law to give rise to a stable SOC state. Overdriving the system near self-organized criticality leads to the excitation of unstable modes characterized by sharp periodic outbreaks in the particle loss current. This instability of the critical state constitutes a fascinating nonlinear system in which SOC and nonlocal properties can appear on an equal footing. The instability cycle is qualitatively similar to the internal kink ('fishbone') mode in tokamaks with high-power beam injection [8]. Relaxations to the critical state are found to be stretched-exponential similar to the typically observed properties of dielectric relaxation in disordered solids and glasses [9,10]. The results presented in this work suggest a type of mixed behavior in which the typical multiscale features due to SOC can coexist along with the global or coherent features. One example of this coexistence is speculated for the coupled solar wind-magnetosphere-ionosphere system. A short account of some of these investigations has been reported recently [11].
The paper is organized as follows. The model, dubbed the dynamic polarization random walk (DPRW) model, is introduced in section 2. Following this is a presentation of analytical theory of SOC dynamics (section 3). A discussion of the various aspects of the DPRW approach to SOC is given in section 4. Instabilities of the SOC state are considered in section 5, where a condition on the strength of the driving is also obtained. We summarize our findings in section 6.

Description of the model
The model is motivated by the well-known problem [12]- [14] that by increasing the probability of site occupancy a lattice is brought to a critical point characterized by fractal geometry of the threshold percolation, i.e. self-similar distribution of arbitrarily large finite clusters, each presenting the same fractal geometry as the infinite cluster. We consider a hypercubic d-dimensional (d 1) lattice confined between two opposite (d − 1)-dimensional hyperplanes, which form a parallel-plate 'capacitor' as shown in figure 1. The plate on the right-hand side is earthed. Free charges are built by external forces on the capacitor's left plate. When a unit free charge is added to the capacitor, the lattice responds by burning a unit 'polarization' charge, which is an occupied site added at random to the lattice. When a unit free charge is removed from the capacitor, a randomly chosen occupied site is converted into a 'hole' site (missing occupied site). A hole will be deleted from the system (converted into an empty site) if/when the corresponding free charge has reached (or been moved to) the ground.
There is a limit, Q max , on the amount of free charges the capacitor can store, and this is defined as Q max = ep c N , where e is the elementary charge (e = −1), p c is the percolation threshold and N is the total number of sites across the lattice. Thus, the ability of the capacitor to store electric charges is limited to the occurrence of the infinite cluster at the percolation point. If, at any time, the above limit is exceeded, a double amount of free charges in excess of Q max will be removed from the capacitor, and will be distributed between the sites of the infinite cluster with equal probability. The implication is that the capacitor leaks electric charges above the percolation point. This property reflects the onset of the dc conduction at the threshold of percolation. The doubling of the amount of electric charges to be removed from the capacitor mimics nonzero inductance of the conduction process.
When a hole appears on the infinite cluster, it causes an activation event with the following consequence: one of the nearest-neighbor occupied sites, which is a random choice, will deliver its charge content to the hole. The hole that has just received the polarization charge becomes an ordinary occupied site, while the donor site becomes a hole. The newborn hole, in turn, will cause a further activation event at the location where it has occurred, thus triggering a chain reaction of redistribution of polarization charges. The chain reaction continues until the hole reaches the earthed plate where it is absorbed (converted into an empty site). When a hole appears on a finite cluster, it causes a chain reaction of activation events in much the same way as on the infinite cluster but with one modification regarding the ending of the activation. The chain reaction stops if (i) similar to the infinite cluster case the hole reaches the earthed plate where it is converted into an empty site, or if (ii) there is no more activity going on in the infinite cluster. In the latter case, the finite cluster freezes in a 'glassy' state with the holes still in it until either a new hole appears on the infinite cluster or one or more occupied sites are added to the lattice by external forces.

Random walk hopping process
Essentially, the holes interchange their position with the nearest-neighbor occupied sites, and it appears reasonable to model this process as an interchange hopping process [15]. In what follows we assume, following [16], that there is a characteristic microscopic hopping time, 5 which is taken to be unity, but more general hopping models can be obtained by introducing a distribution of waiting times between consecutive steps of the hopping motion (i.e. continuous time random walks) [17,18]. With the above assumption that the site acting as donor is a random choice, the transport model is defined as a random walk hopping model. Similarly to the hole case, the free charges are assumed to behave as random walkers after their re-injection into the infinite cluster. They will hop at a constant rate between the nearest-neighbor occupied sites in random directions on the cluster on which they are initially placed until they reach the earthed plate where they sink in the ground circuit (see figure 1). The holes act as conducting sites for the motion of free charges. The charged plate acts as a perfectly reflecting boundary. Hops to empty sites are forbidden. The latter condition limits random walks to fractal geometry of the threshold percolation.

Dynamical geometry of threshold percolation
Overall, one can see that the system responds by chain reactions of random walk hopping processes when it becomes slightly supercritical and is quiescent otherwise. Excess free charges dissipating at the earthed plate provide a feedback mechanism by which the system returns to the percolation point. There will be a slowly (as compared to hopping motions) evolving dynamical geometry of the threshold percolation resulting from competition between the addition of occupied sites to the lattice and the charge-releasing chain reactions. Based on the quantitative analysis below, we identify this state as a SOC state. This general picture based on the idea of a dynamic polarization response with random walk hopping of the charge carriers might be called the dynamic polarization random walk (DPRW) model.
In the DPRW SOC model, chain reactions of the hopping motion acquire the role of 'avalanches' in traditional sandpiles. In this work, we are interested in obtaining the critical exponents of the DPRW model by means of analytical theory. Numerical simulation of the DPRW dynamics is under way for comparison with analytical predictions. At the time of writing, the characteristic signatures of the multi-scale conductivity response of the dynamical system at criticality are being confirmed in the computer simulation model. In figure 2, we illustrate the existence of relaxation events of various sizes due to hole hopping on a 10 × 10 square lattice with random distribution of the conducting nodes and the probability of site occupancy so as to mimic the percolation threshold and conjectured SOC activities.

Dynamics and orderings
Starting from an empty lattice (no potential difference between the plates), by randomly adding occupied sites to it, one builds the fractal geometry of random, or uncorrelated, percolation, characterized by three percolation critical exponents β, ν and µ (connected clusters have fractal dimensionality d f = d − β/ν) [12]- [14]. Note that the infinite percolation cluster, in the true sense of the wording, exists only in the thermodynamic limit when the lattice itself is infinite. This limit arises because of the need to model the system-sized conducting clusters in terms of fractal geometry. In the absence of holes this percolation geometry is static (polarization charges can only move by exchanging their position with a hole), but when the holes appear on the lattice they cause local rearrangements in the distribution of conducting sites. As a consequence, the conducting clusters on which the transport processes concentrate change their shape and their position in real space. In the analysis of this section, we require that the average number density of holes is very small compared to the average number density of polarization charges. The implication is that the system remains near the percolation point despite the slow evolution of the conducting clusters. Note that the lattice rules are such as to preserve the properties of random percolation. In fact, no correlations are introduced in the distribution of the conducting sites at any step of the lattice update.

Frequency-dependent conductivity and power spectral density
Given an input electric driving field E(t, r), the polarization response of the system is defined by where the response function χ (t − t ) is identically zero for t < t as required by causality. We should stress that nonlocal integration over the space variable is not needed here in view of the local (nearest-neighbor) character of the lattice interactions. In a model in which the assumption of locality is relaxed as for instance in models permitting particle Lévy flights, integration over the space variable is expected to produce a physically nontrivial effect. We do not consider such models here. A Fourier transformed χ(t) defines the frequency-dependent complex susceptibility of the system, χ(ω). In a basic theory of polarization response, one also introduces the frequency-dependent complex ac conductivity, σ ac (ω), which is related to χ(ω) by the Kramers-Kronig integral, see equation (6). The dependence of the ac conductivity on frequency, specialized to the random walks on percolation systems [19,20], is given by the 7 scaling relation σ ac (ω) ∝ ω η , where the power exponent η (0 η 1) is expressible in terms of the percolation indices β, ν and µ as η = µ/(2ν + µ − β). A derivation of this result is due to Gefen et al [21]. Note that the scaling σ ac (ω) ∝ ω η incorporates conductivity responses from all clusters at percolation, including those that are finite. In the DPRW model, these implications are matched by the mechanism of the hole conduction permitting the polarization current on both infinite and finite clusters. The general linear-response theory expression [22] for conductivity σ ac (ω) in terms of the mean-square displacement from the origin r 2 (t) is with n d a constant depending on the dimensionality of the lattice and n and e the density and charge of the carriers, respectively. The function D(ω) has the sense of the frequencydependent diffusion coefficient [23]. It is required that the frequency ω be large compared to the characteristic evolution frequency in the distribution of the conducting sites. In what follows, we denote this frequency by ω * . We have ω ω * . Consistently with the above definitions, the inverse frequency 1/ω * is ordered as the characteristic diffusion time on the infinite cluster, T d ξ 2 /D(ω * ). Note that this time will depend on ω * in accordance with equation (2). Hence, the pair connectedness) length, p (0 p 1) is the probability of site occupancy and p c is the percolation threshold. Noting that D(ω * ) ∝ (ω * ) η , one obtains ω * ∝ | p − p c | 2ν/(1−η) . Observe that ω * → 0 for p → p c . Recalling that there is a microscopic hopping time, which is taken to be unity, we assess the Kubo number [24,25] in the vicinity of the SOC state as One sees that Q * → ∞ for p → p c . The Kubo number is a suitable dimensionless parameter that quantifies how the evolution processes in the lattice compare with the microscopic hopping motions. The divergency of the Kubo number at criticality implies that there is a time scale separation: fast hopping motions versus infinitesimal evolution change. For Q * → ∞, we also find that We identify this scaling law with low-frequency 'percolation' scaling [13,26] for particle diffusion on a time evolving fractal structure. By applying the Kramers-Kronig relations Im χ(ω) ∝ σ ac (ω)/ω and it is found that χ (ω) ∝ ω −γ , with γ = 1 − η. A Fourier transformed equation (1) reads as P(ω, r) = χ (ω)E(ω, r). One can see that the power spectral density, S(ω), of the system response to a white-noise perturbation, E(ω, r) = 1, will be proportional to |χ(ω) . Hence, the power spectral density in the DPRW model is given by an inverse power-law distribution, with the α value depending on scaling properties of the ac conductivity response.

Exponent of stretched-exponential relaxation and distribution of relaxation times
Next we obtain the distribution of relaxation times self-consistently. For this, assume that the system is slightly supercritical, then consider a charge density perturbation, δρ(t, r), caused by the presence of either free charges or holes on the conducting clusters. 'Slightly supercritical' means that the dependence of the ac conductivity response on frequency can, with good accuracy, be taken in the power-law form σ ac (ω) ∝ ω 1−γ discussed above. The implication is that on adding δρ(t, r) to the conducting system at percolation we neglect the departure of the system's geometric properties from pure self-similarity. Without loss in generality, we assume that the perturbation δρ(t, r) is created instantaneously at time t = 0. This means that the function δρ(t, r) ≡ 0 for t < 0 for all r. The perturbation δρ(t, r) generates an electric field inhomogeneity δE(t, r) in accordance with Maxwell's equation ∇ · δE(t, r) = 4πδρ(t, r). Consistently with the above discussion, we consider that for t > 0 the decay of δρ(t, r) is due to the spreading of charge-carrying particles (electrons and/or holes) via random walks on the underlying fractal distribution. The polarization response to δE(t, r) is given by where, as usual, The density of relaxation currents is defined as the time derivative of δP(t, r), i.e.
The continuity implies that ∂ ∂t δρ(t, r) Taking ∇· under the integral sign, and then eliminating δE(t, r) by means of Maxwell's equation ∇ · δE(t, r) = 4π δρ(t, r), we find, with the self-consistent charge density, In writing equations (9) and (10), we have also assumed that t > 0. Upon integrating over t, equation (10) becomes Here, the function φ(r) is an arbitrary function of the position vector r, which appears in the derivation as the constant of integration over time. Under the conditions χ(t − t ) ≡ 0 for t < t and δρ(t, r) ≡ 0 for t < 0 for all r, equation (11) reduces to If we allow t → +0, we find that for γ > 0 the integral term on the left-hand side goes to zero (as ∝ t γ ): 9 from which it is clear that φ(r) = lim t→+0 δρ(t, r). We consider this last condition as the initial condition for the relaxation problem. Essentially the same condition holds in the limit γ → 0, provided that lim t→+0 is taken first. A Fourier transformed equation (12) reads as where k is the position vector in reciprocal space and φ(k) is the Fourier image of φ(r). Writing the susceptibility as The quantity T λ has the sense of lifetime of a perturbation with wavelength λ. We expect that T λ ∝ λ z at criticality, where z is a scaling exponent. A derivation of this scaling relation will be given shortly. Separating the variables we write δρ(ω, is the Mittag-Leffler function [27]. The latter is the solution [28,29] of the fractional relaxation equation is the so-called Riemann-Liouville derivative [30]. Partial cases of this derivative are the unity operator for γ → 1, and ∂/∂t for γ → 0. One concludes that the relaxation to SOC of a slightly supercritical state is described by the Mittag-Leffler function E γ [ − (t/T λ ) γ ], and not by a simple exponential function as for standard relaxation. It is noticed, following [29], that the Mittag-Leffler function E γ [ − (t/T λ ) γ ] describes the relaxation toward equilibrium of particles governed by the fractional diffusion equation where (t, r) is the probability density of finding a particle (random walker) at time t at point r, and the Laplacian operator stands for the local (nearest-neighbor) character of the lattice interactions. Being intrinsically nonergodic [31], the fractional diffusion equation, equation (19), aims to incorporate the trapping effect [32] caused by cycles and dead ends of the fractal structure on which the transport processes concentrate (i.e. the 'nodes-links-blobs' model) ( [14] and figure 1 therein). This trapping effect appears in the fact that for γ < 1, the mean squared displacement of charge carriers grows slower-than-linear with time: If we recall that γ = 1 − η, then, with the aid of the expression η = µ/(2ν + µ − β), we find This scaling behavior reiterates one known result [21] for the random walks on percolation systems. A subtle point here is that · · · incorporates the system average over all clusters at percolation including the finite clusters conformally with the implication of the average ac conductivity response, σ ac (ω) ∝ ω η . Note that it is in this 'system-average' sense that the above fractional diffusion and relaxation equations apply to percolation. For times t T λ , the Mittag-Leffler function, equation (16), can be approximated [28] by a stretched-exponential, the so-called Kohlrausch-Williams-Watts (KWW) relaxation function [33] The KWW function, in turn, can conveniently be considered [9] as a weighted average of the ordinary exponential functions, each corresponding to a single relaxation event in the system: The weighting function w γ ( t) is given by equations (51d) and (55) of [9], where one replaces the exponent α with γ , the time constant T with T λ and the variable µ with T λ / t. In our notation where L γ ,−1 is the Lévy distribution function with skewness −1 (e.g. [34]). Assuming a longwavelength perturbation (i.e. the parameter λ being much longer than the microscopic lattice distance: λ 1), and setting T λ / t 1, we can further approximate the Lévy distribution L γ ,−1 by the Pareto inverse-power distribution. This gives L γ , leading to a pure power-law distribution of relaxation times, consistently with the wisdom of SOC: Our conclusion so far is that the relaxations are multi-scale, in accordance with equation (21), and their durations are power-law distributed. The distribution is heavy-tailed in the sense that , the limit T λ → ∞ can only be satisfied for p → p c implying that the system must be sufficiently close to the ideal SOC state. Distributions of the form of equation (23) have been earlier postulated [4] for SOC.

Consistency check
In a basic theory of dielectric relaxation, one writes the frequency-dependent complex dielectric parameter as [9,35] ( where ψ(t) is the relaxation function that describes the decay of polarization after the polarizing electric field has been stepped down or removed instantaneously. In the DPRW model, a stepdown-type electric field occurs as a consequence of re-injection of the free charges to the infinite cluster. The ensuing relaxation dynamics are mimicked by the chain reactions of hole hopping that act to suitably redistribute the polarization charges across the lattice. Based on the above analysis, we identify the relaxation function in equation (24) (24), after simple algebra one obtains where s = ωT λ is dimensionless frequency, and Q(s) and V (s) are the Lèvy definite integrals: In the parameter range of multi-scale relaxation response, T λ / t 1, ωT λ 1, the following series expansions of the Lèvy integrals apply [9]: From equations (28) and (29), one can see that the expansion of Q(s) starts from a term that is proportional to s −(1+γ ) , and so does the expansion of V (s) − 1/s. Hence, up to higher-order terms, (ω) − 1 ∝ s −γ . Given this, one applies the Kramers-Kronig relations sQ(s) ∝ σ ac (ω)/ω and to find the frequency scaling of the ac conduction coefficient to be σ ac (ω) ∝ ω 1−γ . By comparing this with the above expression σ ac (ω) ∝ ω η , one reiterates that consistently with the distribution of durations of relaxation events, equation (23).

Dispersion-relation exponent, Hurst exponent and the τ -exponent
In sandpile SOC models, one is interested in how the lifetime of an activation cluster scales with its size [5]. In the DPRW model, by activation cluster one means a connected cluster of activated sites. An occupied site is said to be 'activated' if it has become a hole or if it contains a free charge. Clearly, activation clusters can only exist above the percolation threshold. Note that activation clusters are subsets of the underlying conducting cluster of polarization charges. The notion of activation cluster is but a visualization of the charge density inhomogeneity δρ(t, r) in terms of a connected distribution of activated sites. Activation clusters decay because the constituent charged particles (holes and/or free charges) diffuse away via the random walks.
Consider an isotropic activation cluster composed of free particles. (The nature of the particles does not matter here-the hole case is similar.) It is assumed for convenience, without loss of generality, that each site of the activation cluster contains only one particle. Thus, the number density of the free particles inside the activation area is equal to one. It steps  [14].

Exponent
Expression down to zero just outside. If the microscopic lattice distance is a (a = 1), then there is a unit density gradient across the boundary of the activation cluster looking inside. Because of this gradient, the activation cluster will be losing particles on average. A particle that has crossed the boundary against the direction of the gradient is considered lost from the cluster. As the particles dissipate, the location of the density pedestal shifts inward with speed u. The local flux density of those particles leaving the activation area per second is just the gradient times the local diffusion coefficient. The latter depends on the frequency of the relaxation process as D(ω) ∝ ω η , in accordance with equation (2). If l is the current size of the cluster, then the corresponding relaxation frequency is ω u/l. Using this, the frequency dependence of the diffusion coefficient can be translated into the corresponding l-dependence, the result being D(l) ∝ l −η . Balancing the rate of decay of the cluster with the outward flux of the particles, we write dl/dt ∝ −l −η . Integrating this simple equation over time from t = 0 to t = T λ and over l from l = λ to l = 0, one finds the dispersion relation T λ ∝ λ z with z = 1 + η = 2 − γ . The persistency [36] of relaxation is measured by the Hurst exponent H , which is related to our z via H = 1/z. Lastly, the τ -exponent, which defines the distribution of particle flows caused by a single chain reaction, is obtained from equation (5) of [4], where one replaces φ with α and the fractal dimension D with the fractal dimension of the infinite percolation cluster, d f = d − β/ν. The end result is τ = 3 − αz/d f . Note that the τ values in [4] and [5] differ by 1. Using known estimates [12]- [14] of the percolation indices β, ν and µ, we could evaluate the critical exponents of the DPRW model in all ambient dimensions d 1. The results of this evaluation, summarized in table 1, are in good agreement with the reported numerical values from the traditional sandpiles (for d = 2, z ≈ 1.29, τ ≈ 2.0; for d = 3, z ≈ 1.7, τ ≈ 2.33) [4] and earlier theoretical predictions (for d = 2, z = 4/3, τ = 2; for d = 3, z = 5/3, τ = 7/3) [5]. We consider this conformity as a manifestation of the universality class of the model. For d = ∞, the model reproduces the exponents of mean-field SOC [3].
Equation (5), when account is taken of the η value at criticality, table 1, yields a lowfrequency percolation scaling for particle diffusion on a time-varying fractal structure for all d 1. As a partial case d = 2, it confirms the behavior D(ω * ) ∝ ω * Q 2/3 * predicted in [37,38] (up to the small deviation between 2/3 and ≈ 0.66). More so, the DPRW model gives a Hurst exponent (for d = 2, H ≈ 0.75; for d = 3, H ≈ 0.6) consistently with the reported narrow range of variation of H as observed in different magnetic confinement systems (Hurst exponent varying between H ≈ 0.62 and 0.75) [39]- [41]. In this respect, we also note that some connection between anomalous transport by plasma turbulence, fractional kinetics and SOC has been addressed on a phenomenological level by Carreras et al [42]. After all, the 13 DPRW model gives a value of the exponent γ (for d = 2, γ ≈ 0.66; for d = 3, γ ≈ 0.4) in good agreement with the typically found values (between 0.3 and 0.8) [9,10,43,44] from dielectric-relaxation phenomena in disordered amorphous materials. This last observation supports the hypothesis [32] that dielectrics exhibiting stretched-exponential relaxations are in a self-organized critical state.

Discussion
Apart from the details of the mathematical formalism, the DPRW model is actually quite simple. The main points are as follows. A lattice site can be either empty or occupied. An occupied site is interpreted as a polarization charge. The equilibrium concentration of the polarization charges depends on the potential difference between the plates. When the potential difference changes, the lattice occupancy parameter adjusts. A dynamical mechanism for this uses holes. The holes are just missing polarization charges. They are important key elements of the model as they provide a mechanism for the polarization current in the system. Beside holes, free charges are introduced. The free charges, too, carry electric current whose very specific role in the model is just to control the potential difference between the plates. The changing amount of free charges in the system has an effect on the lattice occupancy parameter. Nonlinearly, it affects the conductivity of the lattice. This nonlinear twist provides a dynamical feedback by which the system is stabilized at the state of critical percolation. The existence of such feedback proves to be an essential ingredient of the SOC phenomenon [45,46]. In many ways the proposed model is but a simple lattice model for dielectric relaxation in a self-adjusting disordered medium. It is perhaps the simplest model that accounts for the whole set of relaxation processes including hole conduction.
It is worth assessing the advantages and disadvantages of the DPRW approach to SOC. In terms of advantages, the electric nature of the model greatly facilitates analytical theory: not only does it permit quantification of the microscopic lattice rules in terms of the frequencydependent complex ac conductivity, the use of the Kramers-Kronig relation in equation (6) makes it possible to directly obtain the susceptibility function by integrating the conductivity response. As a result, exponents z, γ , α and H are expressible in terms of only one parameter, the exponent of ac conduction η. The latter is obtained as a simple function of the percolation indices β, ν and µ.
With respect to disadvantages, the model is seemingly different from traditional approaches to SOC based on cellular automation (CA), and its integration in the existing family of SOC models might be a matter of debate. Even so, the idea of random walks on a self-organized percolation system as a simplified yet relevant model for SOC has significant appeal. Firstly, it relies on the established mathematical formalism of random walks [19,20,28,47] whose advance on the SOC problem is theoretically very beneficial. Secondly, it offers a clear connection to studies outside the conventional SOC paradigm such as, for instance, to transport of mass and charge in disordered media [15,22,23,48]. Instead, the traditional CA-type models are complicated by a poor analytical description of the microscopic transport mechanisms and their basic physics appreciation is at times uneasy.
It is theoretically important to note that the dielectric context of the considered model, apart from offering a convenient platform for analytical theory, is not however essential for the SOC phenomenon. Indeed the DPRW model could be defined in terms of diffusion processes for neutral particles of different kinds. A formal reason for this is the equivalence [22,23] of 14 the frequency-dependent electrical conductivity problem and the frequency-dependent diffusion problem, specific to hopping conduction. The crucial element of the model is the assumption of random walks, not the nature of the particles.
In the DPRW model, the critical percolation state is made self-organized via the mechanisms of hole hopping by which the system responds to the fluctuating potential difference on the capacitor. We should stress that the holes redistribute the polarization charges so as to preserve the properties of the random percolation. They change the shape and the folding of the percolation clusters in the ambient configuration space, but not the random character in the distribution of the conducting sites. This random character contrasts those models of SOC based on invasion percolation (see [3] for a brief review). Note that invasion percolation [49] has different implications as it connects to phenomena with an extremum property such as for instance to self-organized growth phenomena and diffusion-limited aggregation [20]. We do not consider those models here.
The possible generalizations of the DPRW model correspond to biased random walks of free charges in the direction of the potential drop and/or inclusion of a second critical threshold p cc p c above which the random walk dynamics might change to a biased motion. We consider those generalizations obvious as they mainly intend to modify the value of the exponent η in a certain parameter range, while the basic physical picture of SOC will remain essentially the same.
The observed connection to the fractional diffusion equation, equation (19), may be interpreted, with the aid of the proposed SOC model, in favor of considering SOC as one important case for fractional kinetics [29,50]. It is understood that, in the DPRW model, correlations due to the fractional time derivative 0 D 1−γ t coexist along with the essentially random nature of microscopic particle motions. In many ways, the correlations build up, because the random walks of particles are squeezed to a low-dimensional structure with fractal geometry. Some discussion of these properties can be found in [38], where an approach based on pseudochaotic Hamiltonian dynamics has been suggested.
The final point to be addressed here concerns the issue of universality class. We take notice of the fact that the DPRW SOC model uses the charitable redistribution rule [51] to propagate the activities, similar to the traditional BTW sandpile [2] or others [4,5]. This means that an active site always loses its content to the neighbors. The charitable rule is to be distinguished from the neutral rule, when each of 2d + 1 sites involved in redistribution gets an unbiased random share of the transported quantity. Models using the neutral rule often fall in the universality class of directed percolation (DP) and are characterized by appreciably larger values of the dynamic exponent z (for d = 2, z ≈ 1.73 ± 0.05) [51]. Based on this evidence, we suggest that the DPRW model belongs to the same universality class as the BTW sandpile, and not to the DP universality class, consistent with the values of the critical exponents collected in table 1.

A model for the instability cycle
As the probability of site occupancy p approaches the percolation threshold p c , the pair connectedness length (i.e. the size of the largest cluster) diverges as ξ ∝ | p − p c | −ν . For p > p c , the longest relaxation time in the system is T ξ ∝ ( p − p c ) −zν and the dynamic susceptibility goes to infinity as χ ∝ ( p − p c ) −zνγ . In the DPRW model, p as a function of time is determined dynamically from the competing charge deposition and loss processes. That is, d p/dt = Z + − Z − , where Z + is the net deposition rate of free charges on the capacitor's left plate and Z − is the particle loss rate. The net deposition rate, or the driving rate, is the control parameter of the model: it takes a given value. The particle loss rate is obtained as electric current in the ground circuit, i.e. Z − = I θ( p − p min ), where p min is the lower limit of variation of p. Note that Z − is due to the free particles leaving the system through the earthed plate. The Heaviside θ function indicates that the lattice can release charges only when/if p p min . We expect that p min lies close to, although somewhat lower than, the percolation threshold p c . This is because the conducting cluster can still lose its charge content to the ground circuit even in the absence of a connecting path to the charged plate. The dynamics of I can be estimated from dI /dt = ±I /T ξ , where the upper sign corresponds to the relaxation process in the lattice. Putting all the various pieces together, we write where sign( p − p c ) = +1 (−1) for p > p c ( p < p c ) is the sign-function, and W is a numerical coefficient. Equations (32) and (33) define a simple system of equations for two cross-talking variables, the lattice occupancy per site and the particle loss current. These model equations are perhaps the simplest nonlinear equations describing the generic fueling-storage-release cycle in driven, dissipative, thresholded dynamical systems. An examination of these equations shows that the dynamics are periodic (auto-oscillatory), with the peak value of electric current Here, p max ( p max > p c ) is the upper limit of the p variation. Note that I max → 0 for p max → p c as expected. The auto-oscillatory motions signify that the pure SOC state is destabilized and that the system's phase trajectories enter the supercritical parameter range. When p max → 1, the periodic dynamics acquires a sharp, bursting character. The bursts half-duration (or half-width) equals (1/W )( p max − p c ) −zν . Eliminating the distance to the critical state, one obtains the scaling relation I max ∝ −ζ , where ζ = (zν + 1)/zν. The period between the bursts is found to be b ( p c − p min )/Z + . A pure SOC state with no superimposed periodic bursts arises when b → ∞. This implies that Z + → 0 for p → p c . Thus, criticality requires the vanishing of Z + , in agreement with the result of [3]. The critical state is stable when b T ξ . We have This is satisfied when Z + → 0 faster than This limit exceeded, the system starts to auto-oscillate around the percolation point with a period dictated by the net deposition rate of the polarization charges. The physical origin of this autooscillatory motion lies in the fact that the changing number of free particles provides feedback to the lattice occupancy parameter. It is because of this feedback relation that the DPRW system operates as a self-adjusting, intrinsically nonlinear dynamical system. Whether or not this feedback will excite the instability depends on how the characteristic driving time compares to the characteristic relaxation time. Indeed, focus on the stability condition in equation (34). The system being stable at the percolation point requires that the relaxation time due to the random walks T ξ be short compared to the characteristic driving time, 1/Z + . In this parameter range, any occasional charge density perturbations will dissipate via random walks before input conducting sites are again introduced. When the percolation point is approached, because the time scale T ξ ∝ | p − p c | −zν diverges, it is essential that the system be driven infinitesimally slowly to remain in a pure SOC state. Instability occurs when the relaxation processes operate on a longer time scale than the driving processes. In this regime, the system accumulates the polarization charges 3 , whereas to remain at criticality it should get rid of them. The accumulation of the polarization charges has a direct effect on the conductivity between the plates, which steps up with the lattice overshooting the percolation threshold. When p max → 1, the system can be thought of as facing the typical conditions of electrostatic discharge in the regime of short circuit. It should be emphasized that the feedback mechanism does a twofold job: (i) it stabilizes the system at the state of critical percolation in a regime when the driving rate is infinitesimal and (ii) it excites cross-talk between the conductivity and the lattice occupancy parameters when the driving rate is faster than the relaxation rate.
In the parameter range in which the strength of the driving vanishes, the multi-scale geometry of the critical percolation is dominant in providing the major transport characteristics for the DPRW lattice. The situation changes drastically when the strength of the driving increases above some level. With the system's departure away from the percolation point, the multi-scale features will soon be lost substituted by the bulk-average nonlinearities. The fact that equations (32) and (33) above are formulated in terms of the system-average parameters, p and I , merely reflects that the system is allowed to appreciably depart from the state of marginal stability, or the SOC state (this means that p max can be rather closer to 1 than to p c ), and that the effect of overdriving readily calls for global features to come into play. It is noted that, in general, the multi-scale properties due to SOC can coexist along with the global or coherent features. We illustrate this type of mixed SOC-coherent behavior on an example below: substorm behavior of the dynamic magnetosphere.
The end result of the discussion above is that the strength of the driving plays a crucial role in dictating both linear and nonlinear behavior in the DPRW model. To obtain a pure SOC state, the driving rate should go to zero fast enough as the critical point is approached. The main effect that overdriving has on the DPRW dynamics is to excite unstable modes associated with periodic bursts in the particle loss current. Accordingly, the system auto-oscillates between a subcritical ( p min < p c ) and a supercritical ( p max > p c ) state in response to external forcing. These dynamical properties are summarized in figure 3.
The transition to auto-oscillatory dynamics signifies the increased role of global and nonlinear behavior in the strongly driven DPRW system as compared to a pure SOC system. The borderline between the two regimes corresponds to T ξ ( p c − p min )/Z + . The stability condition T ξ ( p c − p min )/Z + has serious implications for the achievable SOC regimes. It imposes one important restriction on the net deposition rate of free particles against the longest relaxation time on the incipient percolation cluster.

Why fishbone-like instability
To help judge the result obtained, let the critical exponents take their mean-field values: z = 2, ν = 1/2. In this limit, equations (32) and (33) above reproduce, up to a change of variables, equations (13) and (14) of [8]. The latter set of equations appear in the basic theory of Alfvèn instabilities as a simple model for the coupled kink-mode and trapped-particle system in a magnetically confined toroidal plasma where beams of energetic particles are injected at high power. The mode dubbed 'fishbone' is characterized by large-amplitude, periodic bursts of magnetohydrodynamic (MHD) fluctuations, which are found to correlate with significant losses of energetic beam ions [8]. By comparing the two sets of equations, one can see that the lattice occupancy per site p corresponds to the effective resonant beam-particle normalized pressure within the q = 1 surface (q is the familiar safety factor used in tokamak research), p c corresponds to the mode excitation threshold, and the particle loss current I corresponds to the amplitude of fishbone. This direct correspondence between the two models suggests considering the instability in equations (32) and (33) as analog 'fishbone' instability for SOC dynamics.
This correspondence is not really surprising. Mathematically, it stems from the resonant character of fishbone excitation [52], implying that the energetic particle scattering process is directly proportional to the amplitude of fishbone [8]. This resonant property dictates a specific nonlinear twist to the fishbone cycle, differentiating it from other bursting instabilities in magnetically confined plasmas. It is this 'resonant' twist observed in the DPRW model system that identifies the analog 'fishbone' mode for SOC. The instability is global in that it involves oscillations of system-average parameters as a result of the lattice conductivity properties nonlinearly changing with the varying strength of the driving and is accompanied by a transition from weak, subdiffusive to strong, ballistic transport on the lattice, similar to the plasma confinement case [53,54].
We reiterate that the DPRW instability cycle exhibits characteristic aspects of strongly driven MHD instabilities in magnetic confinement devices [52,54], and its occurrence is a reminder of the basic physics of resonant fishbone excitation. Conversely, the nonlinear science of fishbone [55] might be thought of as being connected with those properties generic to strongly driven, interaction-dominated, thresholded dynamical systems such as the DPRW supercritical system, rather than unique to toroidally confined fusion plasmas.

Coexistence of self-organized criticality and coherent features: substorm behavior of the dynamic magnetosphere
The idea of fishbone-like instability in self-organized critical dynamics is very appealing as it addresses a type of behavior in which multi-scale features due to SOC coexist along with global or coherent features. One example of this coexistence can be found in the solar wind−magnetosphere interaction. It has been discussed by a few authors [56]- [58] that the coupled solar wind-magnetosphere-ionosphere system operates as an avalanching system and that there is a significant SOC component in the dynamics of magnetospheric storms and substorms, along with a coherent component [59] that evolves predictably through a sequence of clearly recognizable phases [60]. Here, we advocate a way of thinking [61] in which the SOC component is associated with the properties of self-organization of electric currents and magnetic field fluctuations in the plasma sheet of the Earth's magnetotail [62], whereas the coherent component is subordinate to the global instability of the SOC current system and the phenomenon of tail current disruption [63]. This means that the dynamic magnetosphere survives through a mixed SOC-coherent behavior.
In this spirit, we expect the input power due to magnetic reconnection at the Earth's dayside magnetopause to self-consistently control the system's departure from the state of marginal stability, or the SOC state, with stronger departures favoring coherent features. The magnetotail current system being close to SOC implies that the power spectral density of the magnetic fluctuation field is, in the parameter range of sufficiently low frequencies, given by a power-law distribution S(ω) ∝ ω −α , with α ≈ 1.3 (see table 1). The latter value falls within the typical range of variation of α as found from in situ satellite measurements [63]- [65]. With the input power going above a certain level, the dynamical system at SOC fails to accommodate the increased potential difference between the two flanks of the magnetotail. At this point, a portion of the cross-tail electric current will be redirected to the ionosphere (here considered as the analog 'ground circuit'), thus triggering a magnetospheric disturbance or a substorm. These processes are schematically illustrated in figure 4. The phenomenon of substorm bears signatures [62] enabling it to be considered as a second-order phase transition in the presence of a coexisting nonlocal symmetry [61], consistent with the description in terms of the fractional Ginzburg-Landau equation [66].
We should stress that we associate magnetospheric disturbances with instabilities on the top of the underlying SOC state, and not with the SOC behavior itself, contrary to the implication of [58]. In the above coupled system of equations, equations (32) and (33), we identify the driving rate, Z + , with the magnetic reconnection rate at the Earth's dayside magnetopause; the lattice occupancy parameter, p, with the average normalized energy density of magnetic fluctuation field in the magnetotail current sheet; and the particle loss current, I , with electric current in the ionosphere. Note that the period between the fishbone events (major magnetospheric disturbances) may vary with varying Z + . This is going to be the case in the realistic magnetosphere owing to the solar wind variability.

Summary and final remarks
In this paper, we have proposed a simple model of self-organized criticality, the DPRW model, which addresses the SOC problem as a transport problem for electric charges. The novel concepts of our study are (i) a theory of self-organized criticality based on the analogy with dielectric-relaxation phenomena in self-adjusting random media, and (ii) prediction of a 'resonant' instability of SOC due to the nonlinearities present. The system adjusts itself to remain at the critical point via the mechanisms of hole hopping associated with the random walk-like motion of lattice defects on a self-consistently evolving percolation cluster. The relaxation to SOC of a slightly supercritical state is described by the Mittag-Leffler function E γ [ − (t/T λ ) γ ], this being the solution [28] to the fractional relaxation equation, and not by a simple exponential function as for standard relaxation. The durations of relaxation events are power-law distributed, with diverging upper cut-off relaxation scale. The ideal SOC state requires that the driving rate goes to zero faster than a certain scaling law as the percolation point is approached. The model belongs to the same universality class as the BTW sandpile, and should be distinguished from the DP-like SOC models.
In a narrower sense, the DPRW approach to SOC offers a simple yet relevant lattice model for dielectric-relaxation phenomena in systems with spatial disorder. One by-product of this approach is the case for stretched-exponential, the so-called KWW relaxation function, equation (20), which is often found empirically in various amorphous materials such as for instance in many polymers and glass-like materials near the glass transition temperature (for reviews, see [10] and [67], and references therein). In this connection, we reiterate that the DPRW SOC model gives an exponent of the stretched-exponential relaxation (for d = 2, γ ≈ 0.66; for d = 3, γ ≈ 0.4) in good agreement with the experimentally observed values [9,10,44]. Overdriving the DPRW system near self-organized criticality has a destabilizing effect on the SOC state. The fundamental physics of this instability consists in the following. Because of rapid, in the sense that Z + Z +max ∝ | p − p c | zν , accumulation of the conducting sites, the system departs from the percolation point, and its geometry changes from fractal-like to crystalline-like. This means that the average conductivity inside the capacitor has greatly increased. As the lattice conducts more electricity, losses increase in the ground circuit. However, because the particle loss current has feedback on the lattice occupancy parameter, cross-talk is excited between the system's average conductivity response and the distance to the critical state. We have observed that the instability cycle is qualitatively similar to the excitation of the internal kink ('fishbone') mode in tokamaks with high-power beam injection (the lattice occupancy per site p corresponds to the effective resonant beam-particle normalized pressure within the q = 1 surface, p c corresponds to the mode excitation threshold and the particle loss current I corresponds to the amplitude of fishbone). The instability is 'resonant' in that the particle loss process is directly proportional to I . This resonant property dictates a specific nonlinear twist to the fishbone cycle, differentiating it from other bursting instabilities in magnetically confined plasmas.
We have discussed that the resonant instability, the fishbone, may have serious implications for the functioning of complex systems. The existence of an instability on the top of SOC dynamics conforms with the results of [68] in which the traditional (sandpile) SOC model has been modified by adding diffusivity, giving rise to periodic relaxation-type events as a function of the system drive, while a pure SOC state requires a vanishing drive 4 .
The excitation of fishbone-like instability in SOC systems leads to a type of behavior in which the multi-scale features due to SOC coexist along with the global or coherent features (i.e. mixed SOC-coherent behavior). One example of this coexistence is found in the solar wind-magnetosphere interaction. We expect the concept of mixed SOC-coherent behavior to be a plausible statistical picture for thresholded, dissipative, nonlinear dynamical systems in the parameter range of nonvanishing external forcing. In this respect, we speculate that some of the 'extreme' events, or system-scale responses, observed in complex natural and social systems [7,69] may, in fact, be the fishbone-like instabilities of SOC predicted by the present theory. It will be interesting to investigate if such phenomena as, for instance, the El Niño South Pacific Oscillation or the Atlantic Multi-decadal Oscillation might be envisaged as fishbone-like events in the Earth's SOC climate system as a consequence of ocean-atmosphere coupling [70,71]. In the same spirit, the periodic occurrences of glacial periods on Earth might perhaps be considered as globally induced unstable climate modes stemming from cross-talk between air temperature and dust concentration. Support for this suggestion can be found in Antarctic ice records as discussed in [72].
All in all, the fishbone phenomenon warns of repeating severe events being virtually unavoidable in driven systems.