Strongly Nonlinear Diffusion in Turbulent Environment: A Problem with Inﬁnitely Many Couplings

: The ﬁeld theoretic renormalization group is applied to the strongly nonlinear stochastic advection-diffusion equation. The turbulent advection is modelled by the Kazantsev–Kraichnan “rapid-change” ensemble. As a requirement of the renormalizability, the model necessarily involves inﬁnite number of coupling constants (“charges”). The one-loop counterterm is calculated explicitly. The corresponding renormalization group equation demonstrates existence of a pair of two-dimensional surfaces of ﬁxed points in the inﬁnite-dimensional parameter space. If the surfaces contain infrared attractive regions, the problem allows for the large-scale, long-time scaling behaviour. For the ﬁrst surface (advection is irrelevant), the critical dimensions of the scalar ﬁeld ∆ θ , the response ﬁeld ∆ θ (cid:48) and the frequency ∆ ω are nonuniversal (through the dependence on the effective couplings) but satisfy certain exact identities. For the second surface (advection is relevant), the dimensions are universal and they are found exactly.

In this paper, we study strongly nonlinear diffusion of a scalar field θ = θ(x) = θ(t, x) (temperature, entropy, density of a pollutant) in a turbulent environment described by the velocity field v = {v i (x)}, i = 1, . . . , d, where d is the arbitrary (for generality) dimension of space.
This problem is not related to any kind of phase transition, but it is expected that various correlation and response functions exhibit scaling behaviour in the infrared (IR) range (large distances and time intervals). In particular, for the structure functions one expects: C 2n (t, r) = [θ(t, x) − θ(0, 0)] 2n r −2n∆ θ F 2n (tr ∆ t ), r = |x|, (1) while for the Green (linear response) function: Here, the brackets · denote averaging over the statistical ensemble, θ = θ (x) is an auxiliary response field (see Section 2), ∆ θ , ∆ θ and ∆ t = −∆ ω are the critical dimensions of the corresponding quantities and F 2n (·), F(·) are certain scaling functions. (In a more traditional notations ∆ θ = −χ, ∆ t = −z.) The aim of the study is to establish scaling relations of the type (1) and (2), on the basis of a certain dynamical model, to calculate the critical dimensions and to investigate their universality (that is, dependence on d and other parameters of the model).
In most approaches, the diffusion equation is taken in the self-sufficient linear approximation; see, e.g., [14] and references therein. In a more general situation, it is natural to include a nonlinearity of an appropriate form (such that the diffusion equation retains its meaning as a conservation law). Then the dimensional analysis and more sophisticated renormalization considerations immediately show that infinite number of nonlinear terms are equally relevant and should be all included into the model; see, e.g., [15]. Thus, the nonlinear model necessarily involves infinitely many coupling constants from the very beginning.
In the quantum field theory, such models for a long time faced negative attitudes for being nonrenormalizable and therefore for having no predictive power. Fortunately, now the situation has changed, although a generally accepted interpretation is still not achieved; for a recent discussion see, e.g., [16] and references therein.
In this paper, we reformulate the stochastic nonlinear advection-diffusion equation as a certain field theoretic model with infinitely many couplings. We show that this full-scale model can be considered multiplicatively renormalizable in a wider, than usual, sense of the term. This allows one to derive the renormalization group equations, to investigate possible asymptotic regimes of the problem and to obtain some exact relations and expressions for the critical dimensions in scaling laws.
Experience with nearly-equilibrium nearly-critical systems suggests that they can be drastically affected by the motion of the constituting or surrounding medium. Indeed, the critical scaling behaviour can be destroyed in favour of the mean-field behaviour or a completely new non-equilibrium universality classes [17][18][19][20]. Thus, it is highly desirable to study the effects of the medium motion (which can hardly be excluded in real experimental settings) on the non-equilibrium critical behaviour.
In this paper, the random velocity field v(x) will be described by the Kazantsev-Kraichnan "rapid-change ensemble", a Gaussian distribution δ-correlated in time and power-like in momentum. In the end of the previous century this model attracted enormous attention among turbulent community because of the deep insight it offers into the origin of intermittency and anomalous scaling in fluid turbulence [21]. It is therefore interesting to employ that ensemble in our problem. In our previous study [15], the velocity was governed by a stochastic Navier-Stokes (NS) equation.
The plan of the paper is the following. In Section 2 we give detailed description of the model, including its field theoretic formulation. In Section 3 we show, on the basis of dimensional and symmetry considerations, that the model is multiplicatively renormalizable in the infinite-dimensional space of parameters (coupling constants). The corresponding one-loop counterterm is derived in Section 4 in an explicit closed form, which allows one to write down the full infinite set of the renormalization constants in the leading one-loop approximation. For a multiplicatively renormalizable model, the differential renormalization group (RG) equations can be derived in a standard manner; Section 5. The corresponding RG functions (anomalous dimensions and β functions) are presented in the one-loop order. In Section 6 we show that the renormalization group equations have two attractors in the form of two-dimensional surfaces of fixed points in the infinite-dimensional parameter space. If those surfaces contain infrared attractive regions, the problem allows for the large-scale, long-time scaling behaviour. For the first surface (advection is irrelevant), the critical dimensions of the fields and the frequency are nonuniversal (they depend on the specific choice of a fixed point on the attractor surface) but satisfy certain exact identities. For the more interesting second surface (advection is relevant), the dimensions are shown to be universal (independent of the choice of a fixed point on the sheet) and they are found exactly. Section 7 is reserved for discussion and conclusion.

Formulation of the Model
The general advection-diffusion equation for the scalar field θ = θ(x) = θ(t, x) has the form: is the Lagrangian (Galilean covariant) derivative, J i (θ) is the current density, ∂ t = ∂/∂t and Summation over repeated indices is always implied. Within the field theoretic RG the problem of passively advected scalar field was studied in [22] in the simplest linear approximation J i ∝ ∂ i θ and without the noise f . The authors justified the phenomenological "four-thirds law" by Richardson [14] and calculated the effective turbulent Prandtl number. Subsequent works dealt with the Batchelor constant in the spectrum of the scalar field, chemically decaying admixture, effects of anisotropy, and so on; see, e.g., Sections 3.1-3.3 and 3.5 in [23] and references therein.
The random noise f = f (x) models interaction with the large-scale modes and boundaries that supply kinetic energy to the system. In order for the standard field theoretic renormalization to be applicable, the noise is taken to be Gaussian with zero mean, δ-correlated in time (this is required by Galilean symmetry), with the pair correlation function: Here, the cut off at k = m provides IR regularization, 1/m has the meaning of the outer turbulence scale (the largest scale in the problem). It also ensures vanishing of the correlation function (5) at k = 0, so that the noise f does not violate the conservation law expressed by the Equation (3). Furthermore, 0 < y < 2 is an arbitrary exponent with the most realistic value y → 2: then, with an appropriate choice of the amplitude, the function (5) and (6) can be viewed as a power-like representation of the function δ(k). All these issues are discussed in more detail (also in connection with the stochastic NS equation) in the paper [15], in the review paper [24] (Sections 2.1 and 3.7) and in the monographs [23] (Sections 1.1 and 2.8) and [25] (Section 6.3).
Since we are interested in the IR asymptotic behaviour, only the first-order term should be retained in the expansion J i = ∂ i V(θ) + O(∂ 3 ) of the current density J i (θ) in powers of the gradients ∂. On the contrary, in the expansion of the coefficient function V(θ) in powers of θ the dimensional analysis shows that all the terms should be retained (see Section 3). Thus, the higher powers in θ cannot be neglected in comparison with the lower-order ones in V(θ).
From the more refined renormalization viewpoints, this means that inclusion of any term θ n with n > 1 into the function (7) gives rise to the infinite set of counterterms θ p with all possible p ≥ 1, needed to eliminate the ultraviolet (UV) divergences, as we will see below in Section 3. In order to ensure multiplicative renormalizability, they should be included into the function V(θ) from the very beginning.
As a result, our model involves infinitely many coupling constants "hidden" in the amplitudes λ n0 with n > 1 (more precisely, see Section 3). In this respect, our model appears sharply different from conventional models of equilibrium critical behaviour such as φ 4 , where all the terms involving extra gradients or extra fields are IR irrelevant and should be discarded; see, e.g., the monograph [25] and references therein.
An important exception is provided by the linear case V(θ) ∝ θ. The linearity is preserved by the renormalization procedure and the higher-order terms θ n with n > 1 do not arise as counterterms (although they are allowed by dimensions); see Section 3. Thus, the linear approximation to V(θ) appears "closed with respect to the renormalization" and can therefore be studied as a self-contained model; see, e.g., [22,23].
A similar situation with infinite number of couplings was encountered earlier in a properly extended version [26] of Pavlik's model [27] of kinetic surface roughening. Another example is provided by stochastic models of landscape erosion [28,29]: there, correct RG analysis also requires to extend the models by adding infinite number of interaction terms and corresponding coupling constants [30][31][32][33].
In the Kazantsev-Kraichnan ensemble the velocity field is taken Gaussian, with zero mean and the given pair correlation function (8) where 0 < ξ < 2 is an arbitrary exponent. The specific choice ξ = 4/3 corresponds to the assumptions of the Kolmogorov-Obukhov phenomenological theory of turbulence: independence of the velocity correlations from the inner and outer turbulence scales in the inertial-convective interval; see, e.g., [14]. The transverse projector P ij (k) = δ ij − k i k j /k 2 ensures fulfilment of the incompressibility condition ∂ i v i (x) = 0. Again, the cut off at k = m provides IR regularization. There is no need to distinguish the parameters m here and in (5): the precise form of IR regularization is unessential from the renormalization point of view; see, e.g., [25].
Thus, we arrive at the equation with the Laplace operator ∂ 2 = ∂ i ∂ i , the covariant derivative ∇ t from (4) and the infinite series from (7). The model involves additive random noise f (x) with the statistics described by the pair correlation function (5) and (6) and multiplicative "noise" v(x) with the statistics specified by (8). This completes formulation of the model. According to the general De Dominicis-Janssen theorem (see, e.g., Chapter 5 in [25] and references therein), the original stochastic problem is equivalent to the field theoretic model of an extended set of fields Φ = {θ , θ, v} with the action functional S (with D f from (5) and (6)) is the action for the stochastic problem (9) at fixed v and provides averaging over the Gaussian velocity statistics defined by Equation (8); D −1 ij is the kernel of the integral operation inverse to D ij = v i v j v . In this formulation, the auxiliary Martin-Siggia-Rose response field θ (x) replaces the additive noise f (x).
Here and below, all the needed integrations over the arguments x = {t, x} and summations over repeated vector indices are implied for all terms in the expressions such as those entering (10), for example, The formulation (10) and (11) means that various correlation and response functions Φ . . . Φ f ,v of the original stochastic problem can be represented as functional averages over the full set of fields Φ = {θ , θ, v} with weight exp S(Φ). This allows one to apply to the problem the field theoretic machinery, including Feynman diagrammatic techniques, functional equations of Schwinger-Dyson type, Ward identities for Galilean symmetry, renormalization theory and renormalization group.

Canonical Dimensions and UV Renormalization
The UV divergences are eliminated from the correlation functions by means of the standard renormalization procedure. However, for our dynamic model with interactions containing gradients, this procedure has a number of important peculiarities.
It is well known that analysis of UV divergences is based on analysis of canonical dimensions, see, e.g., [25] (Sections 1.15 and 1.16). In contrast to conventional static models, dynamic ones have two independent scales: a time scale [T] and a spatial scale [L]; see [25] (Sections 1.17 and 5.14). Thus, the canonical dimension of any quantity F (a field or a parameter) is determined by two numbers, the frequency dimension d ω F and the momentum dimension d k F : The dimensions are found from obvious normalization conditions and from the requirement that all terms in the action functional be dimensionless with respect to both the canonical dimensions separately. The total canonical dimension is defined as The coefficient 2 follows from the relation ∂ t ∝ ∂ 2 in the free theory.) In the renormalization routine, d F plays the same part as does the conventional (momentum) dimension in static models; see [25] (Section 5.14).
In order to derive unambiguous values for the dimensions, one has to get rid of "superfluous" parameters, if any. In our model, the amplitude D 0 can be removed from the correlation function (6) by an appropriate rescaling of the fields and the parameters. In other words, one can set D 0 = 1 without loss of generality, which is always implied in the following. Of course, such a rescaling does not affect the exponents in the expressions such as (1) and (2).
Canonical dimensions of all the fields and the parameters are given in Table 1. It also involves renormalized parameters (without subscript "0") and the renormalization mass µ, an additional parameter of the renormalized theory; they all will appear later on. Table 1. Canonical dimensions of the fields and the parameters in the model (10) and (11).
We also denote λ n0 = ν 0 for n = 1, while for n > 1 we introduce the new variables: They are chosen so that d ω g n0 = 0, while their renormalized counterparts g n are completely dimensionless. We also define B 0 = w 0 ν 0 for the amplitude in (8), so that d k w 0 = ξ and d ω w 0 = 0. Its renormalized counterpart w 0 ∼ wµ ξ will be completely dimensionless. The parameters w 0 and g n0 with n > 1 play the part of coupling constants ("charges") in the original unrenormalized model, while their dimensionless analogs w, g n serve as couplings in the renormalized perturbation theory. Since we are interested in the expansion in the number of loops, the consistency relations g n ∼ g (n−1) 2 should be implied. From Table 1, it follows that all the interactions θ ∂ 2 θ n simultaneously become logarithmic (the bare charges g n0 become dimensionless) at y = 0, while the advection term θ (v i ∂ i )θ becomes logarithmic (the bare charge w 0 becomes dimensionless) at ξ = 0. According to the general ideology of renormalization, the exponents y and ξ that "measure" deviation from the logarithmicity, should be treated as formal small parameters of the same order, y ∼ ξ: only then we obtain the full-fledged model in which all the interactions are simultaneously relevant. It should be stressed that, in contrast to the standard φ 4 model of equilibrium critical behaviour, where the deviation from logarithmicity is measured by the deviation ε = 4 − d of the spatial dimension d from its logarithmic value d = 4 [25], in our case the exponents y and ξ are completely unrelated to the dimension d which, therefore, remains an arbitrary parameter. The UV divergences manifest themselves as singularities at y → 0 and ξ → 0 in the correlation functions; in the one-loop approximation they have the form of simple poles.
Superficial UV divergences, whose removal requires introducing the counterterms, can appear in the 1-irreducible correlation functions (Green's functions in the quantum field terminology) for which the formal index of divergence in the logarithmic theory (y = ξ = 0) is a non-negative integer. Here N Φ are the numbers of the fields Φ = {θ , θ, v} entering the Green's function and d Φ are their total canonical dimensions; see, e.g., [23] (Section 1.4) and [25] (Section 5.15).
In addition, when analyzing divergences in the model (10) and (11) one should take into account the following additional considerations.
(i) For any dynamic model of the type (10), all the 1-irreducible functions without the response fields contain closed circuits of retarded propagators and vanish; see, e.g., [25] (Section 5.4). Thus, it is sufficient to consider the functions with N θ > 0.
(ii) In the vertices θ ∂ 2 θ n , the Laplace operator can be moved, using integration by parts, onto the field θ . Thus, in any 1-irreducible diagram each external field θ , attached to one of such vertices, "releases" the square of the corresponding external momentum, and the real index of divergence reduces by the corresponding number of units: δ = δ − 2N θ . Furthermore, the field θ enters the counterterms only in the form of a spatial gradient.
(iii) Owing to the incompressibility condition ∂ i v i = 0, the derivative at the vertex θ (v i ∂ i )θ can also be moved, using integration by parts, onto the field θ . Thus, in any 1-irreducible diagram each external field θ , attached to such vertex, releases an external momentum, and the real index of divergence reduces: δ = δ − N θ . Again, this also means that θ appears in the counterterms only as a gradient.
From (ii) and (iii), it follows that the real index of divergence for a general 1-irreducible diagram satisfies the inequalities δ − 2N θ ≤ δ ≤ δ − N θ .
(iv) The counterterm θ ∂ t θ, allowed by the formal index δ, is in fact forbidden by the items (ii) and (iii): it does not contain a spatial gradient. On the other hand, the Galilean symmetry of the model requires, in particular, that the covariant derivative θ ∇ t θ enter the counterterms as a single unit. This forbids the counterterm θ (v i ∂ i )θ, allowed by the formal index δ and the items (ii) and (iii).
Taking into account all these considerations, one can conclude that in our model superficial UV divergences are present only in the 1-irreducible diagrams of the Green's functions θ θ . . . θ with a single response field θ and n ≥ 1 basic fields θ. For all these functions δ = 2, δ = 0, and the corresponding counterterms necessarily reduce to the forms θ ∂ 2 θ n . All such terms are already present in the action functional (10), so that inclusion of the counterterms can be reproduced by multiplicative renormalization of the model parameters (no renormalization of the fields Φ and the IR scale m is required).
The linear case V(θ) ∝ θ needs some special reservations. Here, the direct analysis of Feynman diagrams shows that for any 1-irreducible Green's function N θ ≥ N θ (to be precise, the difference (N θ − N θ ) is a non-negative even integer; cf., e.g., [34]). Thus, the only superficially divergent 1-irreducible Green's function is θ θ with δ = 2, δ = 0 and the single counterterm θ ∂ 2 θ. As a consequence, the linear model is "closed with respect to the renormalization" and can therefore be studied on its own right; see, e.g., [22,23].
In the general case the renormalized action functional has the form with the renormalized analog of the series (7) Here S vR is the functional S v from (11) expressed in renormalized variables: and the renormalization mass µ is an additional parameter of the renormalized theory. The renormalization constants Z n = Z n (w, g n , d, y, ξ) are chosen to absorb all the UV divergences. In practical calculations, we will employ the minimal subtraction (MS) renormalization scheme, where these constants have the forms Z n = 1+ only pure poles in y, ξ and, in general, their linear combinations; see Section 4. The renormalized action functional (13) is obtained from its unrenormalized analog (10) by the substitution λ n0 = λ n Z n (n ≥ 1), ν 0 = νZ ν , w 0 = wµ ξ Z w , g n0 = g n µ (n−1)y/2 Z g n (n ≥ 2), (16) where the definitions (12) are taken into account. The constants Z n are calculated directly from the diagrams (see Section 4) and the others are found from the relations The first two relations follow from the definitions (12) and (16), the third one is a consequence of the absence of renormalization of S v = S vR , see Equation (15).

Calculation of the One-Loop Counterterm
In this section, we explicitly construct in a closed form the one-loop counterterm in the renormalized action (13), which allows for calculation of all renormalization constants Z n in (14) at the one-loop level. To this end, we adopt the functional techniques developed earlier for various models with infinite number of couplings [15,26,30,31].
Consider the generating functional of the 1-irreducible Green's functions in renormalized theory. Its expansion in the number of loops p has the form: The loopless (tree-like) term is just the (renormalized) action (13), while the one-loop contribution is given by the expression (see, e.g., [25], Section 2.9): where W is a linear operation with the kernel and W 0 is its analog for the free (bilinear in the fields) part of the action. Both W and W 0 are 3 × 3 matrices in the full set of fields Φ = {θ , θ, v}; hence the trace in (19). Within our accuracy we can put Z n = 1 and v = 0 in the one-loop contribution (19) and (20). In the loopless contribution, the constants Z n should be taken in the first nontrivial order in the couplings w and g n , with the additional convention that g n ∼ g (n−1) 2 . In the following, we interpret the quantities such as (7) and (14) as functions of a single variable θ(x), and V (θ), V (θ), etc., as the corresponding derivatives with respect to this variable. Then, the elements of the matrix W can be symbolically represented as (for brevity, we omit the vector indices): with D f = k 2−d−y from (6) and In order to find the constants Z n , we do not need the full expression (19), rather we need its UV divergent part that was previously established to have the form (Section 3): with some function R(θ) analogous to V(θ). This means that (19) is needed only in the first order in its elements W (θθ) , W (θv) and W (vθ) , linear in θ ; see Equation (21). In order to extract this part, we use the well-known formula δ(Tr ln K) = Tr(K −1 δK) for any variation δK. By varying the matrix W with respect to θ we obtain, with the needed accuracy, Tr ln(W/W 0 ) −I 1 + 2I 2 , where and In these expressions D (ΦΦ) are the corresponding elements of the matrix W −1 for zero values of θ and v. By its very meaning, D (θθ) is the propagator θθ 0 with Z n = 1 and with ν∂ 2 substituted by ∂ 2 V in the denominator, that is, the propagator in the external "background" field θ. Similarly, D (θ θ) is the propagator θ θ 0 with the same substitutions, while D (vv) ij is just the correlation function (8). For what follows, it is crucial that, once the two external derivatives are moved to the external fields θ(x) and θ (x) in the expressions (22) and (23), the remaining integrals diverge only nearly logarithmically. This allows one to set all the external momenta and frequencies equal to zero in the calculation of the divergent parts of these "diagrams". Then the IR regularization is provided by the cut off m in Equations (5) and (8). In other words, after the derivatives ∂ 2 θ (x) in I 1 and ∂ i θ(x), ∂ j θ (x ) in I 2 are isolated as extra factors, one can neglect the inhomogeneity of the fields θ(x) in the remaining integral factors and treat them as if they were merely constants. In the Fourier representation, this means that no frequencies or momenta flow out through the outer vertices of these diagrams, while a single integration momentum and a single frequency circulate inside the whole one-loop diagram. Then the integrals are easily calculated by going over to the Fourier (momentum-frequency) representation.
For the integral entering (22) one has: In the following we denote where S d is the surface area of the unit sphere in d-dimensional space, and Γ(·) is Euler's Γ function.
For the integral entering (23) one has: This integral is ill-defined and requires a careful interpretation. Consider the Fourier transform: where Θ(·) is the Heaviside step function. What we need in (26) is the value of (27) at t = t . This ambiguity can be resolved on physical grounds. The δ function in (8) is in fact an idealized model for a very narrow function f (t, t ) concentrated near t = t and rapidly decaying for |t − t | → ∞. The key observation is that this function must be symmetric in t, t because it represents a pair correlation function. Thus, for any physical choice of "regularized" function f (t, t ) we will unambiguously obtain Θ(0) = 1/2 in the limit f (t, t ) → δ(t − t ). In terms of stochastic differential equations theory this means that the multiplicative noise v is interpreted in Stratonovich's sense [35,36]. Then (26) gives: Substitution of the obtained expressions (24) and (28) into the diagrams (22) and (23) yields the desired explicit expression for the divergent part of the one-loop contribution (19): where V(θ) is given by (14) with Z n = 1; also note that the expression (30) has the same canonical dimensions as V(θ). Both terms appear in (29) with the positive sign: the overall (−1/2) factor in (19) changes the sign and the coefficient in the expression for Tr ln(W/W 0 ) −I 1 + 2I 2 , while the integration by parts that moves the derivative from the field θ to θ in (23) provides an additional (−1) factor for the term coming from I 2 .
Of course, it is tempting to proceed with the renormalization analysis in terms of the closed functional representations (29) and (30), but the authors find no way to do so (see the remarks in Section 7). Hopefully, this can be accomplished with the aid of the functional renormalization group [33]. Thus, for lack of a better, here we expand the function F(θ) back in the powers of θ and arrive at the series of the type (7) and (14): with dimensionless coefficients r n . This decomposition with V R (θ) from (14) with the substitution Z n = 1 (see above) determines r n as the functions of g n . The first few coefficients have the forms that illustrate their general structure: and so on; we recall the consistency relation g n ∼ g (n−1) 2 . From the requirement that the poles in y and ξ be cancelled in the sum of first two terms Γ (0) (Φ) + Γ (1) (Φ) of the loop expansion (18) for the renormalized functional Γ R (Φ), the explicit expressions for the renormalization constants are derived. We use the MS scheme, where they have the forms "Z n = 1+ only pure poles in y, ξ", so that the factors such as (µ/m) y and (µ/m) ξ should be replaced by unity. Finally, one obtains: with r n from (31) and (32) and a d from (25).

Renormalization Group Equations, Anomalous Dimensions and β Functions
Since the model in (10), (11) and (13) is multiplicatively renormalizable, the differential renormalization group (RG) equation can be derived in a standard fashion; see, e.g., [25] (Section 1.24). It has the form: Here, G(·) is a certain renormalized (and hence UV finite) Green's function of the model (13), expressed in the full set of renormalized variables e = {w, g n , ν, µ, m}; the ellipsis stands for the other variables such as coordinates/momenta and times/frequencies.
Here and below ∂ x = ∂/∂x and D x = x∂ x for any variable x.
The coefficients in the RG differential operator (34), the anomalous dimensions γ and the β functions are defined as: where D µ is the differential operation D µ at fixed bare (unrenormalized) parameters. From the definitions (35) and the relations (16) and (17) all the RG functions can be expressed in terms of the basic anomalous dimensions γ n as follows: β n = g n [−(n − 1)y/2 − γ g n ] = g n [−(n − 1)y/2 − γ n + (n + 1)γ 1 /2].
The operation D µ , when acting on the functions that depend only on the couplings, takes on the form: Within our accuracy, this expression reduces to We note that, with the consistency assumption g n ∼ g (n−1) 2 , all the one-loop contributions in (33) are of the same order, r n /g n ∼ g 2 2 (one should also assume w ∼ g 2 2 ). As a result, one can show that D g r n = (n + 1) r n , D g (r n /g n ) = 2(r n /g n ).
Then substitution of the one-loop expressions (33) finally gives: r n g n with r n from (31), (32) and a d from (25).

Attractors of the RG Equations, IR Scaling and Critical Dimensions
Asymptotic behaviour of the Green's functions of a multiplicatively renormalizable model is governed by the attractors of the RG equations, determined by the requirement that all its β functions vanish; see, e.g., [25] (Section 1.42). In our case, it translates to β w (w * , g * n ) = 0, β n (w * , g * n ) = 0 (n > 1).
In conventional models with finite number of couplings the solution to (40) is normally given by a finite number of fixed points; the IR attractive ones (see below) govern the largedistance, long-time scaling behaviour.
In both cases, consecutive substitution of the solutions of equations β w = 0, β k = 0 with k ≤ n into the remaining equations β k = 0 with k > n allows one to express all the coordinates g * n with n > 3 in terms of the two parameters g * 2 and g * 3 . Indeed, it is clear from the explicit forms (32) of the coefficients r n that the function β n involves one "new" coupling g n+2 , not encountered in the previous functions β k with k < n, and the "old" couplings g k with k < n + 2, which were already expressed in terms of g * 2 and g * 3 during preceding steps.
We conclude that the attractor for the system (40) consists of a pair of two-dimensional surfaces. Clearly, the first sheet with w * = 0 corresponds to the situation when the turbulent advection is irrelevant for the scaling behaviour (that is, irrelevant in the sense of Wilson 1 ); for the second sheet the advection is important.
In dynamical models, the critical dimension of any quantity (a field or a parameter) is given by the expression (see, e.g., Section 2.1 in [23], Section 3.1 in [24] and Sections 5.16 and 6.7 in [25]): (with the standard normalization convention that ∆ k = −∆ x = 1). We recall that γ * denotes the value of the anomalous dimension γ * at some point of the attractor surface determined by Equation (40). From the data in Table 1 and the relations γ θ = γ θ = γ m = 0 (no renormalization of the fields Φ and the IR scale m is needed; see Section 3) one obtains: where in the one-loop approximation for the sheet with w * = 0 while for the sheet with w * = 0 exactly (that is, in all orders of the double expansion in ξ and y). This is a consequence of the second relation in (36), the expressions for β w in (37) and for ∆ ω in (41). Thus, one can see that for the second case the critical dimensions ∆ F and ∆ ω appear universal in the sense that they do not depend on the coordinates g * 2 and g * 3 , that is, do not depend on the specific choice of the point on the attractor sheet. This result is probably the main outcome of our study.
In both cases, these dimensions are subject to exact relation another exact relation 2∆ v = −ξ + ∆ ω follows simply from the definition (8). Also note that for w * = 0 Equation (42) can be viewed as certain exact universal relations between the nonuniversal dimensions ∆ θ , ∆ θ and ∆ ω . In general, it is difficult to explore the character of attractiveness for the aforementioned surfaces.
According to the general rule, a fixed point (in our case, a fragment of a surface) is IR attractive if all the eigenvalues of the matrix (in the following discussion it is convenient to tentatively denote w = g 1 ): Ω nm = ∂β n /∂g m | g * , n, m = 1, 2, . . . at the fixed point (in our case, at some point on the attractor surface) have strictly positive real parts; see, e.g., [25] (Section 1.42).
One can ascertain that there are regions on the attractor surfaces where all the diagonal elements Ω nn are positive. Indeed, the calculation gives: In these cases, the inequalities Ω nn > 0 can be transformed into the inequality of the form: All three coefficient functions A, B, C are positive, bounded from above, and have finite limits when n → ∞ (A and C tend to zero while B tends to 2). Thus, inequality (46) provides us with some regions on the attractor surfaces where all Ω nn are positive.
If the matrix Ω were symmetric, the inequalities Ω nn > 0 would give us an infinite set of necessary conditions for Ω > 0 (i.e., for positive eigenvalues), but it is not so. Thus, we are left with the single necessary condition Tr Ω = ∑ n Ω nn > 0. Although it is only a necessary and not a sufficient condition of IR attractiveness (and rather weak due to the infinite number of terms), it still leaves the possibility of existence of IR regions on the attractor surfaces.
If such regions indeed exist, the Green's functions of the model will exhibit scaling behaviour of the type (1) and (2) with nonuniversal exponents ∆ θ and ∆ t = −∆ ω from (42) and (43) for the surface with w * = 0 and exactly known universal ones from (42) and (44) for the surface with w * = 0. Due to the exact relation (45), the overall exponent in (2) is known exactly for both surfaces.
As an application, consider the spreading of a cloud of admixed particles in a turbulent medium. The effective radius of the cloud of such particles at time t > 0 which started at t = 0 from the origin x = 0 is given by: Substituting the scaling representation (2) and taking into account relation (45), one readily arrives at the spreading law or, equivalently, where we used expression (44) derived for the case w * = 0. As already mentioned, the special choice ξ = 4/3 in the Kazantsev-Kraichnan ensemble (8) corresponds to the assumptions of the Kolmogorov-Obukhov theory of turbulence. This gives dR 2 /dt ∝ R 4/3 , which is nothing other than the statement of Richardson's "four-thirds" law [14].

Conclusions and Discussion
We applied the field theoretic RG to the problem of strongly nonlinear diffusion in a turbulent medium. The model is described by the advection-diffusion Equation (3) with the nonlinearity (7) and the random stirring noise (5) and (6). The turbulent incompressible velocity field is modelled by the Kazantsev-Kraichnan "rapid-change" ensemble with the correlation function (8).
We showed that the multiplicative renormalizability requires that an infinite number of interaction terms be taken into account in the advection-diffusion equation; as a result, the model involves infinitely many coupling constants. Nevertheless, we were able to construct the one-loop counterterm in an explicit closed form. This allows one to find, in the one-loop approximation, the full infinite set of the β functions entering the RG equations.
The analysis shows that the RG equations have two attractors in the form of a pair of two-dimensional surfaces of fixed points. The first sheet corresponds to irrelevance of turbulent advection while for the second sheet the advection is important. If those attractor surfaces contain some IR stability regions, the model demonstrates IR scaling behaviour of the form (1) and (2).
For the first sheet, the corresponding critical dimensions ∆ θ , ∆ θ and ∆ ω appear nonuniversal: they depend on the specific choice of a fixed point on the attractor surface. However, they satisfy certain exact relations (42) and (45).
For the more interesting second case the dimensions are shown to be universal (independent of the choice of a fixed point on the sheet) and are found exactly in (42) and (44). This is probably the most interesting result of our study.
The exact spreading law for a cloud of impurity particles is derived in (47) and (48). It is interesting to compare our results with those obtained in [15], where the dynamical stochastic NS equation was employed instead of the Kazantsev-Kraichnan ensemble (8). There, the RG equations also have two attractor surfaces (in fact, only the second one was discussed). For the second surface (advection is relevant), the critical dimensions are also universal and also are found exactly. However, there is a notable difference between the two models: the one studied in [15] involves an additional parameter, the ratio u of the viscosity coefficient in the NS equation (absent in the ensemble (8)) and the diffusivity coefficient ν in the advection-diffusion equation. As a result, the surfaces are parametrized by the two free parameters u * and g * 2 , while all g * n with n ≥ 3 are expressed in them. In order to further explore this distinction it would be interesting to consider the "intermediate" Gaussian velocity ensemble with finite correlation time, proposed in [37]. On the one hand, it involves the parameter u and the corresponding additional function β u . On the other, it contains the rapid-change ensemble (8) as a special limiting case. This work is already in progress.
From a more theoretical point of view, it is desirable to formulate the RG equations and to find their attractors directly in terms of the function V(θ), so that, instead of infinitely many β n functions for the couplings, we would have the single functional β(V) with the single functional argument V(θ); cf. the discussion in [38] for a similar case.
As a possible first step towards this goal, one can speculate that for a fixed point in the functional space the counterterm (30) may have the same functional form as the initial "input" function V in (7), so that V ∼ V /V . Possible solutions have the forms V(θ) tg (aθ + b) and V(θ) (1 + c exp(αθ))/(1 − c exp(αθ)) with some constants a, b, c, α (c > 0) and demonstrate pole singularities in θ. Thus, their physical significance is far from clear. Hopefully, application of the nonperturbative functional renormalization group will shed light on this interesting problem, in analogy with infinite-dimensional models of landscape erosion [33]. This work remains for the future.
Author Contributions: All authors contributed equally to this paper. All authors have read and agreed to the published version of the manuscript.
Funding: This research received no external funding. 1 To avoid possible misunderstandings, it should be stressed that the advection term is still present in the model. It can affect stability of the scaling regime, determine corrections to the leading-order IR asymptotic behaviour (1) and (2), contribute to various amplitude factors and so on.