Effects of random environment on a self-organized critical system: Renormalization group analysis of a continuous model

We study effects of random fluid motion on a system in a self-organized critical state. The latter is described by the continuous stochastic model, proposed by Hwa and Kardar [{\it Phys. Rev. Lett.} {\bf 62}: 1813 (1989)]. The advecting velocity field is Gaussian, not correlated in time, with the pair correlation function of the form $\propto \delta(t-t') / k_{\bot}^{d-1+\xi}$, where $k_{\bot}=|{\bf k}_{\bot}|$ and ${\bf k}_{\bot}$ is the component of the wave vector, perpendicular to a certain preferred direction -- the $d$-dimensional generalization of the ensemble introduced by Avellaneda and Majda [{\it Commun. Math. Phys.} {\bf 131}: 381 (1990)]. Using the field theoretic renormalization group we show that, depending on the relation between the exponent $\xi$ and the spatial dimension $d$, the system reveals different types of large-scale, long-time scaling behaviour, associated with the three possible fixed points of the renormalization group equations. They correspond to ordinary diffusion, to passively advected scalar field (the nonlinearity of the Hwa--Kardar model is irrelevant) and to the"pure"Hwa--Kardar model (the advection is irrelevant). For the special choice $\xi=2(4-d)/3$ both the nonlinearity and the advection are important. The corresponding critical exponents are found exactly for all these cases.


Introduction
Numerous physical systems reveal self-similar (scaling) behaviour over the extended ranges of spatial or temporal scales with highly universal exponents. The most prominent example is provided by the critical behaviour of various equilibrium systems near their second-order phase transition points; see, e.g., [1,2] and references therein. An essentially different example is given by the phenomenon of self-organized criticality (SOC), typical of open nonequilibrium systems with dissipative transport; see [3] and references therein. In contrast to equilibrium systems, they do not include a tuning parameter (like the temperature for the second-order transitions) and evolve to the critical state as a result of their intrinsic dynamics. Self-organized critical states are believed to be ubiquitous in the Nature and are observed also in biological, ecological and social systems [4].
Investigation of the effects of various kinds of deterministic or chaotic flows on the behaviour of critical systems (like liquid crystals or binary mixtures near their consolution points) has shown that the flow can destroy the usual critical behaviour: it can change to the mean-field behaviour or, under some conditions, to a more complex behaviour described by new non-equilibrium universality classes [5]- [10].
In this paper, we study effects of turbulent advection on a system in a self-organized critical state by means of the field theoretic renormalization group (RG). As a rule, SOC is studied on the base of discrete models with discrete time steps (like cellular automata). We employ the continuous model proposed in [11], described by the stochastic equation for a smoothed (coarse-grained) height field h(x) = h(t, x) of a certain profile (sand pile or a landscape). Some modifications of the model [11] were proposed in [12,13] and the effects of a flow on the erosion of landscapes were discussed in [13], but the velocity field as an independent dynamical variable was not introduced there.
In this paper, we explicitly introduce coupling with the velocity field. It is modelled by the strongly anisotropic Gaussian ensemble, with vanishing correlation time and prescribed power-like pair correlation function, the d-dimensional generalization of the ensemble introduced and studied in [14] and, at the same time, the anisotropic modification of the popular Kraichnan's rapid-change model; see [15] for the review and the references. It is natural here to consider anisotropic flow, because the original model [11] already involves intrinsic anisotropy, related to the overall tilt of the landscape or the sand pile surface.

The model. Field theoretic formulation and renormalization
Let n be a unit constant vector that determines a certain preferred direction. Then any vector can be decomposed into the components perpendicular and parallel to n, for example, x = x ⊥ + nx with x ⊥ · n = 0. Let ∂ i = ∂/∂x i with i = 1 . . . d be the derivative in the full d-dimensional x space, ∂ ⊥ = ∂/∂x ⊥i with i = 1 . . . d − 1 is the derivative in the subspace orthogonal to n, and ∂ = n · ∂.
The stochastic differential equation for the height field h(x) = h(t, x) is taken in the form [11]: where ∂ t = ∂/∂t, ν ⊥0 and ν 0 are the bare diffusivity coefficients and f (x) is a Gaussian random noise with zero mean and the pair correlation function The velocity field will be taken in the form v = nv(t, x ⊥ ), where v(t, x ⊥ ) is a scalar function independent of x . Then the incompressibility condition is automatically satisfied: For v(t, x ⊥ ) we assume a Gaussian distribution with zero mean and the pair correlation function with k ⊥ = |k ⊥ | and the scalar coefficient functions of the form Here B 0 > 0 is a constant amplitude factor and ξ is an arbitrary exponent, which (along with the conventional ε = 4 − d) will play the role of a formal RG expansion parameter. The infrared (IR) regularization in (2.4) is provided by the cutoff k ⊥ > m.
[Its precise form is unimportant; sharp cutoff is merely the most convenient choice from the calculational viewpoints.] According to the general statement (see, e.g, the monographs [1,2] and the references therein), our stochastic problem is equivalent to the field theoretic model of the extended set of fields Φ = {h ′ , h, v} with action functional (2.6) All the required integrations over x = {t, x} and summations over the vector indices are implied. The last term in (2.6) corresponds to the Gaussian averaging over v with correlator (2.4) and has the form where This allows one to apply the field theoretic renormalization theory and renormalization group to our stochastic problem. The model (2.6) corresponds to a standard Feynman diagrammatic technique with three bare propagators: vv 0 , given by (2.4), (2.5), and the propagators of the scalar fields (in the frequency-momentum representation): The propagator h ′ h ′ 0 vanishes identically for any field theory of the type (2.6). The model also involves two types of vertices corresponding to the interaction terms −h ′ ∂ h 2 /2 and −h ′ (v∂ )h. The corresponding coupling constants g 0 and w 0 are defined by the relations where is ℓ has the order of the smallest length scale in our problem. The analysis of canonical dimensions (see, e.g., [1,2]) along with symmetry considerations shows that our model is multiplicatively renormalizable with the only independent renormalization constant Z ν : all the renormalization constants are equal to 1. In the leading (one-loop) approximation Z ν has the form Z ν = 1 − ag/ε − bw/ξ (the use of the minimal subtraction scheme is implied). Here g and w are renormalized analogues of the bare parameters in (2.10), which play the role of the coupling constants (dimensionless expansion parameters) in renormalized perturbation theory, and a, b are numerical coefficients. Their precise values are unimportant (they can simply be absorbed by the redefinition of g and w), but it is important for the following that they are positive: a, b > 0.
To proceed, we have to find the RG functions: the anomalous dimensions and the β functions, which appear as coefficients in the differential RG equations; detailed discussion can be found in [1,2]. The anomalous dimension γ F for a given renormalization constant Z F is defined as where D µ is the differential operation µ∂ µ for fixed e 0 , the latter being the full set of bare parameters e 0 = {ν 0 , ν ⊥0 , w 0 , g 0 }. The β functions for the two dimensionless couplings g and w are Taking into account the relations between the bare parameters (2.10), absence of renormalization of D 0 , B 0 and the explicit form of the constant Z ν , we obtain: where the ellypsis stands for the higher-order corrections in g and w.

Fixed points and scaling regimes
It is well known that possible scaling regimes of a renormalizable model are associated with IR attractive fixed points of the corresponding RG equations. In the present case, the coordinates g * , w * of the fixed points are found from the equations β g (g * , w * ) = 0, β w (g * , w * ) = 0, (3.1) with the β functions given in (2.13). The type of a fixed point is determined by the matrix Ω = {Ω ij = ∂β i /∂g j }| g * ,w * , where β i is the full set of the β functions and g j = {g, w} is the full set of coupling constants. For IR attractive points the matrix Ω is positive, that is, the real parts of all its eigenvalues Ω i are positive. The β functions (2.13) involve the same anomalous dimension γ ν and thus satisfy the exact relation wβ g − 3gβ w /2 = gw(−ε + 3ξ/2). In turn, this means that for general ε and ξ the equations (3.1) can be satisfied only if at least one of the coordinates g * , w * vanishes identically. There are three such fixed points: (1) Gaussian (free) fixed point with g * = w * = 0 and γ * ν ≡ γ ν (g * = w * = 0) = 0. The eigenvalues of the matrix Ω are equal to its diagonal elements: Ω g = ∂ g β g = −ε, Ω w = ∂ w β w = −ξ. Thus, the point is attractive for ε, ξ < 0.
From the physics viewpoints, in this critical regime the nonlinear term in the equation (2.1) does not affect the leading terms of the IR asymptotic behaviour (it is IR irrelevant in the sense of Wilson).
(3) The point with w * = 0, g * = ε/a and γ * ν = 2ε/3. The latter expression is exact (no corrections of order ε 2 and higher) owing to the exact relation between β g and γ ν in (2.13). Furthermore, the derivative ∂ g β w = w∂ g γ ν vanishes for w * = 0, the matrix Ω is block triangular and its eigenvalues again coincide with the diagonal elements: Ω w = −ξ + 2ε/3 (exact result) and Ω w = ε (with corrections of order ε 2 and higher). We conclude that this point is attractive for ξ < 2ε/3, ε > 0. In this critical regime, the advection is irrelevant and the IR behaviour of the model coincides with that of the original Hwa-Kardar model.
It is also worth noting, that for all fixed points (1)-(3) in the regions of their IR attraction, their coordinates lie in the physical region g * ≥ 0, w * ≥ 0.
For the special choice ξ = 2ε/3 the functions (2.13) become proportional, and equations (3.1) only impose the restriction ag * + bw * = ξ; that is, we have a line of fixed points in the g-w plane. However, the exact value γ * ν = ξ = 2ε/3 is the same for all points on that line. The direct calculation of the elements of the matrix Ω shows that one of the eigenvalues vanishes (this exact result is a consequence of the degeneracy of the fixed point) while the other, equal to 3ag * /2 + bw * , is positive in the physical region g * , w * > 0, as well as the exponent ξ.
In figure 1 we show the regions of IR attraction for all the fixed points in the ε-ξ plane, that is, the regions where the eigenvalues of the matrix Ω are both positive.

Critical scaling and critical dimensions
Existence of IR attractive fixed points of the RG equations leads to appearance of the scaling behaviour of the Green functions in the IR range. In dynamical models with two momentum scales, like our present model, the critical dimension ∆ F of an IR relevant quantity F (a field or a parameter) is given by the relation Here d F are the canonical dimensions (see, e.g., [2]), γ * F is the fixed-point value of the corresponding anomalous dimension, are the critical dimensions of the momentum in the direction parallel to n and of the frequency, respectively, with the normalization ∆ ⊥ = 1 for the momentum in perpendicular subspace. Detailed derivation and discussion of the relations of the type (4.1), (4.2) can be found, e.g., in ref. [9]. In our model these relations simplify owing to the absence of renormalization of all the fields (so that γ * F = 0) and of the parameter ν ⊥0 (so that γ * ν ⊥ = 0); see sec. 2. What is more, γ * ν is known exactly for all the fixed points. This gives the following exact expressions: for the Gaussian fixed point (1), for the fixed point (2) and the line of the fixed points, and for the point (3). In particular, for the pair correlation function of the main field this gives: where r ⊥ = |x ⊥ |, r = x and F is a certain scaling function of critically dimensionless arguments. It remains to note that for the fixed point (3) our results agree, up to the notation, with those obtained in [11] (one has to identify z = ∆ ω /∆ , ζ = 1/∆ , χ = −∆ h /∆ ).

Conclusion
We studied effects of fluid motion, including turbulent mixing, on a system in the state of self-organized criticality. The latter was modelled by a stochastic differential equation for a height field (2.1), (2.2), while the velocity of the fluid was described by a random Gaussian ensemble (2.4). Both the equation and the velocity ensemble are strongly anisotropic.
The full model can be reformulated as a renormalizable field theoretic model with a single independent renormalization constant. The corresponding RG equations possess several IR attractive fixed points, corresponding to different possible types of critical behaviour. They correspond to ordinary diffusion, passively advected scalar field and to the original SOC model without the mixing. Their regions of IR attraction in the plane of the model parameters d and ξ and the critical dimensions of the basic fields and parameters are found exactly. For the special choice ξ = 2(4 − d)/3 an intermediate regime arises, where both the nonlinearity and the advection are important.
In the following, it would be interesting to study other models of SOC and to employ more realistic models for the turbulent velocity field, e.g., models with finite correlation time and non-Gaussianity. This work is in progress.