Symmetries and Fixed Point Stability of Stochastic Differential Equations Modeling Self-Organized Criticality

A stochastic nonlinear partial differential equation is built for two different models exhibiting self-organized criticality, the Bak, Tang, and Wiesenfeld (BTW) sandpile model and the Zhang's model. The dynamic renormalization group (DRG) enables to compute the critical exponents. However, the nontrivial stable fixed point of the DRG transformation is unreachable for the original parameters of the models. We introduce an alternative regularization of the step function involved in the threshold condition, which breaks the symmetry of the BTW model. Although the symmetry properties of the two models are different, it is shown that they both belong to the same universality class. In this case the DRG procedure leads to a symmetric behavior for both models, restoring the broken symmetry, and makes accessible the nontrivial fixed point. This technique could also be applied to other problems with threshold dynamics.


I. INTRODUCTION
In the last decade much attention has been paid to the phenomenon known as self-organized criticality (SOC). Bak, Tang, and Wiesenfeld (BTW) [1] studied a cellular automaton model as a paradigm for the explanation of two widely occurring phenomena in nature: 1/f noise and fractal structures. Both have in common a lack of characteristic scales. The SOC models, although not always show 1/f noise, have no characteristic scales too, and this scale invariance suggests that these systems are critical in analogy with classical equilibrium critical phenomena, but in SOC one deals with dynamical nonequilibrium statistical properties. Moreover, the system evolves naturally to the critical state without any tuning of external parameters, that is, in a self-organized process.
Several cellular automata and coupled map lattices models exhibiting SOC have been reported in the literature. In the original sandpile model of Bak et al. [1] the system is perturbed externally by a random addition of sand grains. Once the slope between neighboring cells has reached a threshold value, sand is transferred between them in a fixed amount. Taking this model as a reference, different dynamical rules have been investigated leading to a wide variety of universality classes. Continuous variables with a full transfer from a cell instead of a fixed discrete amount [2][3][4][5], directed flows [6], threshold condition imposed on the height, on the gradient, or even on the laplacian [7], anisotropy [8], are a few examples. These randomly driven models do not exhibit SOC when the interaction rules are not conservative [9]. Later on, other deterministically driven models have been introduced where conservation is not a necessary condition [10][11][12][13][14][15]. Much more recently, sandpile models with deterministic perturbations but intrinsic randomness in the threshold dynamics have been used to reproduce experimental results on transport properties on rice piles [16]. The close connection between these sandpile models and interface depinning has been establish in Ref. [17].
Some authors have attempted to connect the randomly driven models to stochastic differential equations [18,19]. These continuous descriptions are built according to the symmetry rules obeyed by the discrete models in order to achieve a generic scale invariance condition [20]. Nevertheless none of them neither explicitly nor implicitly includes the threshold condition, which is one of the main characteristics of SOC models. On the other hand anomalous diffusion equations with singularities in the diffusion coefficient have been considered in order to study the deterministic dynamics of the avalanches generated in the critical state [21,22]. Latterly, a different approach has been introduced by Pietronero and co-workers, using a real-space renormalization procedure to determine the dynamical exponent as well as the avalanche size exponent. [23].
In a previous paper [24] one of us studied two nonlinear stochastic differential equations derived from the discrete dynamical rules of two models with different symmetry properties. In principle one would expect, due to this reason, different critical behavior. However, it was shown analytically, by means of the dynamic renormalization group (DRG) [25][26][27], that both models belong to the same universality class. The threshold condition was kept but the step function was regularized in order to allow a power series expansion. In the limit that recovers the threshold it was shown that the coupling constants that distinguish both models become decoupled of the common coupling constants; since the critical exponents only depend on the latter constants, one obtains the same values for both models. Once this equivalence was established, the most symmetric model was considered, showing that an infinite number of coupling constants was relevant below the upper critical dimension d c = 4; by expanding in the number of nonlinearities, the DRG procedure, up to first order in ǫ = 4 − d, gave an estimate of the dynamical exponent close to the value obtained by scaling arguments and in the numerical simulations [2,5]. This value of the dynamical exponent is obtained when the flow in parameter space reaches the nontrivial stable fixed point; nevertheless, when taking into account the physical values of the parameters, they do not lie in the basin of attraction of the fixed point, thus making this computed value in some sense speculative, since it cannot be ensured that the flow in parameter space will be able to reach the attractor.
Our goal in this paper is to complement the previous work in order to check the validity in the calculation of the critical exponents at the stable fixed points and to analyze the role played by symmetries in randomly driven SOC models and, in general, in other models where an infinite hierarchy of nonlinear terms is required. Our procedure also illustrates the effect of symmetry breaking in DRG calculations as a mechanism to make the attractors in parameter space accessible for the physical values of the parameters in the original equations. The continuum equation for the BTW and Zhang's models is introduced in Sect. II, as well as the alternative regularization which breaks the symmetry that distinguish both models. In Sect. III we develop the DRG procedure and show that this symmetry is irrelevant, in view of the fact that the nontrivial fixed point is not modified by this new approach. Moreover, the effect of symmetry breaking allows the flow of the original parameters to reach the nontrivial fixed point, where the critical exponents can now be computed. Finally we present our conclusions in Sect. IV.

II. MODELS AND SYMMETRIES
First, we describe briefly the dynamics of the two SOC models under consideration. The first model was originally proposed by Zhang [2], and it consists of a d-dimensional lattice in which any site can store some continuously distributed variable E. This variable, that we will call energy, can have different physical interpretations [3]. The system is perturbed by adding a random amount of energy δE > 0 at a randomly chosen site. Once a site reaches a value of the energy greater than the threshold value E c , this site becomes active, and transfers all its energy to its nearest neighbors. At this point the input of energy from the outside is turned off. The energy transferred to the neighboring sites can make them active, giving rise to an activation cluster or avalanche, that ends when all the sites have reached a value of the energy smaller than E c . It is only when the avalanche has stopped that energy is added again, otherwise the system remains quiescent. In this way there is a clear time-scale separation in the dynamics. The external noise acts in a slow time scale, whereas the avalanches evolve infinitely fast, in comparison. The second model differs from Zhang's model only in the amount of energy an active site transfers to its neighbors, which is a fixed amount E c , instead of its whole energy E. Therefore it is closer to the original sandpile model of Bak et al. [1], but continuous in E. When δE is not random but fixed this difference becomes irrelevant. For this reason, it will be referred as BTW model. Notice that both models are conservative, in the sense that the added energy (always positive) is only dissipated at the open boundaries.
The microscopic evolution rules can be written from time t (in the fast time scale) to t + 1 for each site i as and for Zhang's and BTW model, respectively. The sum runs over the q nearest neighbors of site i, labeled as nn, and the threshold condition enters through the Heaviside step function Θ, defined as Θ(x < 0) = 0 and Θ(x > 0) = 1. Due to the continuous nature of the models the value Θ(x = 0) is irrelevant and we can keep it undefined, by now. For the external noise ξ i (t), that drives the system only when there are no active sites, one can formally write where δ i,n(t) is the Kronecker delta symbol and n(t) is a random vector pointing the site of the lattice that will be perturbed with a random amount of energy δE (in the original BTW model δE = E c /q).
The product runs over all the lattice sites.
When applying the DRG one deals with infinite systems and then the important effect of dissipation at the open boundaries is not taken into account. However in SOC models a distribution of absorbing defects through the lattice plays the same role as the open (absorbing) boundaries [28], as we have verified through computers simulations [29]. Then, we can redefine our models in an infinite lattice but with a quenched distribution of defects. The results are not modified with this assumption. Another possibility is to consider that each site of an infinite lattice has a small probability of dissipating an amount of energy E c /q when it topples, instead of transferring it to a certain neighbor. This procedure, that represents the assumption of random boundaries, implies that when a site receives a toppling from some neighbor, it has a small probability of not to accept the amount of energy E c /q, that is lost [30]. This dissipation can be included as a new term in the noise, and Eq. (2) has to be replaced by where ζ nn is a dichotomous noise, taking value 0 with large probability and value E c /q (dissipation) with a small one. If ζ nn depends on t, i.e., ζ nn = ζ nn (t) we are dealing with annealed random boundaries, whereas if it only depends on the position, we have quenched random boundaries or absorbing defects. In terms of a rescaled energy E − E c → E, and introducing a parameter Z to unify the description, we have for both models where Z = 1 for Zhang's model and Z = 0 for BTW. Equation (4) defines a stochastic coupled map lattice. Moreover, notice that the deterministic BTW equation displays invariance under a parity transformation of the order parameter, E → −E. This one is the only symmetry that the BTW model does not share with Zhang's model. The common symmetries are invariance under spatial translations, rotations, and reflections, as well as conservation of the order parameter. Equation (4) can be coarse-grained in order to obtain a continuum equation for the effective E( r, t). Then, by using the prescriptions for the temporal derivative and for the Laplace operator: where α is a coefficient that depends on the lattice spacing, the unit time step and the coordination number q. The noise η( r, t) accounts for the effective external noise as well as for the internal noise that appears due to the elimination of microscopic degrees of freedom.
Up to now, Eq. (5) describes truly the coarse-grained evolution of the system, but we have not yet characterized the noise η( r, t), which derives from ξ i (t). The product in Eq. (3) makes η( r, t) to be a multiplicative noise that depends on the whole lattice state, and the problem is intractable. We are going to ignore the restrictions imposed by the step functions in Eq. (5), breaking the time scale separation. Then the noise η( r, t) acts continuously in time and can provoke avalanches to overlap. However for small enough noise this is very unlikely, and one can still identify avalanches in computer simulations. Moreover the dynamical exponent does not change with this assumption [5], due to the fact that the added noise is orders of magnitude smaller than the energy transferred by the avalanche and then its dynamics is not affected. In this case we still have two time scales, although they overlap. In what follows η( r, t) will be considered as an additive random process, including two effects, the external driving, always positive, and the dissipation at the (random) boundaries, always negative. In the statistical stationary state the random input of energy must equal on average the output at the boundaries. Then we assume that < η( r, t) >= 0. (6) In fact this is the same assumption done in all the studies of SOC by means of DRG [18,19], and it is somehow equivalent to the stationary condition used in Ref. [23]. Moreover we are mainly interested in the spatio-temporal propagation of a perturbation through the system, that is, in measuring the value of the dynamical exponent. For this purpose we have to look at the system in the fast time scale, i.e., the scale of the evolution of the avalanches. In Ref. [5] it was argued that in this case one can understand the noise as a quenched Gaussian process uncorrelated in space, and then its correlation function is given by When looking at the system at the slow time scale one cannot use this prescription for the noise, which has to be uncorrelated in time too, i.e., < η( , and this prescription is mainly related to the interface roughness between avalanches [31]. Equations (5), (6), and (7), together with the fact that the noise is a Gaussian process, define completely our model. However, the presence of the step function in Eq. (5) gives rise to a strong nonlinearity. A perturbative expansion of this equation can be performed if one regularizes the step function as and makes a series expansion of f (βE) in powers of E [22,24]. The function f (x) must be monotonously increasing with f (−∞) = 0 and f (∞) = 1. Moreover we choose f (x) − 1/2 as an odd function, so f (0) = 1/2. Several functions of this type have been used in the literature, but that coming from the error function as allows a power expansion that has an infinite radius of convergence, in contrast with previous choices [22,24]. In any case, the relevant results do not depend on the particular form of f (x). The regularization given by Eq. (8) keeps the symmetry of the step function and therefore the invariance under a parity transformation in the BTW model. As an alternative regularization that breaks this invariance we propose the following: with K an arbitrary constant. Although in the limit β → ∞ we recover the step function, we do not recover its symmetry anymore, because Θ(E = 0) = f (K) = 1/2 if K = 0, and this is the reason for the breaking of the symmetry in the BTW model. Now we perform a series expansion of the regularizing function f (βE) in powers of E, obtaining where the coefficients a n (β, K) are given by being f (n) (K) the n-th order derivative of f (x) at x = K. Substituting the expansion (11) into Eq. (5) we can write where the effective diffusion constant D and the coupling constants λ n (that make the equation nonlinear) take different values depending on the model: On the one hand, for K = 0, since all the even derivatives verify f (2n+2) (0) = 0, all even coupling constants vanish for the BTW model, whereas they do not for the Zhang's one. Using Eq. (13) this fact allows to verify the symmetry of the BTW model under the parity transformation of the order parameter E. On the other hand, for K = 0, the even coupling constants do not vanish in any case, and this fact constitutes the symmetry breaking for the BTW model. Then, under this condition, the only difference between both models is that the constants depend on β in a different way; however it is easy to see that in the limit β → ∞ both sets of constants are identical and then the Zhang's model and the broken-symmetry BTW model have to belong to the same universality class. This fact can only be shown for K = 0. Nevertheless, the fact of considering K = 0 only introduces a difference in the value of Θ(0), which is irrelevant in a continuous model, and then one can include the (symmetric) BTW model in this universality class too. At this point it is worth noting that we have transformed a stochastic coupled map lattice, that involves a threshold condition and presents a clear separation of time scales, into a nonlinear stochastic partial differential equation, where the nonlinearity of the threshold is described by an infinite series of powers and the randomness enters via a Gaussian process, with zero mean to account for the dissipation at the boundaries. Along this transformation, and due to the approximations we have performed concerning the noise correlation, we have broken the time scale separation since the noise acts constantly in time. Nevertheless, we expect that such an equation explains the dynamical properties of the system within the fast time scale of the propagation of the avalanches. As we have mention before and it is discussed in [31], to deal with the slow time scale, where the avalanches are instantaneous, another noise correlation is more appropriate.

III. DYNAMIC RENORMALIZATION GROUP PROCEDURE
The model to be studied by the DRG is defined by the nonlinear partial differential equation (13) and the Gaussian noise given by (6) and (7). As a previous step we can check the relevance of the different coupling constants in this equation by naive dimensional analysis: a change of scale b = e l > 1: is performed in Eqs. (13) and (7), being χ the roughness exponent, which is related to the hydrodynamic exponent, and z the dynamical exponent. Then one obtains that the parameters transform as Under this scaling transformation, z and χ are chosen to keep the linear model scale invariant, i.e., the parameters D and Γ have not to be modified. This choice gives z = 2, χ = (4 − d)/2, and Then one can see that when we apply iteratively the transformation (b → ∞) for d > 4 all the nonlinear terms vanish and they are irrelevant. However, all the coupling constants go to infinity for d < 4, and hence all nonlinear terms become relevant; this implies that the upper critical dimension is d c = 4, and nontrivial values of the exponents are expected below it. The relevance of all the terms makes our problem much more complicated than for instance the Kardar-Parisi-Zhang model of interface growth, where only the first nonlinear term is relevant [32]. The appropriate treatment of Eq. (13) would be to renormalize the infinite number of relevant coupling constants that are involved. Of course this is impossible to do in practice. In [24] an expansion in the number of coupling constants for the BTW model was performed with only odd terms, i.e., without symmetry breaking (K = 0). The critical exponents where obtained as a function of the highest coupling constant, up to λ 9 . Fortunately, the dynamical exponent was well behaved and could be extrapolated up to λ ∞ . However, keeping the symmetry of the step function, the nontrivial fixed point of the DRG is unreachable using the parameters given by Eqs. (14) and (15), even for the Zhang's model. We want to show that with the proposed alternative regularization of the step function, that breaks the symmetry of the BTW model and allows the existence of even coupling constants, the DRG fixed points are not changed but now they are accessible to the flow when the parameters take their real values. For this reason, and as a first attempt to justify our hypothesis as well as the conclusions of Ref. [24], we will focus on Eq. (13) with only its first two nonlinear terms, i.e., λ 2 and λ 3 , and see how they behave under a DRG transformation, The DRG procedure consists of the removal of the fast modes (large wavenumber k) in the momentum space, followed by a rescaling of a factor e l in order to recover the original Brillouin zone [25][26][27]. After this transformation, one obtains an equation that is equivalent to the original one but with different (effective or renormalized) coefficients. Successive iterations of this transformation give the flow of the coefficients in the parameter space. If this flow converges towards a fixed point, the system presents "scale invariance" in the hydrodynamic limit (large-distance and long-time behavior). Then, the fluctuations of the order parameter verify the scaling equation where the critical exponents χ and z are those that ensure the existence of the fixed point. However, it is worth mentioning that with this procedure the scaling function F (x) remains unknown [33]. We now outline the DRG calculation. First of all we write Eq. (19) in Fourier space: Here E( k, ω) and η( k, ω) are defined as the Fourier transforms of E( r, t) and η( r, t), i.e., whereas is called the bare propagator. The symbol " * " represents the convolution product, defined as Fig. 1 (a) shows the expression of Eq. (21) in terms of Feynman diagrams. As the intensity of the noise Γ has also to be renormalized by the DRG transformation, we need to consider the equation for the correlation function of the transformed energy < E( k, ω)E( k ′ , ω ′ ) >, which, up to one-loop order, is: where the prime denotes a dependence on k ′ , ω ′ instead of the dependence on k, ω. The diagrammatic representation of this equation is shown in Fig. 1(b). Eqs. (21) and (25), that are the ones that we are going to renormalize, hold for 0 < k < Λ, where Λ is the wavenumber cutoff due to the underlying discrete structure. The transformed noise η( k, ω) turns out to be also a Gaussian process with zero mean, but with correlation  21) and (25), defined in the range 0 < k < Λ. The double bar with the cross × at its end is the order parameter E, the single bar with the cross represents G0η, whereas the single bar alone is G0. A vertex with n branches (n = 2 or 3 in the figure) represents a convolution product of n elements, including a prefactor −λnk 2 /(2π) (n−1)(d+1) . The circles correspond to the average over the noise.  Fig. 1 allows to define new coefficients. Observe that the new averages only affect the outer shell. The notation has been simplified in respect to Fig. 1, suppressing the symbol × at the end of the vertices and also the arrows.
The first step of the DRG transformation consists in splitting the Fourier space in two shells: an inner shell, that contains the slow modes, i.e., 0 < k < e −l Λ, and an outer shell, containing the fast modes, e −l Λ < k < Λ. Both modes are coupled through the convolution products in (21) and (25). We consider the diagrams for the slow modes and perform a perturbative expansion of the fast ones up to the lowest order in the intensity of the noise (see Appendix for more details). Then we integrate out these modes by an average over the noise in the outer shell. After this transformation the resultant equations are shown diagrammatically in Fig. 2. It is clear that we can obtain new equations which are formally equivalent to the initial ones, (21) and (25), defining the new coefficients as the original ones plus the corresponding integrals over the outer shell. With the noise correlation (26) these integrals can be easily computed in the hydrodynamic limit ( k → 0, ω → 0), as it is shown in the Appendix, and then the coefficients transform according to where and S d is the complete solid angle in d dimensions. However, the new equations are only defined in the inner shell 0 < k < e −l Λ. The second step allows to recover the original Brillouin zone 0 < k < Λ by rescaling the equations using transformation (16), which in Fourier space writes The combined effect of both transformations, in the limit l → 0, constitutes an infinitesimal DRG transformation, which gives the flow equations of the parameters in parameter space. In these flow equations instead of λ 2 and λ 3 it is suitable to use the dimensionless coupling constantsλ 2 andλ 3 , given byλ where I (1) We are interested in the invariance of the parameters under DRG transformations. This means that we have to look for the fixed points of the flow equations; if we write Eqs. (31) as dα i /dl = g i (...α j ...), where α j represents any coefficient, then the fixed points verify g i (...α * j ...) = 0, ∀i. Considering D = 0 and Γ = 0 we obtain four algebraic equations with four unknowns, χ, z,λ 2 , andλ 3 ; their solutions will give us the fixed points of the transformation,λ * 2 andλ * 3 , as well as the values of the exponents z and χ that guarantee that the DRG transformation leads to a scale-free behavior. Notice that the particular values of Γ and D play no role in the existence and location of the fixed points. We can also find the stability of the fixed points under small perturbations using a linear stability analysis: the fixed point {α * j } is stable (i.e., an attractor) if all the eigenvalues (or their real parts) of the matrix ∂g i /∂α j evaluated at this fixed point are negative.
The results are the following: for d > 4 one obtains six different fixed points, but the only stable one corresponds to where as usual ǫ is defined as ǫ = d c − d = 4 − d. This is the trivial or Gaussian fixed point, that gives a normal (or Brownian) diffusive behavior because of the vanishing of the coupling constants. The values of the exponents do not correspond with those of the Edwards-Wilkinson model, used in the study of surface growth, because the noise correlation is different [34]. For d < 4 this fixed point becomes unstable and the only stable one is that was unstable for d > 4. In this case the diffusion is anomalous, to be more precise, the fact that z < 2 gives a superdiffusive behavior in the hydrodynamic limit. Note that the one-loop expansion in the intensity of the noise Γ gives a nontrivial fixed point which is expressed as a perturbation of the Gaussian one in a first order ǫ-expansion. Observe also that the breaking of symmetry does not modify the value of the fixed point obtained without taking into account the even coupling constant λ 2 [24]. Moreover, the fact that the nontrivial fixed point is an attractor of the dynamics contrasts with equilibrium critical phenomena, where this point is stable only along one direction. In this fact lies the difference between fine tuning of parameters for equilibrium systems at the critical point and self-organization towards criticality for nonequilibrium processes. Now we know the attractors in the parameter space, but this is not enough in our case; since our stochastic equation (13) is derived directly from the discrete rules of the BTW and Zhang's models, we also need to know the basins of attraction of the stable fixed points and whether our initial conditions, that is, the initial values of the coefficients corresponding to our physical problem, are inside these basins. These values for the dimensionless coupling constants (30) can be calculated from Eqs. (14) and (15), and they are result that also holds for K = 0, where we obtainλ o 2 = 0 even for the Zhang's model. The superscript "o" indicates the initial value of the coefficient, that is, its value before any renormalization. As we have no restriction for Γ o (except that it has to be small), and K can take any arbitrary real value, this implies that the initial dimensionless coupling constants will be defined in the following region: having used for f (x) the explicit form given by Eq. (9). Clearly, the stable fixed point for d < 4 (33) is outside the region of initial conditions defined by Eq. (35). It will be of the maximum interest however to know whether these conditions will drive the system towards the nontrivial fixed point or not. We first consider K = 0, which impliesλ o 2 = 0, corresponding to the case studied in Ref. [24]. For d < 4 one gets a different behavior depending on the sign ofλ o 3 . Fig. 3 shows that whenλ o 3 is positive it flows towards the stable fixed pointλ * 3 = ǫ/27 giving a dynamical exponent z = 2 − ǫ/9. A negativeλ o 3 , which is our case of physical interest, flows away. An exact solution of Eq. (31d) withλ 2 = 0 gives thatλ 3 would reach −∞ in a finite l and then would reappear asλ 3 = ∞, being then under the attraction of the nontrivial fixed point. However, our one-loop calculation forces the flow of the coupling constant along the parameter space to stay of order ǫ, and one cannot sustain the validity of the preceding description. Then it is not possible to predict the renormalization ofλ 3 . It will be either renormalized toλ * 3 or other fixed points will appear along the flow (corresponding to strong coupling and not given by the one-loop ǫ-expansion). Therefore, the fixed point given by Eq. (33), although it is an attractor, is unreachable from our initial conditions (λ o 3 < 0). For that reason the conclusions of Ref. [24] were incomplete. On the other hand, above the upper critical dimension the system evolves towards the trivial fixed pointλ * 3 = 0 giving a diffusive behavior with z = 2 provided thatλ o 3 is not too much negative (see Fig. 3). This behavior of the fixed pointλ * 3 as a function of ǫ corresponds to a transcritical bifurcation. Now, by introducing the alternative regularization (K = 0), we will see the effect of the symmetry breaking. First of all, we insist that the stable fixed points are the same as for K = 0, due to the fact thatλ 2 renormalizes to zero. Moreover, as can be seen in Fig. 4, where we have plotted the flow lines of Eq. (31) obtained by numerical integration, the basin of attraction of the nontrivial fixed point is delimited by the parabolaλ which is also a particular solution of the flow equations, no matter the value of d. This parabola is inside the region defined by Eq. (35), and this fact implies that the new regularization makes possible to reach the attractor for d < 4 starting in the region of physical meaning. Using Eqs. (34) and (36) together with (9) one gets that the condition to converge towards the nontrivial fixed point is K 2 > 7 2 . Then, the parameter that breaks the symmetry in the regularization of the step function, which in principle was arbitrary, determines the behavior of the system in the hydrodynamic limit. In general for any d < 4 the results are qualitatively the same. Dots correspond to the numerical integration of Eqs. (31c) and (31d), and the thin line is Eq. (36), that clearly delimits the basin of attraction of the nontrivial fixed point, as it is seen in the plot. Below the continuous thick line the values of the parameters correspond to our physical situation, Eq. (35). Squares correspond to the fixed points. Observe that forλ2 = 0 we obtain the same results as in Fig. 3.
For d > 4 the flow is more complex, because of the six fixed points, but the result is that convergence towards the Gaussian one also happens for our initial conditions, as Fig. 5 shows. The linear stability analysis of the fixed points gives the same results as the numerical integration shown in the figure. However this linear analysis fails for d = 4, where all the fixed points collapse towards the Gaussian one. It is by means of the numerical integration that we verify that it is an attractor for the region above the parabola given by Eq. (35), but for the region below it is a repeller. This strange behavior appears because in d = 4 we are at the bifurcation point. In Ref. [24] it was shown for the BTW model and K = 0 that λ 2n = 0, whereas for the Zhang's model, although the even coupling constants do not vanish, it was argued that their flow equations became decoupled from the odd ones in the limit β → ∞. This fact enabled to establish the same universality class for both models, and to deal with only odd terms in Eq. (13). Then, an expansion in the number of coupling constants was performed, whose extrapolation compares well with the results of the simulations [2,5]. Note that in the simulations one computes the dynamical exponent relating the characteristic length and lifetime of the avalanches, whereas within the DRG framework one computes the dynamical exponent from the fluctuations of the order parameter [35]. The agreement between these calculations confirms the basic scaling hypothesis that in both cases length and time are related by means of the same exponent. However, the problem of this calculation was that the nontrivial fixed point was unreachable for the original equation.
In our approach, due to the symmetry breaking, we have to consider also the effect of even coupling constants. In the present work we have dealt with a restricted problem with only the lower-order even and odd coupling constants, λ 2 and λ 3 , showing that λ 2 renormalizes to zero, supporting the calculation of Ref. [24]. Then the stable fixed points are not modified by the presence of an even coupling constant in the model, but due to the symmetry breaking that we have introduce, the nontrivial one is an attractor in the parameter space when the parameters corresponding to the real model are taken into account. This behavior should be the same for any even coupling constant; actually, preliminary calculations including λ 4 and λ 5 in Eq. (19) make us to suspect that all even coupling constants renormalize to zero. This fact means that in the hydrodynamic limit the solution of both models has to be symmetric under parity transformations of the order parameter; then, for the BTW model the DRG restores the broken symmetry, whereas for the Zhang's model we conclude that its asymmetric nature is irrelevant in the behavior at large distances and long times. Therefore this validates the extrapolation performed in [24] since now we have showed that the symmetry breaking makes the stable fixed points reachable, when starting in the region of physical interest in the space of parameters. Let us finally mention that in a recent work, Ghaffari and Jensen [36] perform a different extrapolation of the same results which show a better agreement with large-scale simulations and with real-space renormalization calculations for the dynamical exponent [23]. It is noticeable than the same technique has been applied to the study of the effect of dissipation in a uniformly driven BTW model [37].

IV. CONCLUSIONS
We have studied analytically two models which show self-organized criticality. The difference between them is that the second one (BTW) is symmetric under a parity transformation, whereas the first (Zhang's) model is not. From the microscopic rules one writes a effective long wave-length equation involving the threshold condition, which enters into the equation through a step function, making the equation unapproachable under this form. We have introduced a new regularization of the step function that breaks the symmetry of the BTW model. After a power series expansion, the equation is suitable for the application of the dynamic renormalization group, although it contains an infinite number of relevant coupling constants. In consequence one has to truncate at some point the expansion in the coupling constants. The results only have sense if it is possible to extrapolate the values of the exponents up to an infinite number of coupling constants. We obtain the fixed points of the transformation in parameter space and study carefully their stability and basins of attraction. Then we find that with this regularization it is possible to reach the nontrivial fixed point for d < 4, that was unreachable in a previous work, where symmetry was not broken. This means that in the hydrodynamic limit the models display "scale invariance". Moreover, in this limit we obtain a symmetric behavior under parity transformations for both models, and therefore the recovery of the broken symmetry for the BTW model and the irrelevance of this symmetry for Zhang's one. Although we have dealt with a simplified version of the problem, we expect this behavior to be the same for the complete problem, in the sense that all even coupling constants renormalize to zero, validating the calculation of Ref. [24]. The application of this technique should also be useful for other kinds of problems in which one deals with thresholds or with an infinite number of nonlinear terms, for instance interface dynamics. Moreover, the performed DRG calculation is interesting because it provides an example showing how much important is to know not only the stable fixed points of a DRG transformation but also their basins of attraction. It is remarkable the fact that a simple symmetry breaking can solve the problem of the inaccessibility of the attractors in parameter space.
The fact that the parameter that breaks the symmetry determines the behavior in the hydrodynamic limit is difficult to understand and we believe that it is an artificiality introduced in the calculation by the truncation in the coupling-constants expansion. We expect that higher orders in this expansion will give a behavior independent on the K value.
The DRG procedure eliminates the modes of the outer shell, within the same philosophy than the Kadanoff transformation does in real space. Then one is only interested in E < ( k, ω) and < E < ( k, ω)E < ( k ′ , ω ′ ) >, whose equations turn out to be equivalent to (21) and (25), but with additional terms due to the coupling between the two shells, via the convolution products. The fact that E > ( k, ω) appears in the inner-shell equations allows a perturbative expansion in the form E > ( k, ω) = G > 0 ( k, ω)η > ( k, ω) + ... (using the equivalent of Eq. (21) but in the outer shell). Then the noise in the outer shell enters into the equation for E < ( k, ω). A similar perturbative expansion is done for < E < ( k, ω)E < ( k ′ , ω ′ ) >. By averaging over η > ( k, ω), the contribution of the fast modes is eliminated from the inner shell. This is done up to one-loop order in the perturbative expansion, that is, the lowest order in the intensity of the noise Γ, which implies that it has to be small enough. This tedious calculation becomes more appealing using the diagrams of Fig. 1, instead of the corresponding equations. After this process, the relevant diagrams that survive the averaging are shown in Fig. 2. As an example let us consider one of them, shown in Fig. 6 and denoted by V ( k, ω): where the symbol < ... > > stands for an average over the outer shell. Using the noise correlation given by Eq. (26) we can integrate over Ω, Ω ′′ , and q ′′ , and then we have As the bare propagator is a known function, given by Eq. (23), we are also able to perform the integral over q, that is This integral is a function of k, ω, q ′ , and Ω ′ . However, we are going to evaluate it in the hydrodynamic limit, by taking k, q ′ → 0 and ω, Ω ′ → 0. Then It is easy to check that this result is also valid for d = 4. We have used the explicit form of the bare propagator (23) and also that d d q = S d q d−1 dq, with S d the complete solid angle in d dimensions, that is, the area of a unit (d)-sphere. Then, by making use of Eq. (28), we obtain and substituting into Eq. (A3), It is clear from Fig. 2 that after the first step of the DRG we have the same as at the beginning (Fig. 1), but defined only in the inner shell, plus a lot of diagrams of the same type as the one in Fig.  6. These diagrams, which contain integrals over the outer shell, renormalize the other diagrams that are only defined in the inner shell. For instance, if we consider the following diagram and compare it with with Eq. (A7) we observe only an additional term I d (l)λ 2 2 Γ/D 4 that comes from the outer shell integration. So, the diagram shown in Fig. 6 contributes to the renormalization of (A8), that is, renormalizes the coupling constant λ 2 . As Fig 2(a) states, the diagram in Fig. 6 appears eight times in the perturbative expansion, and the new λ 2 , after the first step of the transformation will be modified by In the same way one can perform the outer-shell integrals of the rest of diagrams in Fig. 2(a). A general result for its contribution to the renormalization of any coupling constant λ n or to the diffusion coefficient D (which will be referred here also as −λ 1 ) is given by where v is the number of vertices each diagram has, three for our example (since the dashed line in Fig. 6 forms a triangle), b(m) is the number of branches of the m-th vertex, two for each one in the example, and B is the number of branches of the diagram that is renormalized (two in the example), and fulfills B = v m=1 b(m) − v − 1. Note that the magnitude in Eq. (A10) is dimensionless. Using this equation and Fig. 2(a) the derivation of Eqs. (27b-27d) is then straightforward.
For the renormalization of the intensity of the noise Γ we have only one diagram, the dashed one in Fig. 2(b). It is immediate to see that the integral over the outer shell, no matter is value, is multiplied by a factor k 2 k ′ 2 . Then where A is simply a numeric factor, and hence, in the hydrodynamic limit and up to one-loop order, the intensity of the noise is not renormalized, after the first step of the DRG, as Eq. (27a) states. * Electronic address: alvaro@ulyses.ffn.ub.es * *