Passive advection of percolation process: Two-loop approximation

TL;DR

This paper investigates the directed percolation process under external velocity fields using a field-theoretic approach, analyzing fixed points with a two-loop approximation to understand phase transition behavior.

Contribution

It introduces a two-loop renormalization group analysis of directed percolation in a velocity field, extending previous one-loop studies and incorporating Kraichnan ensemble effects.

Findings

Fixed point structure characterized at two-loop level

Velocity fluctuations influence percolation spreading

Enhanced understanding of phase transition dynamics

Abstract

The paradigmatic model of the directed percolation process is studied near its second order phase transition between an absorbing and an active state. The model is first expressed in a form of Langevin equation and later rewritten into a field-theoretic formulation. The ensuing response functional is then analyzed employing Feynman diagrammatic technique and perturbative renormalization group method. Percolation process is assumed to occur in external velocity field, which has an additional effect on spreading properties. Kraichnan rapid change ensemble is used for generation of velocity fluctuations. The structure of the fixed points structure is obtained within the two-loop approximation.

Tables

Table 1. Table 1: Canonical dimensions of the bare fields and bare parameters for the model

$Q$	$\psi$	$\tilde{\psi }$	$𝒗$	$D_0$	${\tau }_0$	${\lambda }_0$	$g_0$	$u_0$
$d_Q^k$	$d/2$	$d/2$	$-1$	$-2$	$2$	$\epsilon /2$	$\xi$	$\epsilon$
$d_Q^{\omega }$	$0$	$0$	$1$	$1$	$0$	$0$	$0$	$0$
$d_Q$	$d/2$	$d/2$	$1$	$0$	$2$	$\epsilon /2$	$\xi$	$\epsilon$

Equations

$${\partial }_t\psi =D_0\left[{\nabla }^2-{\tau }_0\right]\psi -\frac{{\lambda }_0{D}_0}{2}{\psi }^2+\zeta \sqrt{\psi },$$

$$⟨\zeta (t,\mathbf{x})\zeta (s,\mathbf{y})⟩={\lambda }_0{D}_0\delta (t-s){\delta }^d(\mathbf{x}-\mathbf{y})$$

$${\mathcal{S}}^{DP}=\tilde{\psi }[-{\partial }_t+D_0{\nabla }^2-D_0{\tau }_0]\psi +\frac{D_0{\lambda }_0}{2}\left[{\tilde{\psi }}^2\psi -\tilde{\psi }{\psi }^2\right],$$

$$\psi (t,\mathbf{x})⟷-\tilde{\psi }(-t,\mathbf{x}).$$

$$⟨\psi \tilde{\psi }{⟩}_0=\theta (t)\exp [-D_0(k^2+{\tau }_0)t],$$

$$⟨\psi \tilde{\psi }{⟩}_0=\frac{1}{-i\omega +D_0(k^2+{\tau }_0)}.$$

$$V_N(x_1,\dots ,x_N;\phi )=\frac{{\delta }^N{\mathcal{S}}_{\text{int}}[\phi ]}{\delta \phi (x_1)\dots \delta \phi (x_N)},\phi \in \{\tilde{\psi },\psi ,\mathbf{v}\},$$

$$V_{{\tilde{\psi }}^2\psi }=-V_{\tilde{\psi }{\psi }^2}=D_0{\lambda }_0.$$

$${\partial }_t\to {\nabla }_t={\partial }_t+(\mathbf{v}\cdot \nabla ),$$

$$⟨v_i(t,\mathbf{x}),v_j(t^{\prime },{\mathbf{x}}^{\prime })⟩$$

$$D_{ij}(\mathbf{x}-{\mathbf{x}}^{\prime })$$

$${\mathcal{S}}_{\text{vel}}=\frac{1}{2}\int \mathrm{d}{t}_1\int \mathrm{d}{t}_2\int {\mathrm{d}}^d{x}_1\int {\mathrm{d}}^d{x}_2{\mathbf{v}}_i(t_1,x_1)D_{ij}^{-1}(t_1-t_2,{\mathbf{x}}_1-{\mathbf{x}}_2){\mathbf{v}}_j(t_2,{\mathbf{x}}_2),$$

$$\mathcal{S}=\tilde{\psi }[-{\partial }_t-(\mathbf{v}\cdot \nabla )+D_0{\nabla }^2-D_0{\tau }_0]\psi +\frac{D_0{\lambda }_0}{2}\left[{\tilde{\psi }}^2\psi -\tilde{\psi }{\psi }^2\right]+\frac{1}{2}\mathbf{v}{D}^{-1}\mathbf{v}.$$

$$d_{\Gamma }=d+2-n_{\phi }{d}_{\phi },$$

$$u_0={\lambda }_0^2,$$

$${\mathcal{S}}_R=\tilde{\psi }[-Z_1{\partial }_t+Z_{2}D{\nabla }^2-Z_{3}D\tau ]\psi +\frac{Z_{4}D\lambda {\mu }^{\epsilon /2}}{2}[{\tilde{\psi }}^2\psi -\tilde{\psi }{\psi }^2]-Z_1\tilde{\psi }(\mathbf{v}\cdot \nabla )\psi +\frac{1}{2}\mathbf{v}{D}^{-1}\mathbf{v},$$

$$g_0{D}_0=gD{\mu }^{\xi },$$

$$Z_g{Z}_D=1,Z_g^{-1}=Z_D=Z_2{Z}_1^{-1}.$$

$$\psi \to \psi Z_{\psi },$$

$$D_0=D{Z}_D,$$

$$Z_1=Z_{\psi }{Z}_{\tilde{\psi }},Z_2=Z_D{Z}_{\psi }{Z}_{\tilde{\psi }},Z_3=Z_{\tau }{Z}_D{Z}_{\psi }{Z}_{\tilde{\psi }},Z_4=Z_u^{1/2}{Z}_D{Z}_{\psi }^2{Z}_{\tilde{\psi }}.$$

$$Z_{\psi }=Z_{\tilde{\psi }}=Z_1^{1/2},Z_D=Z_2{Z}_1^{-1},Z_{\tau }=Z_3{Z}_2^{-1},Z_u=Z_4^2{Z}_2^{-2}{Z}_1^{-1}.$$

$${\Gamma }_R\left(e,\mu ,\dots \right)=Z_{\phi }^{n_{\phi }}\Gamma \left(e_0,\dots \right),$$

$$\left[D_{RG}-n_{\phi }{\gamma }_{\phi }\right]{\Gamma }_R=0,$$

$$D_{RG}\equiv \mu {\partial }_{\mu }+{\beta }_u{\partial }_u+{\beta }_g{\partial }_g-\tau {\gamma }_{\tau }{\partial }_{\tau }-D{\gamma }_D{\partial }_D.$$

$${\gamma }_F={\tilde{D}}_{\mu }\ln Z_F,F\in \{\psi ,\tilde{\psi },D,\tau ,u,g\}.$$

$${\beta }_u\equiv {\tilde{D}}_{\mu }u=u\left[-\epsilon -{\gamma }_u\right],{\beta }_g\equiv {\tilde{D}}_{\mu }g=g\left[-\xi -{\gamma }_g\right].$$

$${\gamma }_1=2{\gamma }_{\psi },{\gamma }_2=2{\gamma }_{\psi }+{\gamma }_D,{\gamma }_3=2{\gamma }_{\psi }+{\gamma }_D+{\gamma }_{\tau },{\gamma }_4=3{\gamma }_{\psi }+{\gamma }_D+\frac{1}{2}{\gamma }_u.$$

$${\beta }_g(g^{\ast },u^{\ast })={\beta }_u(g^{\ast },u^{\ast })=0.$$

$$\Omega_{ij}=\frac{\partial\beta_{i}}{\partial g_{j}}\biggl{|}_{*},\quad i,j\in\{g,u\}.$$

$$Z_i=1+\text{ pole terms in }\epsilon ,\xi \text{ and their linear combinations}.$$

Full text

1 Pavol Jozef Šafárik University, Šrobárova 2, 041 54 Košice, Slovakia

2 Joint Institute for Nuclear Research, 141980 Dubna, Russia

3 Institute of Experimental Physics, Slovak Academy of Sciences, Watsonova 47, 040 01 Košice, Slovakia

4 SUPA, The University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9 3FD, United Kingdom

1 Introduction

In almost every realm of everyday life physical systems under non-equilibrium conditions are encountered. Mutual interplay between dissipative and driving forces give rise to a complicated and intriguing macroscopic behavior [1, 2, 3, 4]. Among most interesting, and at the same time very difficult to be tackled theoretically, are systems far from thermal equilibrium. Despite a lot of effort that has been made during last decades, fundamental understanding of non-equilibrium physics is still missing.

Reaction-diffusion problems appear commonly in biological systems and due to its very nature they could not be described by equilibrium statistical physics. Spreading constituents (atoms or more generally agents) interact with each other and thus the number of agents is not conserved. As control parameters are changing, it might happen that a underlying reaction scheme allows an existence of so-called absorbing state. Once system enters this state, it could not leave it. Clearly, this causes the system to be non-ergodic and thus impossible to study by equilibrium methods. It might happen that absorbing state is separated from an active state, i.e., a state with fluctuating (non-zero) constituents, by a critical point. At this point an emergent scaling behavior is observed, quite analogous to critical points in equilibrium systems [5, 6]. It is well-known that at critical point large scale spatio-temporal fluctuations govern the overall statistical properties and the resulting collective behavior can be effectively described by a certain set of continuous fields. The usual manifestation of criticality as a presence of divergences in various correlation functions is expected. A classical example is provided by the directed percolation process (DP), also known as Gribov process in hadron physics [2, 7, 8]. DP is mainly used as a simple model for a description of a population dynamics on the edge of extinction. Other possible applications embody high-energy physics, fluid turbulence, ecology and others [3, 9, 10]. In order for a system to be in corresponding universality class is a fulfillment of four conditions: (i) a unique absorbing state, (ii) short-ranged interactions, (iii) a positive one-component order parameter, (iv) no additional property (symmetry, presence of additional slow variables, etc.) [11, 12]. As a prominent example DP was studied by diverse analytical and numerical methods [2]. Therefore it is natural to consider DP when the main aim is to improve an existing method, what is part of our goal. Invaluable theoretical framework for an analysis of the scaling behavior is the renormalization group (RG) method [5, 6, 13]. In terms of RG flows of effective charges and accompanied existence of fixed points, divergent (power-law) behavior of various quantities can be naturally explained. Moreover, RG enables us with different computational approaches for an approximate estimation of universal quantities in a controllable fashion. Famous scheme consists in dimensional regularization augmented by so-called minimal subtraction (MS) scheme. For DP this yields a perturbative calculation in formally small parameter $\varepsilon$, where $\varepsilon$ is the deviation from the upper critical dimension $d_{c}=4$. Regarding this line of reasoning an existing research has been mostly restricted to the two-loop approximation of the perturbation theory [10, 14]. The main reason is that three-loop calculations put high demands on analytical methods [15, 16] and usually it is not possible to evaluate Feynman diagrams save by some efficient computational procedure.

Among the conditions of the DP universality class the item (iv) is probably most relevant from the experimental point of view. In realistic setups impurities and defects are expected to cause deviations from DP universality class. This is believed to be one of the reasons why there are not so many direct experimental realizations [17, 18, 19, 20] of DP. A study of deviations from the ideal situation could proceed in different routes and this still constitutes a topic of an ongoing debate [2]. A substantial effort has been made in studying a long-range interaction using Lévy flights [21, 22, 23], effects of immunization [10, 14], or in the presence of spatially quenched disorder [24]. Hence, as a further possible application of our methods, we analyze DP model in a presence of external velocity field that adds up to diffusion motion additional stirring effects. In this paper, we focus on DP in the presence of advective velocity fluctuations, which are generated by means of Kraichnan model. Such problem was first proposed in the work [32]. There the model was analyzed using field-theoretic renormalization group to a leading one-loop approximation.

Basic idea of the model is to assume that the velocity field can be imagined as a random Gaussian variable with prescribed statistical properties [25, 26, 27]. Despite obvious simplification in comparison to realistic flows, Kraichnan model is heavily used in a fluid dynamics. His role is especially important in intermittency studies, because it is one of the few models that allows an exact solution [27].

The remainder of the paper is organized as follows. In Sec. 2, we introduce a coarse-grained formulation of the DP problem, and we give a brief description of Kraichnan model for velocity fluctuations. Next we reformulate both models into a field-theoretic language. In Sec. 3, we present main steps of the perturbative RG analysis, and DP in presence of advection velocity field is renormalized to two-loop order. Sec. 4 is saved for a concluding summary.

2 Description of the Model

The stochastic reaction-diffusion equation for a positive coarse-grain density of percolating particles $\psi=\psi(t,{{\bm{x}}})$ has following form [2, 3]

[TABLE]

where $\partial_{t}=\partial/\partial t$, ${\bm{\nabla}}^{2}$ is the Laplace operator, $D_{0}$ is a diffusion constant, $\lambda_{0}$ is a positive coupling constant and $\tau_{0}$ is a deviation from the threshold value of injected probability. It can be interpreted as a formal analog of a deviation from critical temperature in static models [13]. Hereinafter the subscripts [math] will always indicate an unrenormalized (bare) quantity. The random Gaussian variable $\zeta(t,{{\bm{x}}})$ can be chosen in the following form [3]

[TABLE]

with $d$-dimensional version of Dirac $\delta(x)$ function, i.e. $\delta^{d}({{\bm{x}}}-{{\bm{y}}})=\delta(x_{1}-y_{1})\cdots\delta(x_{d}-y_{d})$. The average $\langle\cdots\rangle$ corresponds to a functional averaging over all noise realizations.

Further step consists in an utilization of famous De Dominicis-Janssen formalism [28, 29, 30] that allows us to map the stochastic problem (1)-(2) onto a field-theoretical model. Effectively one gets rid of noise variable, but on the other hand number of fields is doubled. This could be done in a standard fashion [3, 13] and the resulting response functional for the percolation process [2, 8, 10] reads

[TABLE]

where $\widetilde{\psi}$ is an auxiliary Martin-Siggia-Rose response field, and the integration over the spatio-temporal arguments is implicitly assumed. For instance, first term on the right hand side actually stands for the expression $\widetilde{\psi}\partial_{t}\psi=\int{\mathrm{d}}^{d}x\int{\mathrm{d}}t\,\widetilde{\psi}(t,{\bm{x}})\partial_{t}\psi(t,{\bm{x}})$.

The percolation model is manifestly invariant [2] with respect to so-called rapidity reversal symmetry

[TABLE]

This symmetry plays an important role in analysis of statistical quantities and among the practical consequences is a reduction of number of independent critical indices [2].

Main object of interest for the stochastic problem (1)-(2) are statistical quantities, which correspond to mean averages of expressions involving product of arbitrary number of fields. In field-theoretic formulation they are equivalently given as functional averages over the full set of fields with the “weight” functional $\exp(\mathcal{S}^{DP})$ and are known as correlation and response functions. These functions are conveniently represented in the diagrammatic form of Feynman graphs [5, 13].

Field-theoretic response functional (3) is amenable to a standard field-theoretic perturbative analysis. Free part of the response functional yields just one bare propagator $\langle\psi\widetilde{\psi}\rangle_{0}$, which takes the following form

[TABLE]

in time-momentum representation, whereas in frequency-momentum representation it takes the form

[TABLE]

Note that in (5) time flows from $\widetilde{\psi}$ to $\psi$ field. Non-linear terms in the response functional (3) give rise to two cubic vertices, whose vertex factors [13] can be obtained using general formula

[TABLE]

where $\mathcal{S}_{\text{int}}$ is a non-linear part of the response functional. It is easy to verify that from the action (3) we obtain two vertex factors

[TABLE]

In this paper, we analyze DP in a presence of additional velocity fluctuations. Basic underlying assumption is that DP is advected, but does not exert any backward influence on velocity field itself, i.e. it is an passive quantity [26]. Following works [26, 27, 31, 32] the turbulent mixing is described by velocity ensemble with prescribed statistics. Inclusion of velocity field ${{\bm{v}}}(t,{\bm{x}})$ corresponds to a replacement

[TABLE]

where $\nabla_{t}$ is the Lagrangian derivative [33]. The velocity field will be assumed incompressible, i.e. condition ${\bm{\nabla}}\cdot{\bm{v}}=0$ is fulfilled. According to Kraichnan suggestion we assume velocity field ${\bm{v}}$ to be a Gaussian random variable with zero mean and prescribed correlator

[TABLE]

where $P_{ij}(k)=\delta_{ij}-k_{i}k_{j}/k^{2}$ is transversal projection operator, $g_{0}$ is small coupling constant, and the cutoff $k=m$ provides an infrared regularization. Scaling exponent $\xi$ is related to a power-law behavior of energy spectrum [26], and in RG approach plays a role of a formally small expansion parameter. Since the velocity fluctuations are governed by the Gaussian statistics, the corresponding averaging procedure is performed with the quadratic functional

[TABLE]

where $D_{ij}^{-1}$ is the kernel of the inverse linear operation for the function $D_{ij}({{\bm{x}}}-{{\bm{x}}^{\prime}})$ in (11).

The full field theoretic model of the three fields $\varphi=\{\widetilde{\psi},\psi,{{\bm{v}}}\}$ is described by the response functional with the following abbreviated form

[TABLE]

We see that total response functional contains additional propagator $\langle vv\rangle_{0}$ and triple advection vertex $\widetilde{\psi}\psi v$. Its vertex factor (7) is proportional to the momentum $ik_{j}$ of auxiliary field $\widetilde{\psi}$. The graphical representation of this vertex can be found in Fig. 1.

3 Renormalization

A starting point of the perturbation theory is a free part of the response functional given by expression (13). By graphical means, it is represented as lines in the Feynman diagrams, whereas the non-linear terms correspond to vertices connected by these lines.

For the calculation of the RG constants we employ dimensional regularization in the combination with the modified minimal subtraction ($\overline{\text{MS}}$) scheme [5]. It must be borne in mind that now we are dealing with double expansion approach [34]. Therefore poles to the two-loop order that we encounter are of three types: either $1/\varepsilon$, $1/\xi$ and $1/(\varepsilon+\xi)$. This simple picture pertains only to the lowest orders in a perturbation scheme. In higher order terms, poles of general linear combinations in $\varepsilon$ and $\xi$ are expected.

The detailed examination of UV divergences is typically based on the analysis of canonical dimensions [5, 13]. Dynamical models have two scales, i.e. the canonical dimension of some quantity $Q$ is described by two values, the momentum dimension $d_{Q}^{k}$ and frequency dimension $d_{Q}^{\omega}$. First, it is needed to introduce normalization conditions $d^{k}_{{\bm{k}}}=-d^{k}_{{\bm{x}}}=1$, $d^{\omega}_{{\bm{k}}}=-d^{\omega}_{{\bm{x}}}=0$, $d^{k}_{t}=d^{k}_{\omega}=0$, $d^{\omega}_{\omega}=-d^{\omega}_{t}=1$. The dimensions are then found from the requirement that each term of the response functional remains dimensionless (with respect to the momentum and frequency dimensions separately). Further, the total canonical dimension $d_{Q}=d^{k}_{Q}+2d^{\omega}_{Q}$ plays the same role as momentum dimension in static models [13] and all canonical dimensions are given in Tab.1.

Crucial objects in an analysis of translationally invariant theories are $1$-irreducible Green function $\Gamma=\langle\varphi\dots\varphi\rangle_{\text{1-ir}}$, which are derived from connected Green function by an appropriate Legendre transformation [13]. Canonical dimension of $\Gamma$ is given by the relation

[TABLE]

where $n_{\varphi}=\{n_{\widetilde{\psi}},n_{\psi},n_{{\bm{v}}}\}$ represent the number of fields appearing in the function $\Gamma$. Superficial UV divergences can be generated only in those Green functions for which $d_{\Gamma}$ is a nonnegative integer. For the pure DP model [10] UV divergences are present in the following $1$-irreducible functions: $\langle\widetilde{\psi}\psi\rangle_{\text{1-ir}}$, $\langle\widetilde{\psi}\psi\psi\rangle_{\text{1-ir}}$, $\langle\widetilde{\psi}\widetilde{\psi}\psi\rangle_{\text{1-ir}}$ and their corresponding counterterms are already present in response functional (3).

By direct inspection of the Feynman diagrams one can observe that the real expansion parameter in perturbation theory is $\lambda_{0}^{2}$ instead of $\lambda_{0}$. This could be easily seen by a direct examination of Feynman diagrams and can be regarded as a direct consequence of the rapidity-reversal symmetry (4). For this reason it is convenient to introduce a new charge

[TABLE]

where $u_{0}$ has canonical dimension $4-d=\varepsilon$. The perturbative calculation is then made in terms of $u_{0}$.

The total renormalized response functional for DP in a presence of advecting velocity fluctuations takes the following form

[TABLE]

where $\mu$ is renormalization mass [5, 13]. The model is assumed to be in a scaling region, which is obtained for $\tau_{0}$ close enough to its critical value. In order to preserve the Galilean invariance [13] the advection term and the term containing temporal derivative (9) have to be renormalized by the same renormalization constant. In addition, the last term in the action (16) is not renormalized at all due to a passive nature of the advecting scalar field $\psi$. This ensures nonexistence of nontrivial Feynman diagrams for velocity propagator [26]. The amplitude factor is then expressed as

[TABLE]

and the renormalization constants are related as follows

[TABLE]

The renormalized response functional can also be obtained by renormalization of fields and parameters

[TABLE]

where we have symbolically expressed needed renormalization of fields. Relations among the renormalization constants take following form

[TABLE]

On the other hand, the renormalization constants for fields and parameters can be expressed by the inverse formulas

[TABLE]

As has been pointed out, passive nature of the problem ensures that the renormalization constant for velocity field $v$ is simply $Z_{v}=1$.

Next, we briefly show derivation of RG equation [5, 13] needed for an overall analysis of scaling behavior. The basic idea is corroborated by the claim that renormalized Green functions differ from unrenormalized ones by rescaling of the fields and choice of parameters. The fundamental relation $\mathcal{S}_{R}(\varphi,e,\mu)=\mathcal{S}(\varphi,e_{0})$ between response functional (3) and (16) leads directly to the formula

[TABLE]

where $e_{0}=\{D_{0},\tau_{0},u_{0},g_{0}\}$ is full set of bare parameters and $e=\{D,\tau,u,g\}$ are their renormalized counterparts, the ellipsis stands for the other arguments (time/frequency, coordinates/momenta of appearing fields). Further, we denote $\mu$-derivatives at fixed bare parameters by $\widetilde{D}_{\mu}=\mu\partial_{\mu}$. The equation $\widetilde{D}_{\mu}\Gamma=0$ then yields basic RG differential equation for the renormalized Green function $\Gamma_{R}$

[TABLE]

where $D_{RG}$ is the operation $\widetilde{D}_{\mu}$ expressed in the renormalized quantities

[TABLE]

The anomalous dimension $\gamma_{F}$ for any quantity $F$ is given by the relation

[TABLE]

$\beta$-functions for the dimensionless couplings $u$ and $g$ are

[TABLE]

The relations among the anomalous dimensions fields, parameters and $\gamma_{i}$ take the form

[TABLE]

Let us note that in what follows our main aim is to analyze phase structure of the theory. To this end the calculation of $\gamma_{3}$ is not needed and therefore we do not consider it here.

Large-scale (macroscopic) regimes of a given renormalizable field theoretic model are associated with IR attractive fixed points of the corresponding RG equations [5, 13]. A fixed point (FP) is defined as a such point $(g^{*},u^{*})$ for which $\beta$ functions $\beta_{g}$ and $\beta_{u}$ simultaneously vanish, i.e.

[TABLE]

The IR stability of a fixed point is then determined by the matrix of first derivatives of $\beta$ functions

[TABLE]

Asterisk in this equation indicates a corresponding fixed point value. The IR-asymptotic behavior is governed by IR-stable fixed points, for which real parts of all eigenvalues of matrix (28) are positive.

The renormalization constants absorb all divergences at $\varepsilon$, $\xi\rightarrow 0$ and the renormalized functions are finite for $\varepsilon$, $\xi=0$. For the RG calculation the modified minimal subtraction scheme $\overline{\text{MS}}$ has been chosen [35]. Difference between $\overline{\text{MS}}$ and minimal subtraction is that factor $S_{d}/(2\pi)^{d}$ is not expanded in $\varepsilon=4-d$. In the present case the RG constants attain a general form

[TABLE]

Coefficients in this expansion depend on two coupling constant $u$ and $g$. The coupling constants have been rescaled by a convenient constant factor

[TABLE]

Dyson equation for the function $\langle\widetilde{\psi}\psi\rangle$ reads

[TABLE]

where $\sum^{\text{perc}}$ contains only contributions arising solely from the pure DP process, whereas $\sum_{\widetilde{\psi}\psi}^{\text{adv}}$ corresponds to diagrams containing velocity propagator $\langle vv\rangle$. Two-loop 1-irreducible Feynman diagrams with nonzero contribution to Green function (31) are the following

[TABLE]

The counterterms in the $\overline{\text{MS}}$ scheme are polynomials in IR regulators. Furthermore, the RG constants $Z_{i};i=1,2,3,4$ could not depend on the choice of the IR regularization [5, 13]. In the case of the pure DP process $\tau$ is the IR regulator [10]. From the practical point of view, it is advantageous to set $\tau=0$ in the response functional (propagator $\langle\psi\widetilde{\psi}\rangle$ and cut off the momentum integrals at $k=m$ (by dimension $\tau\sim m^{2}$).

The renormalization constants can be expressed as follows

[TABLE]

Let us illustrate a calculation of anomalous dimension $\gamma_{1}$ using (24). As a first step we derive approximate relation

[TABLE]

where the last three terms stand for higher order terms that are neglected in what follows. Next, we need a formula $(\beta_{u}\partial_{u}+\beta_{g}\partial_{g})$, which can be approximated as follows

[TABLE]

where $\gamma^{(1)}_{u}$ and $\gamma^{(1)}_{g}$ are appropriate gamma functions up to the first order in perturbation theory. Finally, we get

[TABLE]

where terms proportional to $ug\xi/(\varepsilon(\varepsilon+\xi))$ and $u^{2}/(\varepsilon+\xi)$ drop out. In a similar fashion anomalous dimensions $\gamma_{2}$ and $\gamma_{4}$ can be calculated. Using them we derive in a straightforward way anomalous dimensions for $\psi$ field and diffusion constant $D$

[TABLE]

where the ratio $\varepsilon/\xi=\zeta$ is a finite quantity [34]. In the models with two regulators such as $\varepsilon$ and $\xi$, it is usually assumed that they are of the same order $\varepsilon=\mathcal{O}(\xi)$. For $g=0$ (DP model without advection interactions) expressions (37)-(38) coincide with known two-loop results [10, 11]. Further, Feynman diagrams contributing to RG constant $Z_{4}$ are the following

[TABLE]

where the factor two in front of the diagrams accounts for an internal symmetry of the graph, i.e. it expresses number of ways for drawing a given topological configuration. Altogether, the final expression for RG constant $Z_{4}$ reads

[TABLE]

Then, using (21) we can calculate anomalous dimension of the charge $u$

[TABLE]

This expression finalizes two-loop perturbative renormalization of the model. Using explicit information about $Z_{1},Z_{2}$ and $Z_{4}$ allows one to determine fixed points’ structure and thus find scaling regimes.

There are no other two-loop diagrams of an order $g^{2}$ due to a fact that in pure Kraichnan model for passive admixture all higher order corrections vanish [26]. Knowledge of anomalous dimensions $\gamma_{g}$ and $\gamma_{u}$ along with $\beta$-functions (25) allows us to perform a full two-loop analysis of fixed points’ structure. First and foremost there is a trivial or Gaussian fixed point FPI with coordinates

[TABLE]

This corresponds to a fixed point (FP) with irrelevant both DP interactions and advection process, and standard perturbation theory is applicable. As expected, this regime is IR stable in the region

[TABLE]

The former condition ensures that we are above the upper critical dimension $d_{c}=4$.

Next, there is a FP point FPII that corresponds to a pure DP process without advection. Its coordinates are

[TABLE]

Condition $g^{*}=0$ ensures that velocity propagator is effectively irrelevant. Eigenvalues of the matrix (28) are

[TABLE]

The first eigenvalue $\lambda_{1}$ agrees with a known two-loop result [10]. From an inspection of second eigenvalue $\lambda_{2}$ we observe that $\xi$ is restricted by a parabolic function of $\varepsilon$.

Coordinates of third FPIII are

[TABLE]

and it is IR stable in the region

[TABLE]

This FP corresponds to the pure advection process for which DP non-linearities are irrelevant.

Last FPIV is most interesting, because both non-linearities are IR relevant. Mutual interplay between DP and advection processes give rise to non-trivial behavior. The coordinates are

[TABLE]

and this FP is stable in region

[TABLE]

where the second inequality is obtained by numerical calculation with error smaller then $10^{-4}$. While the second-loop approximation does not qualitatively change one-loop results [32], we see that now the boundaries between the regions of stability are described not by lines in contrast to the one-loop result, but rather by parabolic curves. For a better visual aid, stability regions in $(d,\xi)-$plane are depicted in Fig. 2.

Note that only boundary between fixed points FPII and FPIV becomes parabolic due to two-loop corrections.

4 Conclusion

In this paper, we have investigated an effect of incompressible velocity fluctuations on the directed percolation process. Main points of field-theoretic formulation with inclusion of the advecting velocity field have been shown together with a renormalization group analysis.

We have established that depending on the values of a spatial dimension $d=4-\varepsilon$ and scaling exponent $\xi$, describing scaling properties of velocity fluctuations, the model exhibits four distinct universality classes. They correspond to: the Gaussian (free) fixed point, a directed percolation without advection, a passive scalar advection, and fully non-trivial regime, in which both percolation and advection interactions are relevant. All relevant quantities, such as fixed points coordinates, regions of stability and anomalous dimensions $\gamma_{\psi},\gamma_{D}$ and $\gamma_{u}$ have been calculated up to two-loop approximation. Despite obvious technical difficulties related to two-loop calculations, main physical consequences are in accordance with the previous one-loop result [32].

The purpose of this paper was two-fold. First, our aim was to improve existing results in non-equilibrium physics, which are mostly restricted to one-loop order. Second, this article may be considered as a first step in more challenging attempt, which would correspond to velocity field generated by some microscopic model such as stochastic Navier-Stokes equation in two-loop approximation.

Acknowledgment

The authors thank to Loran Ts. Adzhemyan, Juha Honkonen and Nikolay M. Gulitskiy for illuminating and fruitful discussions. The work was supported by VEGA grant No. 1/0345/17 of the Ministry of Education, Science, Research and Sport of the Slovak Republic, and the grant of the Slovak Research and Development Agency under the contract No. APVV-16-0186.

Appendix A: Explicit calculation of Feynman diagram

In this section we present main steps of a typical calculation of a divergent part of a Feynman diagram. Let us consider second diagram in an expansion (32) for 1-irreducible Green function $\langle\widetilde{\psi}\psi\rangle_{\text{1-ir}}$. We choose external parameters $p=({\bm{p}},\Omega)$ to enter the diagram from the left. There are internal variables $k=({\bm{k}},\omega_{k})$ and $q=({\bm{q}},\omega_{q})$ over which we have to integrate. Using the standard Feynman diagrammatic technique, we construct the following algebraic expression for the diagram

[TABLE]

where indices $m,n$ denote vector components of a velocity propagator (10). Using Cauchy integral formula integration over frequency variables $\omega_{k}$ and $\omega_{q}$ is readily performed. In addition, a straightforward simplification of the tensor structure is possible and finally we arrive at the following expression

[TABLE]

We are interested in UV divergent parts, which are known to be proportional to external frequency $\Omega$, square of external momentum ${\bm{p}}$ and mass term $\tau$. We expand (52) in a Taylor series and make a following substitution $\cos\theta\rightarrow z$ for a scalar product $({\bm{k}}\cdot{\bm{q}})=kq\cos\theta$ between internal momenta ${\bm{q}}$ and ${\bm{k}}$. In other words $\theta$ is the angle between the vectors ${\bm{k}}$ and ${\bm{q}}$.

[TABLE]

where $S_{d}=2\pi^{d/2}/\Gamma(d/2)$ is the surface of a $d$-dimensional sphere, and for a calculation of ${\bm{p}}^{2}$-term following formula

[TABLE]

was used [13]. Using a definition of charge $u$ from Eq. (15), substitution (30) and a relation $d=4-\varepsilon$, we finally obtain a expression for the UV divergent part of a diagram (51)

[TABLE]